text:
# Evidence for two-photon exchange contributions in electron-proton and positron-proton elastic scattering

## I INTRODUCTION

Measurements of $`G_E/G_M`$, the ratio of the proton electric and magnetic form factors, from Rosenbluth separation and polarization transfer techniques yield significantly different results Arrington (2003a) at large values of $`Q^2`$, the four-momentum transfer squared. The systematic uncertainties of both the Rosenbluth Arrington (2003a) and polarization Punjabi et al. (2003) measurements have been studied in detail, and no explanation for the discrepancy in terms of experimental problems has been found. If the discrepancy is not due to errors in the experiments or analyses, it may indicate a more fundamental problem with one of the techniques. Until this discrepancy is understood, there will be large uncertainties in our knowledge of the proton form factors. Since the polarization transfer measurements have only extracted the ratio $`G_E/G_M`$, cross section measurements are still needed to determine the magnitude of the individual form factors. So even if it is shown that the polarization transfer measurements are correct, and that the problem is due to unaccounted for corrections in the cross section measurements, there will still be uncertainties in the form factors until we fully understand these corrections Arrington (2003b).

Because the discrepancy grows rapidly with $`Q^2`$, it has typically been assumed that it is a problem with the cross section measurements, where a fixed error in the $`\epsilon `$ dependence of the cross sections would yield an error in $`(G_E/G_M)^2`$ that grows approximately linearly with $`Q^2`$. Assuming that the difference is due primarily to missing corrections in the cross section measurements, the discrepancy requires an error in the $`\epsilon `$ dependence of the cross section of approximately 5–8% for $`1<Q^2<6`$ GeV2 Arrington (2003a, b); Guichon and Vanderhaeghen (2003).

In order for a modification to the cross section to change the extracted form factor ratio, it must modify the $`\epsilon `$ dependence of the reduced cross section,

$$\sigma _R\frac{d\sigma }{d\mathrm{\Omega }}\frac{\epsilon (1+\tau )}{\sigma _{\mathrm{Mott}}}=\tau G_M^2(Q^2)+\epsilon G_E^2(Q^2),$$
(1)
where $`\tau =Q^2/4M_p^2`$, and $`\epsilon `$ is the longitudinal polarization of the virtual photon \[$`\epsilon ^1=1+2(1+\tau )\mathrm{tan}^2(\theta /2)`$\]. Several attempts have been made to find effects that might introduce an additional $`\epsilon `$ dependence to the measured cross section, thus modifying the extracted Rosenbluth form factors. Coulomb corrections, when implemented in a simple effective momentum approximation Higinbotham (2003), do modify the $`\epsilon `$ dependence of the cross section, but yield a very small effect compared to the size needed to explain the discrepancy. For the most part, investigations have focussed on the effect of two-photon exchange corrections Guichon and Vanderhaeghen (2003); Blunden et al. (2003); Rekalo and Tomasi-Gustafsson (2003) beyond the limited contributions that are already included in the traditional calculations of radiative corrections Mo and Tsai (1969); Tsai (1971); Meister and Yennie (1963).

While these works have shown that it is possible that a two-photon correction could explain the discrepancy, the only quantitative calculation Blunden et al. (2003) is limited to the elastic part of the two-photon contributions, i.e., the box and crossed-box diagrams considering only the case where the intermediate state is a proton, and yield only a 2% $`\epsilon `$ dependence, less than half the size necessary to explain the discrepancy. In this work, we reexamine positron measurements that were designed to test for two-photon contributions in elastic scattering in light of the possibility that they may be responsible for the discrepancy.

## II EXPERIMENTAL LIMITS ON TWO-PHOTON CONTRIBUTIONS

The effect of two-photon exchange terms can be observed in several processes. The imaginary part of the two-photon amplitude can, in principle, be measured in polarization observables. Measurements of the normal polarization, $`P_N`$, which is zero in the Born approximation, have been made Bizot et al. (1965); Giorgio et al. (1965); Lundquist et al. (1968), but no statistically significant indication of two-photon contributions has been seen. Similarly, the asymmetry $`A_N`$ has been measured for both elastic and inelastic scattering Chen et al. (1968); Powell et al. (1970), again with only null results. Thus far, the only observations of possible two-photon effects are in the asymmetry of scattering of transversely polarized electrons from protons. These asymmetries are extremely small, of order $`10^5`$, and so extremely difficult to measure. However, they have been observed by the SAMPLE experiment at MIT-Bates Wells et al. (2001) and the PVA4 Collaboration at Mainz Maas . However, these polarization observables are related to the imaginary part of the two-photon amplitude, while the cross section measurements are related to the real part. Therefore, while these data can be used to test models of the two-photon exchange, they do not directly constrain the two-photon contributions to the unpolarized cross sections, which might explain the discrepancy.

There are two ways to look for the effects of two-photon exchange corrections in the unpolarized elastic electron-proton cross section. First, one can look for deviations from the linear $`\epsilon `$ dependence in Eq. 1. There are a few Rosenbluth separation measurements Borkowski et al. (1974); Andivahis et al. (1994); Christy et al. (2004) that cover both large and small $`\epsilon `$ values and have small uncertainties (1–2%). These measurements do not show any significant deviations from linearity, but they have limited sensitivity because they have little data below $`\epsilon =0.4`$ and no data below $`\epsilon =0.2`$. Data from different experiments can be combined to expand the $`\epsilon `$ range, but normalization uncertainties between different experiments reduce the significance of such tests, while attempts to normalize across the datasets Walker et al. (1994); Arrington (2003a) rely on the linearity when determining normalization factors. So while these data set limits on non-linearities at large $`\epsilon `$, the limits are less significant at low $`\epsilon `$. In addition, these measurements are insensitive to two-photon correction terms that are constant or vary linearly with $`\epsilon `$. Such corrections would modify the extracted values of $`G_E`$ and $`G_M`$ without spoiling the linearity of the Rosenbluth plot.

The second approach is to compare positron-proton scattering to electron-proton scattering. For positron-proton scattering, the interference term between the one-photon and two-photon amplitudes changes sign, yielding a ratio $`R\sigma (e^+p)/\sigma (e^{}p)1+4\text{Re}(B)/A`$, where $`B`$ is the two-photon amplitude, and $`A`$ is the one-photon amplitude Mar et al. (1968). The modification to the electron cross section is approximately $`12\text{Re}(B)/A`$, and so any change in the electron cross section will yield roughly twice the change in $`R`$, but with the opposite sign. In the simplest approximation, one expects an additional factor of $`\alpha `$ in the two-photon amplitude relative to one-photon amplitude, yielding corrections to the electron cross section of roughly $`2\alpha 1.5\%`$, and to the ratio R of $`3`$%. Additional differences come from Bremsstrahlung corrections where proton recoil is taken into account, but these are included in the usual radiative corrections. An analysis by Mar and collaborators Mar et al. (1968) found these additional differences to be relatively small, typically less than 1–2% for their kinematics, and to be identical to better than 0.3% in different prescriptions Meister and Yennie (1963); Tsai (1961) of the radiative corrections.

Figure 1 shows the existing data Yount and Pine (1962); Browman et al. (1965); Anderson et al. (1966); Bartel et al. (1967); Cassiday et al. (1967); Anderson et al. (1968); Bouquet et al. (1968); Mar et al. (1968) for the ratio of positron-proton to electron-proton elastic cross sections as a function of $`Q^2`$. While there is some hint of a $`Q^2`$ dependence, the large $`Q^2`$ data have large uncertainties, and a fit to the data of the form $`R=a+bQ^2`$ yields $`b=0.0085\pm 0.0063`$, less then 1.5 standard deviations from zero. A fit of the ratios to a constant value yields $`R=1.003\pm 0.005`$, with $`\chi _\nu ^2=0.87`$, which corresponds to a two-photon correction to the electron cross sections of ($`0.15\pm 0.25`$)%. This result has been interpreted to mean that the two-photon corrections must be even smaller than the naive estimate, limiting the effect on the electron-proton cross section to less than one percent.

However, the low intensity of the secondary positron beams used in these experiments makes it difficult to perform precise measurements where the cross section is small. Because of this, the data in Fig. 1 are limited to low $`Q^2`$ values ($`<1`$ GeV2) or small scattering angles ($`\epsilon >0.7`$), where the cross section is large. While the existing data do place tight limits on the size of two-photon corrections in some regions, they do not place any limits on two-photon contributions at low $`\epsilon `$ except at relatively low $`Q^2`$ values ($`<`$1 GeV2). So it is still possible that the discrepancy in the extracted form factors is due to two-photon corrections to the cross sections, if the correction is only large for small $`\epsilon `$ values.

If we assume that the two-photon corrections are responsible for the discrepancy between polarization transfer and Rosenbluth measurements, we can make specific predictions about how these corrections would impact the positron measurements, and use this to examine the existing data more carefully. In order to explain the discrepancy, the effect must increase the slope of the Rosenbluth plot, and so must increase the cross section at large $`\epsilon `$ relative to the low $`\epsilon `$. Based on the size and $`Q^2`$ dependence of the discrepancy, the $`\epsilon `$ dependence in the electron cross section must be 5–8%, depending only weakly on $`Q^2`$, for $`Q^2>2`$ GeV2. It must also be reasonably close to linear in $`\epsilon `$, or else it would introduce visible non-linearities in the Rosenbluth plot. This implies that the ratio $`R`$ should have a 10–15% $`\epsilon `$ dependence, approximately linear in $`\epsilon `$, and of the opposite sign as in the electron cross section, i.e., the positron to electron ratio must either increase at small $`\epsilon `$ or decrease at large $`\epsilon `$.

Unfortunately, there is very little positron data above $`Q^2`$=2, and it covers a very limited $`\epsilon `$ range. The data from Mar, et al. Mar et al. (1968), has four points above $`Q^2=2`$ GeV2, yielding $`R=1.034\pm 0.024`$. These data are all at large $`\epsilon `$ values ($`\epsilon =0.88`$), and so do not exclude significant two-photon corrections at large $`Q^2`$. If the two-photon correction is small at $`\epsilon =0`$, then a 10–15% decrease in $`R`$ at large $`\epsilon `$ would be necessary to explain the discrepancy, and this is clearly ruled out by the high $`\epsilon `$ data. So any two-photon corrections would have to increase $`R`$ (decrease the electron cross section) at low $`\epsilon `$ in order to explain the discrepancy and still be consistent with the positron data.

