text:
# Introduction to Statistical Issues in Particle Physics

## I Particle Physics

### I.1 The Subject

Particle Physics emerged as a discipline in its own right half a century ago. It pioneered ‘big science’; experiments are performed at accelerators of increasing energy and complexity by collaborations of many physicists from many institutes. It has evolved a research methodology within which statistics is of great importance, although it has done so without strong links to the statistics community – a fault that this conference exists to remedy. Thus although a statistician will be familiar with the research methods and statistical issues arising in, say, agricultural field trials or clinical testing, they may be interested in a brief description of how particle physicists do research, and the statistical issues that arise.

Particle physics is also known as High Energy Physics 111The terms are almost equivalent; strictly the phrase ‘High Energy’ means ‘above the threshold for pion production’, i.e. the energy at which a collision between two protons can produce three outgoing particles. and the names are sometimes merged to give High Energy Particle Physics. Whatever it is called, its field of study is all the ‘Elementary’ Particles that have been discovered:

* The 6 quarks ($`u,d,s,c,b,t`$)
* The 6 leptons (e, $`\mu `$, $`\tau `$, $`\nu _e`$, $`\nu _\mu `$, $`\nu _\tau `$)
* The intermediate bosons: $`W`$, $`Z`$, $`\gamma `$, g
* The 100+ hadrons made from two quarks( $`\pi `$, $`K`$, $`D_s(2317)`$… ) or three quarks ($`p`$,$`n`$, $`\mathrm{\Lambda }`$…) or five quarks ($`\mathrm{\Theta }^+`$…)

To this long list must be added the corresponding list of antiparticles. However this is not all: the domain of particle physics also includes all the particles that have not yet been discovered – some of which never will be discovered:

* Higgs boson(s)
* squarks and sleptons
* Winos and Zinos/charginos and neutralinos
* further hadrons
* etcetera etcetera…

This list of proposed particles is limited only by the imagination of the theorists who propose them – which is no limitation at all.

For each species of particle we want to establish:

* Does it exist?
* If it does exist, what are its properties: its mass, its lifetime, its charge, magnetic moment and so on?
* If its lifetime is not infinite, what particles does it decay into? What are the branching fractions to different decay modes? What are the distributions in the parameters (energies and directions) of the particles in the final state? Do they agree with our theoretical models?
* What happens when it collides with another particle? What processes can occur and with what probabilities (expressed as cross sections)? What are the distributions of the parameters of the particles produced? Answers will depend on the target particle and the collision energy.

### I.2 Template for an Experiment

To study some phenomenon $`X`$, which could be any of the above, a particle physics experiment goes through the following stages:

* Arrange for instances of $`X`$

This may involve a beam of particles, directly from an accelerator or through some secondary system, striking a target; the beam and target particles and the energy are chosen as being favourable for $`X`$. It may entail a colliding beam machine like LEP for the $`Z`$ or BABAR for CP violation in the $`B`$ system or the LHC for the Higgs. It may be done by producing particles and then letting them decay, as in the studies of CP violation in the $`K^0`$ system. An extreme example is proton lifetime studies, where one just assembles a large number of ordinary protons (perhaps as hydrogen in water) in suitably low-background conditions deep underground and waits to observe any decays.

For important studies dedicated experiments (even accelerators) are built. For many more, the experimenter utilises data taken with an experiment designed primarily for another purpose but also favourable for $`X`$. An example is the study of charm mesons at BABAR, Belle and CLEO, for which the primary purpose is B physics.
* Record events that might be $`X`$

A detector is built (or an existing detector is utilised). ‘Events’ – interactions or decays – are observed by a whole range of detectors (tracking detectors like drift chambers and silicon detectors, calorimeters that measure deposited energy). Fast logic and/or online computers select the events that look promising, and these are recorded: the phrase ‘written to tape’ is used even though today the recording medium is generally disk storage.
* Reconstruct the measurable quantities of the visible particles.

The electronic signals are combined and interpreted: points are joined to form tracks, and measurement of their curvature in a magnetic field gives the particle momentum. A calorimeter may give the energy, a Cherenkov counter the velocity. From this emerges a reconstructed ‘event’ as a list of the particles produced, their kinematic quantities (energies and directions) and possibly their identity (as pions or kaons or electrons, etc.)
* Select events that could be $`X`$ by applying cuts

Knowing the pattern one is looking for, one can then select the events that contain the phenomenon being studied.

A key point is that this selection (and also the electronic selection described above) is not going to be perfect. There will always be a selection efficiency which is less than 100%.

There is also a chance that some of the events that look like $`X`$ and survive the selection and the cuts are actually from some other process. There will be a background which is greater than zero. Statistical techniques are obviously important for the treatment and understanding of efficiency and background.
* Histogram distributions of interesting variables

Relevant quantities, sensitive to $`X`$, are formed from the kinematic variables of the particles detected and measured. These are typically displayed in a histogram, or histograms. (Joint two-dimensional plots are also common. Sometimes, but rarely, the data at this stage is a single number.)

These distributions are then compared with the theoretical predictions, of which there may be several. One will be the predicted distribution if $`X`$ is not present. Another may be the prediction if $`X`$ is present in the amount, and with the properties, predicted by an expected theory such as the ‘Standard Model’ SM  of Particle Physics. There may also predictions obtained within the framework of a particular model, but with one or more parameters adjusted to fit the data.

An example of such a result is shown in Figure 1 (taken from babar1 ). In the top plot the phenomenon $`X`$ is the decay of the $`B^0`$ meson to $`D\pi `$, in the lower plot the decay to $`D^{}\pi `$. The distributions show the invariant mass, which is the quantity given by

$$M^2c^4=\left(\underset{i}{}E_i\right)^2\left(\underset{i}{}\stackrel{}{p}_ic\right)^2$$
(1)
where the sums run over the two final-state particles. If the two observed particles do indeed come from the decay of a $`B^0`$ particle then this quantity should be $`5.28\mathrm{GeV}/c^2`$, though this is smeared out by experimental resolution. The plots show the predictions of a theory in which this decay does not occur (and all events are background) and also a prediction in which the decay is produced, with a normalisation adjusted to give the best fit to the data. The result of this fit gives the number of signal events, from which the branching ratio can be obtained (though in fact that was not done in this example).

If that looks trivial, a harder example is the decay $`B\pi ^0\pi ^0`$, taken from babar2  and shown in Figure 2. (To be fair, things are not quite as bad as this 1 dimensional plot implies.)

In this confrontation of theory with experiment, one can then ask: is there any evidence for $`X`$ or is the null hypothesis unrefuted? Given that there is $`X`$, what is the best estimate for the normalisation (and perhaps other) parameter(s) involved in our model for $`X`$? Are these results compatible with the standard prediction for $`X`$? These are familiar statistical questions.

### I.3  Statistics in HEP

From the above description we can bring out some features of the way statistics is used in HEP.

Firstly, everything is a counting experiment. To measure a branching ratio or a cross section, one counts the number of events produced and observed. To measure the mass of a particle one uses a histogram where the number of entries in each bin is a random Poisson process. (The data of Figure 1 could be used to fit the mass of the $`B^0`$ meson, were it not already well known.) Poisson statistics is of paramount importance. Even the Gaussian (Normal) distribution plays its main rôle as the large $`N`$ limit of the Poisson. (There are exceptions to this generalisation, but they occur in the details of the reconstruction of particle quantities.)

This unpredictability is not due to any lack of knowledge on our part: not sampling error, or measurement error, or due to unconsidered effects. It is true and absolute randomness, driven by the fundamental nature of quantum mechanics. We know that, for instance, a $`K_s^0`$ particle may decay into two charged pions or two neutral pions, with probabilities of 69% and 31% respectively. That is all we can ever know. A sample of $`K^0`$ particles will decay to $`\pi ^+\pi ^{}`$ and $`\pi ^0\pi ^0`$ in a ratio of roughly 2:1 even if they are prepared absolutely identically – we have no hope of ever being able to say which ones are more likely to ‘choose’ one path rather than another. Likewise the timing of a decay is absolutely random in that the probability that a particle existing at time $`t`$ will decay before time $`t+\delta t`$ is a constant, independent of the value of $`t`$; there is no ‘ageing’ process.

But the Poisson distributions that result are just like any conventional Poisson process. These and other uncertainties, are (almost always) controlled and understood. These distributions have standard deviations known to be $`\sqrt{N}`$. The Gaussian used for the signal distributions in Figure 1 is well established (it has a mean of 5.28 and a standard deviation of $`0.0025\mathrm{GeV}/c^2`$).

So, in common with the other physical sciences, the distributions involved (signal, backgrounds) are given by functions known up to a few parameters – which can be fitted for. The approach to the data is not descriptive (identifying features, looking for trends) but prescriptive: the distribution is taken as having some functional form, and one has a pretty good idea as to what that functional form is, apart (possibly) from a few adjustable parameters.

### I.4 Unused Statistical Methods

A consequence of this knowledge of uncertainty – the fact that we know what it is that we don’t know – is that many techniques commonly used in the broad field of statistics are little used (or not used at all) in particle physics.

Student’s $`t`$ is unknown. This is a technique used to handle small numbers of values from a distribution of unknown mean and unknown standard deviation, but our uncertainties come from known measurement errors. (If a measurement error is not known, a separate large-number determination is made.) The $`F`$ test and ANOVA, tools for studying problems with unknown variances, are similarly of little use. The whole experimental design field – Latin squares and similar techniques used to minimise uncontrollable effects – is not needed as such effects are not a problem.

Another set of neglected techniques are those handling Time Series and Markov chains. Changes with time can be relevant in some studies, but it appears in them as another quantity to be measured and histogrammed. The development with time of a particle is basically smooth, punctuated by radical transformations (such as the decay of a particle to two or more lighter ones) which occur at random times.

Non-parametric Statistics are also barely featured, as all these distributions, which are believed to be true idealisations of what ‘really’ happens, or at least good approximations to them, are parametrised.

The notion of a Parent population is not helpful: a sample of particles is taken, but the randomness is (as stated earlier) inherent in the nature of particle behaviour and not produced by the sampling. If there is a parent distribution, it is an infinite set of particles produced under these conditions – all the events we might have seen.

The point here is not that particle physics has nothing to learn from standard statistical techniques. The Statistician has many implements in their toolbox. Different fields of application will call for different tools; some of those heavily used in other fields are of less relevance in this one.

## II Tools

Having seen that particle physics makes little use of some statistical tools, we take a more detailed look at the ones it does utilise.

### II.1 Monte Carlo Simulation

Theoretical distributions for the quantitites being studied are predicted by quantum mechanics – perhaps with a few unknown parameters – and are often beautiful and simple. Angular distributions may be flat, or described by a few trigonometric terms; masses often follow a Cauchy function (which the particle physicists call the Breit-Wigner), time distributions may be exponential, or exponential with a sinusoidal oscillation.

These beautiful and simple forms are generally modified by unbeautiful and complicated effects (higher-order calculations in perturbation theory, or the fragmentation of quarks into other particles). Furthermore the measurement and reconstruction process that the detector does for the particles is not completely accurate or completely efficient.

