text:
# Materials Science and Protein Crystallography Using the MX Beamline Control Toolkit

## I Introduction

The MX beamline control toolkit is a general package for control of synchrotron radiation beamlines and laboratory X-ray generator systems that is being jointly developed by the MR-CATMRCATwebsite MRCAT , IMCA-CATIMCACATwebsite IMCACAT1 IMCACAT2 , and SER-CATSERCATwebsite  sectors at the Advanced Photon Source. It has been described at previous NOBUGS meetings LavenderNOBUGS1997  LavenderNOBUGS2000  and elsewhereLavenderSRI99 . This paper will focus on changes and additions to MX made since the last NOBUGS meeting in 2000.

## II Review of MX

MX is a portable beamline control toolkit available from the web site http://www.imca.aps.anl.gov/mx/. Current users of MX are listed in Table 1. The primary goal of MX is to make it possible to write beamline control system software that can be reused in almost any environment. This goal was conceived based on the author’s frustration with several previous control systems he had worked with that paid little or no attention to transportability of the code.

MX attempts to achieve this goal in several ways:

* MX is designed to be easily ported to new operating systems. It already runs on Linux, Microsoft Win32, MacOS X, Solaris, Irix, HP/UX, and Cygwin with legacy support for SunOS 4, and MSDOS.
* MX is written in ANSI C to maximize its ability to be invoked from other packages. Most applications and scripting languages that support an external calling interface can easily invoke code written in C.
* MX localizes dependencies on particular network protocols to relatively small areas of the code. Thus, while MX comes with its own network protocol, it is designed to make it easy to use other network protocols instead. For example, it is easy to use MX on a beamline that uses only EPICSEPICS  network protocols.
* The core of MX is designed as middleware that can be inserted into other systems as necessary. This is done by providing a high level application programming interface (API) to applications or servers which tries to avoid making assumptions about the nature of those applications and servers. At the low end, a device driver API is provided that tries to avoid making assumptions about the nature of the hardware underneath.

MX has been under development since 1995 and has developed a relatively large number of device drivers. Currently there are 34 drivers for various motor controllers, devices that can be treated as motors, and for access to motors controlled by other control systems such as EPICS. There are also 17 pseudomotor drivers for items like X-ray energy, monochromator control and goniostat table control. In addition, there are a variety of drivers for other devices such as scalers, timers, MCAs, MCSs, and so forth with over 250 drivers in all.

While MX can be viewed as a set of components, it can also be used as a complete beamline control system. It provides an MX server that can be used to manage beamline hardware and a set of graphical user interface tools. At present, MX supports developing GUI programs via its Python and Tcl/Tk scripting interfaces. The most important of these include:

* Imcagui \- A primary user interface for crystallography seen in Figure 1. It provides support for monochromator, slit, and filter control. It also provides a way of executing fluorescence and absorption scans for multiwavelength anomalous diffraction (MAD) experiments as seen in Figure 2 which are then processed by Chooch to compute $`f^{}`$ and $`f^{\prime \prime }`$ values. It was originally developed for IMCA-CAT but is now used at other beamlines as well.
* Optimize and optimize\_auto \- Automatic beamline intensity optimization programs that give beamline users a simplified interface for maximizing X-ray intensity on their own initiative.
* Mxgui \- A more staff-oriented GUI for performing motor moves and executing arbitrary scans.

There are also a number of command line programs and utilities for MX. The most important of these is motor which is routinely used at MR-CAT for materials science experiments.

Most of our existing GUI programs are written in Tcl/Tk. However, we now plan to write most new GUIs in Python instead. Our experience has been that the syntax of Python seems more natural to scientists than that of Tcl. There also appears to be a groundswell of support for Python at synchrotron radiation facilities as well, given Python’s extensive support for scientific and numerical programming. With that in mind, we have now created an MX interface to Python called MP and have begun writing applications using it.

## III XIA Multichannel Analyzers

A major task has been the development of software to interface with the new high performance multichannel analyzers (MCAs) from X-Ray Instrumentation Associates in Newark, California. Our efforts have focused on the DXP-2X and the Saturn, formerly called the DXP-X10P.

The DXP-2X is a CAMAC based system that contains four MCA channels per CAMAC module. The existing MX interface uses the Xerxes library provided by XIA to control the DXP-2X. The Xerxes library provides direct access to many of the low level features of the MCA. Since the DXP-2X is a very powerful and complicated system, it took a substantial amount of effort to correctly program the MX software interface. XIA now provides a new high-level software library called Handel, but this was not available at the time of the original development of our software.

There were a number of problems encountered in getting the DXP-2X to run reliably, but these were partially due to the fact that the DXP-2X firmware was still under development at the time. Most of the problems we encountered have been solved now and the DXP-2X owned by MR-CAT has been reliably used with a 13-element Ge detector from Canberra for several experimental runs now. It has proved to be capable of handling much higher input count rates than the other competing systems that were available at the time of purchase. The DXP-2X was recently used at MR-CAT to acquire data from a sample of Pu(VI) adsorbed on MgO at input count rates of up to $`1.510^6`$ counts per second per MCA channel.

MR-CAT currently controls the DXP-2X via an MX server running on a Windows 98 computers. Since the XIA staff currently do all of their development and testing on Windows platforms, it was regarded as desirable to use the vendor library as provided to maximize XIA’s ability to diagnose problems. Most of the rest of the beamline runs under Linux, so we anticipate moving control of the XIA MCAs over to Linux at some point. If it is moved to the Linux user interface computer, this will also improve performance by eliminating a network transfer that must otherwise be done.

Another XIA multichannel analyzer that we are interested in is the Saturn, which was formerly known as the DXP-X10P. This is a single channel MCA that is controlled via a PC parallel port. The Xerxes library used with the DXP-2X also supports the Saturn. Thus, after a few initial difficulties, we were able to get the MX driver to work with the Saturn as well. An optional feature of the Saturn that is not available for the DXP-2X is 16 TTL output connectors that are associated with the 16 regions of interest (ROIs) of the Saturn. In this configuration, when a detector pulse falls within the boundary of a region of interest, a TTL pulse is synthesized on the corresponding output channel. Thus, in this mode one can think of the Saturn as a set of 16 programmable single channel analyzers whose output can be fed into conventional pulse-counting scaler/timer systems. This is the mode that DND-CAT currently uses.

## IV Quick Scans

Another major effort over the last two years has been the development of software to support quick scans, also known as slew scans or fast scans. This was under initial development at the time of the last NOBUGS meeting and has been routinely used now for two years at MR-CAT.

Quick scans are designed as a direct replacement for many types of step scans. The problem with step scans is that there is a large delay associated with slowing down the motor to a stop to take the data point and then accelerating back to speed. Quick scans avoid this issue by acquiring data continuously while the motors are in motion. Quick scans are available for most motor drivers and for many of the pseudomotor drivers as well. For example, it is possible to quick scan monochromator energy pseudomotors.

Of course, quick scans need a place to store the data as it is acquired. This role is normally played by a multichannel scaler (MCS). A multichannel scaler counts pulses from a set of inputs and then, on receipt of a channel advance pulse, writes the counter values to a buffer and then starts acquiring again. This is then repeated for each data point and at the end of the scan, the buffer containing all of the measurements is read out.

For the data to be useful, we must also record the motor positions for each point in the scan. This is most easily done by recording some signal that depends on the motor position. The most obvious candidate is a quadrature encoder signal from the motor controller that is slaved to the position of the motor. At MR-CAT, this is handled by taking the encoder signals into a small electronics module that converts the signals into two pulse trains. One pulse train corresponds to motion in the positive direction, while the second pulse train corresponds to motion in the negative direction. If the motor controller starts the move a few measurements after the MCS has started acquiring data, then the absolute position of the motor can be determined by taking the difference between the pulse count in the positive motion channel and the pulse count in the negative motion channel.

This mode of quick scanning is currently only possible for motors that have associated encoder outputs. At MR-CAT, only two motors have this capability, namely the monochromator Bragg angle and the 8-circle diffractometer’s primary two theta angle. However, SER-CAT plans to make this possible for almost all motors on the beamline. For SER-CAT, this is practical because they are using the Turbo PMAC line of motor controllers from Delta Tau. The Turbo PMAC is an up to 32 axis motor controller capable of driving both stepper and servo motors, although SER-CAT is using only servo motors. Turbo PMAC controllers are able to slave one of its motor axes to any of the other 31 motors controlled by it. The slave axis can then be commanded to generate a step and direction output which can be converted by external electronics into the same kind of positive and negative pulse trains as are used by MR-CAT. Since the association of master and slave axes can be changed at run time, this makes it possible to quick scan any of the other 31 motors controlled by this PMAC. While SER-CAT does not yet routinely operate in this mode, a prototype of this configuration has already successfully been tested.

At present, most of the support for multichannel scalers in MX is for the Struck SIS3801 MCS. There is a small amount of support for alternate MCS systems, but the main focus has been on the SIS3801. The SIS3801 is a VME-based MCS that can have up to 32 input channels. The SIS3801 acts as an MCS by transferring its measurements into a FIFO that can hold up to 128K scaler values upon receipt of a channel advance pulse. If all 32 channels are in use, this allows up to 4K measurements per channel. However, the SIS3801 can be programmed to ignore some or even most of its input channels, which then leaves more space for the remaining channels. For example, if only 4 channels are used, then up to 32K measurements per channel can be taken.

Currently, MX has two drivers for the SIS3801. The first, called epics\_mcs, is for SIS3801s controlled by the EPICS MCA record written by Mark Rivers of the University of Chicago EpicsMCARecord . The second, called sis3801, talks directly to the module via MX VME I/O. This configuration is intended for use with VME crates that are controlled via PCI-to-VME bus interfaces or by MX servers running directly on a VME crate controller. The existence of two drivers allows MX to support both VME crates controlled by EPICS and VME crates that are not controlled by EPICS. At a higher level, the two MX drivers behave similarly. The primary difference is that the epics\_mcs driver is currently restricted to a maximum of 4000 measurements per channel due to a restriction in the maximum length of EPICS Channel Access messages.

The SIS3801 can be triggered to move to the next measurement by either an internal or an external channel advance. Use of the built-in 10 MHz clock with the internal channel advance allows the SIS3801 to run in a stand alone fashion. However, sometimes it is desirable to synchronize measurements by an SIS3801 with external devices. For this reason, MX also supports the use of the external channel advance signal. Consequently, MR-CAT and SER-CAT have both purchased Struck SIS3807 pulse generator modules to serve as master clocks for the beamline. Theoretically, one could synchronize the operation of multiple SIS3801s in this fashion.

A more interesting use of the external channel advance support is for synchronization with the MultiSCA firmware for the XIA DXP-2X MCA. With this special firmware, the DXP-2X acts similarly to a multichannel scaler in that each channel buffers a new set of 16 ROI integral measurements upon receipt of an external channel advance signal. With a common external clock such as the SIS3807, it will be possible to implement quick scans that also record MCA measurements for each point of the quick scan. Although some preliminary work has been done, this mode of operation has not yet been fully implemented by MX yet. However, we hope to be able to provide this mode of operation at MR-CAT by the end of 2003.

## V Performance

The most important current task in MX development is improvement in the performance of pieces of MX that are perceived to be slow. When examined in detail, most of these issues revolve around performance problems with the default MX network protocol. We are currently addressing this issue in more than one way.

First of all, one possible solution is to use a different network protocol. In an APS environment, the most obvious alternative is the EPICS Channel Access protocol. The author originally was concerned about its dependence on VME crates and an expensive real time operating system. However, recent changes in EPICS licensing and the porting of the server side of EPICS IocCore to Linux, Solaris, Win32, and RTEMS have gone a long way toward mitigating these concerns. MX can already be used as an EPICS client via MX drivers such as epics\_motor, epics\_scaler, and epics\_mcs. We have also begun work on making the MX drivers and toolkit available from the EPICS IocCore server which is described further in Section VII.