The positron data with significant $`\epsilon `$ range is limited to $`Q^2<2`$ GeV2. Figure 2 shows these data as a function of $`\epsilon `$, and a significant $`\epsilon `$ dependence can be seen. A linear fit, neglecting any $`Q^2`$ dependence, yields an slope of $`(5.7\pm 1.8)`$%, with $`\chi ^2=11.1`$ for 22 degrees of freedom. The extremely low $`\chi ^2`$ indicates that the uncertainties in the data have most likely been overestimated, and that the effect may be more significant than indicated by the fit uncertainty.

The observed increase in the positron to electron ratio at low $`\epsilon `$ corresponds to an increase of 2.8% in the observed slope in the Rosenbluth extraction. This implies that the value of $`G_E`$ extracted from the Rosenbluth separation will be larger than the true value, while the extracted $`G_M`$ value will be smaller. While the 2.8% $`\epsilon `$ dependence is only half the size necessary to explain the discrepancy at large $`Q^2`$ value, data covering a wide range of $`\epsilon `$ values are only available at low $`Q^2`$. The average $`Q^2`$ value of the data in Fig. 2 is only 0.5 GeV2, and the lower $`Q^2`$ data are generally more precise, making the weighted average $`Q^2`$ value less than 0.4 GeV2.

We can estimate the $`\epsilon `$ dependence necessary to explain the discrepancy in the form factors at large $`Q^2`$, but at these low $`Q^2`$ values the polarization transfer and Rosenbluth form factors are not precise enough to determine if there is an inconsistency, and so cannot be used to estimate the size of the two-photon corrections. A decrease in the size of the $`\epsilon `$ dependence at low $`Q^2`$ could easily yield the slope observed in Fig. 2, yet still be large enough to fully explain the discrepancy between polarization and Rosenbluth extractions at larger $`Q^2`$ values. At larger $`Q^2`$ values, where the size of the corrections can be estimated from the discrepancy, the effect decreases somewhat as $`Q^2`$ decreases, and is approximately 5% for $`Q^2=1`$–2 GeV2. In addition, the correction must become smaller for very small $`Q^2`$ values ($`0.010.1`$ GeV2), or the decrease in the low-$`\epsilon `$ cross sections would lead to significant reductions in the extracted values of $`G_M`$. The extractions of $`G_M`$ are not precise enough to conclude that the corrections must go to zero, but they must be significantly smaller than the 5% corrections observed at larger $`Q^2`$ values. Thus, we expect this slope in the positron to electron comparison to be less than 10%.

A global analysis of the cross section and polarization transfer data was used to try and estimate the low $`Q^2`$ behavior. In Ref. Arrington (2003b), a global analysis of the cross section and polarization transfer data, assuming a fixed 6% $`\epsilon `$-dependent correction to the cross section, was used to extract the ‘Polarization form factors’. A modified version of this global analysis was performed, but rather than extracting $`G_E`$ and $`G_M`$ with a fixed two-photon correction, we extract $`G_E`$, $`G_M`$, and the $`Q^2`$ dependence of the slope of the linear $`\epsilon `$-corrections. Several different functional forms were tried, and a range of curves, which all gave good fits, are shown in Fig. 3. While the fits were not constrained to go to zero, they all yield a much smaller value as $`Q^20`$. When these curves are used to estimate the $`\epsilon `$ dependence for the correction to the electron cross section at $`Q^2=0.4`$ GeV2, they yield slopes of (1.8–3.3)%, implying a slope in the positron to electron ratio of $``$(3.7–6.8)%, in agreement with the observed $`5.7`$% slope.

The $`\epsilon `$ dependence extracted above assumed no $`Q^2`$ dependence to the size of the correction, and a simple linear $`\epsilon `$ dependence. While it is in agreement with the estimated $`\epsilon `$ dependence from the form factor discrepancy, the estimate relied on an extrapolation to lower $`Q^2`$ values. However, while the above analysis had to make some assumptions about the $`Q^2`$ dependence of the two-photon effects, we can also make some significant model-independent statements from this data.

For the low $`\epsilon `$ data ($`\epsilon <0.5`$) $`Q^2=0.5`$ GeV2 (weighted average) and $`R=1.049\pm 0.014`$, clearly demonstrating that two-photon effects decrease the electron cross section at low $`\epsilon `$ and low $`Q^2`$. Note that only one point below $`\epsilon =0.5`$ is above $`Q^2=1`$ GeV2, and it has a positron to electron ratio of $`1.079\pm 0.046`$. This is consistent with the 10–12% increase that would explain the discrepancy in form factor measurements for $`Q^2=`$1–2 GeV2, but also only two standard deviations from $`R=1`$.

In addition to the observation of a significant $`\epsilon `$ dependence at low $`Q^2`$ values, the data also sets significant limits on possible two-photon exchange corrections at large $`Q^2`$ for $`\epsilon >0.8`$. For $`Q^2>1`$ GeV2, $`R=1.020\pm .015`$. So the 95% confidence region for $`R`$ at large $`\epsilon `$ is 0.99–1.05, yields limits on the electron cross section modification of -2.5% to +0.5%. To increase the slope of the Rosenbluth plot, they would have to increase the high-$`\epsilon `$ electron cross section, and such an enhancement is limited to $`<0.5`$%.

## III IMPLICATIONS OF THE POSITRON MEASUREMENTS

The $`\epsilon `$ dependence of the two-photon effects seen here is consistent with the calculations of Refs. Blunden et al. (2003); Afanasev , which have small corrections at large $`\epsilon `$, and a significant decrease of the electron cross section at low $`\epsilon `$. However, it rules out the form of Ref. Rekalo and Tomasi-Gustafsson (2003) as an explanation for the discrepancy between polarization and Rosenbluth extractions of $`G_E/G_M`$. In Ref. Rekalo and Tomasi-Gustafsson (2003), the authors predict a specific $`\epsilon `$ dependence, which is zero at $`\epsilon =0`$, and grows rapidly at large $`\epsilon `$. If the size of the correction is made small enough to be consistent with the constraints from the positron measurements at large $`\epsilon `$, then the two-photon effects at smaller $`\epsilon `$ values will be negligible.

The positron measurements are also inconsistent with the corrections obtained in Ref. Guichon and Vanderhaeghen (2003), at least given the specific approximations that the authors use to obtain the two-photon effects from the discrepancy in form factor measurements. They do not assume single photon exchange, but instead write a more general expression for the cross section in terms of two generalized form factors, $`\stackrel{~}{G}_E`$ and $`\stackrel{~}{G}_M`$, along with a third term, $`\stackrel{~}{F}_3`$, which is zero in the Born approximation. By assuming that the two-photon contributions are negligible in $`\stackrel{~}{G}_E`$ and $`\stackrel{~}{G}_M`$, they extract values of $`Y_{2\gamma }`$, a dimensionless parameter related to the size of $`\stackrel{~}{F}_3`$, such that the effect of $`Y_{2\gamma }`$ on the cross section and polarization transfer data resolve the discrepancy. Under this assumption, the two-photon effects on the cross section measurements are approximately proportional to $`\epsilon `$, and are $`>`$5% at $`\epsilon =1`$ for $`Q^2>1`$ GeV2. This would yield $`R0.9`$ for the large $`Q^2`$ positron measurements, which is clearly ruled out by the data (Fig. 1). It is possible that the two-photon effects in their $`Y_{2\gamma }`$ terms could be cancelled by two-photon effects that modify $`\stackrel{~}{G}_E`$ and $`\stackrel{~}{G}_M`$. If this is the case, the formalism may still allow a connection between the two-photon effects in polarization transfer and cross section measurements, but it is no longer possible to determine the two-photon terms directly from the extraction of $`Y_{2\gamma }`$, because of the sizable two-photon contributions to $`\stackrel{~}{G}_E`$ and $`\stackrel{~}{G}_M`$.

This observation of the form of the two-photon effects can also be used to assist in the extraction of form factors from Rosenbluth and polarization transfer data. To have a consistent extraction of the form factors from the cross section and polarization data, we have to assume something about the nature of the discrepancy. An $`\epsilon `$ dependence of the cross section, of the form observed in the positron data, is consistent with the assumption used in Ref. Arrington (2003b). In this case, it was assumed that the cross sections were modified by two-photon exchange terms that were zero at $`\epsilon =1`$, linear in $`\epsilon `$, and large enough to explain the discrepancy. This assumed modification to the cross sections was removed to correct for the two-photon effects, with no correction at $`\epsilon =1`$, and a 6% increase in the cross sections at small $`\epsilon `$ ($`\sigma _{c2}`$ of Ref. Arrington (2003b)). The size of this correction is such that the Rosenbluth data approximately reproduce the polarization transfer values of $`G_E/G_M`$, and the 6% increase in the $`\epsilon =0`$ cross section yields a value of $`G_M`$ that is approximately 3% higher than a direct Rosenbluth extraction from the unmodified cross sections (e.g. the parametrization of Ref. Bosted (1994) or the ‘Rosenbluth form factors’ of Ref. Arrington (2003b)).

A similar combined extraction of form factors in Ref. Brash et al. (2002) used the polarization transfer values for $`G_E/G_M`$ to fix the slope of the reduced cross section, and used the uncorrected cross sections to extract the magnitude of $`G_E`$ and $`G_M`$. For a data point at $`\epsilon =1`$, the change in the assumed slope of Eq. 1 yields an increase in the extrapolation to $`\epsilon =1`$ of 5–8%, the size implied by the discrepancy, and so gives similar results to the extraction of Ref. Arrington (2003b). For a measurement at very low $`\epsilon `$, the change in slope leaves the $`\epsilon =0`$ extrapolation unchanged. So depending on the mean $`\epsilon `$ value of the data in a given $`Q^2`$ range, the extracted value of $`G_M`$ will be 0–4% lower in this analysis, compared to Ref. Arrington (2003b), with a typical difference of 1–2%, as there are more data at large $`\epsilon `$ values. Note that both of these combined extractions rely on the assumption that the two-photon exchange terms have little or no effect on the polarization transfer results.

## IV CONCLUSIONS