The translation from knowing the distributions in principle to knowing them in practice is done by Monte Carlo simulation. Particles are generated according to the original simple distributions, and then put through repeated random processes to describe the theoretical complications and then the passage of particles through the detector, including probabilities for colliding with nuclei in the beam pipe, slipping through cracks in the acceptance, or other eventualities. A complete software representation of all the experimental hardware has to be coded. The effects of the particles on the detector elements is simulated and the information used to reconstuct the kinematic quantities using the same programs that are run on the real data. This provides the full theoretical distribution function that the data is predicted to follow, albeit as a histogram rather than a smooth curve.

These programs are large and slow to run. Significant resources (both people and machines) are put into them. The generation of ‘Monte Carlo data’ is a significant issue for all experiments. Cases are known where data has been taken and analysed but results delayed because of lack of the correct Monte Carlo data MC .

### II.2 The Likelihood

Having the parametrised theoretical description of the distribution means the likelihood function is always known, and it assumes an overwhelmingly important position. Writing this function – where the $`x_i`$ are the data and $`\theta `$ the unknown parameter(s)

$$L(x_1,x_2\mathrm{}x_N|\theta )=P(x_i|\theta )\mathrm{ln}L=\underset{i}{}\mathrm{ln}P(x_i|\theta )$$

the form $`p(x|\theta )`$ is totally known, and $`L`$ (or $`\mathrm{ln}L`$) follows.

Having the likelihood function, the Maximum Likelihood estimator is then easy to implement, and is very widely used. Even estimators like least-squares are, at least by some, ‘justified’ as being derivable from Maximum Likelihood. Its (asymptotic) efficiency, and its invariance properties are desirable and useful.

In some cases the ML estimate leads to an algebraic solution but in general, and in complex analyses, the physicist just maps out $`\mathrm{ln}L`$ for their dataset as a function of $`\theta `$ and reads off the ML estimator from the peak, as can be done in Figure 3. This also produces an interval estimate as part of the minimisation process. Following the value of $`\mathrm{ln}L`$ until it falls off by $`\frac{1}{2}`$ from its maximum gives the 68% central confidence interval. Strictly speaking this is valid only for large $`N`$, but this restriction is generally disregarded. Perhaps we should not be so cavalier about doing so.

Maximum Likelihood methods can also be used for functions with several parameters, as illustrated in Figure 4. Confidence regions are mapped out by reading off the likelihood contours. This is done in many analyses and the MINUIT program MINUIT  is widely used in exploring the likelihood and parameter space.

### II.3 Fitting Data

Fitting the parametrised curve to the experimental data is done by several techniques.

1) $`\chi ^2`$ using $`\sigma ^2=n`$ i.e. minimising $`\chi ^2=_i\frac{(y_if(x_i|\theta ))^2}{n}`$ has the advantage that the minimisation can be done by differentiating and solving the normal equation, which is especially simple if $`f`$ is linear in $`\theta `$. However the use of the observed number rather than the predicted number in the denominator is recognised to lead to bias (downward fluctuations get an undue weight) and this cannot safely be used if $`n`$ is small. (Actually in many cases what happens is that one of the bins has $`n=0`$, and the physicist gets divide-by-zero messages and then starts to worry.)

2) $`\chi ^2`$ using $`\sigma ^2=f`$ i.e. the predicted value rather than the actual number, avoids the bias (and the divide-by-zero problem) but gives nonlinear equations. It still suffers from using a Gaussian probability as an approximation to a Poisson distribution and is thus not the ‘real’ maximum likelihood estimator.

3) ‘Binned Maximum Likelihood’ uses the Poisson likelihood in each bin rather than the $`\chi ^2`$. It is therefore a proper Maximum Likelihood estimator. Efficiency is lost (only) if the bins are wider than the structure of the data.

4) Full maximum likelihood does not use binning at all. It can be useful for very small event samples. For large samples it becomes computationally intensive (as there is a sum over events rather than a sum over bins) though with today’s computers this is hardly important. Perhaps a more significant factor for physicists is that it does not have the readily interpretable graphic image given by a histogram and fitted curve.

### II.4 Goodness of Fit

Having found a fit, one has to judge whether to believe it. Whether the question is ‘Does the curve really describe the data?’ or ‘Do the data really fit the curve’ depends on one’s point of view.

The likelihood value does not contain the answer to this question. This appears counter-intuitive and many people have wrestled (unsuccessfully) to produce ways that the likelihood can be used to say something about the quality of the fit.

The $`\chi ^2=_i\left(\frac{y_if(x_i|\theta )}{\sigma _i}\right)^2`$ certainly does give a goodness of fit number. It is heavily used for GoF and 2-sample tests: researchers may quote $`\chi ^2`$ or $`\chi ^2/N_D`$ or the probability of exceeding this $`\chi ^2`$.

Alternative measures of goodness of fit have never really caught on. The Kolmogorov-Smirnov test is occasionally used – generally misleadingly, in my opinion. This is a totally robust test but pays the price for that by being weak. If you know anything about the data, e.g. that the numerical value of the parameter means something, then a more powerful test should be available. The KS test is being used to certify that distributions are in agreement when a more powerful approach would show up a difference.

### II.5 Toy Monte Carlo

The ‘Toy Monte Carlo’ has emerged as a technique made possible by modern computing resources. Having obtained a result, it may be hard or impossible to obtain significance levels or confidence regions in the traditional analytic way, for instance if the likelihood function one is studying does not even plausibly resemble a distorted parabola, but instead some shape with multiple maxima.

As an alternative approach, starting with an estimate $`\widehat{\theta }_{exp}`$ from the data, say $`\{x_1\mathrm{}x_N\}`$, how can one establish a confidence region? Consider any particular $`\theta `$. Use the known $`L(x|\theta )`$ to generate a set of $`N`$ values of $`x`$ – an “experiment’. Use this in your estimator (whatever that is) to find a corresponding $`\widehat{\theta }`$. Repeating many times gives the probability that this $`\theta `$ will give an estimate below (or above) the experimental one. This is just what the Neyman construction uses. To find a particular confidence region one has to explore the parameter space until one finds the limits one wants.

## III Topics

Having explained the basic and generally agreed techniques used, there are a number of topics where advances are being made, or which are the subject of heated dicussion and argument, or both.

### III.1 Bayesian Probability

The religious war which has been waged over the past few years has now cooled – although some isolated zealots remain on both sides. The ‘frequentists’ have come to accept that the use of Bayesian techniques can be illuminating and helpful, and sometimes provide more useful information than a frequentist confidence level, especially for measurements of bounded parameters (e.g. masses). The ‘Bayesians’ are recognising that Bayesian confidence levels will not totally replace the use of frequentist levels, and that they do have to take on board the issue of robustness (or otherwise) under changes of prior.

A real benefit of this debate has been to bring the subject out into the open. The classic statistics texts Orear ; Frodesen , from which many particle physicists first learned the subject, slide swiftly between the frequentist and Bayesian concepts of probability, never really acknowledging that they are using two very different quantitites.

### III.2 Small Signals and Confidence Regions

The ‘Energy Frontier’ is a cutting edge of particle physics: new, more powerful, accelerators open up new areas for investigation and new particles are discovered. Another cutting edge is the ‘Luminosity Frontier’: the discovery that processes hitherto thought to be impossible do actually occur, albeit very rarely. The discovery of CP violation CP : that the probability of the decay $`K_L^0\pi ^+\pi ^{}`$ was not zero but $`0.2\%`$, was enormously important despite the smallness of the figure. Many of today’s experiments are looking for phenomena which are known to be exceedingly rare, at the parts-per-million level at best.

Although the implications can be spelt out quite simply and dramatically – ‘If the AMS experiment sees even one $`\overline{{}_{}{}^{12}C}`$ nucleus, our entire view of the universe will change.’– in practice things are not so clear-cut because of the presence of background. Also one has to be able to handle not just the dramatic discoveries, but the much more frequent useful analyses that make no discovery but push back the limits and the region in which any discovery may be made.

An experiment that sees no events will note the standard result from Poisson statistics that an observed number of zero translates to a limit on the true value of less than 3 events, with 95% confidence. This can then be converted (using the figures for this particular experiment) into a limit on the branching ratio or cross section for the process concerned, and then possibly into a limit on a mass or coupling constant. If there is an expected background for this process equivalent to, say, 0.2 events, then the amount for the branching ratio limit is reduced to 2.8. But this clearly has problems: suppose the predicted background were 3.1 and no events were observed (unlikely but not impossible), what can one then say about the limit?

There has been a lot of activity and discussion recently in this area. Indeed it sparked off the workshop Limitsconf  of which this conference is the successor. The standard frequentist (Neyman) construction may result in statements about results in the non-physical region (here, a negative number of signal events) which, though statistically correct, appear nonsensical. Bayesian methods avoid this problem, as does the frequentist technique proposed by Feldman Cousins FeldmanCousins  which switches smoothly and automatically between quoting central and one-sided confidence regions.

### III.3 When to Claim a Discovery?

Another area of discussion is over the form of reporting non-zero signals. When the number of signal events is much larger than the expected background, or a fitted parameter is significantly different from the theoretical prediction, then clearly the experiment can claim a discovery. If the numbers or parameter values are compatible, the experiment quotes an upper limit. But there is an area in between where the probability of the null hypothesis giving the result is small enough to be interesting, but not so small as to be completely negligible. The experiment must not be rash, phoning the New York Times with a discovery which turns out to be a statistical fluctuation, nor must it be too cautious or the subject can never progress. Such results are bound to occur – the probability that an experiment will produce a value in this region is by definition small-but-not-negligible, or better. Given the large number of busy experiments reporting results, this is a real problem.

Some experiments have policies such as 4$`\sigma `$ for ‘evidence for’, 5$`\sigma `$ for ‘discovery of’ – significance levels are often presented in terms of the equivalent discrepancy in standard deviations. Is it possible to report a two-sided result (as the Feldman Cousins technique will sometimes produce) and yet not claim a discovery? ‘We report with 95% confidence that the branching ratio lies in the range $`(2.3\mathrm{to}3.4)\mathrm{\hspace{0.17em}10}^6`$ but we’re not actually claiming to have seen it.’ Such ‘discoveries’ are reported in a way which must be affected by the prior (subjective) probability, in exactly the way the Bayesians describe. Statistically identical data on the decays $`B^+\pi ^+\pi ^0`$ and $`B^+\pi ^{}\pi ^0`$ would be reported completely differently.

### III.4 Blind Analysis

In recent years particle physicists have become aware of practitioner bias. This has been fuelled particularly by reports from the Particle Data Group, which has the job of reporting and combining the measurements of particle properties PDG , who show how some values change significantly over time, but never by more than one standard deviation. Another source of disquiet was the Electroweak measurements from LEP and the SLC which agree with each other and with the Standard Model far too well with a $`\chi ^2`$ per degree of freedom well below 1 LEP .

This practitioner bias is against claiming differences from the null hypohesis. The experiment template presented in section I.2 often continues

* Extract result, usually by fitting parametrised distribution(s) to data.
* Compare your result with that of accepted theory and/or other experiments.
* If it disagrees, look for a bug in your analysis. You will probably find one. Keep searching and fixing until the agreement is acceptable.

The mistake in method is that the experimenter stops looking for bugs when they have agreement, not when they honestly believe that all (substantial) biases are accounted for. To guard against this the data can be ‘blinded’. There are two techniques used, covering two types of situation

* In the extraction of a result, this can be encoded by some unknown offset.
* Choosing the cuts which select the data is done on Monte Carlo data, or on real data in sidebands – regions close to but not actually including the region where the signal is expected. Otherwise the temptation to nudge a cut slightly to include a few more events is too great.