We are also pursuing improvements to the existing MX network protocol. It has been suggested that much of the overhead may be due to the formatting and parsing of data sent over the network and due to the lack of caching of the identity of remote records. The author will be focusing his efforts on this and similar ideas for the rest of 2002 and on into 2003. Some of the testing will revolve around improving the performance specifically of network XIA MCA readout, since this makes particularly heavy use of the network.

Perceived responsiveness is also important for user satisfaction GUIBloopers . Early versions of Imcagui, the main user GUI for crystallography, had a beamline status update routine that could take around a second to run. This meant that if the user typed or clicked something on the screen, it might take up to a second for the GUI to acknowledge this. When the status update routine was modified to force frequent updates of the screen, the delay in responding to user input dropped down to a small fraction of a second. The users were much happier with this version of the GUI even though it was really doing work at about the same rate as before.

## VI New and Upgraded Beamlines

A significant amount of the recent work in MX has been in porting the software to new beamlines. Table 1 lists the beamlines where MX is currently in routine use and highlights of the equipment used there. One thing that should be apparent from Table 1 is that MX is able to support a heterogeneous mix of equipment.

Much of the new work has been in support of the SER-CAT sector (Sector 22) at the Advanced Photon Source. SER-CAT uses almost exclusively Delta Tau Turbo PMACs for motion control which are controlled via RS-422 by several Linux-based MX servers. Struck SIS3801 multichannel scalers are used for counter/timer support. They are controlled from a National Instruments PCI-to-VME bus interface via MX drivers.

SER-CAT beamline specific work has included the development of pseudomotors to handle A-frame detector supports of the type originally developed for SBC-CAT as shown in Figure 3. These pseudomotors transform the vertical and horizontal motions of the detector into detector two-theta angle, detector distance, and detector offset. Other work has included the development of software to control the PMAC-based beamline goniostat from the MarCCD and Bruker Proteus CCD systems via MX.

At IMCA-CAT, a recent development has been support for the Rigaku/MSC sample changing robot derived from the original Abbott design. This task and support of the Bruker Proteus system at SER-CAT has motivated the development of a simplified server controlled via ASCII commands. The idea here is to make it as easy as possible for vendors to interface their systems to the beamline.

A recent addition has been the new DND-CAT protein crystallography beamline at APS Sector 5-BM. This beamline implements low level hardware control using their own protocol called SCIPE SCIPEreference . MX applications are supported on this beamline via a set of MX devices drivers such as scipe\_motor, scipe\_scaler, and so forth that are analogous to the ones used to communicate with EPICS controlled devices. Another new user of MX is the Gulf Coast Protein Crystallography Consortium at CAMD. They are operating their beamline with Compumotor 6K motor controllers, Ortec 974 counter/timers, and XIA Saturn MCAs.

One task that will become important in the near future will be more complete automation of crystallography data collection. The ideal will be to make it so that crystals can be shipped to the beamline in standardized containers. These will then be automatically mounted and aligned, followed by CCD data acquisition and processing. The plans for how this will be implemented are not complete yet, but SER-CAT plans to begin working on this very soon. In fact, some early steps have already been taken such as the purchase of x-y sample positioning stages from Oceaneering Space Systems that are to be used as part of automatic crystal alignment.

The first phase for SER-CAT will involve development of a standard interface for interaction with the remote experimenter, including a method for transmitting images of diffraction patterns for remote evaluation. At first, the staff will do the work, and then transmit the results. Later, more real-time interaction will be developed. This will probably need to be in the form of a CGI to allow the interaction and authentication. Control will come later, and is actually easier. The other push will be the development of robotics for automation of alignment and sample changing, which will require motion control programming.

Another aspect of the automation is the development of automatic beamline alignment and optimization scripts, further automation of the fluorescence measurement, automation of data collection strategy for multiaxis diffractometry, and automated data reduction.

## VII MX as Device Support for EPICS

A recent major development in the EPICS world has been the porting of the EPICS server software, called IocCore, to Linux, Solaris, Win32, and RTEMS. This is in beta releases of what is to become EPICS release 3.14. Previous versions of IocCore only ran on VxWorks. However, the device drivers for previous versions of IocCore were generally also VxWorks specific and most have not been ported. This means that the workstation version of IocCore has relatively few device drivers at this point. One of the design goals for MX has been to make it easily embeddable in other packages. Thus, an interesting idea would be to embed MX inside of EPICS IocCore. The EPICS motor record now has a partially finished port to EPICS 3.14, so it is an ideal test case for this. Correspondingly, I am now working with Ron Sluiter to develop such an interface.

Since MX presents the same abstract high level API for all of its motor drivers, my plan is to essentially treat MX as an EPICS driver support package. Then, all that remains to be done is to create a device support module for the EPICS motor record that translates requests from EPICS into a form that MX can understand and then translate MX’s responses back into a form that EPICS can understand. I believe that there is a good chance that this can be completed quickly and hope to have a working version by the end of 2002.

The benefits of implementing such an interface are:

* The EPICS MEDM program is a good static interface builder. This would allow MEDM to be used with MX-controlled motors.
* This would allow Spec to operate MX-controlled motors.
* MX has a fairly large motor driver library which currently contains 34 motor drivers and 17 pseudomotor drivers. These would all then be available more or less immediately to the rest of EPICS.

The initial development work for this will be done on Linux. Once it works well, it will then be ported to Solaris, Win32, and maybe even RTEMS.

## VIII Conclusions

In summary, much progress has been made in the last couple of years. MX is now a full featured beamline control system for materials science, protein crystallography, and other experiments. It is expected that many new features will be added as needed to support better beamline automation and new experimental techniques.

The author would like to thank MR-CAT, IMCA-CAT, and SER-CAT for the support they have provided during the development of this software. MR-CAT related development was supported in part by funding from the U. S. Department of Energy under grant number DEFG0200ER45811. IMCA-CAT related development was supported by the companies of the Industrial Macromolecular Crystallography Association through a contract with Illinois Institute of Technology (IIT), executed through the IIT’s Center for Synchrotron Radiation Research and Instrumentation. SER-CAT related development was supported by the institutes of the South East Regional Collaborative Access Team through a contract with the University of Georgia. Use of the Advanced Photon Source was supported by the U. S. Department of Energy, Basic Energy Sciences, Office of Science, under Contract No. W-31-109-Eng-38.


format: text

text:
# Analytical solution of 1D lattice gas model with infinite number of multiatom interactions

## I Introduction

In recent years there has been growing interest in studies of monatomic chains obtained in the processes of heteroepitaxial growth Mo et al. (1994); Yang et al. (2001); Segovia et al. (1999) at various substrates. A major goal of these studies is to obtain the so-called quantum wires for potential use in microelectronic applications. In these experiments one usually aims at obtaining wires of infinite length for subsequent use in the studies of the Luttinger liquid Segovia et al. (1999). In practice, however, one frequently encounters the monatomic wires of finite length. These may be useful in more practical applications, such as the microelectronic circuitry Orlov et al. (1996). Such finite 1D clusters were observed, e. g., in Refs. Gambardella et al. (2000a, b). Current theoretical interest is in magnetic properties of finite monatomic clusters Dorantes-Dávila and Pastor (1998) in view of their potential use in magnetic memory devices. The latter application, however, will require the development of the techniques of mass production of such objects. In this context it would be interesting to study the possibility of their self-assembly and self-organization similar to analogous processes in 2D heteroepitaxial systems Eaglesham and Cerullo (1990); Notzel et al. (1994).

In the present paper we consider an analytic solution for the cluster size distribution in the framework of a 1D lattice gas model with an arbitrary number of cluster interactions within contiguous atomic chains plus the next nearest neighbor pair interaction. Such model can be justified in the framework of the Frenkel-Kontorova model of strained epitaxy Tokar and Dreyssé (2002) but also can be useful in other cases, e. g., in ab initio approaches, where the cluster interactions appear because of the many-body nature of the electron interactions which cannot be reduced to the pair interatomic potential.

## II The model

Let us consider a lattice gas model defined by the function $`E(l)`$ describing the dependence of the energy of 1D atomic clusters (or chains) on their length $`l`$. Such a model was already used for the description of the self-assembly phenomena in Ref. Priester and Lannoo (1995). We note, however, that in 2D this approach can be only phenomenological because the cluster size does not characterize a 2D structure uniquely. In contrast, in 1D the cluster is defined unambiguously by its length and the function $`E(l)`$ is also uniquely defined provided the clusters do not interact.

To apply standard tools of statistical physics we need to find the Hamiltonian corresponding to our model. To this end we consider the following trial expression

$$H=\underset{i;l=2}{\overset{\mathrm{}}{}}V_ln_in_{i+1}\mathrm{}n_{i+(l1)},$$
(1)
where $`V_l`$ are numerical coefficients to be fitted to reproduce the chain energies $`E(l)`$ and $`n_i`$ the occupation numbers taking values 0 or 1 depending on whether site $`i`$ is empty or is occupied by an atom. Assuming we have already fitted $`V_{l1},V_{l2},\mathrm{},V_2`$ to the energies $`E(l1)E_{l1},E_{l2},\mathrm{},E_2`$ let us consider the chain of length $`l`$. From Eq. (1) one finds

$$E(l)=E(l1)+V_l+V_{l1}+V_{l2}+\mathrm{}+V_2,$$
(2)
where the $`V`$ terms on the right hand side (r. h. s.) account for the interactions of the newly added atom with the rest. Similarly

$$E(l1)=E(l2)+V_{l1}+V_{l2}+\mathrm{}+V_2.$$
(3)
Subtracting Eq. (3) from Eq. (2) we get

$$E(l)E(l1)=E(l1)+V_lE(l2)$$
(4)
from where it follows that the Hamiltonian of our model is Eq. (1) with the coefficients given by the recursion relation

$$V_l=E_l2E_{l1}+E_{l2}$$
(5)
initialized by $`E_0=E_1=0`$. For the system to be well defined in the thermodynamic limit the chain energy $`E_l`$ cannot grow quicker then linearly when $`l\mathrm{}`$. This means that the the cluster interactions $`V_l`$ tend to zero at large $`l`$ because according to Eq. (5) they are equal to the discrete second derivative of $`E_l`$ with respect to $`l`$.

In the above model atoms interact only when they belong to the same contiguous chain, so that separate chains are not coupled to each other. To make the model more realistic we allow for the interchain coupling by adding to the Hamiltonian (1) the next nearest neighbor (NNN) pair interaction term

$$H=V_{NNN}\underset{i}{}n_in_{i+2}+\underset{i;l=2}{\overset{\mathrm{}}{}}V_ln_in_{i+1}\mathrm{}n_{i+(l1)}.$$
(6)
## III The canonical ensemble solution

We consider a 1D lattice with a fixed number of atoms $`N`$ which can occupy $`I>N`$ lattice sites. Our major goal in this paper is to find the equilibrium distribution of the atomic chain sizes at finite temperature $`T`$ in the model defined by Hamiltonian (6). But first we will consider the above problem using a simpler Hamiltonian (1), i. e., by neglecting the NNN interaction which will be accounted for later. We will seek the canonical ensemble solution of our problem by generalizing the approach of Ref. Vavro (2001). In this approach one first have to compute the energy of the system in terms of cluster variables $`m_l`$—the number of clusters of length $`l`$ and their energy $`E(l)`$. Because the $`l\mathrm{}`$ limit is not completely trivial, we calculate the configurational energy by assuming that all $`V_l`$ with $`l>L`$ are equal to zero. In this case the total energy can be calculated as

$$=\underset{l=1}{\overset{L2}{}}E_lm_l+E_L^{}N_>m_>\underset{l=2}{\overset{L}{}}(l1)V_l,$$
(7)
where $`E_L^{}=_{l=2}^LV_l`$ is the first (discrete) derivative of $`E_l`$, $`m_>`$ is the total number of clusters of length exceeding $`L2`$ which below we will call the big clusters, and $`N_>=\left(N_{l=1}^{L2}lm_l\right)`$ is the number of atoms contained in the big clusters. The last two terms on the r. h. s. of Eq. (7) can be verified by first checking their correctness for $`m_>=1`$ and then noting that cutting a big cluster into two ones amounts to cutting $`l1`$ couplings of type $`V_l`$. To complete the calculation of the total free energy $`F_{tot}=E_{tot}TS`$ it remains to compute the entropy $`S`$. The number of atomic configurations can be expressed through the cluster variables as