This is the first direct experimental evidence for large two-photon corrections in the unpolarized elastic electron-proton cross section. The effect is only observed for low $`Q^2`$ values, and so cannot be directly compared to the two-photon contributions necessary to explain the discrepancy between Rosenbluth and polarization transfer measurements of the proton form factors. However, the size and $`\epsilon `$ dependence of these effects are consistent with simple estimates based on the observed discrepancy, and so this observation supports the idea that two-photon contributions may significantly modify the Rosenbluth extraction of nucleon form factors.

Additional comparisons of positron to electron scattering over a range in $`Q^2`$ and $`\epsilon `$ would provide the most direct extraction of these two-photon corrections. With precise measurements over an adequate range in $`\epsilon `$ and $`Q^2`$, we could determine if the two-photon effects can fully explain the difference between polarization transfer and Rosenbluth measurements of the form factors, and could also provide significant data with which to constrain models of the real part of the two-photon amplitude. However, at the present time it is unclear where such a program could be carried out over the necessary kinematic range, and with the precision needed to map out these corrections. In the meantime, the existing data can be used to test calculations of the two-photon effects, and already are sufficient to rule out approaches with large electron cross section enhancements at large $`\epsilon `$ values. These data also provide information about the size and $`\epsilon `$ dependence of the two-photon effects, important when attempting to extract the proton form factors from a combined analysis of Rosenbluth and polarization transfer data.

###### Acknowledgements.

This work was supported by the U. S. Department of Energy, Nuclear Physics Division, under contract W-31-109-ENG-38.


format: text

text:
# Adaptive quantum design of atomic clusters

## I Introduction

In the near future, it will likely be possible to control the precise spatial positions of atoms and molecules using the experimental techniques now being developed by nanoscience. To complement these emerging capabilities it is clear that a new set of theoretical tools has to be developed to assist in the exploration of a potentially vast number of atom configurations and a corresponding enormous range of physical properties. In this paper, we outline an approach, which we call “adaptive quantum design”, that sets out to address this challenging task. In contrast to classical systems, atomic scale devices exhibit quantum fluctuations and collective quantum phenomena caused by particle interactions. Besides offering an excellent testing ground for models of correlated electrons, they also force us to reconsider conventional paradigms of condensed matter physics, such as crystal symmetries that are imposed by nature. In some instances, such symmetries need to be explicitly broken in order to enable or optimize a desired system response. Consider for example the quasiparticle density of states in tight-binding systems, which is the subject of the present study. In translationally invariant structures, i.e. crystals, it is well known that the spectral response function exhibits van-Hove singularities at positions of low dispersion, such as the band edges in a one-dimensional chain or the band center in a two-dimensional square lattice. These enhancements of the density of states can be very useful in amplifying system responses such as optical conductivity at specific incident energies. It is therefore important to be able to control the positions and shapes of such features by adaptive design techniques applied to models which capture the essential degrees of freedom of interacting atomic clusters.

Traditional ad-hoc methods for the design of nanoscale devices will likely miss many possible configurations. At the same time, it is unrealistic to expect individuals to manually explore the vast phase space of possibilities for a particular device function. The proposed solution to this difficult design problem is to employ machine-based searches of configuration space that enable user-defined target functions. Adaptive quantum design solves the inverse problem by numerically identifying the best broken-symmetry spatial configuration of atoms and molecules that produce a desired target function response.

The two major ingredients of adaptive quantum design are the physical model, which in this work evaluates the electronic density of states for a particular spatial arrangement of the atoms, and the search algorithm that finds the global minimum in the parameter space of all possible configurations. This problem is typically highly underdetermined in the sense that there can be several atomic configurations that yield a system response very close to the desired target function. Often, the associated landscape of solutions is shallow and has many nearly degenerate local minima.

Adaptive quantum design therefore relies heavily on efficient algorithms that accurately model the interacting system and find its optimal configuration, i.e. the global minimum in the available parameter space that yields the best match to the desired target function. This target response may be a specific angular transmission of a photonic crystal, a fine-tuned Josephson current along a superconducting junction array, or an energy-dependent density of states profile of an atomic cluster - the case on which we focus here. In this work, various numerical search tools are applied, including guided random walk, simulated and triggered annealing, and genetic algorithms. In particular the last set of techniques is most efficiently implemented on parallel computers. It is observed that in practice hybrids of these methods yield the best results, i.e. fast convergence towards a global minimum.

The particular example of a long-range tight-binding model is chosen for this study because it captures essential features of correlated quantum mechanical systems, and yet permits fast numerical diagonalization of relatively large clusters with broken translational symmetry. In a typical optimization run, these “function calls” occur 100-10,000 times. This model is therefore suitable for developing and testing adaptive design algorithms, and should be viewed as an initial step towards the design of atomic clusters. The goal of this work is to test the applicability and limits of adaptive design techniques on a simple but non-trivial quantum system.

## II Long-range tight-binding model

The tight-binding approach is an effective tool to describe the band structure of electronic systems. It is commonly used to model the relevant bands close to the Fermi level, obtained from complex density functional theory calculations. Since in this work symmetry-breaking, non-periodic configurations are considered, a long-range variant of the tight-binding model with overlap integrals depending on the variable inter-atomic distance has to be used. Its Hamiltonian is given by

$`H={\displaystyle \underset{i,j}{}}t_{i,j}\left(c_i^{}c_j+c_ic_j^{}\right),`$ (1)
where $`c_i^{}`$ and $`c_i`$ are creation and annihilation operators at a site $`\stackrel{}{r}_i`$, and the sum over pairs of atoms is restricted to avoid double-counting. Here the spatial decay of the overlap integral $`t_{i,j}`$ is parametrized by a power-law,

$`t_{i,j}={\displaystyle \frac{t}{|\stackrel{}{r}_i\stackrel{}{r}_j|^\alpha }},`$ (2)
where the exponent $`\alpha `$ is taken to be 3.0 throughout the paper unless mentioned otherwise. harrison  This parametrization reflects an algebraic variation of the overlap integral with inter-atomic separation. The choice of sign in the Hamiltonian follows the convention for s-orbitals. However, this simple implementation of the tight-binding model does not account for the orbital directionality of realistic Au, Ag, or Pd atomic wave functions. More sophisticated and numerically expensive techniques, such as the local density approximation, would be needed to make quantitative predictions for these systems.

The Hamiltonian matrix of the long-range tight-binding model in the basis of single-particle states is non-sparse, and only its diagonal matrix elements vanish. In order to obtain the spectrum, the matrix is diagonalized numerically for finite clusters. In Fig. 1, the resulting densities of states of translationally invariant chains and square lattices are shown for the nearest-neighbor tight-binding model with $`t_{i,j}=t\delta _{i,j}`$ in Figs. 1(a) and (c), and for the case of long-range overlap integrals in Figs. 1(b) and (d). Characteristic van-Hove singularities are observed, in one dimension at the band edges, and in two dimensions at the band center. For the long-range model, the particle-hole symmetry is broken because of frustration introduced by the competing overlaps, leading to asymmetries in the density of states. While the system sizes in Fig. 1 are chosen to be rather large in order to make contact with the familiar thermodynamic limit, there are still some visible finite-size remnants, i.e. a faint pole structure due to the discreteness of the system. These features become much more pronounced for the few-atom clusters that are studied in the next sections of this paper.

## III Target functions and convergence criterion

While “nature” gives us densities of states that are constrained by the dimensionality and symmetry of the underlying lattice, our objective is to engineer spectral responses with specific shapes that are useful in designing nanoscale devices. For example, we may wish to produce a quasi-two-dimensional spectrum in a one-dimensional system or to concentrate spectral weight in particular energy windows. These goals are achieved by placing the atomic constituents into optimized symmetry-breaking configurations which are determined by numerical searches.

The specific target functions that are studied in this work are shown in Fig. 2. They are (a) a flat top hat density of states centered at $`E=0`$,

$`N(E)_{target}=\theta (EE_c)\theta (E+E_c),`$ (3)
where $`\theta (x)`$ is the Heavyside function, and $`E_c`$ is an energy cut-off, (b) a symmetric two-peak function, centered at $`E=0`$,

$`N(E)_{target}`$ $`=`$ $`\theta (EE_{c2})\theta (E+E_{c1})`$ (4)
$`+`$ $`\theta (EE_{c1})\theta (E+E_{c2}),`$
and (c) a particle-hole-symmetry-breaking function with two unequal peaks, i.e. more spectral weight on the quasi-hole side ($`E<0`$) than on the quasi-electron side ($`E>0`$) of the spectrum,

$`N(E)_{target}`$ $`=`$ $`\theta (EE_{c2})\theta (E+E_{c1})`$ (5)
$`+`$ $`\theta (EE_{c4})\theta (E+E_{c3}).`$
In systems with finite numbers of tight-binding atoms, these continuous shapes are approximated by equally spaced poles within the energy windows where $`N(E)_{target}`$ is a non-vanishing constant. Here, the delta-functions are given a finite width of 0.02t. As more atoms are added to the system, these peaks merge together, approaching the bulk result. For other non-flat target functions the quasiparticle peak spacing can be varied, e.g. following a Gaussian or Lorentzian shape. Naturally, not all targets can be achieved equally well. Factors that influence the achievable match to a target include the number of available number of atoms and, as will be shown, a continuous or discrete number of accessible spatial positions. The dimensionality of the system poses additional constraints. In particular, there are less available configurations in lower dimensions. This study focuses on three prototype spectral responses which are targeted by numerically optimizing configurations of clusters with up to 48 atoms in two spatial dimensions.

The optimization algorithms seek to minimize the deviation from a given target density of states, defined by the error function

$`\mathrm{\Delta }={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E[N(E)N(E)_{target}]^2`$ (6)
which is the least-square difference between the system response for a given configuration and the target response function. We have explored a number of numerical techniques, including the Newton-Raphson steepest descent method, simple downhill random walk, simulated and triggered annealing, and genetic algorithms. The advantages and disadvantages of these techniques are briefly discussed in the appendix. In general, it is found that hybrids of these methods tend to work best. In the next section, we focus on adaptive design of atomic clusters in continuous configuration space without any restrictions to underlying discrete lattices. In this case, the search space is infinite which generally allows better convergence to a given target response than in finite configuration space. However, some experimental realizations of such structures require deposition of atoms on substrates with discrete lattice structures.nazin ; nilius ; wallis  Therefore, this case is addressed separately in the subsequent section.

## IV Atoms-up design of tight-binding clusters in continuous configuration space