### III.5  Systematic Errors

In the early days of particle physics, the 50s and 60s, a typical experiment would get handfuls of events – a few hundred if lucky – from painstaking analysis of bubble chamber pictures. Statistical errors were thus $`10\%`$ and were so large that the effect of systematic uncertainites was generally small.

In the 70s and 80s, the development of counter experiments led to event samples in the tens of thousands. Statistical errors were now at the per cent level, and systematic errors began to be more important.

The current generation of experiments – the $`Z`$ factory at LEP, the $`B`$ factories, Deep Inelastic Scattering at HERA – deal with millions of events. Statistical errors are at the level of $`0.1\%`$ and we have learned how to talk about ‘parts per mille’.

Systematic errors (uncertainties in factors systematically applied in the analysis) can no longer be fudged. The word ‘conservative’ has been grossly overused in this context. It sounds safe and reassuring; in practice it is usually a sign of laziness or cowardice. The experiment perhaps cannot be bothered to evalute an uncertainty and makes a guess, and then it inflates that guess to cover the possibility that they’ll be caught out, and calls it a ‘conservative’ estimate of the systematic error.

Particle physicists also confuse the evaluation of systematic errors with overall consistency checks. There is bad practice being spread to and between graduate students. They will identify all the calibration constants and parameters that contribute to the final result and vary those by their appropriate error, and fold the resultant variation into the systematic error. This is correct procedure. But they will also vary quantities like cut values, which should not in principle affect the result, by some arbitary amount and then solemnly fold those resulting variations into the systematic error. This is nonsense. Looking at what happens when you change a cut value is a good and sensible thing: a (say) looser cut will give a higher efficiency and a higher background and thus more observed events, but after correcting for the new efficiency and background the result should be compatible with the original. This is a useful check that one understands what’s going on and that the analysis is consistent. But it does not feed into a numerical uncertainty.

### III.6 Unfolding

Measurements of the properties of particles in events are made with finite resolution, so the plots of these quantities, and functions of these quantities, are ‘smeared out’. Events move between histogram bins. Sharp peaks become broad, edges become slopes.

The recovery of the original sharp distribution from the observed one is known as ‘unfolding’. This is an alternative use of the Monte Carlo simulation process: rather than compare the data with a theoretical prediction smeared by Monte Carlo simulation, one compares the original theory with the de-smeared data. Clearly this is preferable, if it can be done, as the unfolding process depends only on the experiment and not on the original theory, and so once unfolded the data can be compared with any prediction.

It looks a simple problem: given an original distribution as a histogram, the probability of migration from any bin $`i`$ to any bin $`j`$, $`P_{ji}`$, can be estimated from a Monte Carlo sample (this includes the probability that it may not be accepted: $`_jP_{ji}1`$). The matrix is inverted, and then applied to the data histogram to give the reconstructed original.

Unfortunately it is not at all simple Blobel . In the matrix inversion the errors on the $`P_{ji}`$ from finite statistics have devastating consequences and produce unrealistic results. There is a lot of activity in handling this in a sensible way, and in investigating other approaches, such as Maximum Entropy techniques.

### III.7 Combining Results

The combination of compatible measurements with different errors is straightforward. However results are sometimes incompatible, or marginally compatible. But something must be done with the results, as the community needs a way of using the combined number. Indeed it is the responsibility of the Particle Data Group PDG  to combine measurements and form ‘world average’ results in a meaningful way.

There is also a problem in combining limits. If two experiments report 95% confidence level upper limits of, say, 0.012 and 0.013, how can one combine these two measurements? This question was put forcefully by the Higgs searches at the end of the LEP run. The four experiments reported results separately compatible and possibly marginally suggestive of a signal from a Higgs boson of mass around $`114\mathrm{GeV}/c^2`$. Did four possibles make a probable? The answer to that statistics question determined whether or not LEP would run another year, at a cost of millions not only in power bills but in its impact on the construction schedule for the LHC. The CERN management decided that the answer in this case was ‘no’. History will be their judge.

In combining experiments the likelihood function contains much more information than a simple limit, or value and error. There is a suggestion that these should be routinely published, and we are probably going to see that happening a lot in the future.

### III.8 Multivariate Classification

The classification of events (usually ‘signal’ and ‘background’) and particles (pion, kaon…) by means of a cut on a discriminator variable is a basic hypothesis testing problem. However there may be several variables, each containing useful information, and the best choice will be made by combining these in some way.

The Fisher Discriminant has been re-discovered as a technique which is good if the means of distributions differ between the two samples. The Neural Network (feed-forward ‘perceptron’ configuration) has become a standard item in the toolbox which can handle more general differences, and there are many developments going on in this area.

The use of cuts is deeply engrained. In many cases it is simple and appropriate. However in cases where there are no clean boundaries it may be better to consider all events, weighting them according to their signal-like or background-like nature.

## IV Conclusions

I have given several talks on ‘Statistics for Particle Physicists’ but ‘Particle Physics for Statisticians’ has been a new and interesting experience. This has been a very broad view. Particular topics will be considered in detail in the subsequent talks in this conference, in both plenary and parallel sessions. Hopefully the account here will provide you with a map which will help you place them in context.

###### Acknowledgements.

The author gratefully acknowledges the support of the Fulbright Foundation.


format: text

text:
# Symmetry breaking and structural distortions in charged XH4( X = C, Si, Ge, Sn and Pb ) molecules

## I  Introduction

Although neutral XH4( X= C, Si, Ge, Sn and Pb ) molecules have tetrahedral symmetry, charging these molecules breaks this symmetry. The symmetry breaking in general is associated with a large structural distortion and in many cases unusual H-H bonding arrangement. Some of these clusters have been investigated earlier and the symmetry breaking has been attributed to the Jahn-Teller effect. In the Jahn-Teller effectBersukerCM01v101 ; solidv20 ; Englman , electrons occupying a degenerate level favor a low symmetric structure by lifting the degeneracy. For example, the highest occupied molecular orbital(HOMO) of tetrahedral SiH4 is triply degenerate. When an electron is removed from SiH4, its structure gets distorted drastically due to the Jahn-Teller effectFreyJCP88v89 ; ProftCPL96v262 ; KudoCP88v122 ; CaballolCPL86v130 ; PorterJCP01v114 .

The structure of SiH4+ has been investigated by a large number of workersFreyJCP88v89 ; ProftCPL96v262 ; KudoCP88v122 ; CaballolCPL86v130 ; PorterJCP01v114  and it took almost a decade to settle its ground state structure. Now it is generally agreed that the ground state structure of SiH4+ has Cs symmetry. Similarly, the structure of CH4+ has been studied extensivelyPorterJCP01v114 ; WetmoreJCP99v110 ; VagerPRL86v57 ; SignorellJCP99v110 ; VegaJMOL85v123 ; BoydJCP91v94 ; SignorellJES00v108 ; TakeshitaJCP87v86 ; FreyJCP88v88 , though there is still ambiguity regarding its ground state structure. The calculations done by the density functional theory( DFT ) methods show that the ground state structure of CH4+ has D2d symmetryPorterJCP01v114 ; WetmoreJCP99v110 , while the Hartree-Fock( HF ) and post HF calculations show the ground state structure as C2vWetmoreJCP99v110 ; FreyJCP88v88 . The lowering of tetrahedral symmetry of CH4+ and SiH4+ has been attributed to the Jahn-Teller effect. However, the structure of CH4+ is quite different from that of SiH4+. It is difficult to understand the difference between the CH4+ and SiH4+ distortions only from the Jahn-Teller effect. This is because the Jahn-Teller theorem only suggests the existence of distortion and says nothing regarding the nature and magnitude of the distortion.

In this paper we focus on the ground state structures of charged XH4 molecules using the first-principles electronic structure methods and provide an insight into the symmetry breaking mechanism. We find that upon charging most of the molecules undergo large structural distortion which can be attributed to the Jahn-Teller effect. Further, the nature of distortions of positively charged molecules is different from that of negatively charged molecules. Interestingly, most of the positively charged molecules have an unusual H-H bond and appear like transition states of unstable XH2+ and H2 complex. We try to understand the distortions, unusual H-H bonding and stability of these molecules. Vibrational analysis shows that these molecules are not transition states. The calculated fragmentation energy indicates that the stability of XH4+ molecules decreases from CH4+ to PbH4+.

Using various exchange-correlation functionals such as local density approximation(LDA), local spin density approximation(LSDA), gradient corrected approximation(GGA) and spin-polarized GGA(GGA-SP) we have investigated the effect of different occupations in charged tetrahedral XH4 molecules. Our calculations show that the configuration corresponding to integer occupancy and asymmetrical charge density leads to the lowest energy. We provide simple arguments in terms of charge asymmetry and electrostatic interaction between atoms which explain the nature of distortions in charged XH4 molecules. For example, we are able to explain why the nature of distortion of SiH4+ is different from CH4+ and why the positively charged molecules distort differently from the negatively charged molecules. We study the structural distortions of SiH4+ in detail using conjugate gradient method and show that the evolution of the structure is governed by the electrostatic interaction between the atoms. Our investigations reveal that the distribution of charges on H atoms plays an important role in distorting the charged molecules. If the charges on H atoms are not equal, the structure gets distorted due to different electrostatic forces on the atoms. We also discuss how the charge asymmetry is created in these molecules by the directionality and occupation of the HOMO.

The plan of the paper is as follows. In section II we give computational details of the present work. Section III contains results and discussion and finally in section IV we give our conclusions.

## II  Computational details

We have investigated the structure of neutral and charged XH4 using Vienna *ab-initio* simulation package( VASP )KressePRB96v54  and Gaussian98Gaussian98  package. We have cross checked our calculations using both packages and find that the results obtained by both methods are in good agreement. A number of calculations were done using the unrestricted Hartree-Fock( UHF ) and DFT methods using the Gaussian98 package. In the DFT calculations we have used several exchange-correlation functionals such as Slater’s exchange with the local correlation functional of Perdew( SPL ), Becke’s exchange which includes gradient density with gradient corrected correlation functional( BPW91) and Becke’s 3 parameter hybrid functional( B3LYP )Gaussian98 ; Foresman . In these calculations we have employed 6-311g\** basis set for molecules having C and Si atoms and LanL2DZ basis set with effective core potential for other molecules. The plane wave ultrasoft pseudopotential calculationsKressePRB96v54  were done with the local density approximation( LDA )CeperleyPRL80v45 ; PerdewPRB81v23  employing a simple cubic supercell of 20$`\AA `$. The number of plane waves used in the calculation is set by energy cut-off( Ecut). We have used 287 eV Ecut for neutral and charged CH4 molecules and 150 eV Ecut for other molecules. The use of supercell technique for charged systems causes electrostatic interaction within the images of supercellsNeugebauerPRB92v46 ; MakovPRB95v51 . Calculation of charged systems are efficiently handled in VASP by applying a background charge to maintain the charge neutrality and by adding dipole and quadrupole moment correctionsKressePRB96v54 . The structure and ionization potential of SiH4+ predicted by using this procedure are in close agreement with those of earlier calculations based on the HF, MP2 and LSDA methodsKudoCP88v122 ; PorterJCP01v114 .

## III  Results and Discussion

### III.1  Ground state structures