$$\mathrm{\Omega }_{at}=\frac{m_{tot}!}{m_1!m_2!\mathrm{}m_>!}\left(\begin{array}{c}N_>(L2)m_>1\\ m_>1\end{array}\right)$$
(8)
In this equality the first multiplier on the r. h. s. counts all inequivalent permutations between the clusters with all big clusters being considered as equivalent while the second factor accounts for the number of ways to divide $`N_>`$ atoms into $`m_>`$ big clusters including their permutations. We find this factor by observing that the division of $`N_>`$ atoms into clusters of sizes exceeding $`L2`$ is equivalent to dividing $`n=N_>(L2)m_>`$ atoms into $`k=m_>`$ clusters. The latter quantity was calculated in Ref. Vavro (2001) to be equal to the binomial coefficient

$$C_{k1}^{n1}=\left(\begin{array}{c}n1\\ k1\end{array}\right)$$

provided $`nk1`$. Otherwise it is equal to zero. The meaning of this formula is simple. In a contiguous chain of $`n`$ atoms there is $`n1`$ places to cut the chain into pieces. $`k`$ pieces can be obtained with $`k1`$ cuts. Hence, the number of the possible cuts of $`n`$ atoms into $`k`$ pieces is equal exactly to the above combinatorial coefficient. To complete the calculation of the total number of configurations we have to multiply $`\mathrm{\Omega }_{at}`$ by the corresponding factor responsible for the configurations of vacant sites Vavro (2001). But before doing this we repair Eq. (7) by accounting for the omitted NNN pair interaction. To this end we note that as long as the contiguous chains are concerned, $`V_{NNN}`$ interaction amounts simply to renormalization of $`V_3`$. Additionally this term introduces the interchain coupling between the chains separated by a single vacant site (a vacancy). Therefore, the term $`V_{NNN}k_1`$, where $`k_1`$ the number of vacancies, should be added to the r. h. s. of Eq. (7) (here and above we do not pay attention to the boundary conditions Vavro (2001) because in this paper we are interested only in the thermodynamic limit). So when counting the configurations of empty sites we should separate the configurations with different numbers of vacancies. This is achieved by using Eq. (8) with the total number of empty sites $`IN`$ instead of $`N`$, the big clusters in this case should exceed $`L2=1`$. Denoting the total number of clusters as $`k`$ from Eq. (7) we get

$$\mathrm{\Omega }_{vac}=\frac{k!}{k_1!(kk_1)!}\left(\begin{array}{c}ILk1\\ kk_11\end{array}\right).$$
(9)
Our final expression for $`S`$ is obtained from Eqs. (8) and (9) with the use of the Stirling formula as

$`{\displaystyle \frac{S}{k_BI}}=2c\mathrm{ln}c{\displaystyle \underset{l=1}{\overset{L2}{}}}c_l\mathrm{ln}c_lv\mathrm{ln}v2(cv)\mathrm{ln}(cv)`$
$`+(1\theta c)\mathrm{ln}(1\theta c)(1\theta 2c+v)\mathrm{ln}(1\theta 2c+v)`$ (10)
$`2c_>\mathrm{ln}c_>+[\theta _>(L2)c_>]\mathrm{ln}[\theta _>(L2)c_>][\theta _>(L1)c_>]\mathrm{ln}[\theta _>(L1)c_>],`$
where $`\theta =N/I`$ is the total coverage, $`c_l=m_l/I`$ is the concentration of clusters of length $`l`$, $`c=_lc_l`$ the total cluster concentration, $`v=k_1/I`$ the vacancy concentration, $`\theta _>=\theta _{l=1}^{L1}lc_l=N_>/I`$. This expression for reduced entropy together with the reduced energy obtained from Eq. (7) are sufficient to obtain the reduced free energy as a function of the unknown quantities $`\{c_l\}`$, $`c_>,`$ and $`v`$. Minimizing it with respect to these variables we arrive in the limit $`L\mathrm{}`$ at the following set of equations

$`c_l={\displaystyle \frac{c^2(1\theta 2c+v)^2}{(cv)^2(1\theta c)}}\mathrm{exp}\left({\displaystyle \frac{\mu lE_l}{k_BT}}\right)`$ (11)
$`(cv)^2=v(1\theta 2c+v)\mathrm{exp}\left({\displaystyle \frac{V_{NNN}}{k_BT}}\right).`$ (12)
To find the above limit it is necessary to know the asymptotic behavior of $`c_l`$ at large $`l`$. The limit was taken by assuming the exponential behavior $`c_l\mathrm{exp}(\lambda l)`$ which is consistent with the resulting equation (11). We note that Eq. (12) is of the second order in $`v`$, so one can exclude this variable from the first equation to obtain a closed equation for $`\{c_l\}`$. From the above set, however, it is easier to see that if the NNN interaction is negative, than $`vc`$ as $`T0`$, i. e., the vacancy concentration becomes equal to the concentration of clusters which means that the system becomes ordered: the clusters self-organize into chains separated by monovacancies. This is further confirmed by the fact that the entropy Eq. (III) tends to zero at zero temperature provided $`V_{NNN}<0`$ and the clusters are size calibrated. Indeed, in the limit $`L\mathrm{}`$ only the first two lines of Eq. (III) survive

$`{\displaystyle \frac{s}{k_B}}`$ $`=`$ $`2c\mathrm{ln}c{\displaystyle \underset{l=1}{\overset{L2}{}}}c_l\mathrm{ln}c_lv\mathrm{ln}v2(cv)\mathrm{ln}(cv)`$ (13)
$`+(1\theta c)\mathrm{ln}(1\theta c)(1\theta 2c+v)\mathrm{ln}(1\theta 2c+v),`$
where $`s=S/I`$ is the entropy per site.

In connection with this expression it is worth noting that Eqs. (11) and (12) can formally be obtained in the thermodynamic limit from the variation with respect to $`c_l`$ of the expression for the free energy density

$$f=\underset{l=1}{\overset{\mathrm{}}{}}c_lE_lk_BTs\mu (\underset{l=1}{\overset{\mathrm{}}{}}lc_l\theta ),$$

where $`\mu `$ is the Lagrange multiplier. This derivation, however, may cause doubts in the case of clusters of sizes $`l=O(N)`$ when the terms $`lc_l=m_lO(N)/N`$ in the above equation acquire discrete values and so are not suitable for the variational treatment. The derivation presented above is more rigorous and besides can be modified to be applicable to finite systems Vavro (2001); Hodak and Girifalco (2001).

At $`T=0`$ and $`c=v`$ the entropy density (13) takes the form

$$\frac{s}{k_B}=c\mathrm{ln}c\underset{l=1}{\overset{L2}{}}c_l\mathrm{ln}c_l.$$
(14)
Thus, if the clusters at $`T=0`$ are size calibrated, i. e., if for some value of $`l=l_0`$ $`c_{l_0}=c`$ while $`c_l=0`$ for all other $`l`$, then from the above equation it follows that $`s=0`$. For $`V_{NNN}0`$ the entropy is positive even at $`T=0`$ meaning disordered state. The issue of size calibration will be considered in more detail in section V.

## IV The Ising model