The first target density of states we would like to study in detail is the particle-hole symmetric top hat function, centered at energy E=0, shown in Fig. 2(a). This function represents a constant density of states for a bulk solid between the energy cut-offs $`\pm E_c`$, giving a bandwidth $`W=2E_c`$. Here, we choose $`E_c=3t`$. However, it should be noted that with the adaptive quantum design approach we are not restricted to this choice, and target densities of states with quite different bandwidths can be matched, although often to a lesser degree of accuracy.

In Fig. 3, the solution for a system with 16 atoms confined to a box with periodic boundary conditions is shown. The guided random walk method is applied to optimize the configuration of atoms by iterative local updates of their positions in order to match the top hat density of states. As observed in Figs. 3(a) and (c), good convergence to the target function can be achieved for this case after less than 2000 updates. A contour plot of the potential $`_{ij}t_{ij}/|\stackrel{}{r}_i\stackrel{}{r}_j|^\alpha `$ for the resulting spatial configuration is shown in Fig. 3(b). Here, equipotential lines are used to denote the positions and overlapping wave functions of the atomic constituents in the system. For this target function one discovers the formation of dimers with a wide range of inter- and intra-dimer spacings. These self-assembled building blocks have variable directional orientation, and are closely packed.

To explore the generic aspects of this result, let us now turn to the two-peak target function shown in Fig. 2(b). This density of states has a gap in the center of the spectrum, as may be desired for the construction of two-state systems such as q-bits. As shown in Figs. 4(a) and (c), the convergence to this target function is not quite as good as for the top hat function, but saturation occurs already at half the number of iterations compared to the previous example. Interestingly, the optimized spatial configuration that is found in Fig. 4(b) also displays a preference for dimer formation.

In the long-range interacting systems we are studying it is, at least at first sight, not obvious why these target functions prefer dimer building blocks. Let us address this question by examining the individual spectra of the molecular building blocks shown in Fig 5. Each dimer molecule contributes to the density of states a positive and negative pole with energies $`E=\pm t_{12}`$, where $`t_{12}`$ is the hopping integral between the two participating atoms. The zero-energy quasiparticle peaks of two isolated atoms are split into bonding and anti-bonding combinations once they form dimer molecules. Therefore, isolated dimers are ideal building blocks for particle-hole symmetric densities of states, such as the top hat and the two-peak target function. The intra-dimer spacing determines the positions of the $`E=\pm t_{12}`$ poles via Eq. 2. With an appropriate distribution of these distances the full target spectrum of the top hat function can be covered. The poles close to the band center ($`E=0`$) are provided by the less tightly bound dimers. For the two-peak target function the dimer building blocks required to realize the target spectrum are more tightly bound, and the intra-dimer spacings need to vary less to achieve this target.

The idea that dimers can be used as building blocks for the particle-hole symmetric target functions strictly applies only to isolated dimers, i.e. the dilute limit, or when potential gradients across dimers pairs from the presence of adjacent atoms do not break particle-hole symmetry. The absence of such gradients, even for relatively high atom densities in a long-range interacting system, accounts for the success of dimers in satisfying the target function.

Lower symmetry building blocks, such as trimers and quadrumers shown in Fig. 5, can achieve more complex target functions, in particular those with broken particle-hole symmetry. Due to frustration, trimer molecules intrinsically have asymmetric densities of states with unequal spectral weights on the electron and the hole side of their spectra. Quadrumers have a symmetric spectral response in the absence of longer-range frustrating interactions. As an example of how these building blocks enable more complex target functions, let us consider the asymmetric two-peak density of states given by Eq. 5 and shown in Fig. 2(c), which has a narrow upper peak and a wider lower peak, separated by a gap. In Fig. 6 it is observed that an approximate match can be achieved by the adaptive method. As expected, the building blocks for this particle-hole-asymmetric target function are combinations of dimers, trimers, and quadrumers, which partially recombine into larger clusters. Obviously, in order to achieve asymmetric target functions more complex building blocks are required. Especially for systems which contain only a small number of atoms it is important whether the required building blocks are available, and whether there are non-participating unbound atoms that may deteriorate convergence and match to the target function.

Let us finish this section by addressing the convergence properties of the three cases that were discussed. So far, only the simplest guided random walk optimization method was considered in which every “downhill step” is accepted. These “downhill steps” are local updates of individual atomic position that lead to a better match of the spectrum to the target density of states. Especially for the most symmetric target function, this algorithm converges efficiently, with a small remaining error. For the symmetric two-peak function, convergence occurs even faster because the required dimer building blocks are more uniform than for the first case. However, the remaining error is slightly larger, indicating that this may be a metastable solution which could be improved by global updates in which whole subclusters are simultaneously updated. Finally, the convergence plot for the asymmetric two-peak target function (Fig. 6(c)) shows several plateaus, indicating metastable configurations that exhibit a high resistance against local updates. For this case, a much larger number of local updates is required to achieve an acceptable match. Hence, more complex target functions clearly call for more sophisticated numerical search tools, including annealing steps, parallelization, and global updating schemes when available.

## V Adaptive design in discrete configuration space

In the previous section local updates were considered which allow atomic positions to change continuously within a given radius. However, for the case of atoms deposited on a substrate with a given lattice structure, the set of available positions is usually discrete, although it may be very large. This has significant consequences for the adaptive design approach. The search space of solutions is finite in this case, which makes it feasible to study more atoms with similar computational effort compared to the continuous case. At the same time, the discreteness of the lattice can prohibit favorable configurations that are available in the continuous case, thus deteriorating convergence properties of the optimization procedure. In practice, we find that the feasibility of computations for larger numbers of atoms in the discrete case helps to achieve better matches, as long as the lattice spacings remain sufficiently small. In this section, we explore the effects of an underlying grid on the adaptive design procedure. Also, more advanced techniques are implemented for the numerical optimization, including hybrids of the genetic algorithm, simulated annealing, and the guided random walk method.

In Fig. 7, adaptive design results are shown for systems with 48 atoms on a square lattice with spacing 0.01. The length scale is set by the linear size of the two-dimensional box (48 $`\times `$ 48) to which the particles are confined. These systems are in the dilute limit, and hence the convergence to the target functions, chosen to be the same as in the previous section, is good. Also, because of the larger number of particles, there are less finite size effects. For the top hat target function one obtains a small matching error of $`\mathrm{\Delta }=`$ 0.000643, for the symmetric two-peak function one finds $`\mathrm{\Delta }=`$ 0.000605, and for the asymmetric two-peak function the error is $`\mathrm{\Delta }=`$ 0.150347. Thus, analogous to the continuous case, the more symmetric targets are easier to achieve. Again, clustering into dimers is observed for the particle-hole symmetric target functions, whereas trimers are the preferred building blocks for the asymmetric target. Since atomic densities in these examples were chosen to be in the dilute limit, boundary effects and interactions between the building blocks are small. The resulting configurations have the character of liquids, governed by weak interactions between the molecular building blocks, and relatively strong confining forces that lead to the formation of dimers and trimers.

Next, let us explore the dependence of these solutions on the coarseness of the underlying lattice. In Fig. 8, the grid spacing is varied over two orders of magnitude from 0.01 up to 1.0. As expected, the convergence to the target top hat function deteriorates dramatically as the substrate is made coarser. For the smallest spacing of 0.01 one converges to a final error of $`\mathrm{\Delta }=`$ 0.28457, for a spacing of 0.3 the error is $`\mathrm{\Delta }=`$ 2.45047, for a spacing of 0.7 the error becomes $`\mathrm{\Delta }=`$ 3.81779, and for the coarsest case with spacing 1.0 the error is $`\mathrm{\Delta }=`$ 36.3721, indicating failure of convergence. The corresponding configurations in Fig. 8 show a strong dependence of the clustering sizes as the system is trying to cope with less available phase space to match the target function. For the finest grid spacing (Fig. 8(a)) one observes almost entirely isolated dimers. As the spacing is increased (Fig. 8(b)), a few small strings and groups are formed. At grid spacing 0.7 (Fig. 8(c)), the solution is made up mostly of long strings and dimers in close proximity to each other. Ultimately, for the coarsest grid spacing of 1.0, (Fig. 8(d)), the final configuration consists of square and rectangular blocks of atoms. This result demonstrates that there is a hierarchy of building blocks. The most suitable building blocks for the given target function are dimers. When these become less available due to lattice constraints, the adaptive method selects higher order solutions, i.e. larger size clusters, in order to cope with the more restricted phase space of possible configurations. Hence, for each level of coarseness, adaptive design discovers solutions that enable - up to a given degree of accuracy - a targeted system response.

By widening the grid spacings and keeping the box size constant, the density of atoms in the system is effectively increased, and simultaneously the available energy levels are spaced further apart. Let us examine these dependencies by varying the power-law governing the atomic overlaps (Eq. 2), and by increasing the number of atoms, i.e. the filling fraction, on a fixed lattice. Results for the achieved convergence are shown in Fig. 9. These demonstrate that excellent target matches can be achieved in the dilute limit with filling densities of a few percent. For larger coverages, the numerical search becomes exponentially less effective, indicating increasing frustration effects that need to be addressed by global updating schemes. Higher power-laws imply effective shorter-range atomic overlaps, thus rendering the system more dilute. This is reflected in Fig. 9, where the departure from the regime of negligible matching errors is pushed toward higher filling percentages as $`\alpha `$ in Eq. 2 is increased.

Some aspects of the long-range tight-binding model on a substrate have already been confirmed experimentally using scanning tunneling microscopy (STM) to precisely position gold atoms on the surface of a nickel-aluminum crystal. In a recent study performed at UCI by Dr. Ho’s group,nilius  STM measurements show that the splitting in the value of eigenenergies for Au dimers on NiAl depends inversely on Au atom separation corresponding to $`\alpha =1`$. Here, the grid spacings are dictated by the 0.29 nm lattice periodicity of available add-atom sites on the NiAl substrate. Remarkably, the power-law dependence of the effective overlaps $`t_{ij}`$ between the deposited atoms takes into account interactions with the substrate which are typically difficult to model by first principle computations.

## VI Conclusions