Fig. 1 shows the ground state structures of XH4 molecules under neutral and charged states using the Gaussian98 package with BPW91 functional. We first discuss the ground state structures of the positively charged molecules. We see from the figure that the tetrahedral symmetry is broken when XH4 molecules are ionized. The structure of CH4+ has D2d symmetry, while other positively charged molecules have Cs symmetry. We have investigated various possible symmetries of XH4+ molecules and our results are summarized in Table. I. We see from the table that C2v and D2d structures of CH4+ are very close in energy and differ only by 0.11 eV. The other positively charged molecules also have C2v as the next higher energy structure. We could not find local minima of C3v and Cs symmetries for CH4+ and C2v symmetry for SiH4+ at the BPW91 level. These missing structures of CH4+ and SiH4+ molecules in our calculations have been studied earlier using different exchange-correlation potentials and different methodsFreyJCP88v89 ; ProftCPL96v262 ; KudoCP88v122 ; WetmoreJCP99v110 ; TakeshitaJCP87v86 ; FreyJCP88v88  which show that these are not the ground state structures.

It is interesting to note that in the Cs symmetry, the positively charged molecules have a structure which is severely distorted from the tetrahedral structure. Furthermore, there is an unusual bonding between two hydrogen atoms. This bonding is confirmed by the charge density calculation as shown in Fig. 2 which shows large electron density between the two hydrogen atoms of SiH4+. To ascertain that this is not an artifact of the BPW91 calculation, we have re-optimized the ground state structures using the HF and DFT methods with different exchange correlation functionals such as SPL and B3LYP. The H-H bond lengths obtained from various methods are essentially the same. This kind of over-coordination of hydrogen has also been observed in CH5+MarxSci99v284 ; MarxAng97v36 , SiH5+BooJCP95v103  and small hydrogenated silicon clustersBalaPRB01v64 . In Fig. 3, we have shown the total energy curve of SiH4+, H2 and H2+ as a function of H-H distance. The total energy of SiH4+ is calculated by keeping Si-H bond length fixed but varying H-H bond distance. From the figure we see that the curvature of SiH4+ curve is smaller than H2 and larger than H2+ curves. This shows H-H bond in SiH4+ is stronger than H2+ but weaker than H2. Vibrational analysis also indicates that the H-H bond is sufficiently strong in these molecules. We see from Table. II that the force constant for H-H stretch mode of H-H pair is smaller than H2 and larger than H2+ which implies that the strength of H-H bond is weaker than H2 and stronger than H2+.

For the ground state structures of the negatively charged molecules, we note from Fig. 1 that the distortions of negatively charged molecules are different from those of positively charged molecules. XH4- molecules have C2v symmetry except CH4- which has tetrahedral symmetry. The two hydrogen atoms of XH4- molecules are pushed away from each other which is in contrast to the H-H bond of XH4+ molecules with Cs symmetry.

### III.2  Stability

One may argue that the structure of XH4+ molecules with Cs symmetry appears more like XH2+ and H2 complex and hence unstable. Frey and DavidsonFreyJCP88v89  studied SiH4+ using CI singles root calculation and showed that Cs structure is more stable than C2v against fragmentation into SiH2+ and H2. We have calculated the fragmentation energy which is the energy required to fragment a XH4+ molecule into its binary products and has been calculated by taking the total energy difference between the XH4+ and the possible binary products. Of all the possible channels, the fragmentation of XH4+ into XH2+ and H2 requires least energy. We give the fragmentation energies of XH4+ molecules for this process in Table. III. We see from the table that the energy required for the fragmentation of XH4+ into XH2+ and H2 decreases for X=C to Pb. This means that the stability of XH4+ molecules decreases from C to Pb molecules.

To examine the stability further we have done vibrational analysis of these molecules. In Table. IV we present the vibrational frequencies of charged and neutral XH4 molecules. Our calculated vibrational frequencies 880, 950, 2175 and 2192cm-1 of SiH4 are in good agreement with those of experimentally obtained valuesCoatsJMol94v320  of 914, 953, 2189 and 2267cm-1. We see from the table that PbH4+ has a negative frequency and other molecules have only positive frequencies. Further, we see that one frequency is very low for Cs symmetric XH4+ molecules and decreases from SiH4+ to GeH4+ and become negative for PbH4+.

The positive frequencies of a molecule indicate that the structure corresponds to an energy minimum. If one frequency is negative, the structure may be a transition state of the moleculeForesman . Since the frequency calculation is highly dependent on the method, the small negative frequency of PbH4+ may be an artifact of the method. To clarify this we have done vibrational analysis on re-optimized structures of Cs symmetric XH4+ molecules using various methods. Calculation with the SPL exchange-correlation functional shows that one frequency of PbH4+ which was negative in the calculation with BPW91 functional is positive but near zero. Similarly, calculation with the UHF method also shows that all the frequencies of PbH4+ are positive. Irrespective of the sign of calculated frequencies, all these methods show that one vibrational mode of Cs symmetric XH4+ molecules has a low frequency. Also we find that the force constant corresponding to this mode is nearly zero, which implies that the H atoms in these mode are moving on nearly flat potential energy surfaces.

We see from Fig. 1 that the mirror symmetry plane of Cs structure is formed by the H-H pair and X atom. The vibrational mode with low frequency corresponds to the vibration of H-H pair about the mirror symmetry plane of Cs structure. When the H-H pair vibrates about the plane, one atom of H-H pair moves in and the other moves out of the mirror plane without change in the H-H distance. Since the force constant is nearly zero for this vibration, the motion of H-H pair is more like an internal rotation about an axis which is perpendicular to the H-H bond and passing through the X atom. This internal rotation of H-H pair in XH4+ is similar to the internal rotation of O-H in CH3OH moleculeAyalaJCP98v108 ; QuadeJMS98v188 ; DuanCP02v280 . Thus the low frequency mode in the present case does not indicate a transition state.

### III.3  Effect of different occupations

Before proceeding further, we would like to resolve an interesting question regarding the occupation of degenerate states in charged tetrahedral XH4 clusters. Let us ionize a charged XH4 cluster keeping its tetrahedral structure intact. In this structure the HOMO is triply degenerate and is occupied by only 5 electrons. Hence, the HOMO of XH4+ is not completely filled. In such case of incompletely filled degenerate systems, there are many possible ways to fill the levels. The question is what is the right way of filling the levels and which distribution leads to the lowest energy. The usual practice in such situation is to fill all levels with fractional number of electrons, thereby keeping the symmetry intact.

To resolve the issue of occupation we have performed paramagnetic and spin-polarized calculations on tetrahedral SiH4+ and CH4+ using VASP package with local and non-local functionals for various possible electronic distributions. In paramagnetic calculations we have chosen two possible electronic distributions. The configuration with integer( 2,2,1 ) occupation corresponds to asymmetric charge density and the other with fractional( $`\frac{5}{3}`$, $`\frac{5}{3}`$, $`\frac{5}{3}`$ ) occupation corresponds to symmetric charge density. Note that the spin multiplicity of SiH4+ and CH4+ is unity. The multiplicity is preserved in charge asymmetric configuration whereas it is zero for charge symmetric configuration. In spin-polarized calculations we have chosen three possible electronic configurations. The configuration with $`\frac{5}{3}`$$``$, $`\frac{5}{3}`$$``$, $`\frac{5}{3}`$$``$ ( $`\frac{5}{3}`$$``$, $`\frac{5}{3}`$$``$, $`\frac{5}{3}`$$``$ ) occupation corresponds to charge symmetric density and multiplicity zero. Another configuration with 1$``$, 1$``$, 1$``$( $`\frac{2}{3}`$$``$, $`\frac{2}{3}`$$``$, $`\frac{2}{3}`$$``$) occupation also corresponds to symmetric charge density but multiplicity one. The other configuration with 1$``$, 1$``$, 1$``$( 1$``$, 1$``$, 0$``$) occupation corresponds to asymmetric charge density and multiplicity one. We note that in all these configurations, the Kohn-Sham wave function is a single slater determinant. Hence the energy of the multiplet can be calculated directly using von Barth schemeBarthPRA79v20 . The calculated total energy of tetrahedral SiH4+ and CH4+ molecules are summarized in Table. V.

From the table we see that paramagnetic calculations with local( LDA ) and non-local( GGA ) functionals do not give much difference between total energies of the molecules for symmetric and asymmetric charge density. On the other hand, the spin polarized calculations with local functional( LSD ) favors charge symmetry for SiH4+ and charge asymmetry for CH4+. However, for both the molecules LSD calculations show that the symmetric charge density which has multiplicity of zero is highly unfavored. We note that in spin polarized calculations the exchange interaction is taken care of more accurately than in paramagnetic calculations. Hence the unphysical zero multiplicity is strictly not favored by spin polarized calculations. But LSD does not give consistent results for multiplicity preserved symmetric and asymmetric charge density of SiH4+ and CH4+ molecules. This is a well known shortcoming of the LSD in the context of atomic calculationsJanakPRB81v23 . This deficiency of LSD is removed by the gradient corrected functionalsKutzlerPRL87v59 . We have therefore done spin polarized GGA( GGA-SP) calculations on SiH4+ and CH4+ molecules. Our calculations with GGA-SP functional give the consistent result that the asymmetric charge density is the ground state configuration for both SiH4+ and CH4+ molecules in tetrahedral structure.

Apart from the total energy calculation, one more important issue is the calculation of forces in such systems. We find that enforcing charge symmetry with fractional $`\frac{5}{3}`$, $`\frac{5}{3}`$, $`\frac{5}{3}`$ occupation in the degenerate level of XH4+ in LDA calculations does not distort the tetrahedral structure. This is in contrast with the structural relaxation of tetrahedral XH4+ with 2,2,1 occupation of degenerate level where the molecules results in distorted structure. Although LDA does not distinguish clearly the charge symmetric and asymmetric configurations, the calculated forces are different for the two configurations. This indicates that by enforcing symmetric charge distribution, the forces on atoms of XH4+ molecules are incorrectly estimated. It was also noted by Janak and WilliamsJanakPRB81v23  that enforcing spherical charge distribution on atoms can give rise to errors in Hellman-Feynman forces. Therefore, our results clearly indicate that the electron charge density of tetrahedral XH4+ is asymmetric.

### III.4  Structural evolution

To gain insight into the mechanism of structural distortion in SiH4+ we look at the evolution of SiH4+ structure by using the VASP package with LDA functional and the conjugate gradient(CG) methodSteep . We first ionize SiH4 by holding atoms to obtain tetrahedral SiH4+ and then relax the structure. We find that with 64 steps the system rolls down from the initial tetrahedral to a well converged ground state structure. Each CG step corresponds to a particular geometrical configuration. The bond length between silicon atom and the distance between the hydrogen atoms as a function of CG steps are shown in Fig. 4(a). Initially one H atom moves away from the silicon atom and attains maximum Si-H bond length of 1.92$`\AA `$ at 8th step. Then this repelled H atom moves closer to the nearest hydrogen from 13th step which can be seen from the H-H distance shown in Fig. 4(b). The two hydrogen atoms come very close to each other to a distance of 0.84$`\AA `$ and form an H-H pair. At this stage the structure of SiH4+ is C2v. After 23rd step, this H-H pair gets rotated towards the other two hydrogen atoms without changing the H-H pair distance which finally results in Cs symmetric SiH4+.