To check Eqs. (11) and (12) we will apply them to the known exactly solvable problem—the 1D Ising model. As is known, it is equivalent to the lattice gas model with pair interatomic interaction. This can be easily shown (see Ducastelle (1991) Ch. II) by substituting $`n_i=(1\sigma _i)/2`$ ($`\sigma _i=\pm 1`$—the Ising spin variable) into the LGM Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}V_{ij}n_in_j={\displaystyle \frac{1}{8}}{\displaystyle \underset{ij}{}}V_{ij}(1\sigma _i)(1\sigma _j)`$ (15)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}(V_{ij}/4)\sigma _i\sigma _j+{\displaystyle \underset{j}{}}(V_{ij}/4)(12\sigma _j)`$
Thus, the Ising model with NN interaction can be solved in the canonical ensemble formalism with the use of our formulas. In this section we will obtain the appropriate formulas in closed analytical form. Also, we will compare our results obtained in the canonical ensemble with more familiar results obtained in the grand canonical formalism via the transfer matrix method. The latter solution can be found e. g., in Ch. V of Ref. Ducastelle (1991), where the 1D Ising model Hamiltonian was considered which in our notation can be written as (cf. Eq. (15)

$$H_{IM}=(V_{NN}/4)\underset{i}{}\sigma _i\sigma _{i+1}\frac{1}{2}h\underset{i}{}(\sigma _i+\sigma _{i+1})$$
(16)
was considered. The reduced free energy calculated in the grand canonical ensemble was found to be

$$\frac{F_{GC}}{Nk_BT}=\mathrm{ln}\left(e^{K/4}\mathrm{cosh}\overline{h}+\sqrt{e^{K/2}\mathrm{sinh}^2\overline{h}+e^{K/2}}\right),$$
(17)
where $`K=V_{NN}/k_BT`$ and $`\overline{h}=h/k_BT`$ is the external field (divided by $`k_BT`$) which fixes the magnetization $`M`$. The latter is connected to the coverage $`\theta `$ by the relation

$$\theta =(1M)/2$$
(18)
and can be calculated as

$$M=F_{GC}/\overline{h}.$$
(19)
To obtain the canonical ensemble free energy we have to perform the Legendre transform (see Eq. (16))

$$\frac{F_C}{Nk_BT}=F_{GC}+\overline{h}\frac{F_{GC}}{\overline{h}},$$
(20)
where the partial derivative is calculated from Eq. (19) as

$$M=\frac{e^{K/4}\mathrm{sinh}\overline{h}(1+e^{K/4}\mathrm{cosh}\overline{h}/\sqrt{e^{K/2}\mathrm{sinh}^2\overline{h}+e^{K/2}})}{e^{K/4}\mathrm{cosh}\overline{h}+\sqrt{e^{K/2}\mathrm{sinh}^2\overline{h}+e^{K/2}}}.$$
(21)
Finally, to obtain the LGM free energy $`F`$ in the canonical ensemble formalism we have to add the statistical average of the last term on the r. h. s. of Eq. (15):

$$\frac{F}{Nk_BT}=\frac{F_C}{Nk_BT}+(K/4)(12M).$$
(22)
In our formalism the free energy $`F`$ for the NN LGM can be calculated as follows. According to Eq. (12), in the case $`V_{NNN}=0`$

$$v=c^2/(1\theta ).$$
(23)
Substituting this into Eq. (11) we get in the case of the Ising model

$$c_l=(1\theta c)\mathrm{exp}[\overline{\mu }+(\overline{\mu }K)(l1)],$$
(24)
where $`\overline{\mu }=\mu /k_BT`$. From Eq. (24) we easily find

$$\{\begin{array}{ccc}\hfill c& =& \underset{l=1}{\overset{\mathrm{}}{}}c_l=(1\theta c)e^{\overline{\mu }}/(1e^{\overline{\mu }K})\hfill \\ \hfill \theta & =& \underset{l=1}{\overset{\mathrm{}}{}}lc_l=c+c^2e^K/(1e^{\overline{\mu }K}).\hfill \end{array}$$

These two equations are sufficient to express the unknown quantities $`c`$ and $`\overline{\mu }`$ through the coverage $`\theta `$ and the interaction parameter $`K`$—the independent variables in the canonical ensemble. From Eq. (IV) we have after some algebra

$$\{\begin{array}{ccc}\hfill \overline{\mu }& =& \mathrm{ln}\left[\frac{1}{2\theta }x+\sqrt{\left(\frac{1}{2\theta }x\right)^2+x(1x)}\right]\hfill \\ \hfill c& =& \frac{1}{2x}\left[1\sqrt{14\theta (1\theta )x}\right],\hfill \end{array}$$

where $`x=(1e^K)`$. These expressions can be used to obtain analytic expressions for all quantities of interest.

In the canonical ensemble case the free energy is

$$F=k_BT\left[\overline{\mu }\theta (1\theta )\mathrm{ln}\left(1+\frac{e^{\overline{\mu }}}{1e^{\overline{\mu }K}}\right)\right].$$
(25)
As we see, it is quite different from Eq. (22). We were unable to compare the two solutions analytically and used the numerical procedure. By fixing $`\overline{h}`$ and $`K`$ we computed $`M`$ with the use of Eq. (21) and $`F/Nk_BT`$ from Eq. (22). Then $`\theta `$ was calculated according to Eq. (18) and $`F/Nk_BT`$ was computed from Eq. (25). Both values coincided within 13-15 significant digits depending on how singular the values of $`F/Nk_BT`$ are at the chosen values of $`\overline{h}`$ and $`K`$.

Thus we have solved the 1D NN LGM in the canonical ensemble approach and expressed all quantities of interest in terms of the atomic density (or coverage) $`\theta `$ which is more natural and easier to measure quantity then the chemical potential variable of the grand canonical formalism.

## V The size calibration

According to Ref. Priester and Lannoo (1995), the system exhibits the size calibration of the atomic clusters if the energy density per atom $`E(l)/l`$ has a minimum. The size calibration is easiest to see in the case $`V_{NNN}=0`$ when the distribution Eq. (11) takes the form

$$c_l=(1\theta c)\mathrm{exp}\beta (\mu lE_l).$$

Except for the excluded volume prefactor this expression coincides with the formula proposed in Ref. Priester and Lannoo (1995) where a detailed discussion of the size calibration issue is given and which fully applies to our case provided the pair attraction $`V_{NN}`$ is not too strong. In Ref. Tokar and Dreyssé (2002) we have shown that in the framework of the Frenkel-Kontorova model the effective chain energy (in fact, the free energy—see Ref. Tokar and Dreyssé (2002) for details) has the form

$$E_l=V_{NN}(l1)\underset{\text{chains}}{}(W_lTS_l)$$
(26)
where the last two terms represent the relaxation free energy. The length dependence of the relaxation energy $`W_l`$ is governed by the dimensionless parameter $`\alpha =k_s/k_p`$, where the spring constants $`k_s`$ and $`k_p`$ are the second derivatives of the potentials which bind the atom to the substrate ($`k_s`$) and of the pair interatomic interaction ($`k_p`$), respectively. In the illustrative calculations below the energy unit was chosen to be equal to the energy of the unrelaxed misfit strain $`k_pf^2`$, where the interatomic spring constant $`k_p`$ was defined earlier and $`f`$ is the misfit as defined in Ref. Tokar and Dreyssé (2002). The relaxation entropy $`S_l`$ is practically linear in $`l`$ and so according to Eq. (5) essentially contributes only into the NN interaction.

In our calculations we chosen $`\alpha =10^5`$ to be quite small in order to visualize all qualitative details of the self-assembly behavior. Physically this would correspond to very weak binding to the substrate. In the above energy units $`V_{NN}`$ was chosen to be equal to -0.25 so that there was a minimum in $`E(l)/l`$ at $`l100`$ (see Fig.1). The calculation of the size distribution for this case is shown in Fig. 2. It qualitatively agrees with the Monte Carlo simulations of Ref. Meixner et al. (2001) (cf. Fig. 2 of that reference). The most notable feature of the above calculations is the change of the position of the maximum of the cluster size distribution with lowering temperature—the feature not clearly understood from the theory of Ref. Priester and Lannoo (1995).

Although the entropic contribution into $`E_l`$ is rather small due to the small temperatures considered, the shallowness of the minimum in $`E(l)/l`$ makes it possible that the above shift of the maximum is due to the shift in the position of the minimum in $`E_l=W_lTS_l`$ with lowering temperature. To check this possibility we repeated calculations with $`E_l`$ replaced by temperature independent part $`W_l`$. The results of this calculation is shown in Fig. 3. The only qualitative difference with Fig. 2 is that the bimodal distribution is never observed in the case of $`T`$-independent case. Presumably, this is due to the 1D character of our model, because in 2D case where also $`T`$-independent interactions were used we found a strongly bimodal distribution Tokar and Dreyssé (2002).

## VI Conclusion

In this paper we considered an analytical solution of a 1D lattice gas model with an arbitrary number of cluster interatomic interactions within contiguous atomic chains. We illustrated the model by applying it to the description of the equilibrium distribution of atomic clusters in strained epitaxy in the thermodynamic limit. However, the model can be used also in other cases when the multiatom interactions are of the above type, as well as applied to finite systems, like those discussed in Refs. Vavro (2001); Hodak and Girifalco (2001). Our analytical solution confirms the observation made in Ref. Meixner et al. (2001) in the framework of a phenomenological approach that the maximum in the equilibrium distribution of the cluster sizes moves to higher values with lowering temperature. This can lead to an interesting kinetic phenomenon consisting in the saturated ripening when the system is quenched to lower temperature.

###### Acknowledgements.

One of us (VT) expresses his gratitude to University Louis Pasteur de Strasbourg and IPCMS for their hospitality.


text:
# Pressure-induced structural transitions in MgH2

## Abstract

The stability of MgH2 has been studied up to 20 GPa using density-functional total-energy calculations. At ambient pressure $`\alpha `$-MgH2 takes a TiO2-rutile-type structure. $`\alpha `$-MgH2 is predicted to transform into $`\gamma `$-MgH2 at 0.39 GPa. The calculated structural data for $`\alpha `$\- and $`\gamma `$-MgH2 are in very good agreement with experimental values. At equilibrium the energy difference between these modifications is very small, and as a result both phases coexist in a certain volume and pressure field. Above 3.84 GPa $`\gamma `$-MgH2 transforms into $`\beta `$-MgH2; consistent with experimental findings. Two further transformations have been identified at still higher pressure: i) $`\beta `$\- to $`\delta `$-MgH2 at 6.73 GPa and (ii) $`\delta `$\- to $`ϵ`$-MgH2 at 10.26 GPa.

preprint: vajeeston et al

The utilization of high-pressure technology, which can reach the range of giga Pascal (GPa), has made considerable progress in experimental studies of hydrogen storage materials. It has been demonstrated for a number of metal$``$hydrogen systems that application of high-pressure is an effective tool to produce vacancies in the host metallic matrix, which in turn leads to various novel properties fukai94 . The Mg-based alloys possess many advantageous functional properties, such as heat resistance, vibration absorption, recycling ability, etc. orimo01 . In recent years, therefore, much attention has been paid to investigations on specific material properties of Mg alloys. Magnesium is an attractive material for hydrogen-storage applications because of its light weight, low manufacture cost and high hydrogen-storage capacity (7.66 wt. %). On the other hand, owing to its high operation temperature (pressure plateau of 1 bar at 525 K) and slow absorption kinetics, practical applications of magnesium-based alloys have been limited. However, it has recently been established that improved hydrogen absorption kinetics can be achieved by reduced particle size and/or addition of transition metals to magnesium and magnesium hydrides chen95 . $`\gamma `$-MgH2 often occurs as a by product in high-pressure syntheses of technologically important metal hydrides like Mg2NiH4. Hence, a complete characterization of $`\gamma `$-MgH2, in particular, knowledge about its stability at high pressure is highly desirable. Present theoretical investigation assumes importance as the high-pressure x-ray diffraction and neutron diffraction studies are unable to identify the exact position of hydrogen atoms owing to its very-low scattering cross-section. This contribution represents the first theoretical report on the high-pressure properties of MgH2 as evident from ab initio density functional calculations.

$`\alpha `$-MgH2 crystallizes with TiO2-r-type (r = rutile) structure at ambient pressure and low temperature bastide80 ; zachariasen63 . At higher temperatures and pressures tetragonal $`\alpha `$-MgH2 transforms into orthorhombic $`\gamma `$-MgH2. Recently Bortz et al. bortz99  solved the crystal structure of $`\gamma `$-MgH2 ($`\alpha `$-PbO2 type) on the basis of powder neutron diffraction data collected at 2 GPa. The $`\alpha `$\- to $`\gamma `$-MgH2 transition pressure is not yet known. Owing to the difficulties challenged in establishing hydrogen position in a metal matrix by x-ray diffraction high-pressure information on MgH2 is scarce. The aim of the present investigation is to remedy this situation by examining MgH2 theoretically at high pressure in eleven closely related structural configurations.

We used the ab initio generalized gradient approximation (GGA) perdew  to obtain accurate exchange and correlation energies for a particular structural arrangement. The structures are fully relaxed for all the volumes considered in the present calculations using force as well as stress minimization. Experimentally established structural data were used as input for the calculations when available. The full-potential linear muffin-tin orbital (FP-LMTO)wills  and the projected augmented wave (PAW) implementation of Vienna ab initio simulation package (VASP) vasp1  were used for the total-energy calculations to establish phase stability and transition pressures. Both methods yielded nearly the same result, e.g., the transition pressure from $`\alpha `$\- to $`\gamma `$-MgH2 came out as 0.385 and 0.387 GPa from FP-LMTO and VASP, respectively. Similarly almost identical values were obtained for ground-state properties like bulk modulus and equilibrium volume. Hence, for rest of the calculations we used only the VASP code because of its computational efficiency. In order to avoid ambiguities regarding the free-energy results we have always used the same energy cutoff and a similar k-grid density for convergence. In all calculations 500 k points in the whole Brillouin zone were used for $`\alpha `$-MgH2 and a similar density of k points for the other structural arrangements. The PAW pseudo potentials paw  were used for all our VASP calculations. A criterion of at least 0.01 meV/atom was placed on the self-consistent convergence of the total energy, and the calculations reported here used a plane wave cutoff of 400 eV.

In addition to the experimentally identified $`\alpha `$\- and $`\gamma `$-modifications of MgH2, we have carried out calculations for several other possible lowther99  types of structural arrangements for MgH2 (see Table 1 and Figs. 1 and 2 ). The calculated total energy vs. volume relation for all these alternatives are shown in the Fig. 1. In order to get a clear picture about the structural transition points, in Fig. 2 we have displayed the Gibbs free energy difference between the pertinent crystallographic structures of MgH2 with reference to $`\alpha `$-MgH2 as a function of pressure. The equilibrium volumes, (30.64 and 30.14 Å3/f.u. for $`\alpha `$\- and $`\gamma `$-MgH2, respectively) are within 1 % of the experimental values indicating that the theoretical calculations are reliable. The calculated positional parameters (Table 1) are also in excellent agreement with the experimental data. The calculated values for the bulk modulus ($`B_0`$, see Table 1) vary between 44 to 65 GPa for the various structural arrangements. Among the eleven possibilities considered the Ag2Te-type leads to the smallest $`B_0`$ value and the AlAu2-type to the highest (Table 1).

The calculated transition pressure for the $`\alpha `$\- to $`\gamma `$-MgH2 conversion is 0.387 GPa (Fig. 3) and as the free energy of the two modification is nearly the same at the equilibrium volume (Fig. 1), it is only natural that these phases coexist in a certain volume range bortz99 . Structurally $`\alpha `$\- and $`\gamma `$-MgH2 are also closely related, both comprising Mg in an octahedral coordination of 6 H which in turn are linked by edge sharing in one direction and by corner sharing in the two other directions. The chains are linear in tetragonal $`\alpha `$-MgH2 and run along the four-fold axis of its TiO2-r-type structure, while they are zig-zag shaped in $`\gamma `$-MgH2 and run along a two-fold screw axis of its orthorhombic $`\alpha `$-PbO2-type structure. The octahedra in $`\gamma `$-MgH2 are strongly distorted. The pressure induced $`\alpha `$-to-$`\gamma `$ transition involves reconstructive (viz. bonds are broken and re-established) rearrangements of the cation and anion sublattices. The occurrence of a similar phase transition in PbO2 suggests that thermal activation or enhanced shear is of importance haines96 .

The subsequent phase transition from $`\gamma `$\- to $`\beta `$-MgH2 occur at 3.84 GPa. Bortz et al. bortz99  found no evidence for the formation of such a $`\beta `$ modification up to 2 GPa, whereas Bastide et al. bastide80  found that at higher pressure (4 GPa) and temperature (923 K) there occurs a new phase (viz. $`\beta `$ in accordance with our findings). However, the H positions are not yet determined experimentally and hence our optimized structural parameters for $`\beta `$-MgH2 should be of particular interest. The volume shrinkage at the transition point is 1.45 Å3/f.u. On further increase of the pressure to around 6.7 GPa, $`\beta `$-MgH2 is predicted to transform into $`\delta `$-MgH2 (orthorhombic $`Pbc2_1`$) with a volume shrinkage of 1.1 Å3/f.u. In the pressure range 6.7$``$10.2 GPa the structural arrangements $`\beta `$-MgH2, AuSn2 and $`\delta `$-MgH2 are seen to lie within a narrow energy range of some 10 meV (Fig. 1). This closeness in energy suggests that the relative appearance of these modification will be quite sensitive to, and easily affected by, other external factors like temperature and remnant lattice stresses. A transformation to $`ϵ`$-MgH2 (AlAu2-type structure also called cotunnite-type) is clearly evident from Fig. 1 with a volume change of 1.7 Å3/f.u. at the transition point (see also Figs. 2 and 3). It is interesting to note that similar structural transition sequences are reported for transition-metal oxides like HfO2 lowther99 . Recently, the cotunnite-type structure of TiO2 (synthesized at pressure above 60 GPa and at high temperatures) has been shown dubrovinsky01  to exhibit an extremely high bulk modulus (431 GPa) and hardness (38 GPa). Our present study predicts that, it is possible to stabilize MgH2 in the AlAu2-type structure as a soft material ($`B_0`$ = 65 GPa) above 10.26 GPa. Hence, the present prediction of stabilization of a contunnite-type structure for MgH2 at high pressure should be of considerable interest.

The calculated unit-cell volume difference between the involved phases at their equilibrium volume relative to the $`\alpha `$-phase is 0.497 (for $`\gamma `$), 2.94 (for $`\beta `$), 3.60 (for $`\delta `$), and 5.96 (for $`ϵ`$) Å3/f.u., which is approximately 1.6, 9.5, 11.8, and 19.5 %, respectively, smaller than the equilibrium volume of the $`\alpha `$-phase. The energy difference between $`\alpha `$-phase and the $`\gamma `$, $`\beta `$, $`\delta `$, and $`ϵ`$ phases at their equilibrium volume is 0.81, 43, 66, and 197 meV/f.u., respectively, which is considerable smaller than for similar types of phase transitions in transition metal dioxides lowther99 . The known stabilization of the high-pressure phase of TiO2 by using high-pressure$``$high-temperature synthesis dubrovinsky01  indicates that it may be possible to stabilize the high-pressure phases of MgH2 due to their small energy differences and (partly) reconstructive phase transitions. Such stabilization of MgH2 high-pressure phases at room temperature would reduce the volume considerably, and in the extreme case the expected volume reduction would be ca. 19.5 % for $`ϵ`$-MgH2 compared with $`\alpha `$-MgH2. This would imply an increased volumetric storage capacity of hydrogen for MgH2. Of interest would be to explore the possibility of stabilizing the high-pressure phases by chemical means. At the same time this may possibly open up for improved kinetics with respect to reversible hydrogen absorption/desorption. Hence, the observation of high-pressure phases in MgH2 may have technological significance.

The calculated energy gap (from the top of the valence band (VB) to bottom of the conduction band) is 4.2 and 4.3 eV for $`\alpha `$\- and $`\gamma `$-MgH2, respectively. According to Fig. 4 the width of VB is 7.5 and 8.3 eV in $`\alpha `$\- and $`\gamma `$-MgH2, respectively. The increased width of VB in $`\gamma `$-MgH2 is due to the reduction in the Mg$``$H separation (1.93 Å in $`\alpha `$-MgH2 vs. 1.91 Å in $`\gamma `$-MgH2). An experimental UV absorption study prommer94  of $`\alpha `$-MgH2 gave an absorption edge at 5.16 eV, whereas earlier theoretical electronic structure investigation yu88 ; prommer94  reported a band gap around 3.4 eV. Hence, our calculated energy gap value is in better agreement with the experimentally measured value than the earlier values. In fact, the present underestimation of the band gap by about 17 % is typical for the accuracy obtained by first principle methods for semiconductors and insulators. Such distinctions probably originates from the use of GGA-based exchange correlation functionals.

In conclusion, the calculated structural parameters for $`\alpha `$\- and $`\gamma `$-MgH2 are in excellent agreement with the experimental findings. We found that the ground state of $`\alpha `$-MgH2 becomes unstable at higher pressure and the calculated transition pressures for $`\alpha `$-to-$`\gamma `$\- and $`\gamma `$-to-$`\beta `$-MgH2 are found to be in very good agreement with the experimental findings. We predict that further compression of $`\beta `$-MgH2 will lead to a phase transition to $`\delta `$-MgH2 (orthorhombic; $`Pbc2_1`$) at 6.73 GPa and that the AlAu2-type structure of a hypothetical $`ϵ`$ modification stabilizes above 10.26 GPa. Our total energy study suggested that the formation of the high-pressure $`ϵ`$-phase could also be possible at the ambient pressure/temperature by choosing an appropriate preparation method. This may reduce the volume requirement by ca. 19.5 % compared with $`\alpha `$-MgH2. In agreement with experiments, $`\alpha `$-MgH2 turns out to be an insulator, and the theoretical treatment shows that the high-pressure modifications also should exhibit insulating behavior.

The authors greatfully acknowledge the Research Council of Norway for financial support and for the computer time at the Norwegian supercomputer facilities and R. Vidya for critical reading of this manuscript.


text:
# SU(3) Predictions of 𝐵→𝑃𝑃 Decays in the Standard Model

## 1 SU(3) Parameters for $`BPP`$

In Standard Model (SM) the decay amplitudes for $`B`$ meson to two pseudoscalar SU(3) octet mesons, $`BPP`$, have both tree and penguin contributions and can be written as

$`A(BPP)=<PP|H_{eff}^q|B>`$
$`={\displaystyle \frac{G_F}{\sqrt{2}}}[V_{ub}V_{uq}^{}T(q)+V_{tb}V_{tq}^{}P(q)].`$ (1)
As far as the SU(3) structure is concerned, the effective Hamiltonian contains $`\overline{3}`$, $`6`$, and $`\overline{15}`$ representations of SU(3) which define several invariant amplitudes,

$`T=A_{\overline{3}}^TB_i(\overline{3})^i(M_l^kM_k^l)+C_{\overline{3}}^TB_iM_k^iM_j^k(\overline{3})^j`$
$`+A_6^TB_i(6)_k^{ij}M_j^lM_l^k+C_6^TB_iM_j^i(6)_l^{jk}M_k^l`$
$`+A_{\overline{15}}^TB_i(\overline{15})_k^{ij}M_j^lM_l^k+C_{\overline{15}}^TB_iM_j^i(\overline{15})_l^{jk}M_k^l,`$
where $`B_i=(B_u,B_d,B_s)`$ and $`M`$ is the SU(3) pseudoscalar octet. $`C_6A_6`$ always appear together. We use $`C_6`$ only.

In general there are both tree and penguin amplitudes $`C(A)_{\overline{3},6,\overline{15}}^{T,P}`$. In the SM the amplitudes $`C(A)_{6,\overline{15}}^T`$ and $`C(A)_{6,\overline{15}}^P`$ are related by the Wilson coefficients $`c_i`$ of the relevant operators,

$`{\displaystyle \frac{C_6^P}{C_6^T}}`$ $`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{c_9^{tc}c_{10}^{tc}}{c_1c_23(c_9^{uc}c_{10}^{uc})/2}},`$
$`{\displaystyle \frac{C_{\overline{15}}^P(A_{\overline{15}}^P)}{C_{\overline{15}}^T(A_{\overline{15}}^T)}}`$ $`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{c_9^{tc}+c_{10}^{tc}}{c_1+c_23(c_9^{uc}+c_{10}^{uc})/2}}.`$ (2)
Without loss of generality, one can set $`C_{\overline{3}}^P`$ to be real. One only needs to use 13 real independent parameters to describe $`BPP`$ in the SM which we choose to be

$`C_{\overline{3}}^P,C_{\overline{3}}^Te^{i\delta _{\overline{3}}},C_6^Te^{i\delta _6},C_{\overline{15}}^Te^{i\delta _{\overline{15}}},`$
$`A_{\overline{3}}^Te^{i\delta _{A_{\overline{3}}^T}},A_{\overline{3}}^Pe^{i\delta _{A_{\overline{3}}^P}},A_{\overline{15}}^Te^{i\delta _{A_{\overline{15}}^T}}.`$ (3)
The phases are defined such that all $`C_i^{T,P}`$ are real positive numbers. The amplitudes $`A_i`$ correspond to annihilation contributions and are expected to be small. If they are neglected, there are only 7 independent hadronic parameters.

Using available data on, $`BK\pi ,\pi \pi ,KK`$, one can obtain information about the hadronic parameters and predict other decay branching ratios and CP violating asymmetries. We carried out a $`\chi ^2`$ analysis assuming SU(3) symmetry with the $`BPP`$ data shown in Table 1. In our numerical analysis the KM matrix elements $`V_{us}=\lambda `$, $`V_{cb}=A\lambda ^2`$, $`V_{ub}=|V_{ub}|exp(i\gamma )`$ fixed by $`\lambda =0.2196`$, $`A=0.835`$ and $`|V_{ub}|=0.09|V_{cb}|`$ and $`\gamma =59^{}`$ as determined from other data. We study both the cases with and without annihilation contributions and compare the results.

## 2 Without Annihilation Contributions

For the case without annihilation contributions, there are only 7 hadronic parameters. The results for the best fit values and their 1$`\sigma `$ ranges are shown in Table 2.

Using the above determined hadronic parameters, one can easily obtain the branching ratios and CP asymmetries for other $`BPP`$. We used the following definition for the CP violating rate asymmetry,

$`A_{PP}^{B_i}={\displaystyle \frac{\mathrm{\Gamma }(B_iPP)\mathrm{\Gamma }(\overline{B}_i\overline{P}\overline{P})}{\mathrm{\Gamma }(B_iPP)+\mathrm{\Gamma }(\overline{B}_i\overline{P}\overline{P})}}.`$ (4)
In general $`P`$ can be any one of the SU(3) pseudoscalar octet mesons, $`\pi `$, $`K`$ and $`\eta _8`$. Here we will limit our study to $`P=\pi ,K`$ to avoid complications associated with $`\eta _1`$ and $`\eta _8`$ mixings. In this case there are total 16 decay modes. Among them the decay amplitudes for $`B_dK^{}K^+`$, $`B_s\pi ^{}\pi ^+,\pi ^0\pi ^0`$ only receive annihilation contributions. Since we have neglected annihilation contributions they would have vanishing branching ratios. At present none of them have been measured experimentally. The present bound on $`B_dK^{}K^+`$ is consistent with this prediction.

The predictions for the branching ratios of $`B_sK^+\pi ^{},K^0\pi ^0,K^{}K^+,K^0\overline{K}^0`$ decays are shown in Table 3. These decay modes are predicted to be large and can be measured at hadron colliders. The standard model and SU(3) flavor symmetry can be tested.

The CP asymmetries for some of the decays are shown in Table 4. In the SU(3) limit there are some relations between rate differences defined as, $`\mathrm{\Delta }_{PP}^{B_i}=\mathrm{\Gamma }(B_iPP)\mathrm{\Gamma }(\overline{B}_i\overline{P}\overline{P})`$, between $`\mathrm{\Delta }S=0`$ and $`\mathrm{\Delta }S=1`$ modes due to a unique feature of the SM in the KM matrix element that $`Im(V_{ub}V_{ud}^{}V_{tb}^{}V_{td})=Im(V_{ub}V_{us}^{}V_{tb}^{}V_{ts})`$. One has the following relations,

$`\mathrm{\Delta }_{\pi ^+\pi ^{}}^{B_d}=\mathrm{\Delta }_{K^+\pi ^{}}^{B_s}=\mathrm{\Delta }_{\pi ^+K^{}}^{B_d},`$
$`\mathrm{\Delta }_{\pi ^0\pi ^0}^{B_d}=\mathrm{\Delta }_{K^{}\pi ^0}^{B_s}=\mathrm{\Delta }_{\pi ^0\overline{K}^0}^{B_d}.`$ (5)
As can be seen from Table 4 that the best fit values for $`A_{PP}^{B_i}`$ can be large with several of them to be more than 30%, such as the asymmetries for $`B_sK^+\pi ^{},K^0\pi ^0`$, and $`B_d\pi ^0\pi ^0,\pi ^+\pi ^{}`$. The size of $`A_{PP}^{B_i}`$ for these modes are large, but can be easily understood from the fact that they all have relatively small branching ratios. Using the above relations, one would obtain

$`A_{\pi ^+\pi ^{}}^{B_d}=A_{K^+\pi ^{}}^{B_s}`$
$`=A_{\pi ^+K^{}}^{B_d}{\displaystyle \frac{Br(B_d\pi ^+K^{})}{Br(B_d\pi ^+\pi ^{})}},`$
$`A_{\pi ^0\pi ^0}^{B_d}=A_{K^0\pi ^0}^{B_s}`$
$`=A_{\pi ^0\overline{K}^0}^{B_d}{\displaystyle \frac{Br(B_d\pi ^0\overline{K}^0)}{Br(B_d\pi ^0\pi ^0)}}.`$ (6)
In all the above cases the ratios of the branching ratios are larger than one, a small $`A_{PP}^{B_i}`$ of the decay mode on the right hand side can induce a large $`A_{PP}^{B_i}`$ for the decay modes on the left hand side. These predictions can provide interesting tests for the SM.

## 3 With Annihilation Contributions

In the analyses of the previous sections we have neglected annihilation contributions to $`BPP`$ decays. In this section we study the effects of the annihilation terms on $`BPP`$ decays. In this case we would have total 13 parameters. From Table 1 we see that there are 15 experimental data points. In principle, the 13 hadronic parameters under consideration can be determined. In Tables 5, 6 and 7 we show the results on the hadronic parameters, and some of the $`BPP`$ branching ratios and CP asymmetries.

From Table 5 we see that the size of the best fit annihilation parameters $`A_i`$ are small compared the non-annihilation terms $`C_{3,\overline{15}}`$. This confirms the expectation that annihilation contributions are small. The allowed ranges are large and therefore can not rule out the possibility of having significant annihilation contributions. We have to wait improved experiments to obtain more precise information.

The branching ratios for $`B_dK^{}K^+`$, $`B_s\pi ^+\pi ^{},\pi ^0\pi ^0`$ which only receive contribution from annihilation are not vanishing any more. The branching ratios are expected to be small. From Table 6, we indeed find that these branching ratios are smaller than others.

It is interesting to note that although the annihilation amplitudes are small, in certain decay modes, such as $`B_sK^+K^{}`$ and $`B_sK^0\overline{K}^0`$, the effects can be significant, the branching ratios almost doubled. This is because that although $`A_{\overline{3}}^P`$ is small compared with $`C_{3,\overline{15}}`$, it is comparable with $`C_6^T`$, but enhanced by a KM factor $`|V_{tb}V_{ts}^{}/V_{ub}V_{us}^{}|`$. These modes provide good places to study the annihilation contributions.

## 4 Discussions and Conclusions

In this talk we have presented SU(3) predictions of branching ratios and CP asymmetries for some $`BPP`$ decays in the Standard Model. There can be large CP violation in $`B_d\pi ^0\pi ^0,\pi ^+\pi ^{}`$ and $`B_sK^0\pi ^0,K^{}\pi ^{}`$. Also several $`B_s`$ decays can have large branching ratios. We presented results obtained with and without annihilation contributions. The results indicate that the annihilation contributions are in general small, but can still have large effect in several $`B_sK^+K^{},K^0\overline{K}^0`$ decay modes.

In SM predictions for branching ratios and CP asymmetries are possible is beacuse that SU(3) symmetry relates different decay modes. This symmetry is expected to be broken. In that case more parameters are needed to describe the decays. The effects of SU(3) breaking have to be further studied. We however expect that the general feature will not be be altered dramtically, and results obtained here can still provide some guidance in searching for large CP asymmetries and branching ratios in $`BPP`$ decays. Future experiments on $`BPP`$ will provide valuble information about the Standard Model.


text:
# Dynamically enhanced 𝑩→𝑲^∗𝝅 decays in perturbative QCD

## Abstract

We investigate the $`BK^{}\pi `$ decays, one of hardly understandable processes among charmless B-meson decays, within the perturbative QCD method. Owing to the dynamically enhanced mechanism in PQCD, we obtain large branching ratios and large direct CP asymmetries: $`Br(B^0K^\pm \pi ^{})=(9.1_{3.90.2}^{+4.9+0.3})\times 10^6`$, $`Br(B^\pm K^0\pi ^\pm )=(10.0_{3.5}^{+5.3}\pm 0.0)\times 10^6`$; $`Acp(B^0K^\pm \pi ^{})=(19.2_{1.7}^{+0.5})\%`$, and $`Acp(B^\pm K^\pm \pi ^0)=(43.7_{4.2}^{+4.0})\%`$. The branching ratios are consistent with experimental data and large direct CP violation effects will be tested by near future experimental measurements in Asymmetric B-factory.

The predictive power of the perturbative QCD (pQCD) approach has been demonstrated successfully in exclusive 2-body B-meson decays, especially charmless B-meson processes LiYu ; KLS01 ; LUY ; LuYa ; phiK ; phiKst ; KK ; keum  which is based on $`k_T`$ factorization theorem BS . This is a modified version of the pQCD theory for exclusive processesBL . The idea is to separate hard scattering kernels from a high-energy QCD process, which are calculable in a perturbative way . Nonperturbative parts are organized into universal hadron distribution amplitudes, which can be determined from experimental data. By introducing parton transverse momenta $`k_{}`$, we can generate naturally the Sudakov suppression effect due to the resummation of large double logarithms $`Exp[\frac{\alpha _sC_F}{4\pi }\mathrm{ln}^2(\frac{Q^2}{k_{}^2})]`$, which suppress the long-distance contributions in the small $`k_{}`$ region and give a sizable average $`<k_{}^2>\overline{\mathrm{\Lambda }}M_B`$. This can resolve the end point singularity problem and allow the applicability of pQCD to exclusive B-meson decays. We found that almost all of the contributions to the matrix element come from the integration region where $`\alpha _s/\pi <0.3`$ and the pertubative treatment can be justified.

In the pQCD approach, we can predict the contribution of non-factorizable term and annihilation diagram on the same basis as the factorizable one. A folklore for annihilation contributions is that they are negligible compared to W-emission diagrams due to helicity suppression. However the operators $`O_{5,6}`$ with helicity structure $`(SP)(S+P)`$ are not suppressed and give dominant imaginary values, which is the main source of strong phase in the pQCD approach.

An alternative method to exclusive $`B`$ meson decays is QCD-factorization approach (QCDF)BBNS , which is based on collinear factorization theorem.

For some modes, such as the $`BK\pi `$ decays, the difference between the pQCD and QCDF approaches may not be signigicant. As explained in ref.KLS01 , the typical hard scattering scale is about 1.5 GeV. Since the RG evolution of the Wilson coefficients $`C_{4,6}(t)`$ increase drastically as $`t<M_B/2`$, while that of $`C_{1,2}(t)`$ remain almost constant, we can get a large enhancement effects from both wilson coefficents and matrix elements in pQCD.

In general the amplitude can be expressed as

$$Amp[a_{1,2}\pm a_4\pm m_0^{P,V}(\mu )a_6]<K\pi |O|B>$$
(1)
with the chiral factors $`m_0^P(\mu )=m_P^2/[m_1(\mu )+m_2(\mu )]`$ for pseudoscalar meson and $`m_0^V=m_V`$ for vector meson. To accommodate the $`BK\pi `$ data in the factorization and QCD-factorization approaches, one relies on the chiral enhancement by increasing the mass $`m_0`$ to as large values about 3 GeV at $`\mu =m_b`$ scale. So two methods accomodate large branching ratios of $`BK\pi `$ and it is difficult for us to distinguish two different methods in $`BPP`$ decays. In addition, the direct CP asymmetries in $`BK\pi `$ decays are not large enough to distinguish the two approaches after taking into account the theoretical uncertainties.

However the difference can be detected in the direct CP asymmetry of $`B^0\pi ^\pm \pi ^{}`$ process because of the different power counting rules and the branching ratios of in $`BPV`$ modes since there is no chiral enhanced factor in LCDAs of the vector meson. Due to the different power counting rules of the QCDF and pQCD approaches, based on collinear and $`k_T`$ factorizations, respectively, the vertex correction is the leading source of the strong phase in the former, and the annihilation diagram is in the latter. The strong phases derived from the above two sources are opposite in sign, and the latter has a large magnitude. This is the reason QCDF prefers a small and positive CP asymmetry $`C_{\pi \pi }`$ Ben , while pQCD prefers a large and negative $`C_{\pi \pi }30\%`$ KLS01 ; keum .

We can test whether dynamical enhancement or chiral enhancement is responsible for the large $`BK\pi `$ branching ratios by measuring the $`B\varphi K^{()}`$ modes. In these modes penguin contributions dominate, such that their branching ratios are insensitive to the variation of the unitarity angle $`\varphi _3`$. According to recent works within QCDFCY , the branching ratio of $`B\varphi K`$ is $`(27)\times 10^6`$ including $`30\%`$ annihilation contributions in real part of amplitudes within QCD-factorization approach. However pQCD predicts $`10\times 10^6`$ phiK  with mostly pure imaginary annihilation contributions. For $`B\varphi K^{}`$ decays, QCDF gets about $`9\times 10^6`$CKY , but pQCD have $`15\times 10^6`$phiKst . Because of these relatively small branching ratios for $`BPV`$ and $`VV`$ decays in QCD-factorization approach, they can not globally fit the experimental data for $`BPP,VP`$ and $`VV`$ modes simultaneously with same sets of free parameters $`(\rho _H,\varphi _H)`$ and $`(\rho _A,\varphi _A)`$ zhu . To expalin large branching ratios of $`BVP`$ modes, they have to break the universality of free parameter sets with $`\rho _{H,A}1`$ and finally lost the predictive power.

In this letter we investigate the more complicated $`BK^{}\pi `$ processes, which contain both tree and penguin contributions, while $`B\varphi K^{()}`$ is a pure penguin process. It is well known that it is very difficult to explain the observed $`BK^{}\pi `$ branching ratios using the factorization assumption (FA) BSW  and QCDFzhu : the experimantal measurements are much larger than the theoretical predictions.

The reason is as follows. The measured $`Br(BK^{}\pi )`$ are roughly the same as $`Br(BK\pi )2\times 10^5`$. However, the penguin operators $`O_{5,6}`$ contribute to the latter, but not to the former. Due to the loss of this important piece of contributions, the predicted $`Br(BK^{}\pi )`$ become a quarter of the predicted $`Br(BK\pi )`$. The same difficulty has been encountered in the $`B\rho \pi `$ decays, where the vector meson is replaced by a $`\rho `$ meson LuYa . This controversy remains, no matter how the angle $`\varphi _3`$ is varied WS . Hence, the $`BK^{}\pi `$ decays is worth of an intensive study.

We shall evaluate the branching ratios of the following modes,

$`B^\pm K^0\pi ^\pm ,B_d^0K^\pm \pi ^{},`$
$`B^\pm K^\pm \pi ^0,B_d^0K^0\pi ^0,`$ (2)
and the CP asymmetries, for instance, $`f=K^+\pi ^{}`$

$$A_{CP}=\frac{\mathrm{Br}(\overline{B}_d^0\overline{f})\mathrm{Br}(B_d^0f)}{\mathrm{Br}(\overline{B}_d^0\overline{f})+\mathrm{Br}(B_d^0f)},$$
(3)
as functions of the unitarity angle $`\varphi _3`$. It will be shown that penguin and annihilation amplitudes in the $`BK^{}\pi `$ decays are greatly enhanced by Wilson evolutin effects. There is also small enhancement from the $`K^{}`$ meson wave functions, which are more asymmetric than the kaon wave funcitons, and from the $`K^{}`$ meson decay constant $`f_K^{}`$, which is larger than the kaon decay constant $`f_K`$. It turns out that these enhancements compensate the loss of the $`O_{5,6}`$ contributions, and that PQCD predictions are in agreement with the data.

The decay rates of $`BK^{}\pi `$ have the expressions,

$$\mathrm{\Gamma }=\frac{G_F^2M_B^3}{128\pi }|𝒜|^2.$$
(4)
The decay amplitudes $`𝒜`$ for the different modes are written as

$`𝒜(K^0\pi ^+)`$ $`=`$ $`f_K^{}V_t^{}F_e^{P(d)}+V_t^{}_e^{P(d)}+f_BV_t^{}F_a^{P(u)}`$ (5)
$`+V_t^{}_a^{P(u)}f_BV_u^{}F_a^TV_u^{}_a^T,`$
$`𝒜(K^+\pi ^{})`$ $`=`$ $`f_K^{}V_t^{}F_e^{P(u)}+V_t^{}_e^{P(u)}+f_BV_t^{}F_a^{P(d)}`$ (6)
$`+V_t^{}_a^{P(d)}f_K^{}V_u^{}F_e^TV_u^{}_e^T,`$
$`\sqrt{2}𝒜(K^+\pi ^0)`$ $`=`$ $`f_K^{}V_t^{}F_e^{P(u)}+V_t^{}_e^{P(u)}+f_BV_t^{}F_a^{P(u)}`$ (7)
$`+V_t^{}_a^{P(u)}+f_\pi V_t^{}F_{eK}^P+V_t^{}_{eK}^P`$
$`f_K^{}V_u^{}F_e^TV_u^{}_e^Tf_BV_u^{}F_a^T`$
$`V_u^{}_a^Tf_\pi V_u^{}F_{eK}^TV_u^{}_{eK}^T,`$
$`\sqrt{2}𝒜(K^0\pi ^0)`$ $`=`$ $`f_K^{}V_t^{}F_e^{P(d)}+V_t^{}_e^{P(d)}+f_BV_t^{}F_a^{P(d)}`$ (8)
$`+V_t^{}_a^{P(d)}+f_\pi V_t^{}F_{eK}^P+V_t^{}_{eK}^P`$
$`f_\pi V_u^{}F_{eK}^TV_u^{}_{eK}^T,`$
to which $`𝒜(K^0\pi ^{})`$, $`𝒜(K^{}\pi ^+)`$, $`\sqrt{2}𝒜(K^{}\pi ^0)`$, and $`\sqrt{2}𝒜(\overline{K}^0\pi ^0)`$ are identical, respectively, but with the the product $`V_t^{}`$ ($`V_u^{}`$) of the CKM matrix elements replaced by $`V_t`$ ($`V_u`$).

The detail expression of analytic formulas for all amplitudes ($`F_i^{T,P}`$ and $`M_i^{T,P}`$) will be presented elsewhereKL . In the above expressions $`f_B`$ is the $`B`$ meson decay constant. The notations $`F`$ represent factorizable contributions (form factors), and $``$ represent nonfactorizable (color-suppressed) contributions. The subscripts $`a`$ and $`e`$ denote the annihilation and W-emission topology, respectively. The superscript $`P(T)`$ denotes contributions from the penguin (Tree) operators. $`F_{a(e)}^T`$, associated with the time-like $`K^{}`$-$`\pi `$ form factor ($`B\pi `$ form factor), and $`_{a(e)}^T`$ are from the operators $`O_{1,2}^{(u)}`$. The factorizable contribution $`F_{eK}^P`$ ($`F_{eK}^T`$) is associated with the $`BK^{}`$ form factor from the penguin (tree) operators, and $`_{eK}^P`$ ($`_{eK}^T`$) is the corresponding nonfactorizable contribution.

In our numerical analysis, we use $`G_F=1.16639\times 10^5`$ GeV<sup>-2</sup>, the Wolfenstein parameters $`\lambda =0.2196`$, $`A=0.819`$, and $`R_b=0.38_{0.06}^{+0.10}`$ for the CKM matrix elements, the masses $`M_B=5.28`$ GeV and $`M_K^{}=0.892`$ GeV, and $`\overline{B}_d^0`$ ($`B^{}`$) meson lifetime $`\tau _{B^0}=1.55`$ ps ($`\tau _B^{}=1.65`$ ps) LEP . The angle $`\varphi _3=80^o`$ was extracted from the data of the $`BK\pi `$ and $`\pi \pi `$ decays KLS01 ; keum . With all the meson wave functions given in our previous works, we calculate the contributions from all the topologies as shown in Figs. 2 and 3 in $`BK\pi `$ paperKLS01 . The allowed range of B-meson shape parameter, $`0.36GeV<\omega _B<0.44GeV`$, and chiral factor for pion, $`1.2GeV<m_0^\pi <1.6GeV`$ is determined from the reasonable $`B\pi `$ and $`BK`$ transition form factors.

For example, all amplitudes for the $`B^+K^+\pi ^0`$ modes are listed in Table 1, whose values are mostly the same magitude for other decay channels, because the difference comes only from electroweak penguin contributions. We show in Table 2 the enhancing effect by comparing the decay amplitudes evaluated at the characteristic hard scales $`t`$ in PQCD and $`m_b/2`$ in QCDF. It is also found that the annihilation contributions are sizable in two-body charmless $`B`$ meson decays for the heavy-meson mass around 5 GeV KLS01  and in fact contributed about 60% fraction of the branching ratios, since factorized annihilation penguin contribution has large imaginary part and also the same order of magnitudes in real part as one of the factorized penguin contribution. As expected, the dominant factorizable penguin amplitudes are enhanced by about 30% due to the Wilson evolution, more than the factorizable tree amplitudes are.

The PQCD predictions for the $`BK^{}\pi `$ branching ratios, presented in Table 3, are consistent with present experimental data. Here the first uncertainty comes from the allowed ranges of both $`\omega _B`$ and $`m_0^\pi `$, and the second one comes from the uncertainty of $`R_b=|V_{ub}/V_{cb}|/\lambda `$.

Our predictions for the CP asymmetries in $`BK^{}\pi `$ decays are given in Table 4, which have the same sign as of those in $`BK\pi `$ decays. For $`\varphi _3=80^o`$, CP asymmetries of $`B^0K^\pm \pi ^{}`$ and $`B^\pm K^0\pi ^\pm `$ become large due to the important imaginary penguin annihilation amplitudes.

At last, the dependence of the branching ratios of $`B^0K^\pm \pi ^{}`$ and $`B^\pm K^0\pi ^\pm `$ on the angle $`\varphi _3`$ is shown in Fig.1 and large direct CP asymmetries of the $`B^0K^\pm \pi ^{}`$ and $`B^\pm K^\pm \pi ^0`$ on the angle $`\varphi _3`$ is exhibited in Fig.2. The branching ratios of $`K^\pm \pi ^{}`$ and $`K^\pm \pi ^0`$ increase with $`\varphi _3`$ rapidly, while the $`K^0\pi ^\pm `$ and $`K^0\pi ^0`$ modes are insensitive to the variation of $`\varphi _3`$. The increase with $`\varphi _3`$ is mainly a consequence of the inteference between the penguin contribution $`F_e^P`$ and the tree contribution $`F_e^T`$. The dependence on both the shape parameter $`w_B`$ for the $`B`$ meson wave function and the chiral factor $`m_0^\pi `$ is also shown, which is strong in the $`K^{}\pi ^\pm `$ and $`K^{}\pi ^0`$ modes, and weak in the other two. The sensitivity is attributed to the fact that the former contain both $`F_e^P`$ and $`F_e^T`$, which involve the $`B`$ meson wave function, while the latter contain only $`F_e^P`$.

In this letter in order to explain one of the hardly understandable processes in charmless B-decays, we have investigated the dynamical enhancement effect in the $`BK^{}\pi `$ decays within pQCD method. Owing to the dynamical enhancement of penguin contributions at $`t1.5`$ GeV, pQCD predictions for all the $`BK^{}\pi `$ modes are consistent with the present experimental data, which is a crucial decay process to distinguish pQCD from other approachs. We also predicted large direct CP asymmetry for $`B^0K^\pm \pi ^{}`$ about -20% and for $`B^\pm K^0\pi ^\pm `$ about -40%, which can be tested in near future measurements. More detail works will appear elsewhereKL .

We wishes to thank S.J. Brodsky, H.Y. Cheng, H.-n. Li, A.I. Sanda and G. Zhu for helpful discussions. This work was supported by Grant-in-Aid for Scientific Research from Ministry of Education, Science and Culture of Japan.


text:
# REFERENCES

Eggert, Gustafsson, and Rommer Reply to the “Comment on ‘Phase diagram of an impurity in the spin-1/2 chain: two channel Kondo effect versus Curie law’ ”
The Comment introduces the particular point at $`J_1\mathrm{}`$ as a new fixed point in the phase diagram of the $`J_1`$-$`J_2`$-impurity model which was analyzed in our recent Letter. The point at $`J_1\mathrm{}`$ is certainly worth a separate discussion which we will give here, but our analysis comes to rather different conclusions than Zvyagin.

At $`J_1\mathrm{}`$ the three spins $`\stackrel{}{S}_0`$, $`\stackrel{}{S}_1`$, and $`\stackrel{}{S}_N`$ are strongly coupled and form a complex of total spin $`s=1/2`$, which is characterized by a triplet between the spins $`\stackrel{}{S}_1`$ and $`\stackrel{}{S}_N`$, which in turn is antiferromagnetically correlated with $`\stackrel{}{S}_0`$. This three-spin complex is therefore not decoupled from the rest of the chain, but the effective coupling is given by $`2J/3`$ instead. The point at $`J_1\mathrm{}`$ therefore does not represent a fixed point at all since it does not correspond to a simple boundary condition in the spin chain. This is in contrast to the fixed point $`\mathrm{O}_{\mathrm{N}2}\frac{1}{2}`$ at $`J_2\mathrm{}`$ where the three-spin complex is indeed decoupled from the chain since the spins $`\stackrel{}{S}_1`$ and $`\stackrel{}{S}_N`$ form a singlet. The effective coupling there is given by $`J^2J_1/4J_2^2`$ from a third order perturbation expansion, which is always irrelevant because this coupling is suppressed by one additional power of the cutoff due to the virtual excitations.

For simplicity we will label the point at $`J_1\mathrm{}`$ by $`\mathrm{P}_{\mathrm{N}1}`$ because three spins are effectively replaced by one spin-1/2 complex which remains coupled to the chain. This point $`\mathrm{P}_{\mathrm{N}1}`$ is again characterized by a logarithmically diverging impurity susceptibility and a ferromagnetic correlation $`\stackrel{}{S}_1\stackrel{}{S}_N>0`$, just like the fixed point $`\mathrm{P}_{\mathrm{N}+1}`$. We therefore do not expect any phase transitions or any discontinuities in the order parameter between the two points.

As far as the field theory description near the fixed point $`\mathrm{P}_{\mathrm{N}+1}`$ is concerned, we must emphasize again that the only leading irrelevant operator is given by $`_x\mathrm{trg}`$. The operators mentioned in the Comment are not present at $`\mathrm{P}_{\mathrm{N}+1}`$, since the impurity spin $`\stackrel{}{S}_0`$ has been absorbed in the chain and $`\stackrel{}{S}_0`$ can therefore not be used as an independent degree of freedom to construct operators.

In conclusion we find that the point $`\mathrm{P}_{\mathrm{N}1}`$ is not a fixed point and it appears to be in the same phase as $`\mathrm{P}_{\mathrm{N}+1}`$. The phase diagram as shown in Fig. 1 is therefore complete and correct.

The comparison to the integrable impurity model made in the Comment is certainly interesting. However, obviously the integrable impurity model lives in a different parameter space, so that a direct comparison of the impurity susceptibilities may not be very meaningful. An impurity in a free Fermion model (xx-model) also gives only limited insight since the scaling dimensions are fundamentally different in interacting systems. However, it would be interesting to examine the integrable model in more detail with field theory methods in an expanded parameter space, which, however, is not as trivial as indicated in the Comment. This would answer the question if it may belong into one of the phases that were discussed in our Letter or if it may correspond to a non-generic unstable multi-critical point as has been found for a related integrable impurity model.
Sebastian Eggert, David P. Gustafsson, and Stefan Rommer


text:
# Status report on TAUOLA, its environment, and its applications

## 1 INTRODUCTION

The package TAUOLA  for the simulation of $`\tau `$-lepton decays and PHOTOS  for the simulation of radiative corrections in decays, have a rather long history. Written and maintained by well defined authors, they nonetheless migrated into a wide range of applications, where they became ingredients of complicated simulation chains. As a consequence, a large number of different versions are currently in use. From the algorithmic point of view, they often differ only in a few small details, but incorporate substantial amounts of specific results from distinct $`\tau `$-lepton measurements. Such versions were mainly maintained by the experiments taking precision data on $`\tau `$-leptons. On the other hand, many new applications were developed recently, often requiring program different versions because of interfaces to other packages.

In the following, I will concentrate on those topics where changes with respect to the status presented at the Victoria $`\tau `$ conference two years ago  were introduced. Since that time, there were no changes introduced into the PHOTOS Monte Carlo functionality, and also the TAUOLA interfaces to KORALB , KORALZ , and KKMC  remain unchanged. On the contrary, the universal interface of TAUOLA evolved, and new applications, in particular observables for the measurement of the Higgs-boson parity, became possible. The new program MC-TESTER  instrumental in the development of future versions of TAUOLA, was developed. A new choice of hadronic currents for $`\tau 4\pi \nu `$ decay modes became available. It is based on Novosibirsk data.

Let me concentrate, in the following three sections, on these topics, and close my contribution with a summary.

## 2 THE $`\tau `$-LEPTON FOR HIGGS BOSON PARITY

In many applications, the precision of the control of the $`\tau `$ spin effects is not crucial. This is the case in searches of new particles yet to be discovered, or in applications where $`\tau `$-lepton decays contribute as a final state of some rare decays of known particles. Nonetheless, in such cases, spin effects can be of some use as well. As these, there is no motivation to develop sophisticated spin algorithms for every individual case: less precise, but universal solutions are welcomed.

The universal solution of , based on the HEPEVT common block of FORTRAN77 is now distributed with TAUOLA. The basic idea was to calculate the spin (helicity) state of the decaying $`\tau `$ from the kinematical information available in HEPEVT and some very simple assumptions on the production mechanisms. The program was checking if the production was through $`f\overline{f}Z/\gamma \tau ^+\tau ^{}`$, $`W`$ or Higgs-boson intermediate states/processes. The main properties of the algorithm were presented in Victoria already 2 years ago; modifications necessary for the program to work in the case of some $`\tau `$ production processes involving new, to-be-discovered particles, have gradually been introduced, but one can see  that not in all cases is the effort completed.

In the meantime algorithm functionality was also extended: full spin effects were introduced in the case of Higgs-boson decay . The case of the Higgs boson is exceptionally easy, since the full density matrix of the $`\tau `$-lepton pair produced from a Higgs-boson is fully defined by boson parity and $`\tau `$ leptons four-momenta.

This algorithm helped us to design an observable  for the possible measurement of Higgs boson parity in future accelerators such as Linear Colliders or the LHC. Let us recall the main principle of this technique, which from the point of view of spin analysis is quite involved, as it requires a study of full spin correlation spanning over three levels of decay cascade: $`h/A\tau ^+\tau ^{}`$, $`\tau ^\pm \rho ^\pm \nu `$ and $`\rho ^\pm \pi ^\pm \pi ^0`$. The presence of two non-observable neutrinos as well as our inability to reconstruct the Higgs-boson rest frame sufficiently well (with a precision comparable to the $`\tau `$-lepton mass) complicates the picture even further.

We started from the observation that the distribution of the acoplanarity angle of two planes spanned on decay products of $`\rho ^+\pi ^+\pi ^0`$ and $`\rho ^{}\pi ^{}\pi ^0`$ (defined in the rest frame of a $`\rho ^+\rho ^{}`$ pair) is quite sensitive, see fig. 1,

to the parity of Higgs boson. One has to select sub-samples of particles, with the same sign (opposite sign) of the energy difference between $`\pi ^+`$ and $`\pi ^0`$ in the $`\tau ^+`$ rest frame and that between $`\pi ^{}`$ and $`\pi ^0`$ in the $`\tau ^{}`$ rest frame. This condition is unfortunately impossible to realize in practice, as $`\tau `$-lepton momenta cannot be reconstructed.

In the next step of our study we have weakened the requirement, and we have used a replacement for the $`\tau `$-lepton momentum; the $`\rho ^+`$ and $`\rho ^{}`$ directions (in the rest frame of a $`\rho ^+\rho ^{}`$ pair) replace the unobservable directions of $`\tau `$ flights and the Higgs mass constraint provides the $`\tau `$-lepton energy (also in the same frame). Indeed, with such a technique, and after including some assumption on detector smearings, the sensitivity of our observable diminished significantly (see fig. 2),

but nonetheless remained sizeable. Enough to guarantee the meaurement of Higgs boson parity to a confidence level greater than 95%, if typical assumptions on linear collider luminosity, and (120 GeV) Higgs-boson production mechanism are made. Our method proved quite stable, with a degrading assumption on angular and energy resolution of the detector. Even a reduction by a factor of 2–3 with respect to typical assumptions used in linear collider designs did not destroy the method. To preserve its sensitivity, the angular resolution must be $`𝒪(m_\rho /E_\rho )`$.

## 3  MC-TESTER

In the previous case the HEPEVT data structure was used to search for possible hard processes to calculate the spin state of the decaying $`\tau `$. One can search over an event record such as HEPEVT for other information as well. An essential step before attempting any sizeable rebuilding of the TAUOLA library is to automate at least some of its tests (such rebuilding may soon become necessary because of improving $`\tau `$-lepton data). Such a test-program may be useful for other applications, not necessarily related to generators of $`\tau `$-lepton decay, as well. That is why we have developed an individual program for that purpose. The idea behind MC-TESTER  is quite simple: the user loads an extra library and, after the generation of every event with his generator, calls MC-TESTER, specifying the identifier for the particle to be searched for over the whole event record and studied. A data file from such a run is formed.

In the second (analysis) step, data files from two runs of different Monte Carlo programs can be compared. The output is given in the form of a file, which includes a table of all decay modes found(see fig. 3).

The table is not only a list of these modes, but also, of the branching ratios calculated from the two runs. Number quantifying maximum of shape differences is also given. This number is calculated as a maximum of the shape difference parameters (with a part of the code easy to identify) for all invariant mass distributions of the decay channel under study. We choose to take into consideration invariant mass of every subset of the decay products for the decay channel. The table is followed by a booklet of these invariant mass distributions (see fig. 4 for an example of an individual plot), which are grouped into separate chapters for every decay mode.

## 4  TAUOLA

As was already discussed in Victoria , many options of physics initialization are now available for TAUOLA (they can be constructed from the files included in the program distribution package). They often differ by the models used for fits to the experimental data, and/or by the data themselves. Let me stress that, in general, the best parametrizations will always be in the hands of the leading $`\tau `$-lepton experiments of the time. Our original version of TAUOLA initialization was indeed meant to be used only after re-initialization of its physics content by the users. Nonetheless for experiments in the design phase or for which $`\tau `$-leptons physical processes help in other studies, such a solution is not good. We are gradually including new options into TAUOLA.

In the last two years, parametrizations of the hadronic form factors  based on the Novosibirsk low energy $`e^+e^{}`$ annihilation data  were introduced into TAUOLA as a possible new option for $`\tau ^{}\pi ^{}\pi ^{}\pi ^+\pi ^0\nu `$ and $`\tau ^{}\pi ^{}\pi ^0\pi ^0\pi ^0\nu `$ decay modes . This is the only published contribution documenting upgrades on the physics content of TAUOLA over the last two years.

### 4.1 Comment on one of the old tests

In ref.  (formulae (37) and (38)) we have reported that, for $`\tau `$ decay into 4$`\pi `$, when the chiral limit is taken (the case we used for some tests of phase-space generation), the results of TAUOLA agree with those of ref. , only up to an overall factor111 Let us note that problems with normalization (or/and with misprints) in  were again addressed . The authors of that paper worry that it may be of grave consequences for TAUOLA as well. In practice, even if the normalization of the chiral current was wrong (by the same factor in the Monte Carlo and in the testing function), it would affect neither technical tests nor actual results of TAUOLA generation, since a different current is used then. Note that in cases such as CLEO initialization , the current is normalized directly from the $`\tau `$ data collected by the experiments; the normalization originating from the model disappears completely. This is why the issue is only of marginal interest for TAUOLA as it is now. $`\frac{243}{2n!}\mathrm{cos}^2\theta _c`$. In the case of $`\tau \pi ^{}\pi ^{}\pi ^+\pi ^0\nu `$, ($`n=2`$) the test formula

$$\frac{\mathrm{\Gamma }(2\pi ^{}\pi ^+\pi ^0)}{\mathrm{\Gamma }_e}=\frac{\mathrm{cos}^2\theta _c}{15}\left(\frac{m_\tau }{2\pi f_\pi }\right)^4\frac{1}{128}C_N$$
(1)
returns 0.0271364 with the numerical coefficient $`C_N=1261/120\pi ^2`$ (for compatibility, we use as in  $`m_\tau =1.7842`$, $`f_\pi =0.0933`$ and $`\mathrm{cos}\theta _c=0.975`$). A later calculation (formula 33 in ) finds $`C_N=1009/96\pi ^2`$ and result is 0.0272249. With the computer power available at the time of , we could not distinguish between the two options. At present, such an exercise is straightforward and we find that TAUOLA returns $`0.0272060\pm 0.000011`$. We can see, that the value of $`C_N`$ from ref.  is statistically favoured over the one from ref. .

## 5  SUMMARY

Let me recall finally the main developments of the TAUOLA package in the last two years: as stand-alone generator of $`\tau `$-lepton decays, and for its interfaces to other programs. The progress on $`\tau `$-lepton decay alone consists of a new parametrization of form factors for the $`4\pi `$ decay channels. This parametrization, mainly based on the Novosibirsk data, form a step toward TAUOLA as a framework, where different data and models can be compared. In such multi parametrization approach the appropriate tools are necessary. The MC-TESTER, helpful in different types of works for programs such as TAUOLA became available. It can also be useful, for example, to check if the interfacing of TAUOLA to other programs did not corrupt its initialization. The family of TAUOLA interfaces was enriched with the new option for its universal interface. It now provides complete spin correlations for Higgs boson decay into pair of $`\tau `$ leptons. This option was instrumental in designing an observable that is potentially useful to measure Higgs-boson parity at future experiments.

Let me finally acknowledge people and collaborations who contributed to the present shape of TAUOLA and related projects: A. Bondar, S. Eidelman, P. Golonka, S. Jadach, M. Jeżabek, J.H. Kühn, A. Milstein, T. Pierzchala, E. Richter-Was, N. Root, B.F.L. Ward, M. Worek, and ALEPH, CLEO (in particular A. Weinstein), Delphi, Opal, L3 collaborations, the Karlsruhe theory group, and others. I would like to thank A. Pich for discussion.