In this work, adaptive quantum design techniques were applied to tailor the quasiparticle density of states of atomic clusters, modeled by the long-range tight-binding Hamiltonian. Broken-symmetry spatial configurations of atoms were optimized to match target spectra. By applying adaptive search algorithms, it was shown that matches to target responses can be achieved by forming hierarchies of molecular building blocks that depend on system constraints. For example, symmetric top hat and two-peak target densities of states can be achieved by forming lattices of weakly interacting dimers. While these are the elementary building blocks for particle-hole symmetric case target functions, more complex molecules, such as trimer and quadrumers, are found to dominate the solutions for asymmetric target functions. Implementation of this approach on discrete lattices, corresponding to finite configuration space, introduces frustration effects that distabilize the elementary building blocks in the limit of coarse grid spacings.

The core task of adaptive quantum design is the numerical search for global minima in typically shallow landscapes of configurations with many local minima. Since this procedure typically requires many function calls, an efficient implementation on parallel computers is necessary. For the purposes of this study, the complexity of the physical model was minimized in order to limit the computational expense. While the long-range tight-binding model can be viewed as a semi-realistic testing ground for adaptive quantum design techniques, it is crucial to apply these algorithm to more sophisticated models that include, among other ingredients, orbital directionality, spin degrees of freedom, and electronic correlations. Furthermore, adaptive design is applicable to related areas in nanotechnology, including the design of electro-optic componentsyu  and RF systemsgheorma .

Our enthusiasm for adaptive quantum design is, in part, driven by the potential for scientific and technological discovery built into the methodology. We find it appealing to perform machine-based searches for new configurations, new components, and new sub-systems, as we believe these will inevitably create new understanding, heuristics, and intuition. The fact that searches are performed in an ultra-large configuration space virtually guarantees the discovery of completely new designs and a much more thorough exploration of the possibilities and capabilities of nanotechnology.

## VII Acknowledgements

We are grateful to Ioan Gheorma, Weifei Li, Peter Littlewood, Omid Nohadani, Tommaso Roscilde, and Rong Yu for discussions. We acknowledge financial support by DARPA and the Department of Energy, Grant No. DE-FG03-01ER45908. Computational support was provided by the USC Center for High Performance Computing and Communications and by the NERSCC.

## VIII Appendix: Optimization and search algorithms

The optimization algorithms used in this work are based on local updates of atomic positions in order to minimize the error $`\mathrm{\Delta }`$ defined in Eq 6. Each atom in the system is visited periodically, and a trial change of its position is attempted. Depending on the response in $`\mathrm{\Delta }`$ and on the specific algorithm, this trial step is either accepted or rejected. For the case of continuous configuration spaces, these local updates are random shifts of positions within a given radius. A hard core constraint is implemented which forbids atoms to be placed on top of each other. For discrete configuration space, a stochastic distribution function is used to decide which sites in the neighborhood of the original position of an atom should be visited in a trial step. While this function is naturally peaked at nearest-neighbor sites, it has to include a finite probability of longer range updates in order to avoid getting stuck in local minima of search space. The particular results discussed in this appendix are obtained for the continuous case.

The Newton-Raphson and Broydn methods are multi-dimensional generalizations of the one-dimensional secant methodpress . Unfortunately, their global convergence is rather poor for more than 20 variable parameters. Therefore, they are only applicable to the smallest system sizes, and are not useful for the non-linear multi-parameter searches required for adaptive quantum design.

In the guided random walk or “downhill method”, each random step that results in a smaller target error $`\mathrm{\Delta }`$ is accepted.yu  Random steps are trial spatial variations about a particle’s position. This guided random walk technique is a quickly implemented power horse. However, especially for shallow landscapes of solutions it gets easily stuck in local minima.

Simulated annealing uses an effective temperature representing the likelihood of accepting a step that does not minimize the function. This temperature is lowered slowly - at a rate of 10% - with each iteration. The initial temperature is taken as $`T_{init}=5t`$. This method is better at avoiding local minima than the previous techniques. However, it takes a relatively long time to converge.

The triggered annealing method is a hybrid of the downhill and the simulated annealing method. It implements the downhill method until minimizing steps become hard to find, at which point the simulated annealing method is used to escape from local minima. Parameters are chosen the same as for the simulated annealing method. Triggered annealing tends to converge relatively quickly if there are only a few local minima.

The particle replacement method uses the simple downhill method for guided random updates. In addition, it identifies particles that have not been updated for an extended period, because of being stuck in a local minimum, and assigns them to a new random position within the lattice boundaries. In our implementation, a particle is replaced if there are 10 idle iterations without a successful downhill update for that particular particle. Note that this particle replacement update is different from random step updates because it is independent of the previous position of the particle.

In genetic algorithms a population of possible solutions is created. Those that best minimize the function are allowed to take part in creating a new generation of possible solutions.whitley  These methods are generally good at avoiding local minima, and are also easily implemented on parallel computers. They typically require more function calls than other search algorithms.

In order to illustrate the efficiency of these various approaches each method is used to match a flat top hat target function on a one-dimensional lattice with 24 tight-binding atoms and a box size of 96. Particles exiting the box on one end, enter it from the other side via periodic boundary conditions. The target density of states is chosen to be the symmetric top-hat function, with poles evenly spaced between $`E=4t`$ and $`4t`$. All of these benchmark runs are started with atoms placed randomly along a chain. The computer is a Pentium III 1 GHz with 2GB pc133 memory, and the additional seven nodes used by the genetic algorithm method are Pentium III 850 MHz processors with 1GB memory. Each method is run for 600 seconds.

As shown in Fig. 10, all methods, with the exception of the genetic algorithm, converge rapidly within the first minute of run time, and show only relatively small corrections afterwards. The two techniques with the fastest convergence (inset of Fig. 10) are the downhill and the simulated anealing methods. However, their asymptotic error functions remain relatively large, indicating that they get easily stuck in local minima. In contrast, the triggered annealing and particle replacement methods yield much better matches to the target function, while still converging relatively fast. As shown in the inset, the annealing methods sometimes accept trial steps in the “wrong” direction in order to avoid local minima. Finally, the genetic algorithm takes a relatively long time to converge. However, it yields by far the best match to the target density of states after about 7 minutes run time. In order to ensure best matches to the target, this last method is therefore used whenever the computational effort allows it.


text:
# Spectral Geometry of Operator Polynomials and Applications to QFT

## 1 Introduction

Consider a class of eigenvalue problems which have the following form:

$$P(\omega )\varphi _\omega (P_n+\omega P_{n1}+\mathrm{}+\omega ^{n1}P_1+\omega ^nP_0)\varphi _\omega =0,$$
(1.1)
where $`\omega `$ is a complex spectral parameter and $`\varphi _\omega `$ is an eigenfunction. $`P_k`$ ($`k=0,\mathrm{},n`$) are partial differential operators of $`k`$-th order acting on a Hilbert space $``$. $`P(\omega )`$ represents an example of what is called the operator polynomial or polynomial operator pencil . In general $`P_k`$ and $`P_l`$ do not commute if $`kl`$ and (1.1) is essentially non-linear problem with respect to the eigenvalue $`\omega `$.

Problems like (1.1) appear in different applications (oscillations in a viscous fluid, Schrödinger equation with energy dependent potential and etc), see  for the introduction and references. The focus of the present work is on applications related to quantum field theory.

Consider equation of motion for a free field $`\varphi (x,x^0)`$ ($`x`$ and $`x^0`$ denote spatial and time coordinates, respectively). It can be always written as

$$P(i_0)\varphi (x,x^0)(P_2+P_1(i_0)+P_0(i_0)^2)\varphi (x,x^0)=0,$$
(1.2)
where $`_0=/_0`$. Operators $`P_k`$ depend only on spatial derivatives and may have additional indexes related to a spin-tensor structure. Field $`\varphi `$ propagates in an external classical background. If the background is stationary the variables can be separated by choosing $`\varphi (x,x^0)=e^{i\omega x^0}\varphi _\omega (x)`$. From (1.2) one gets

$$P(\omega )\varphi _\omega (x)=(P_2+\omega P_1+\omega ^2P_0)\varphi _\omega (x)=0.$$
(1.3)
This is a quadratic operator polynomial, a particular case of problem (1.1).

Eigenvalues $`\omega `$ determine the spectrum of energies of single-particle excitations which is used to compute a number of important physical quantities, like the vacuum (Casimir) energy $`E=\pm \frac{1}{2}|\omega |`$, free energy and etc.

Consider for instance a charged Dirac field $`\psi `$ in an external static electromagnetic field with the potential $`A_\mu =A_\mu (x)`$. The equation is

$$[\gamma ^\mu D_\mu +m]\psi =0,$$
(1.4)
where $`D_\mu =_\mu ieA_\mu `$ and $`e`$ is the charge. To get for the single-particle spectrum problem like (1.3) one has to act by the operator $`[\gamma ^\mu D_\mu m]`$ on the left hand side (l.h.s.) of (1.4) and make substitution $`\psi (x,x^0)=e^{i\omega x^0}\psi _\omega (x)`$. This yields

$$[D_kD^km^2+e^2A_0^2\frac{e}{2}(2i\gamma ^0\gamma ^k_kA_0)+\omega 2eA_0+\omega ^2]\psi _\omega (x)=0,$$
(1.5)
where $`=(_jA_k_kA_j)\gamma ^j\gamma ^k`$ and $`k,j=1,2,3.`$ When the time component $`A_0`$ of the potential is a function of spatial coordinates (1.5) cannot be reduced to a standard eigenvalue problem because operator $`P_1=2eA_0`$ does not commute with $`P_2`$.

Sec. 2 of this contribution contains a summary of results obtained in . The aim is to show that the notion of spectral geometry can be extended to a class of quadratic operator polynomials which appear in problems of quantum field theory. Applications to finite temperature theories are discussed in Sec. 3 where, as an example, the obtained results are used to describe the Debye screening in hot plasma. Finally, some important features of operator polynomials for scalar, spinor and gauge fields are analyzed in Sec. 4.

## 2 Spectral Geometry of Quadratic Polynomials

To describe the results we rewrite (1.3) in the form

$$\left[\omega ^2L(\omega )\right]\varphi _\omega (x)=0,$$
(2.1)
$$L(\omega )=L_2+\omega L_1+\omega ^2L_0,$$
(2.2)
where $`L(\omega )=\omega ^2P(\omega )`$ and $`L_k`$ are partial differential operators of the $`k`$-th order. Suppose for a moment that $`\omega `$ is a complex parameter and consider a one-parameter family of operators $`L(\omega )`$. Let $`\mathrm{\Lambda }_k(\omega )`$ be eigenvalues of $`L(\omega )`$ for a given $`\omega `$,