The summation of occupied Kohn-Sham orbital energies as a function of CG step is shown in the Fig. 4(c). From the figure we can see that the sum increases initially up to the 8th step and then decreases enormously between 13th and 23rd step. Then it fluctuates between 24 and 40th step and eventually converges to a constant value. The figure shows that the structural changes occurring from 13 and 23 are due to lowering in the orbital energies. Note that there is a minimum at 23rd step. As mentioned earlier at this CG step the structure of SiH4+ is C2v with the H-H bond. In a sense, the sum of orbital energies shown in the figure is similar to the Walsh diagramWalsh  where one takes into account only the sum over orbital energies in predicting the structure. In Walsh diagram, each orbital energy is plotted as function of bond angle of a molecule. Then the minimum of the sum helps in predicting the structure. Thus in the present case the Walsh diagram would predict the symmetry of SiH4+ as C2v symmetry and not Cs, which is the correct symmetry.

### III.5  Role of electrostatic interaction

We now show that the distortion and formation of H-H bond can be easily understood from arguments based on electrostatic interaction between atoms of the molecule. We first try to understand the structural evolution of SiH4+ as discussed in the previous subsection. Let us first consider neutral SiH4 molecule. Since H is more electronegativeRoblesJACS84v106 than Si, the H atoms in SiH4 has a small negative charge while Si will have a small positive charge. This is also shown by population analysisGaussian98 ; Foresman ; Pop  according to which each hydrogen atom holds a charge of -0.12e and silicon +0.52e. When the molecule is ionized to SiH4+ holding its atoms fixed, population analysis shows that the hydrogen atoms have -0.12e, -0.02e, +0.06e and +0.30e charge and silicon +0.77e charge. Interestingly, hydrogen atoms of tetrahedral SiH4+ do not have equal charge although the structure is symmetric. One hydrogen atom has more positive charge than the others and gets repelled by silicon atom which is also positively charged. This hydrogen moves away from the silicon to reduce the electrostatic repulsion and is attracted by the nearby hydrogen atom which has more electrons and forms an H-H bond. This H-H pair is positively charged and hence attracted by the nearby two hydrogen atoms which results in rotation of the pair towards the other two hydrogen atoms.

Since the electronegativities of Ge, Sn and Pb atoms are lower than that of hydrogenRoblesJACS84v106 , GeH4+, SnH4+ and PbH4+ have distortions and bonds similar to that of SiH4+. However, as noted earlier, the structural distortion of CH4+ is different from that of SiH4+; two hydrogen atoms of CH4+ move away from each other which is in contrast with the formation of H-H bond in SiH4+. We now show that the same arguments based on electronegativity and electrostatic interaction also explain the structural distortion of CH4+. Since carbon has larger electronegativity than H atom, carbon in CH4 will have negative charge while hydrogen will have positive charge. This is confirmed by population analysis which shows that a charge of -0.82e on carbon and +0.20e on each hydrogen atom. Ionizing CH4 in the tetrahedral symmetric structure results in +0.20e, +0.23, +0.42e and +0.55e charge on H atoms and -0.40e charge on carbon. We see that the charge distribution is asymmetric and two hydrogen atoms are more positively charged than others. The two H atoms repel each other and open H-C-H bond angle as seen in Fig. 1.

Now we discuss the distortion of negatively charged molecules and explain it on the basis of the above arguments. We see from Fig. 1 that CH4- has the tetrahedral symmetry while other negatively charged molecules have C2v symmetry. Putting an electron on CH4 results in -1.00e on carbon and nearly zero electronic charge on H atoms. Because of this symmetrical charge distribution on hydrogen atoms, the forces on H atoms are radial and the structure of CH4- remain tetrahedral. In contrast, in other molecules the charge on H atoms are not equal. For example, in SiH4- population analysis shows a charge distribution of -0.14e, -0.15e, -0.16e and -0.37e on H atoms while silicon holds a charge of -0.18e. We can see that the H atom which is more negative gets repelled by other atoms. As a result, one of the H-Si-H bond angle opens up with the corresponding expansions in Si-H bond lengths. The trend in the electron distribution of GeH4-, SnH4- and PbH4- is similar to that of SiH4-. Therefore, the distortion of these molecules are similar to that of SiH4- and can be explained using the same reasoning.

### III.6  Charge asymmetry

It is interesting to note that in most of the charged molecules the symmetry is broken. We find that the symmetry breaking is associated with the creation of asymmetry in the charge distribution on H atoms upon charging. Removing or putting an electron on these molecules without changing the tetrahedral geometry may result in charge imbalance on H atoms or equal charges on H atoms as in the case of CH4-. This charge imbalance or balance on H atoms depends on the electronic structure of the molecule. In particular, on the nature and occupation of the HOMO.

Using molecular orbital theoryLevin  we discuss the nature of HOMO in these molecules and show how it contributes to charge distribution among H atoms. In Fig. 5, we show the MO construction of XH4 from linear combinations of atomic orbitals. As shown in the figure, there are eight MOs, the lower four are bonding states and the upper four are antibonding states. The bonding and antibonding states have either A1 symmetry or T2 symmetry. A1 level is non-degenerate while T2 level is triply degenerate. While the bonding energy level of A1 symmetry is always lower in energy than T2, the order of antibonding levels depends on the molecule. For CH4 antibonding A1 level is lower in energy than T2 while for SiH4, GeH4, SnH4 and PbH4 antibonding T2 level is lower than A1. From the figure we can see that the MO with A1 symmetry is a linear combination of valance atomic s orbitals of X and H atoms. The MO with T2 symmetry is a linear combination of valance p orbitals of X and s orbitals of H atoms. Because of mixing between p and s character, T2 orbitals have strong directionality along H atoms. If T2 orbital is fully occupied, it has equal weight along H atoms and results in a uniform charge distribution on H atoms as in the case of neutral XH4 molecules. If T2 orbital is incompletely filled with unequal electrons, it has unequal weight along H atoms and results in a non-uniform charge distribution on H atoms.

Removing an electron from tetrahedral XH4 results with 5 electrons and the T2 orbital is incompletely filled. Hence the charges on H atoms of tetrahedral XH4+ molecules are not equal. This can be seen from Fig. 6(a) and (b) which show the charge densities of SiH4 and SiH4+ in tetrahedral structure. We see from the figure that while the charge densities on H atoms of SiH4 are same, the charge densities on H atoms of SiH4+ are different. In the case of negatively charged molecules the HOMO can be of T2 or A1 symmetry. The HOMO of CH4- has A1 symmetry. Since A1 is a linear combination of valance s atomic orbitals with equal weights from hydrogen atoms, the charge distribution is equal on H atoms of CH4- as shown in Fig. 6(c). The higher electron density around the carbon atom as seen in the figure can be attributed to the higher electronegativity of carbon compared to hydrogen. The HOMO of other negatively charged molecules has T2 symmetry and is singly occupied. This leads to asymmetric charge distribution on H atoms. Note that since A1 level in CH4- is non-degenerate, CH4- does not show the Jahn-Teller effect, and also no charge asymmetry. On the other hand, if the HOMO has T2 symmetry and is not completely occupied as in the case of SiH4+, the molecule will show the Jahn-Teller effect and the charge asymmetry. This implies that the creation of the charge asymmetry and the Jahn-Teller effect are intimately connected.

## IV  Conclusions

We have investigated the ground state structures of neutral and charged XH4 molecules using first principle electronic structure methods. We find that charging XH4 molecules can break the symmetry and change their structures drastically. Negatively charged molecules are distorted by pushing two hydrogen atoms away from each other while the positively charged molecules get distorted with an unusual H-H bond formation except for X=C. The calculations done by various methods and the charge density analysis confirm the existence of H-H bond. Furthermore, vibrational analysis shows that the strength of H-H bond in the positively charged molecule is stronger than H2+ and weaker than H2. We find that one vibrational mode of Cs symmetric XH4+ molecules has low frequency and corresponds to internal rotation of the H-H pair. Our vibrational analysis indicates that the molecules are not transition states. The fragmentation behavior of XH4+ into XH2 and H2 shows that the stability decreases from CH4+ to PbH4+.

From the total energy calculations using LDA, LSD, GGA and GGA-SP functionals we have shown that the ground state electronic structure of a charged tetrahedral molecule has asymmetric charge density corresponding to integer occupation. We have provided an insight into the symmetry breaking and shown that the distortion and unusual H-H bond formation can be understood easily using simple arguments based on charge asymmetry and electrostatic interaction between the atoms. Unequal charges on hydrogen atoms give rise to different electrostatic forces on different atoms and lead to distortion. On the other hand if the charges are equal on hydrogen atoms, the structural symmetry of the molecule will be preserved. We have shown how the charge asymmetry arises due to character of the HOMO and eventually leads to the symmetry breaking.

###### Acknowledgements.

It is a pleasure to thank Drs. R. Ramaswamy, D. D. Sarma, S. C. Agarwal and N. Sathyamurthy for helpful discussions and comments. This work was supported by the Department of Science and Technology, New Delhi via project No. SP/S2/M-51/96.

FIGURE CAPTIONS

Fig. 1. Ground state structures of neutral and charged XH4 molecules.
Fig. 2 Electron density contour of SiH4+ in the plane of Si and the two bonded H atoms.
Fig. 3 The total energy of SiH4+, H2 and H2+ molecules as a function of H-H distance. The total energy of the molecule is given with respect its ground state energy.
Fig. 4. (a) The Si-H bond lengths and (b) H-H distances in SiH4+ as a function of CG steps. The H atoms have the same label as those of SiH4+ in Fig. 1. (c) The sum of energy levels of SiH4+ as a function of CG steps.
Fig. 5. A schematic energy-level diagram showing the construction of molecular orbitals of XH4 molecule from atomic orbitals. The levels on the extreme left side are the valance states of X atoms and the level on the extreme right is the valance state of H atom. The levels which are in the center are the molecular states of XH4.
Fig. 6. The isosurface contour of valance electrons of (a) SiH4 (b) SiH4+ and (c) CH4- obtained by using the Gaussian package. The value of the contour is chosen such that one able to distinguish the difference between the charge density on the H atoms. The densities shown in figure have the contour value of 0.20 for SiH4 and SiH4+ molecules and 0.3 for CH4-.

Table. I. Relative energies ( in eV) of XH4+ molecules for various symmetries with respect to the ground state structures. The results are obtained using the Gaussian98 with the BPW91 exchange-correlation functional.

Table. II Properties of H-H bond calculated using the Gaussian98 package with the BPW91 exchange-correlation functional.

Table. III. The energy required to fragment XH4+ into XH2+ and H2. The results are obtained using the Gaussian98 package with the BPW91 exchange-correlation functional.

Table. IV The calculated vibrational frequencies $`\nu `$( in cm-1 ) and force constant k( in eV/$`\AA `$ ) of neutral and charged XH4 molecules using the Gaussian98 package with the BPW91 exchange-correlation potential.

Table. V Relative energies( in eV ) of SiH4+ and CH4+ in tetrahedral structure for various electronic configuration with respect to charge asymmetric configuration, using different exchange-correlation functionals.


text:
# Belle 𝐽/𝜓+𝜂_𝑐 anomaly and a very light scalar boson

## I Introduction