$$L(\omega )\varphi _{\mathrm{\Lambda }_k}^{(\omega )}=\mathrm{\Lambda }_k(\omega )\varphi _{\mathrm{\Lambda }_k}^{(\omega )}.$$
(2.3)
If the spectrum $`\mathrm{\Lambda }_k(\omega )`$ is known for any $`\omega `$, the spectrum of problem (2.1) is determined by roots of equation

$$\chi (\omega ,\mathrm{\Lambda }_k)=0,$$
(2.4)
$$\chi (\omega ,\mathrm{\Lambda }_k)=\omega ^2\mathrm{\Lambda }_k(\omega ).$$
(2.5)
Denote these roots by $`\omega _{k,i}`$ and define the following function on the spectrum

$$\chi ^{}(\omega _{k,i})=_\omega \chi (\omega _{k,i},\mathrm{\Lambda }_k)\text{where}\omega _{k,i}^2=\mathrm{\Lambda }_k(\omega _{k,i}).$$
(2.6)
It is assumed that at first the derivative in $`\chi ^{}(\omega _{k,i})`$ over $`\omega `$ is taken for the fixed branch of eigenvalues $`\mathrm{\Lambda }_k(\omega )`$ and after that the result is considered on the corresponding root $`\omega _{k,i}`$. In what follows we write $`\omega `$ instead of $`\omega _{k,i}`$.

We make three assumptions: i) $`L(\omega )`$ is a Laplace type operator of the form:

$$L(\omega )=(_k+iA_k+i\omega a_k)(^k+iA^k+i\omega a^k)+\omega B+V$$
(2.7)
which acts on sections to vector bundles over a $`d`$-dimensional compact Riemannian manifold $`_d`$; ii) the spectrum of $`L(0)`$ is strictly positive; iii) the function $`\chi ^{}(\omega )`$, see (2.6), is positive (negative) on positive (negative) eigenvalues $`\omega `$.

First two conditions are technical (see details in ). The origin of condition (iii) is explained in Sec. 4. At this point we note that in the limit when the polynomial is trivial, $`L(\omega )=L_2`$, this condition is satisfied because $`\chi ^{}(\omega )=2\omega `$.

To describe the spectrum of (2.1) at large $`|\omega |`$ we introduce the pseudo-trace

$$K(t)=\frac{1}{2}\underset{\omega }{}e^{t\omega ^2},t>0,$$
(2.8)
where summation goes over the real eigenvalues $`\omega `$. The asymptotics of $`K(t)`$ at small $`t`$ is related to the distribution of large $`|\omega |`$. Given conditions (i)–(iii) the following asymptotic series takes place at small $`t`$

$$K(t)\frac{1}{(4\pi t)^{d/2}}\underset{n=0}{\overset{\mathrm{}}{}}\left[a_nt^n+b_nt^{n+1/2}\right].$$
(2.9)
$`a_n`$ and $`b_n`$ can be computed by using expansion formula for the heat kernel

$$K_\omega (t)=\text{Tr}e^{tL(\omega )}\frac{1}{(4\pi t)^{d/2}}\underset{n=0}{\overset{\mathrm{}}{}}\left[a_n(\omega )t^n+b_n(\omega )t^{n+1/2}\right]$$
(2.10)
which is well-known . For operators defined in (2.7),

$$a_n(\omega )=\underset{m=0}{\overset{n}{}}a_{m,n}\omega ^m,b_n(\omega )=\underset{m=0}{\overset{n}{}}b_{m,n}\omega ^m$$
(2.11)
and it can be shown that

$$a_n=\underset{m=n}{\overset{2n}{}}(1)^{nm}\frac{\mathrm{\Gamma }\left(\frac{d}{2}+m\right)}{\mathrm{\Gamma }\left(\frac{d}{2}+n\right)}a_{2(mn),m},$$
(2.12)
$$b_n=\underset{m=n}{\overset{2n}{}}(1)^{nm}\frac{\mathrm{\Gamma }\left(\frac{d1}{2}+m\right)}{\mathrm{\Gamma }\left(\frac{d1}{2}+n\right)}b_{2(mn),m}.$$
(2.13)
Coefficient $`a_0`$ is proportional to the volume of $`^d`$, other $`a_n`$, $`b_n`$ are integrals of local functionals of background fields and $`b_n`$ are non-trivial when $`^d`$ has a boundary. The notion of spectral geometry does extend to a class of quadratic operator polynomials.

The proof of (2.9), (2.12), (2.13) is based on relation

$$K(t)=\frac{1}{4\pi }_0^{\mathrm{}}𝑑\omega _C𝑑ze^{\omega ^2(tiz)}\left(2\omega +\frac{1}{iz}_\omega \right)(K_\omega (iz)+K_\omega (iz)),$$
(2.14)
where a contour $`C`$ goes in the complex plane just below the real axis $`\text{Im}z=0`$. Thus $`K(t)`$ is determined by part of $`K_\omega (t)`$ which is symmetric with respect to the transformation $`\omega `$ to $`\omega `$. This property is a direct consequence of condition (iii) which we imposed on the spectrum of the operator polynomial. This also explains why coefficients $`a_n`$ and $`b_n`$ in (2.12), (2.13) are determined by the symmetric part of coefficients $`a_n(\omega )`$ and $`b_n(\omega )`$.

## 3 Applications to Finite-Temperature Theories

Spectral asymptotics of quadratic operator polynomials can be used in a number of physical applications. One of them is studying the behaviour of the free energy of non-interacting quanta at high temperatures. The free energy for Fermi particles at temperature $`T=\beta ^1`$ is defined as

$$F(\beta )=\beta ^1\underset{k}{}\mathrm{ln}\left(1+e^{\beta E_k}\right).$$
(3.1)
Summation goes over single-particle energies $`E_k=|\omega _k|`$ which are determined by the real eigenvalues $`\omega _k`$ of the corresponding wave equations.

It is convenient to start with the case when the single-particle spectrum is determined by the standard linear eigenvalue problem, i.e., when operator $`L(\omega )=L_2`$ and it does not depend on $`\omega `$. This case is well studied. It is known that in three dimensions $`d=3`$ the free energy in the high temperature limit has the following asymptotic expansion

$$F(\beta ,A)\frac{7\pi ^2}{720}T^4a_0\frac{1}{48}T^2a_1\frac{1}{16\pi ^2}\mathrm{ln}(T/\varrho )a_2+O(T^2)$$
(3.2)
(we neglect boundary terms for simplicity). Here $`a_n`$ are the coefficients of the asymptotic expansion of the trace $`\text{Tr}e^{tL_2}`$, $`\varrho `$ is a dimensional parameter.

In the high temperature limit $`F(\beta )`$ is determined by behaviour of the single-particle spectrum at large $`|\omega |`$. For polynomial (2.1) this behaviour is described by asymptotic formula (2.9) which has exactly the same form as expansion of $`\text{Tr}e^{tL_2}`$. Therefore, for systems whose spectrum is determined by (2.1) formula (3.2) can be used if $`a_n`$ are calculated with the help of (2.12) (analogously for boundary terms if they are present).

Let us use this method to derive the free energy of electron-positron plasma in an external static electromagnetic field $`A_\mu (x)`$. (Finite-temperature theories in gravitational backgrounds are discussed in ). Coefficients $`a_n`$ are determined by (2.12) for problem (1.5). The corresponding operator is

$$L(\omega )=𝒟_k𝒟^k+m^2e^2A_0^2+\frac{e}{2}(2i\gamma ^0\gamma ^k_kA_0)2eA_0\omega .$$
(3.3)
By using standard results  one easily finds coefficients of $`L(\omega )`$ and obtains

$$a_0=a_{0,0}=4_^3d^3x,$$
(3.4)
$$a_1=a_{0,1}+\frac{1}{2}a_{2,2}=4_^3d^3x(m^2+2e^2A_0^2).$$
(3.5)
$$a_2=a_{0,2}\frac{1}{2}a_{2,3}+\frac{3}{4}a_{4,4}=\frac{2e^2}{3}_^3d^3xF_{\mu \nu }F^{\mu \nu },$$
(3.6)
$`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu `$. Consider the effective thermodynamic potential of the system at high temperature by neglecting the zero-point energy

$$\mathrm{\Omega }(\beta ,A)=\mathrm{L}(A)+F(\beta ,A)=\frac{1}{4}_^3d^3xF_{\mu \nu }F^{\mu \nu }+F(\beta ,A),$$
(3.7)
where $`\mathrm{L}(A)`$ is the classical Lagrangian of the background potential $`A_\mu `$. Eqs. (3.2), (3.4)–(3.6) yield

$$\mathrm{\Omega }(\beta ,A)=_^3d^3x\left(\frac{c(T)}{4}F_{\mu \nu }F^{\mu \nu }+\frac{1}{2}m^2(T)A_0^2\right),$$
(3.8)
$$c(T)=1\frac{e^2}{24\pi ^2}\mathrm{ln}(T/\rho ),m^2(T)=\frac{1}{3}e^2T^2.$$
(3.9)
These formulas reproduce the result which is well-known in plasma physics. The last term in (3.8) is responsible for screening effects in hot plasma and $`m(T)`$ is the correct expression for the Debye mass.

Note that conventional derivation of (3.8) is based on computation of one-loop diagrams at finite temperature in the limit of vanishing gauge fields. Our calculation is done for arbitrary static gauge potential.

Coefficient $`a_2`$ in (3.6) has a Lorentz invariant structure. For this reason the logarithmic term in expansion (3.2) results in a finite renormalization of the background potential (see the first term in r.h.s. of (3.8)). The Lorentz invariant structure of $`a_2`$ is not a mere coincidence. It reflects a more general property established in : coefficient $`a_2`$ of the pseudo-trace expansion coincides with the coefficient $`A_2`$ of the asymptotic expansion of the heat kernel of the covariant four-dimensional operator $`P(i_0)=_0^2L(i_0)`$.

## 4 Operator Polynomials for Different Field Theories

Spectral asymptotics given in Sec. 2 are obtained by using certain conditions. In particular, the structure of coefficients $`a_n`$ and $`b_n`$ in (2.12), (2.13) is determined by condition (iii). Let us analyze the origin of this condition for models of non-interacting fields.