The best understood QCD system is perhaps the heavy quarkonium, made up of heavy quark and antiquark. Known systems include charmonium, bottomonium, and $`B_c`$ mesons. We thought we understand these systems, but surprises often come up as more and more data become available. A well-known incident was the surplus of $`J/\psi `$ produced at high transverse momentum at the Tevatron cdf . The color-singlet model singlet  missed the cross section by more than an order of magnitude. The color-octet mechanism octet  and subsequently the nonrelativistic QCD model nrqcd  were formulated to explain the surplus.

Recently, the Belle Collaboration observed double charmonium production, namely, $`J/\psi +\eta _c`$ at $`\sqrt{s}=10.6`$ GeV belle . The production rate is anomalous belle ; pakhlov ; talk

$`\sigma (e^+e^{}J/\psi \eta _c)`$ $`\times `$ $`B(\eta _c2\mathrm{charged})`$
$`=(0.033_{0.006}^{+0.007}\pm 0.009)\mathrm{pb},`$
which is roughly an order of magnitude larger than any QCD predictions kaoru . It poses a strong threat to our understanding of quarkonium. However, one has to note that the $`J/\psi `$ is detected via its leptonic decay while $`\eta _c`$ is inferred from the recoil mass spectrum of the $`J/\psi `$. Therefore, the $`\eta _c`$ is actually not identified, and perhaps it could be a $`J/\psi `$, or something else. Based on this argument, Bodwin, Lee, and Braaten bbl  argued that Belle in fact observed double $`J/\psi `$’s, which are produced via double virtual photons. However, the further analysis by Belle does not seem to support this idea pakhlov . Thus, the anomalous production remains one of the challenging problems in QCD, and arouses speculations. For example, a heavy glueball of a mass very near $`m_{\eta _c}`$ with a rather large matrix amplitude has been proposedBrodsky:2003hv  to account for the excess.

In this Letter, we study a possible scenario of the existence of a light scalar boson $`H`$, which has an enhanced coupling to the charm quarks. It contributes to both the $`J/\psi +\eta _c`$ production, and $`J/\psi `$ plus open charm pair production. We show the allowed parameter space region of the mass of this scalar and the coupling to charm, fitted to the Belle data. We find that a coupling constant $`g_H1`$ and $`m_H5`$ GeV can fit the data well. An immediate test for this light scalar boson is to look for a peak around 5 GeV in the recoil mass spectrum of $`J/\psi `$, because this scalar can also be produced directly in association with a $`J/\psi `$. So far, the Belle collaboration has only shown the recoil mass spectrum with very fine bin size up to 4 GeV belle ; pakhlov ; talk , where $`\eta _c,\chi _{c0}`$, and $`\eta _c(2S)`$ can be seen. However, the bin size in the region above 4 GeV is too large. Therefore, we urge the Belle collaboration to scan with very fine steps in the region above 4 GeV. (Even with a large bin size, we notice that there may be something peculiar at around 5.2–5.4 GeV bins belle ; pakhlov .) Finally, we point out that another process $`e^+e^{}H\gamma `$, which is negligible within the SM, becomes sizable in the our scenario. This is another interesting test for the scenario.

## II $`e^+e^{}J/\psi \eta _c`$

Suppose the interaction of the scalar boson $`H`$ with the charm-quark pair is given by

$$=g_H\overline{c}cH,$$
(1)
with a coupling constant $`g_H`$ that will be determined later.

Assuming that the scalar decays into a $`c\overline{c}`$ pair only, the width of $`H`$ is given by

$$\mathrm{\Gamma }_H=\frac{3}{8\pi }g_H^2m_H\left[1\frac{4m_c^2}{m_H^2}\right]^{3/2}.$$
(2)
For $`m_H5`$ GeV and $`g_H1`$, the $`\mathrm{\Gamma }_H0.06m_H`$, which is narrow compared with its mass but much wider than the widths of $`J/\psi `$ and $`\eta _c`$. In the following, we use the Breit-Wigner prescription for the scalar boson propagator.

One of the contributing Feynman diagrams for $`e^+e^{}J/\psi \eta _c`$ is shown in Fig. 1. The sum of the amplitudes of the four contributing Feynman diagrams for $`e^{}(p_1)e^+(p_2)J/\psi (k_1)\eta _c(k_2)`$ is given by

$`i`$ $`=`$ $`4g_H^2e^2Q_eQ_c{\displaystyle \frac{R_\psi (0)R_{\eta _c}(0)}{\pi }}{\displaystyle \frac{1}{s^2}}{\displaystyle \frac{1}{s/4m_H^2+i\mathrm{\Gamma }_Hm_H}}`$ (3)
$`\times `$ $`\overline{v}(p_2)\gamma ^\mu u(p_1)ϵ_{\mu \rho \lambda \sigma }k_2^\rho k_1^\sigma ϵ^\lambda (k_1)`$
where $`ϵ(k_1)`$ is the polarization 4-vector of the $`J/\psi `$, and $`s`$ is the square of the center-of-mass energy of the collision. $`R_\psi (0)`$ and $`R_{\eta _c}(0)`$ are the wavefunctions at the origin, with the normalization $`_0^{\mathrm{}}|R(r)|^2r^2𝑑r=1`$. Here the center-of-mass energy $`\sqrt{s}`$ of the collision is at $`10.6`$ GeV, and so we simply only take the intermediate photon approximation (the intermediate $`Z`$ contribution is suppressed by $`s/m_Z^20.013`$.)

The angular distribution of the process is given by

$`{\displaystyle \frac{d\sigma }{d\mathrm{cos}\theta }}`$ $`=`$ $`\alpha ^2g_H^4Q_e^2Q_c^2{\displaystyle \frac{|R_\psi (0)|^2|R_{\eta _c}(0)|^2}{\pi }}{\displaystyle \frac{1}{s^2}}`$ (4)
$`\times {\displaystyle \frac{1}{(s/4m_H^2)^2+\mathrm{\Gamma }_H^2m_H^2}}\beta ^3(1+\mathrm{cos}^2\theta ),`$
where

$$\beta =\sqrt{\left(1\frac{m_\psi ^2}{s}\frac{m_{\eta _c}^2}{s}\right)^24\frac{m_\psi ^2}{s}\frac{m_{\eta _c}^2}{s}}.$$

Integrating the angle we obtain the total cross section

$`\sigma `$ $`=`$ $`{\displaystyle \frac{8}{3}}\alpha ^2g_H^4Q_e^2Q_c^2{\displaystyle \frac{|R_\psi (0)|^2|R_\eta (0)|^2}{\pi }}{\displaystyle \frac{1}{s^2}}`$ (5)
$`\times {\displaystyle \frac{1}{(s/4m_H^2)^2+\mathrm{\Gamma }_H^2m_H^2}}\beta ^3`$
Let us examine the allowed parameter space of $`m_H`$ and $`g_H`$ that can fit the Belle cross section. Combining the systematic and statistical errors of the Belle data the $`1\sigma `$ range is $`0.033\pm 0.011`$ pb. We show the contours of the central value, $`\pm 1\sigma `$, and $`\pm 2\sigma `$ in Fig. 2. Taking $`g_H=1`$ the allowed central value ($`2\sigma `$ range) for $`m_H`$ is

$$m_H=5.10(4.875.20)\mathrm{GeV}\mathrm{or}\mathrm{\hspace{0.33em}\hspace{0.33em}5.45}(5.365.66)\mathrm{GeV}.$$
(6)
Here we used $`|R_\psi (0)|^2=|R_\eta (0)|^2=0.8\mathrm{GeV}^3`$.

It is straightforward to adapt our scenario to explain other peaks in the recoil mass spectrum between $`33.8`$ GeV, i.e, $`\chi _{c0}`$ and $`\eta _c(2S)`$ belle ; pakhlov . Since our scenario only requires the $`m_H`$ at around $`\sqrt{s}/2`$, it would also explain the enhancement to the production of $`J/\psi +(\eta _c(1S),\chi _{c0},\eta _c(2S))`$. Note that our scenario would not give enhancement to $`J/\psi +(J/\psi ,\chi _{c1},\chi _{c2},\psi (2S))`$ because of the spin-parity-charge-conjugation of the photon propagator. This is in agreement with the observation by Belle pakhlov .

Note that if we choose a smaller $`g_H`$, the corresponding central value $`m_H`$ of the fit to the Belle data would become closer to $`\sqrt{s}/2=5.3`$ GeV, and the $`2\sigma `$ range would also become smaller, as shown in Fig. 2. The smaller the value of $`g_H`$, the more fine-tuned value for $`m_H`$ that the fit gives. So in the following predictions of our scenario we simply take $`g_H=1`$ and the results just scale as some powers of $`g_H`$.

If we replace the scalar $`H`$ by a pseudoscalar boson $`A^0`$, the diagram similar to Fig. 1 with $`H`$ replaced by $`A^0`$ gives vanishing amplitude in the static approximation. There is another diagram with the virtual $`A^0`$ bremsstrahlung off the charm quark turning into the $`\eta _c`$ state. However, such a contribution is important only when the mass $`m_{A^0}`$ is around $`34`$ GeV, which is unfavorable based on the constraint of the Belle search.

The static approximation of using the wavefunction at the origin may not be so valid when the Higgs pole is involved, however our simple calculation serves the purpose of estimating of the size of the production rate due to the light Higgs boson $`H`$.

## III $`e^+e^{}J/\psi c\overline{c}`$

This process is very similar to the original $`J/\psi +\eta _c`$ production, except that the $`c\overline{c}`$ pair does not form a bound state but goes into an open pair. They will hadronize into $`D`$ mesons. The calculation is straightforward, but however the result is not simple enough to be put here. We show the cross section verse $`m_H`$ at $`\sqrt{s}=10.6`$ GeV in Fig. 3. The cross section scales as $`g_H^4`$. At $`m_H5`$ GeV, $`\sigma (e^+e^{}J/\psi c\overline{c})0.15`$ pb for $`g_H=1`$. Therefore, it contributes substantially to $`J/\psi +`$ open charm production.

## IV $`e^+e^{}\psi H`$

The sum of the two Feynman diagrams of the process $`e^{}(p_1)e^+(p_2)\psi (k_1)H(k_2)`$ is given by

$`i`$ $`=`$ $`2\sqrt{3}e^2Q_eQ_cg_H{\displaystyle \frac{R_\psi (0)}{\sqrt{\pi m_\psi }}}{\displaystyle \frac{1}{s}}{\displaystyle \frac{1}{sm_\psi ^2+m_H^2}}`$
$`\times \overline{v}(p_2)\gamma _\mu u(p_1)\left[(sm_H^2+m_\psi ^2)ϵ^\mu (k_1)2k_2ϵ(k_1)k_1^\mu \right]`$
The differential cross section is

$`{\displaystyle \frac{d\sigma }{d\mathrm{cos}\theta }}`$ $`=`$ $`6\alpha ^2g_H^2Q_e^2Q_c^2{\displaystyle \frac{|R_\psi (0)|^2}{m_\psi }}{\displaystyle \frac{1}{(sm_\psi ^2+m_H^2)^2}}\beta `$ (7)
$`\times \left[4{\displaystyle \frac{m_\psi ^2}{s}}+{\displaystyle \frac{1}{2}}\beta ^2(1+\mathrm{cos}^2\theta )\right]`$
where

$$\beta =\sqrt{\left(1\frac{m_\psi ^2}{s}\frac{m_H^2}{s}\right)^24\frac{m_\psi ^2}{s}\frac{m_H^2}{s}}$$

The integrated cross section is