Consider a pair $`f_1`$, $`f_2`$ of solutions to field Eq. (1.2) (or (1.4), for Dirac case) and introduce the relativistic product, a Hermitian bilinear form $`f_1,f_2=f_2,f_1^{}`$ defined on a space-like hypersurface $`\mathrm{\Sigma }`$ in a hyperbolic spacetime $``$ as

$$f_1,f_2=_\mathrm{\Sigma }𝑑\mathrm{\Sigma }^\mu j_\mu (f_1,f_2).$$
(4.1)
The ”current” $`j_\mu (f_1,f_2)`$ is divergence free, $`^\mu j_\mu =0`$, and (4.1) does not depend on $`\mathrm{\Sigma }`$.

If the background is stationary one can introduce a conserved canonical energy $`H[f]`$ for a given solution $`f`$. It can be shown that for free scalar, spinor and gauge fields the energy can be written in the form

$$H[f]=\frac{i}{2}f,_0f+c.c.,$$
(4.2)
where ”$`c.c.`$” is a term obtained from the first term in the r.h.s. of (4.2) by complex conjugation. It follows from (4.2) that for a solution $`f_\omega `$ with certain frequency $`\omega `$ ($`i_0f_\omega =\omega f_\omega `$)

$$H[f_\omega ]=\omega f_\omega ,f_\omega .$$
(4.3)
In general the system may have excitations with complex frequencies $`\omega `$. It is easy to see that they have vanishing relativistic norm $`f_\omega ,f_\omega `$ and, hence, the vanishing energy. Such modes along with zero-frequency modes are not quantized and should be considered separately from real-frequency solutions. This is the reason why pseudo-trace (2.8) is defined for real $`\omega `$.

### 4.1 Scalar Fields

Consider integer spin fields. To have a well-defined theory we require that canonical energy is non-negative,

$$H[f_\omega ]0,$$
(4.4)
for any $`\omega `$. If there are single-particle excitations with negative energy $`H`$ the Fock space in the corresponding quantized theory contains states with negative energies unbounded from below.

Condition (iii) of Sec. 2 for scalar fields $`\varphi `$ follows from (4.4). Indeed it is possible to show  that in this case

$$\varphi _\omega ,\varphi _\omega =\chi ^{}(\omega )(\varphi _\omega ,\varphi _\omega ),$$
(4.5)
where $`\chi ^{}(\omega )`$ was introduced in (2.6) and $`(\varphi _\omega ,\varphi _\omega )=𝑑V|\varphi _\omega |^2`$ is a positive-definite norm in a Hilbert space of single-particle wave functions. (4.4) together with relations (4.3), (4.5) implies that $`\chi ^{}(\omega )=\epsilon (\omega )|\chi ^{}(\omega )|`$ where $`\epsilon (\omega )`$ is the sign function.

### 4.2 Gauge Fields

We discuss now theory of non-Abelian gauge fields with the gauge group $`SU(N)`$. Let $`B_\mu `$ be a solution to Yang-Mills equations, $`𝒟^\mu F_{\mu \nu }^{(B)}=0`$, where $`F_{\mu \nu }^{(B)}`$ is the strength tensor of $`B_\mu `$ (we work in the adjoint representation). In the linear order, perturbations $`A_\mu `$ of the Yang-Mills field on the background $`B_\mu `$ obey the equation

$$𝒟^\nu 𝒟_\nu A_\mu 𝒟_\mu 𝒟^\nu A_\nu +2iF_{\mu \nu }^{(B)}A^\nu =0,$$
(4.6)
where $`𝒟_\mu =_\mu +iB_\mu `$. It is convenient to work in the background gauge $`𝒟_\mu A^\mu =0`$ where equation transforms to

$$𝒟^\nu 𝒟_\nu A_\mu +2iF_{\mu \nu }^{(B)}A^\nu =0.$$
(4.7)
Fixing this gauge leaves a freedom in the gauge transformations $`\delta A_\mu =𝒟_\mu \lambda `$ where $`\lambda `$ is a solution to

$$𝒟_\mu 𝒟^\mu \lambda =0.$$
(4.8)
Suppose now that $`B_\mu `$ is a stationary field. Then Eqs. (4.6)–(4.8) result in operator polynomials of the form (1.3). Let $`\omega `$ be the physical spectrum of (4.6), $`\omega _{(1)}`$ the spectrum related to vector Eq. (4.7) and $`\omega _{(0)}`$ the spectrum of (4.8). It is clear that $`\omega `$’s represent a subset among eigenvalues $`\omega _{(1)}`$. The pseudo-trace for the physical spectrum can be represented as

$$K(t)=\underset{\omega >0}{}e^{t\omega ^2}=K_{(1)}(t)2K_{(0)}(t)\underset{\omega _{(1)}>0}{}e^{t\omega _{(1)}^2}2\underset{\omega _{(0)}>0}{}e^{t\omega _{(0)}^2}.$$
(4.9)
(For $`SU(N)`$ group Eqs. (4.6)–(4.8) are invariant with respect to the charge conjugation and their spectrum is symmetric with respect to change $`\omega `$ to $`\omega `$.) The double subtraction of the sum over $`\omega _{(0)}`$ in (4.9) has the following explanation. One subtraction eliminates solutions of (4.7) which do not respect the gauge conditions, $`𝒟_\mu A^\mu =\varphi 0`$. It is easy to see that such modes have frequencies $`\omega _{(0)}`$ because $`\varphi `$ obeys scalar Eq. (4.8). Additional subtraction of $`K_{(0)}(t)`$ in (4.9) eliminates pure gauge solutions $`A_\mu =𝒟_\mu \lambda `$. The last term in r.h.s of (4.9) can be interpreted as contribution of ghosts.

There is an analog of formula (4.5) for transverse solutions ($`𝒟_\mu A^\mu =0`$) to (4.6). The relativistic norm for the solutions with energy $`\omega `$ can be written as

$$A_\omega ,A_\omega =\chi ^{}(\omega )(A_\omega ,A_\omega ),$$
(4.10)
where $`(A_\omega ,A_\omega )=𝑑V(A_\omega ^+)^\mu (A_\omega )_\mu `$ and quantity $`\chi ^{}(\omega )`$ is defined for the polynomial corresponding to vector problem (4.7). The residual gauge freedom can be used to impose condition $`A_0=0`$. For such solutions the norm $`(A_\omega ,A_\omega )`$ becomes positive-definite.

Suppose that the canonical energy of perturbations $`A_\mu `$ on the background $`B_\mu `$ is non-negative, $`H_{YM}[A]0`$. Require also the same property for canonical energy of the scalar field described by (4.8), $`H[\lambda ]0`$. Then one can show that operator polynomials associated with vector, (4.7), and scalar, (4.8), problems obey condition (iii) of Sec. 2. The proof is based on relations (4.5), (4.10).

### 4.3 Dirac Field

As is known, the classical energy for the Dirac field is not positive-definite. Thus, the proof of condition (iii) in this case should be different. Consider the Dirac field in Minkowsky space-time in the presence of a static gauge potential $`A_\mu `$. The relativistic norm of a single-particle wave function is $`\psi ,\psi =𝑑V\psi ^+\psi `$. If the field has a mass, $`m0`$, one can prove the relation

$$\psi _\omega ,\psi _\omega =\frac{\chi ^{}(\omega _1)}{2m}(\psi _\omega ,\psi _\omega ),$$
(4.11)
where $`(\psi _\omega ,\psi _\omega )=𝑑V\psi _\omega ^+i\gamma _0\psi _\omega `$ and $`\gamma _0`$ is anti-Hermitian. It follows from (4.11) and positivity of the norm $`\psi ,\psi `$ that $`\chi ^{}(\omega )`$ does not vanish for massive fields. Now, if $`A_\mu =0`$, one has $`\chi ^{}(\omega )=2\omega `$ and condition (iii) is satisfied. The validity of (iii) for $`A_\mu 0`$ can be proved by continuity. Because $`\chi ^{}(\omega )0`$ quantity $`\chi ^{}(\omega )`$ cannot change the sign when $`A_\mu `$ is switched on from zero to some value.

## Acknowledgments

The author is grateful to the Organizing Committee of the 6th Workshop on ”Quantum Field Theory under the Influence of External Conditions” for the financial support.


text:
# Bias in the temperature of helium nanodroplets measured by an embedded rotor

## I Introduction

A hallmark of ro–vibrational spectroscopy of molecules in liquid 4He nanodroplets has been the observation of rotational structure, as in the gas phase Hartmann95 ; Callegari01 . The distribution of intensity of transitions in this structure has allowed the determination of the rotational temperature of these molecules, and it has been assumed that this provides a measure of the temperature of the droplets themselves Hartmann95 ; Grebenev98 . This assumption has been supported by the close similarity of the temperatures determined for a wide range of molecules and by the assumption that on the long time scales between pickup of solutes and spectroscopic interrogation, the entire system will come to equilibrium as the droplets cool by evaporation.

Adriaan Dockter and I have modeled the evaporative cooling of pure and doped helium nanodroplets using an angular momentum conserving statistical rate model Dokter03 . These calculations have predict that the droplets cool to a temperature close to those inferred from experiments and previous evaporative cooling calculations, but with a vastly broader distribution of energy and total angular momentum than for a canonical distribution at the same temperature. We also found, much to our initial surprise, that while the populations of the rotational levels of a solvated rotor follow a thermal distribution, the fitted temperature of the rotor populations is higher than that of the droplet itself. While this would seem to violate a fundamental principle of thermodynamics, as will be demonstrated below, this divergence of the two temperatures is in fact a consequence of the constraints imposed by conservation of angular momentum. The temperature of the helium is defined by the inverse of the derivative of the log of the density of states with respect to energy at fixed total angular momentum. However, when a rotor is excited to a state with rotor angular momentum quantum number $`j`$, the internal degrees of freedom of the droplet can have angular momentum values between $`|Jj|`$ and $`J+j`$, where $`J`$ is the total angular momentum quantum number, which is treated as a conserved quantity. Thus the derivatives of the density of droplet states with respect to $`J`$ contributes to the rotor population.

## II State Distributions

At the temperature of helium nanodroplets, $`0.38`$K, the only intrinsic droplet excitations are surface ripplon modes Brink90 . The density of states and other thermodynamic quantities can be written as a function of a reduced energy, $`\overline{E}=E/E_\mathrm{R}`$ where $`E_\mathrm{R}=3.77k_\text{B}\text{K}/\sqrt{N}`$ and $`N`$ is the number of helium atoms in the droplet Brink90 . We denote the density of states as a function of reduced energy and total angular momentum as $`\rho _{\mathrm{R}J}(\overline{E},J)`$, and the density states summed over angular momentum states as $`\rho _\mathrm{R}(\overline{E})=_J(2J+1)\rho _{\mathrm{R}J}(\overline{E},J)`$. Using methods recently described Lehmann03 , these can be well approximated by the analytical expressions

$`\rho _\mathrm{R}(\overline{E})`$ $`=`$ $`\mathrm{exp}\left(a\overline{E}^{4/7}+b\overline{E}^{1/7}\right)`$ (1)
$`\rho _{\mathrm{R}J}(\overline{E},J)`$ $`=`$ $`\rho _\mathrm{R}(\overline{E})(2J+1)\sqrt{{\displaystyle \frac{\beta (\overline{E})^3}{\pi }}}\mathrm{exp}\left(\beta (\overline{E})(J+1/2)^2\right)`$ (2)
$`\beta (\overline{E})`$ $`=`$ $`c\overline{E}^{8/7}+d\overline{E}^{13/7}`$ (3)
if the values a = 2.5088, b = -4.3180, c = 0.8642 and d = -0.3524 are used Note1 . Note that for fixed $`J`$, the density is the number of distinct states, i.e. does not include the $`(2J+1)`$ spatial quantization degeneracy of each such state. We define the ripplon microcanonical temperature, $`T_\mathrm{R}(\overline{E})`$, and with fixed total energy and angular momentum, $`T_{\mathrm{R}J}(\overline{E},J)`$, by

$`{\displaystyle \frac{1}{T_\mathrm{R}(\overline{E})}}`$ $`=`$ $`{\displaystyle \frac{d\mathrm{ln}\left(\rho _\mathrm{R}(\overline{E})\right)}{d\overline{E}}}={\displaystyle \frac{4}{7}}a\overline{E}^{3/7}+{\displaystyle \frac{1}{7}}b\overline{E}^{6/7}`$ (4)
$`{\displaystyle \frac{1}{T_{\mathrm{R}J}(\overline{E},J)}}`$ $`=`$ $`{\displaystyle \frac{d\mathrm{ln}\left(\rho _{\mathrm{R}J}(\overline{E})\right)}{d\overline{E}}}={\displaystyle \frac{1}{T_\mathrm{R}(\overline{E})}}+\left({\displaystyle \frac{3}{2\beta (\overline{E})}}(J+1/2)^2\right)\left({\displaystyle \frac{d\beta }{d\overline{E}}}\right)`$ (5)
where in both cases, temperature is measured in units of reduced ripplon energy, $`E_\mathrm{R}`$, divided by Boltzmann’s constant.

We now consider a droplet that has a solvated rigid linear rotor, with effective rotational constant $`B`$, given again in units of $`E_\mathrm{R}`$. Let $`j`$ be the rotational quantum number of the rotor, $`J_\mathrm{R}`$ the rotational quantum number of the ripplons, and $`J`$ the total rotational quantum number. The total density of states for fixed values of $`\overline{E}`$, $`J`$ is given by

$$\rho (\overline{E},J)=\underset{j}{}\underset{J_\mathrm{R}=|Jj|}{\overset{J+j}{}}\rho _{\mathrm{RJ}}(\overline{E}Bj(j+1),J_\mathrm{R})$$
(6)
$`P(j)`$, the fraction of states with rotor quantum number $`j`$, is given by

$$P(j)=\frac{_{J_\mathrm{R}=|Jj|}^{J+j}\rho _{\mathrm{RJ}}(\overline{E}Bj(j+1),J_\mathrm{R})}{\rho (\overline{E},J)}$$
(7)
The temperature, $`T(\overline{E},J)`$, of the coupled system, at fixed $`\overline{E}`$ and $`J`$ is given by

$$\frac{1}{T(\overline{E},J)}=\frac{d\mathrm{ln}\left(\rho (\overline{E},J)\right)}{d\overline{E}}=\frac{_j_{J_\mathrm{R}=|Jj|}^{J+j}\rho _{\mathrm{RJ}}(\overline{E}Bj(j+1),J_\mathrm{R})\frac{1}{T_{\mathrm{R}J}(\overline{E})}}{\rho (\overline{E},J)}$$
(8)
A rotor temperature, $`T_j(\overline{E},J)`$ can be defined from a Boltzmann fit to the rotor level populations $`P(j)`$. This latter corresponds to the “droplet temperature” that has been reported in numerous experiments. This assignment of the rotor temperature to that of entire system, i.e. that $`T_j(\overline{E},J)=T(\overline{E},J)`$ is natural assumption for a system in microcanonical equilibrium, i.e. by equating $`P(j)`$ to the ratio of density of states. If we ignore angular momentum conservation and sum over all values of $`J`$, then the above identification would be exact as long as the mean energy in the rotor, $`T_r`$, if much less than the total energy in the droplet, $`\overline{E}`$, so that the “heat bath” of ripplon states does not change temperature significantly over the range of significantly populated rotor states. Equivalently, we expect $`T_j(\overline{E},J)=T(\overline{E},J)`$ if that $`\rho _\mathrm{R}(\overline{E})`$ can be approximated by $`\rho _\mathrm{R}(\overline{E}Bj(j+1))=\rho _\mathrm{R}(\overline{E})\mathrm{exp}\left(Bj(j+1)/T_\mathrm{R}(\overline{E})\right)`$. Note that if the finite heat capacity of the droplet heat bath is considered, then we would expect the population in higher rotational levels to fall off faster than predicted by the Boltzmann distribution, and thus the effective temperature determined by a fit to the rotor populations to be lower than the droplet temperature.

When we look at the $`J`$ conserving ensemble, we must consider the fact that the average value of $`J_\mathrm{R}(J_\mathrm{R}+1)`$ equals to $`J(J+1)j(j+1)`$ for fixed $`J,j`$. Thus the population of different rotor levels will be influenced not only by $`T_\mathrm{R}`$, but also by the $`J`$ dependent factor in $`\rho _{\mathrm{R}J}`$, which has the same form as the rotational distribution function of a spherical top in a canonical ensemble. In order to test the size of the expected bias in the estimate of $`T(\overline{E})`$ by $`T_j`$, we will examine numerical results for droplet of $`N=\mathrm{3\hspace{0.17em}000}`$ and $`\mathrm{10\hspace{0.17em}000}`$ helium atoms, for which $`E_\mathrm{R}=0.069`$ and $`0.0377`$ K respectively. Based upon our evaporative cooling calculations Dokter03 , we will assume that the droplets cool to a final temperature of 0.35 K, but have a wide range of final angular momentum values, $`J`$. This condition allows calculation of an isothermal curve $`\overline{E}(J)`$ for each droplet size. We further assume that $`B=1\text{GHz}`$, which equals 0.70 (1.27) in reduced units for droplets of 3 000 (10 000) helium atoms. This is the effective $`B`$ constant for SF6 dissolved in helium, from which the temperature of helium nanodroplets was first measured Hartmann95 . It is found that $`T(\overline{E}(J),J)=T_{\mathrm{R}J}(\overline{E}(J),J)`$ to better than 1 % for all values of $`J`$ over the range from 0–4000. It is found that the values of $`P(j)`$ fit a Boltzmann distribution to high accuracy, so that values of $`T_r`$ are well defined. However, over this range of $`J`$, we find a small be systematic difference in the values of between $`T_j(\overline{E}(J),J)`$ and $`T(\overline{E}(J),J)`$, as demonstrated in figure 1. The ratio of $`T_j(\overline{E}(J),J)`$ to $`T(\overline{E}(J),J)`$ is found to be approximately proportional to $`B^1`$, i.e. is linearly proportional to the effective rotational moment of inertia of the rotor in helium (including the contribution due to motion of the helium). As is evident from the figure, the size of the bias decreases for larger droplets, though less than linearly with number of helium atoms in the droplet.

The present results demonstrate that for helium droplets with significant trapped angular momentum, the population of rotational levels of rotors will be biased from that predicted by a Boltzmann distribution at the temperature of either the droplet excitation modes or the temperature of the entire system. The size of the bias is on the order that should be detectable by careful measurement of the droplet rotational constant as the temperature of the pickup gas (and thus the initial collisional angular momentum) is increased. Previous experiments have reported a range of rotational temperatures, but up to know these do not appear to have been considered physically significant. For example, in a study of the cluster size dependence of the spectrum of SF6, the temperature was found to be constant to within statistical noise ($`0.360.40`$ K) for droplets of more than a few thousand He atoms, but to systematically increase for smaller droplets, down to a value of 0.48(20) K for droplets of a few hundred He atoms Hartmann99 . The authors concluded that the temperature of smaller droplets was higher, but the present result suggest that another possibility is that the effect they observed could be, at least in part, due to the bias effect predicted by the present work. Nauta and Miller found that for the spectrum of HF dimer, the rotational temperature was 0.34(1) K Nauta00b , below the excepted droplet temperature of 0.38(2) K. It should be pointed out that the rotational constant of HF dimer in helium is three times larger than that of SF6, and that the pickup of two HF molecules should deposit considerably less angular momentum than the pickup of SF6. It should be noted that the size of the bias is a function only of the total rotational quantum number, and thus the effect should be essentially unchanged for a symmetric top or spherical top molecule. For an asymmetric top, this effect should lead to a failure of the transition intensities to follow that expected for a gas in thermal equilibrium if this bias effect is ignored. Stimulated by the present analysis, Roger Miller and collaborators have made an attempt to look in a systematic way at the temperature extracted from fits to rotor populations for different molecules. Unfortunately, the asymmetries of the lineshapes often observed for the lowest rovibrational transitions lead to model dependent biases in the fitted rotational temperatures that could not be disentangled from the expected effects RMiller\_pc . This give yet one more reason for the need to address this problem of lineshapes of ro-vibrational transitions in helium nanodroples beyond the one previous attempt Lehmann99b . That previous work assumed a canonical ensemble of internal states which the recent evaporative cooling work has now shown to be a poor assumption.

This work was supported by a grant from the National Science Foundation. The author would like to acknowledge the discussions with Adriaan Dokter, Roman Schmied, and Giacinto Scoles.