$$\sigma =6\alpha ^2g_H^2Q_e^2Q_c^2\frac{|R_\psi (0)|^2}{m_\psi }\frac{1}{(sm_\psi ^2+m_H^2)^2}\beta \left[8\frac{m_\psi ^2}{s}+\frac{4}{3}\beta ^2\right]$$
(8)
In Fig. 4, we show the production cross section of $`\sigma (\psi H)`$ with $`g_H=1`$. The cross section scales as $`g_H^2`$. The scalar boson $`H`$ so produced is heavy enough to decay into an open charm-quark pair. Thus, it will give rise to open charm pair production, with the pair around 5 GeV. If the Belle could scan the recoil mass spectrum of $`J/\psi `$ further up to $`56`$ GeV, they should be able to see if this scalar boson $`H`$ exists. This is the most immediate test for the scenario.

## V $`e^+e^{}H\gamma `$

In the SM, the leading contribution to the process $`e^+e^{}H\gamma `$ comes from the one-loop diagrams with the $`W`$ boson and heavy fermions running in the loop. In our scenario with $`H`$ coupling strongly to charm, this channel could be sizable for detection. The SM calculation was completed in Ref. leve . We use its formulas with the charm-quark contribution only, and modify the SM coupling to our enhanced coupling as follows:

$$\sqrt{2}G_Fm_c^2g_H^2.$$

We show the ratio of $`\sigma (e^+e^{}H\gamma )/\sigma (e^+e^{}\mu ^+\mu ^{})`$ at $`\sqrt{s}=10.6`$ GeV for $`m_H=39`$ GeV with $`g_H=1`$ in Fig. 5. The ratio scales as $`g_H^2`$. Note that with the SM coupling strength $`(gm_c/2m_W)`$, this ratio is very small, of order $`O(10^9)`$. With the enhanced charm coupling, the ratio is of order $`5\times 10^5`$ for $`m_H=5`$ GeV. Thus, the cross section of $`e^+e^{}H\gamma `$ for $`m_H=5`$ GeV at $`\sqrt{s}=10.6`$ GeV is about 40 fb, which will give thousands of events at Belle or BaBar with an integrated luminosity of 100 fb<sup>-1</sup>.

## VI Discussion

One may concern about the high energy analog of $`e^+e^{}c\overline{c}H`$ production. In high energy collisions, we would expect extra $`c\overline{c}c\overline{c}`$ in the final state, with one of the $`c\overline{c}`$ pairs clustering at around 5 GeV. To our knowledge there have not been any measurement on four-charm production, nor any searches explicitly looking for anomalous four-charm events at LEP or at the Tevatron. There were perhaps some mysterious 4-jet events recorded by ALEPH aleph , but it had no explicit identification of charm in those events. Also, the search for four-charm events suffers severely from the weak charm identification, which is, to some extent, overlapped with the b-quark identification. Therefore, there is no explicit constraints on four-charm production coming from high energy experiments.

Another possible source of constraints comes from Upsilon $`\mathrm{{\rm Y}}(1S)`$ decays. From the Particle Data Book PDG , the branching ratio $`B(\mathrm{{\rm Y}}(1S)J/\psi +X)=1.1\pm 0.4\%`$. Given such a large error it is impossible to rule out our scenario, because $`J/\psi +\eta _c`$, $`J/\psi +H`$, or $`J/\psi +`$ open charm could only give a branching ratio below $`0.4\%`$. There is also a limit on $`B(\mathrm{{\rm Y}}\gamma +X)<3\times 10^5`$, where $`X`$ is a pseudoscalar boson. In our scenario, we have a scalar boson and so the limit is not directly applicable. Even so, $`B(\mathrm{{\rm Y}}H\gamma )/B(\mathrm{{\rm Y}}\mu ^+\mu ^{})`$ would be very similar to $`\sigma (e^+e^{}H\gamma )/\sigma (e^+e^{}\mu ^+\mu ^{})5\times 10^5`$ (from Fig. 5) without an explicit calculation. It implies $`B(\mathrm{{\rm Y}}\gamma +H)10^6`$, which is well below the existing limit that we stated above.

The anomalous double charmonium $`J/\psi +\eta _c`$ production measured by Belle posed a very strong challenge to our understanding of QCD. We have investigated the scenario of a light scalar boson of mass about $`5`$ GeV with a large coupling (of order the weak coupling or even larger) to the charm quark. We found that it could explain the $`J/\psi +\eta _c`$ anomaly, and provided the much-needed $`J/\psi `$ plus open charm production. It would also explain the enhancement to the production of $`J/\psi +(\eta _c(1S),\chi _{c0},\eta _c(2S))`$, but not to $`J/\psi +(J/\psi ,\chi _{c1},\chi _{c2},\psi (2S))`$ because of the spin-parity-charge-conjugation of the photon propagator. This is in agreement with the observation by Belle. An immediate test for this light scalar boson is to scan the recoil mass spectrum of $`J/\psi `$ and look for a peak around 5 GeV, because this scalar can also be produced directly in association with a $`J/\psi `$. So far, the Belle collaboration has only shown the recoil mass spectrum with very fine bin size up to 4 GeV belle ; pakhlov , where $`\eta _c,\chi _{c0}`$, and $`\eta _c(2S)`$ can be seen. However, the bin size in the region above 4 GeV is too large. We, therefore, urge the Belle collaboration to scan the recoil mass spectrum with very fine bin size above 4 GeV. (Even with a large bin size, we notice that there may be something peculiar at around 5.2–5.4 GeV bins belle ; pakhlov , but this is only our wild guess.) Finally, we have pointed out that another process $`e^+e^{}H\gamma `$, which is negligible within the SM, becomes sizable in the new scenario of a light scalar boson, which couples strongly to charm. This is another test for the scenario.

One may ask where this scalar boson comes from. Such a light Higgs may come from the general form of the multi-Higgs doublet sector. The required enhanced coupling to the charm quark implies dedicated fine tuning. It also requires additional fine tuning to kill the unwanted flavor changing neutral processes. If we stay within the scope of the two-Higgs doublet sectors with the quark couplings like those in the minimal supersymmetric standard model, the coupling strength of $`H`$ to charm is $`gm_c\mathrm{cot}\beta /2m_W`$. Therefore, in order to have an effective coupling $`g_Hg`$ we need a $`\mathrm{cot}\beta 2m_W/m_cO(100)`$. Such a large $`\mathrm{cot}\beta `$ is unnatural and perhaps causes a big problem to the top Yukawa coupling. We do not offer a good answer to the structure of the Higgs sector in this paper.

We thank Dr. Jungil Lee for discussion.


text:
# Metal-Insulator Transition in Two- and Three-Dimensional Logarithmic Plasmas

(Received )

## Abstract

We consider the scaling of the mean square dipole moment in a plasma with logarithmic interactions in a two- and three-dimensional system. In both cases, we establish the existence of a low-temperature regime where the mean square dipole moment does not scale with system size and a high-temperature regime where it does scale with system size. Thus, there is a nonanalytic change in the polarizability of the system as a function of temperature, and hence a metal-insulator transition in both cases. The relevance of this transition in three dimensions to quantum phase transitions in $`2+1`$-dimensional systems is briefly discussed.

Gauge theories where matter fields are coupled to compact gauge-fields in $`2+1`$ dimensions have recently been much studied as effective theories for Mott insulators with competing orders Sachdev ; Moessner ; Senthil . Compact $`U(1)`$ gauge theories in $`2+1`$ dimensions are capable of sustaining topological defects in the gauge fields in the form of space-time instantons Polyakov . Such theories generically offer the possibility of featuring confinement-deconfinement transitions associated with a proliferation of topological defects in the gauge sector. It is hoped that phase transitions in matter coupled such theories FradkinShenker ; EinhornSavit  may be connected to difficult problems in condensed matter physics, such as the breakdown of Landau Fermi liquid theory and possibly also spin-charge separation in strongly correlated systems at zero temperature in two spatial dimensions SenthilFisher . One such model considered recently is the compact abelian Higgs model where the gauge field is coupled to matter fields with the fundamental charge KleiNogSu . By itself, the compact $`U(1)`$ gauge sector may be mapped onto the $`2+1`$-dimensional Coulomb gas with $`1/R`$-interactions between “point charges”, i.e., the instantonic topological defects of the theory Polyakov . As shown by Polyakov, the latter system is always in a metallic phase (when the instantons are regarded as electric charges) Polyakov . More recently, it was shown that when matter carrying the fundamental charge was coupled to the compact gauge sector, and critical matter field fluctuations were integrated out, the system could be mapped onto a gas of point charges interacting with a potential $`\mathrm{ln}(R)`$ in $`2+1`$ space-time dimensions, instead of $`1/R`$ KleiNogSu . Using RG arguments, it was demonstrated that the system with the $`2+1`$-dimensional logarithmic interaction may undergo a finite-temperature phase transition driven by the unbinding of dipole configurations KleiNogSu , from a low-temperature dielectric regime to a high-temperature metallic regime(again when the instantons are regarded as electric charges). Here, we demonstrate this using Monte Carlo (MC) simulations.

Whether or not there exists a low-temperature dielectric regime separated from a high-temperature metallic regime in a three-dimensional system of point charges interacting with logarithmic interactions and overall charge-neutrality, is presently a matter of debate in both the condensed matter and lattice gauge theory literature KleiNogSu ; Herbut ; Chernodub . The theorem that dipoles cannot screen the Coulomb potential (in any dimension) Frohlich  does not apply to the case of a logarithmic pair potential between point charges in three dimensions. It is due to this inability of the dipoles to screen a Coulomb potential that the Kosterlitz-Thouless (KT) transition is possible: the low temperature dielectric phase is always critical below the KT transition temperature. In three dimensions, on the other hand, Debye-Hückel theory is essentially exact. Thus, starting from a system of charges, screening is such that the screening length cannot become infinite, and the system is always in the metallic phase. On the other hand, if one starts from a three-dimensional dipole system, the system always stays in the dielectric phase, since dipoles cannot screen. Once more, there is no phase transition. If some departure from the Coulomb potential occurs in three dimensions, dipoles could conceivably be able to screen. A $`3D`$ logarithmic interaction could in principle be screened to a potential that decays with separation between the charges Herbut . If so, no metal-insulator transition occurs at any finite temperature, analogous to the situation in a three-dimensional Coulomb plasma Polyakov . If one thinks of the gas of point charges as a gas of instantonic defects in a compact gauge field, the lack of a metal-insulator transition due to non-standard screening by dipoles would correspond to permanent confinement in the compact gauge theory.

We therefore consider a logarithmic plasma in two and three dimensions on the lattice using MC simulations, and study the polarizability of the system as a function of temperature. The screening properties of these systems determine whether they are insulators or metals and are governed by the dielectric constant $`\epsilon `$, which in a low-density approximation is given by the polarizability $`p`$ of the system, $`\epsilon =1+n_d\mathrm{\Omega }_dp`$. Here $`n_d`$ is the dipole density, and $`\mathrm{\Omega }_d`$ is the solid angle in a $`d`$-dimensional system. Since the polarizability of the system is proportional to the mean square separation $`s^2|\stackrel{}{r}_i\stackrel{}{r}_j|^2`$ between the charges constituting the dipoles, it is natural to focus on $`s^2`$ in order to investigate whether the system is a dielectric or a metal. We find a low-temperature regime where dipoles are tightly bound, separated from a high-temperature regime where they are unbound. This implies the existence of a low-temperature insulating regime and a high-temperature metallic regime, separated by a genuine phase transition. For comparison and benchmark purposes, we compute the same quantities for the two-dimensional logarithmic plasma, where a metal-insulator transition in the form of a KT phase transition is known to exist KosterlitzThouless73 ; MinnhagenReview .

We consider the Hamiltonian

$$=\frac{1}{2}\underset{i,j}{}q_iV(\stackrel{}{r}_i\stackrel{}{r}_j)q_j,$$
(1)
with

$$V(\stackrel{}{r})=\frac{4\pi ^2}{Na^d}\underset{\stackrel{}{k}}{}\frac{[e^{i\stackrel{}{k}\stackrel{}{r}}1]}{[2(d_{\alpha =1}^d\mathrm{cos}k_\alpha a)]^{d/2}},$$
(2)
where the sum is over all pairs of sites of a two- or three-dimensional periodic lattice and $`q_i`$ is the charge at site $`i`$, $`N=L^d`$ is the total number of sites, and $`a`$ is the lattice spacing. We have subtracted $`V(\stackrel{}{r}=0)`$ from the potential, since only neutral configurations will contribute to the partition function, and we work at zero chemical potential for the instantons LeeTeitel . This has the usual Coulomb form for $`d=2`$, but differs from it in the $`d=3`$ case: instead of the usual $`1/R`$ potential of the three-dimensional Coulomb gas ($`3D`$CG), the power $`3/2`$ results in a logarithmic potential. Taking the continuum limit of Eq. (2) for $`d=3`$, we find

$`V(\stackrel{}{r})={\displaystyle \frac{d^3k}{(2\pi )^3}\frac{[e^{i\stackrel{}{k}\stackrel{}{r}}1]}{(k^2+m^2)^{3/2}}}`$
$`={\displaystyle \frac{1}{2\pi ^2}}\left[K_0(mr)+{\displaystyle \frac{\mathrm{\Lambda }}{\sqrt{\mathrm{\Lambda }^2+m^2}}}\mathrm{arcsinh}\left({\displaystyle \frac{\mathrm{\Lambda }}{m}}\right)\right],`$ (3)
where $`K_0(x)`$ is a Bessel function, $`\mathrm{\Lambda }=2\pi /a`$, and the Debye-Hückel screening length $`m^1`$ was introduced as an infrared regulator. As $`m0`$ we have

$$V(r)=\frac{1}{2\pi ^2}[1\gamma \mathrm{ln}(\mathrm{\Lambda }r)]+𝒪\left(\frac{m}{\mathrm{\Lambda }}\right)^2,$$
(4)
where $`\gamma `$ is the Euler constant. Therefore, Eq. (2) behaves exactly in the same way as the two-dimensional potential.

Let us first turn to the $`2D`$CG which we use as a benchmark on our method of probing the metal-insulator transition, before proceeding to the three-dimensional logarithmic gas ($`3D`$LG). At low temperatures this model consists of tightly bound dipoles with the two opposite charges separated by a distance $`|\stackrel{}{r}_i\stackrel{}{r}_j|s`$ of the order the lattice constant.

In a $`2D`$ system, it is well-known that dipoles begin to unbind at a critical temperature $`T_{KT}`$ KosterlitzThouless73  and at high temperatures the $`2D`$CG is a fully ionized metallic plasma, separated from the low-temperature dielectric insulating phase by a genuine KT phase-transition. Thus, we expect no finite-size scaling of $`s^2`$ below the KT-transition, whereas we should expect $`s^2L^{\alpha (T)}`$ with $`\alpha (T)2`$ at higher temperatures. Using an intuitive low density argument, neglecting screening effects KosterlitzThouless73 , we can calculate the behaviour of $`s^2`$ to leading order in $`L`$,

$$s^2\{\begin{array}{cc}\mathrm{Const}.\hfill & ;T<T_{KT}\hfill \\ L^{(TT_{KT})/T}\hfill & ;T_{KT}<T<2T_{KT}\hfill \\ L^2\hfill & ;\mathrm{\hspace{0.33em}\hspace{0.33em}2}T_{KT}<T.\hfill \end{array}$$
(5)
Hence, $`\alpha (T)`$ is zero for low temperatures and a monotonically increasing function of temperature just above $`T_{KT}`$.

Including screening effects in $`2D`$ shows that this conclusion still holds, however the temperature at which it occurs is determined by screening. This is ultimately related to the fact that dipoles contribute only through a correction to the coefficient of the bare logarithmic pair-potential. Hence, in $`2D`$ the functional form of the renormalized potential is unaltered, only prefactors are changed. In three dimensions, it is far from obvious that dipoles do not have a much more disruptive effect on a bare logarithmic pair-potential.

To adress the issue of scaling of $`s^2`$, we will use MC simulations and finite-size scaling to study $`s^2`$ both in the $`2D`$CG as well as in the $`3D`$LG. In our simulations, the particle number is not conserved. However, during the simulations the system is maintained electrically neutral. The MC moves involve creation and annihilation of charges, applying the Metropolis algorithm in this process. Starting at some randomly chosen lattice site, an attempt to insert a negative or positive charge at random at this site is made, with an opposite charge at a nearest-neighbour site. The move is accepted with probability $`\mathrm{exp}(\mathrm{\Delta }E/T)=\mathrm{exp}[(_{\mathrm{new}}_{\mathrm{old}})/T]`$, and this is done in all $`d`$ directions, before we move to the next site. It is clear that placing a charge on top of an opposite one corresponds to the annihilation of the existing one. In order to measure $`s`$ we have to keep track of which two charges belong to each other in a dipole, since there is no physical link defining the dipole. Two charges inserted into the lattice at the same MC move is chosen as a dipole, and if one of these charges subsequently is annihilated, the effect of this annihilation is to diffuse one of the charges in the dipole. Hence, such a diffusion may increase or decrease $`s`$ for this dipole footnote . One sweep is defined as going through all the sites in the lattice once, and at every tenth sweep we sample $`s^2`$ averaged over the system. The sample is, however, rejected if there are no charges in the system at this MC time, since there is no information on $`s^2`$ in such a configuration. The thermal average obtained at the end of the simulation is thus taken only over non-empty configurations.

The results from the simulations on the $`2D`$CG are presented in Fig. 1. A total of $`10^5`$ MC sweeps at each temperature are used to produce the first plot. From the $`s^2`$-data of system sizes $`L=`$ 8, 12, 16, 20, 24, 28, 32, 36, 40 and 48, we extract $`\alpha `$ for temperatures in the interval $`(1.0,1.52)`$, which is shown in the second plot. It is evident that there are two distinct regimes of temperatures, one in which the charges of almost all dipoles are bound as tightly as possible, the separation of the charges correspond to the lattice constant. In the high-temperature regime the dipoles have started to separate, reflected by a scaling of $`s^2L^{\alpha (T)}`$ with the system size. The two regimes are necessarily separated by a phase transition, since in the low-temperature regime $`\alpha (T)=0`$ while in the high-temperature regime $`\alpha (T)0`$. This necessarily implies a non-analytic behavior of $`\alpha (T)`$. An attempt to determine the transition temperature from these plots yields approximately $`1.32`$. This is slightly less than the early results of Saito and Müller-Krumbhaar of $`1.35`$ in our units SMK , but is in excellent agreement with much more recent simulations by Olsson Olsson . It provides confidence in the method of locating the critical point by monitoring the quantity $`s^2`$, even when the system is subjected to periodic boundary conditions footnote .

Exactly the same simulation technique is applied to the $`3D`$LG and the results are shown in Fig. 2. The finite-size analysis of $`s^2`$ in order to extract $`\alpha `$ is here done on the basis of system sizes $`L=`$8, 12, 16, 20, 24, 28, 32, 36 and 40 and up to $`2.010^5`$ MC sweeps are used. As in the $`2D`$CG we see that the system exhibits two distinct regimes, one insulating regime consisting of tightly bound dipoles and one metallic. There is no scaling of $`s^2`$ with $`L`$ below $`T0.32`$ with the system size. Note that the change in scaling of the mean square dipole moment occurs at a significantly lower temperature than in the $`2D`$CG. This is to be expected, since there is more configurational entropy available in $`3D`$ than in $`2D`$. It is worthwhile comparing this result with the one obtained in Ref. KleiNogSu  where the coupling $`4\pi ^2/t`$ corresponds to the temperature $`T`$ here. There the critical value $`t_c=12\pi ^2`$ was obtained, corresponding to $`T_c=1/3`$ in our case and agreeing well with our numerical result.

The main result is that in the $`3D`$LG and the $`2D`$CG, a low-temperature regime exists where positive and negative charges are bound in tight dipole pairs. This regime is separated from a high-temperature regime where at least a finite fraction of charges are free. The results obtained in three dimensions have the same features as those found in $`2D`$. The scaling exponent $`\alpha (T)`$ for $`s^2`$ is zero in the low-temperature phase and positive in the high-temperature phase. Such a change in $`\alpha (T)`$ cannot be analytical, and therefore the two scaling regimes must be separated from each other via a phase transition. Since $`s^2`$ is a measure of the polarizability of the system, we conclude that the above demonstrates a non-analytic change in the polarizability of the system, i.e., a non-analytic change in the dielectric function of the system as a function of temperature. Hence, a system of point charges with overall charge neutrality interacting with a bare logarithmic interaction undergoes a metal-insulator transition in both $`2D`$ and $`3D`$.

The $`3D`$LG can be shown to be equivalent to an anomalous sine-Gordon theory in $`d=3`$ where the usual $`k^2`$ dispersion is replaced by a $`|k|^3`$ one KleiNogSu , which is non-analytic. This leads to technical difficulties in a standard renormalization group (RG) calculation, where in a perturbative treatment only analytic singularities are generated. Therefore, in a standard RG analysis an analytic correction $`k^2`$ to the dispersion can be generated. On the other hand, since $`|k|^3`$ is non-analytic, it cannot be renormalized within a standard RG analysis. It is conceivable that this is one of the reasons why a recent RG analysis of this theory did not find evidence for a dielectric phase Herbut .

We emphasize that although it is known that a metal-insulator transition occurs in $`2D`$ via a KT transition KosterlitzThouless73 , our numerics by themselves do not demonstrate this. To establish the KT nature of the transition on purely numerical grounds requires convincing numerical evidence that there is a universal jump in the inverse dielectric constant at the transition. The main point of the present simulations is that they settle the difficult question of whether a low-temperature dielectric regime exists at all in a $`3D`$LG, which is a system not subject to standard electrostatics, and for which the usual theorems on screening of charges by dipoles do not apply. The answer is in the affirmative, and we have applied our method also to the $`2D`$ case for comparison and as a benchmark on the correctness of method of monitoring $`s^2`$ to establish the existence of two scaling regimes separated by a phase transition.

The authors acknowledge stimulating discussions and correspondence with I. F. Herbut and S. Sachdev. This work received financial support from the Norwegian University of Science and Technology through the Norwegian High Performance Computing Program (NOTUR), from the Research Council of Norway, Grants No. 158518/431, 158547/431 (NANOMAT), 157798/432 (A.S), and from DFG Priority Program SPP 1116 (F.S.N.). One of us (S.K.) thanks the Norwegian University of Science and Technology for a PhD Fellowship.