text:
# Creep of superconducting vortices in the limit of vanishing temperature: A fingerprint of off-equilibrium dynamics.

\[

## Abstract

We theoretically study the creep of vortex matter in superconductors. The low temperatures experimental phenomenology, previously interpreted in terms of “quantum tunnelling” of vortices, is reproduced by Monte Carlo simulations of a purely “classical” vortex model. We demonstrate that a non-zero creep rate in the limit of vanishing temperature is to be expected in systems with slow relaxations as a consequence of their off-equilibrium evolution in a complex free energy landscape.
PACS numbers: 74.60.Ge 47.32.Cc 74.50.+r 75.45.+j

\]

There exist an abundance of experimental evidence that the relaxation rate of the magnetisation in type II superconductors (ranging from conventional, to organic, heavy fermions and a variety of $`HT_c`$ superconductors ) does not vanish as the temperature $`T`$ is lowered towards zero (see eg. ). This may seem surprising if one assumes that the mechanism allowing the magnetisation to relax is thermal activation over a characteristic energy barrier $`\mathrm{\Delta }E`$. Namely, when $`T\mathrm{\Delta }E`$ the relaxation rate should vanish since the Arrhenius factor for thermal activation, $`\mathrm{exp}(\mathrm{\Delta }E/T)`$, goes to zero. The question then arises how relaxation can take place at a finite rate while the thermal activation factor exponentially approaches zero.

A number of researchers have suggested that such a phenomenon is caused by quantum tunnelling of vortices through the barriers in the random pinning potential (for a theoretical review see ). The above “quantum” explanation is very intriguing and in good agreement with some experimental results in compounds such as YBCO  or BSCCO . Also other even more exotic materials, such as organic superconductors , give good correspondences. However, in other non conventional systems such as heavy fermion superconductors, the theory of quantum creep is totally unable to describe the observed low $`T`$ relaxation . In fact, strong discrepancies are found in many other systems ranging from PCeCO crystals to YBCO/PBCO multilayers or YBCO and BSCCO films and crystals . One problem is  that the length of the tunnelling vortex segment $`L_c(0)`$ needed to fit the creep rate data, $`S`$, can be orders of magnitude larger that the one theoretically predicted by quantum creep theory . Also the experimental temperature dependence of the creep rate, $`S`$, is often very different from the one predicted by quantum theory (see ).

The above contradictory results suggests to look for additional descriptions of the anomalous low $`T`$ magnetic relaxation. It is worth to stress that the observation of a non-vanishing constant creep rate in the limit $`T0`$ is found under very general circumstances: it does not crucially depend on the thickness of the sample  (i.e. on its dimensionality), nor on whether the pinning is caused by columnar defects or random point pins . Thus, the mechanism behind the low temperature creep seems to be of a fundamental and basic nature.

We demonstrate below that also in a “classical” system (i.e., not “quantum”) logarithmically slow glassy dynamics can naturally persist even at vanishing temperatures and can lead to the experimentally observed phenomenology. This is possible because the low $`T`$ off-equilibrium dynamics consists of searching, among a very large number, for a few “downhill” or “flat” directions in the free energy landscape. The number of these directions decreases as relaxation proceeds though there always remain some. They can be found only by collective cooperative rearrangements of the system, resulting in a slowing down of relaxation .

The model – We study a statistical mechanics model for vortex matter called a Restricted Occupancy Model (ROM) . In the limit of zero temperature and infinite upper critical field, it reduces to a cellular automaton introduced in  to study vortex avalanches. We use Monte Carlo dynamics which enabled the ROM to depict a unified picture of creep and transport phenomena in vortex physics, ranging from magnetisation loops with “anomalous” second peak, logarithmic relaxation, Bean profiles, to history dependent behaviours in vortex flow and I-V characteristics, to the reentrant nature of the equilibrium phase diagram . The model also predicts the existence of a “glassy region” at low temperature with strong “aging” effects .

Here we use the ROM to study the magnetic relaxation rate, $`S`$, in the very low temperature limit. Interestingly, the model reproduces the experimental “anomalous” relaxation and the observed behaviour of $`S`$ .

The ROM model is described in full details in Ref.. A system of straight parallel vortex lines is coarse grained in the $`xy`$-plane by introducing a square grid of lattice spacing $`l_0`$ of the order of the London penetration length, $`\lambda `$ . The number of vortices on the $`i`$-th coarse grained plaquette is denoted by $`n_i`$. The occupancy of each plaquette is a number larger than zero and, importantly, smaller than $`N_{c2}=B_{c2}l_0^2/\varphi _0`$, where $`B_{c2}`$ is the upper critical magnetic field and $`\varphi _0=hc/2e`$ is the magnetic flux quantum. The ROM model is thus defined by the following coarse grained vortex interaction Hamiltonian : $`=\frac{1}{2}_{ij}n_iA_{ij}n_j\frac{1}{2}_iA_{ii}n_i_iA_i^pn_i`$. The first two terms describe the repulsion between the vortices and their self energy. On-site and nearest neighbour interactions are included, i.e., the interaction between vortex lines with a separation greater than the London screening length (which by definition is close to $`l_0`$) is ignored. We choose $`A_{ii}=A_0=1`$; $`A_{ij}=A_1`$ if $`i`$ and $`j`$ are nearest neighbours and $`A_{ij}=0`$ otherwise. The last term in $``$ describes a delta-distributed random pinning $`P(A^p)=(1p)\delta (A^p)+p\delta (A^pA_0^p)`$. Interestingly, the general scenario of creep phenomena we describe below does not depend on the details of pinning in the system (a fact in correspondence with experimental results ). In our model $`A_0`$ sets the energy scale. Below we choose $`A_0=1.0`$; $`A_0^p=0.3`$; $`N_{c2}=27`$; $`p=1/2`$; $`\kappa ^{}A_1/A_0[0.26,0.3]`$.

The relaxation of the model is simulated by use of Monte Carlo dynamics on a square lattice in presence of a thermal bath of temperature $`T`$. The system is periodic in the $`y`$-direction. The two edges parallel to the $`y`$-direction are in contact with a vortex reservoir. Particles can enter and exit the system only through the reservoir, which plays the role of the external magnetic field. Hence the reservoir density, $`N_{ext}`$, is used as the external control parameter. We perform the following zero field cooled experiment: at a low temperature $`T`$ we increase at constant rate $`\gamma =\mathrm{\Delta }N_0/\tau `$ the reservoir density from zero to a working value $`N_{ext}`$. We keep $`N_{ext}`$ fixed while we monitor the time dependence of the magnetisation $`M=N_{in}N_{ext}`$. Here $`N_{in}=_in_i/L^2`$ is the vortex density inside the system (of size $`L^2`$ ). Time is measured in units of single attempted updates per degrees of freedom of the lattice (see Ref. ). The data presented below are averaged over 128 realizations of the pinning background.

In particular, we investigate the creep rate

$$S=\left|\frac{\mathrm{ln}(M)}{\mathrm{ln}(t)}\right|$$
(1)
as function of $`T`$, $`N_{ext}`$ and $`\gamma `$. In typical experiments the nature of the $`t`$ dependence of $`M`$ is such that $`S`$ decreases in time. So usually, one deals with an average creep rate, $`S_a`$, in some given temporal window . As shown in the upper inset of Fig.2, in our model dynamics $`M(t)`$ at low temperatures behaves according to the known logarithmic interpolation formula (see Ref.) found in experiments , namely: $`M(t)M(0)\mathrm{\Delta }M_{\mathrm{}}\{1[1+\frac{\mu T}{U_c}\mathrm{ln}(\frac{t+t_0}{t_0})]^{1/\mu }\}`$. Here, $`M(0)`$ is the magnetisation at the time of preparation of the sample ($`t=0`$), $`\mathrm{\Delta }M_{\mathrm{}}`$ its overall variation, the exponent $`\mu `$ is consistent with 1, $`\frac{\mu T}{U_c}`$ and $`t_0`$ are fit parameters .

Consistently, we define $`S_a`$ as the average value of $`S`$ in the last time decade of our measures (i.e., for $`t[10^5,10^6]`$). The present choice, analogous to those made for experimental data, is the most natural one and the general results presented below do not depend on it.

For the reasons explained in the introduction, a very important physical quantity is the distribution, $`P(\mathrm{\Delta }E)`$, of the energy barriers, $`\mathrm{\Delta }E`$, that vortices segments meet during their motion. Since at low $`T`$ the system is typically off-equilibrium, $`P(\mathrm{\Delta }E)`$ is itself a (logarithmically slow) function of $`t`$ and we consider its average over the last time decade of our measurements.

Results – When the temperature is very low, the model exhibits the same kind of “anomalous” creep found in the experiments on superconductors. In Fig.1, we plot the creep rate, $`S_a`$, as a function of $`T`$ in a broad temperature range. For comparison we present equivalent experimental measurements in a BSCCO single crystal (from Ref.) as inset. The numerical values found for $`S_a`$ at low $`T`$ in our model and in real samples are interestingly very similar. The temperature scales of the simulations and of real experiments can be compared by considering that the ratio $`T/A_0`$ in our model is of the same order of magnitude as $`T/T_c`$ in a real superconductor. This is seen from a comparison of the $`(T,N_{ext})`$ equilibrium phase diagram of our model with the equilibrium temperature-magnetic field, $`(T,H)`$, phase diagram of, say, a BSCCO superconductor (see Ref.).

In both experiments and simulations, $`S_a`$ approaches a finite value, $`S_a^0`$, when $`T0`$. In particular, we find that a linear fit of $`S_a(T)`$ in the low $`T`$ regime is very satisfactory (see Fig.1):

$$S_a(T)=S_a^0+\sigma T$$
(2)
where both $`S_a^0`$ and $`\sigma `$ are a function of the applied field $`N_{ext}`$. We also note that in the present model $`S_a(T)`$ is non monotonous in $`T`$: in the higher $`T`$ region it starts decreasing. This is also a known experimental fact , we discuss it later on. The maximum in $`S_a(T)`$ is just above a characteristic crossover “glassy” temperature, $`T_g`$, defined in .

In our model by varying the applied field we find a range of values for $`S_a^0`$ very similar to experimental ones  (see Fig.2). In particular, $`S_a^0`$ seems to decrease on average by increasing the field $`N_{ext}`$. The overall behaviour can be approximately interpolated with a power law: $`S_a^0(N_{ext})(N_{ext}/N_0)^x`$, where, for $`\kappa ^{}=0.28`$, $`N_00.01`$ and $`x0.6`$. As shown in Fig.2, the presence of a small exponent $`x`$ implies that sensible variations in $`S_0^a`$ can be seen only by changing $`N_{ext}`$ of orders of magnitude. Note that the dips in the $`S_a(0)`$ versus $`N_{ext}`$ data in Fig. 2 at certain values of $`N_{ext}`$ (namely around 3, 13, and 18) are statistically significant. They are related to the low field order-disorder transition, the 2nd peak transition and the reentrant high field order-disorder transition respectively .

In the lower inset of Fig.2, we show that $`S_a`$ is essentially independent of the ramping rate, $`\gamma `$ (the values shown are for $`N_{ext}=10`$ and $`T=10^4`$). This is an other typical experimental observation . However, a very small decrease of $`S_a`$ with increasing $`\gamma `$ cannot be excluded: we show a fit to the form $`S_a^0S_1+s_2\mathrm{ln}(\gamma )`$, with $`S_1=0.11`$ and $`s_2=210^4`$. The fact that $`S_a^0`$ is practically independent on $`\gamma `$, far from being a proof of the presence of equilibrium in the system, is due to the fact that at very low temperatures the characteristic equilibration time, $`\tau _{eq}`$, is enormous (see Ref.). So whenever the driving rate, $`\gamma `$, is much larger than $`\tau _{eq}^1`$ the off-equilibrium state and dynamics of the system are essentially independent of $`\gamma `$. Stronger $`\gamma `$ effects have to be expected when $`\gamma `$ gets closer to $`\tau _{eq}^1`$. In fact, it is experimentally well known that at higher temperatures the systems exhibit strong $`\gamma `$ dependent “memory” effects , the signature of off-equilibrium dynamics. Actually, in the present model at low $`T`$, it is possible to show that $`\tau _{eq}(T)`$ diverges exponentially : $`\tau _{eq}(T)\mathrm{exp}\left(\frac{E_0}{T}\right)`$. In that region, the typical observation time windows, $`t_{obs}`$, are such that $`t_{obs}/\tau _{eq}1`$, and the system is in the early stage of its off-equilibrium relaxation from its initial state. This is schematically the origin of the flattening of $`S_a`$ at very low $`T`$. Notice that, if one could observe the system for an exponentially long time, i.e., if $`t_{obs}/\tau _{eq}1`$, then the creep rate, $`S_a`$, would indeed go to zero.

The above scenario is clarified by the analysis of the energy barrier distribution function, $`P(\mathrm{\Delta }E)`$, recorded during the system evolution at very low $`T`$. Such a quantity also clearly shows the simple mechanical origin of the anomalous creep found at very low temperature in the present model. The function $`P(\mathrm{\Delta }E)`$ (where $`\mathrm{\Delta }E`$ is in units of $`A_0`$), recorded at $`T=10^4`$, is plotted in Fig.3 for two values of the applied field, $`N_{ext}`$. We always find that $`P(\mathrm{\Delta }E)`$ has support also on the negative axis. This is the mark of the off-equilibrium nature of the evolution on the observed time scales. The presence of a $`P(\mathrm{\Delta }E)`$ which extends down to negative values also explains the presence of the recorded relaxation at low $`T`$: in the configuration space the system can still find directions where no positive barriers have to be crossed. The insert in Fig. 3 clarifies the mechanism behind the relaxation. Here, we plot the signal $`A(t)`$ defined, for each single Monte Carlo (MC) step $`t`$, in the following way: $`A(t)=0`$ if the MC trial is rejected; $`A(t)=1`$ if the trial is accepted and the energy reduced, i.e. $`\mathrm{\Delta }E0`$; and finally $`A(t)=2`$ when a trial is accepted with $`\mathrm{\Delta }E>0`$. We plot two sequences of trials. One for $`0t500`$ was measured at the early stage of the relaxation, the second sequence, placed at the interval $`500t1000`$, is measured at the late stage of the relaxation. Most trials are rejected ($`A(t)=0`$) and only once in a while the system does manage to find a route pointing downhill in the energy landscape. $`A(t)=2`$ never occurs. As time proceeds fewer and fewer “negative channels” are available to the relaxation and a decrease in the density of the spikes in $`A(t)`$ is observed.

As the temperature is increased thermal activation over positive energy barriers will become possible as the Arrhenius factor $`\mathrm{exp}(\mathrm{\Delta }E/T)`$ assumes a non-vanishing value for an appreciable range of barrier values $`\mathrm{\Delta }E<T`$. When this happens the relaxation will occur sufficiently fast to allow one, within the experimental time window, to closely approach the equilibrium configurations where the vortex density profile is more or less flat and relaxation ceases, hence $`S_a`$ goes down, as seen in experiments  and in Fig. 1.

Finally, we stress that slow off-equilibrium relaxations at very low temperatures are also observed in glass forming liquids . In that cases too, no activation over barriers occurs and the system simply wanders in its very high dimensional phase space through the few channels where no energy increase is required.

Discussion \- We have above demonstrated that the phenomenological behaviour of the creep rate at low temperatures can be understood in terms of the off-equilibrium nature of the inhomogeneous vortex density profile produced in magnetic creep experiments. In fact, cooperative mechanical rearrangements, possible even at very low $`T`$ (where thermal activation over positive barriers can be negligible), dominate the phenomenon . In this perspective, it is very important to stress that the system’s equilibration time at very low temperature is much larger than any experimentally accessible time window . Accordingly, we can say that experimental findings do not enforce the interpretation in terms of macroscopic quantum tunnelling of vortices.

Relaxation due to quantum tunnelling might be present along with the mechanical relaxation discussed above. It is then important to ask how compelling the quantum tunnelling interpretation is. The theory of quantum tunnelling assumes a London picture and treats the position of the vortex core as the variable that is able to tunnel. It is not entirely clear if this is the right level at which to introduce quantum fluctuations, but, more importantly, the quantum tunnelling description assume the existence of a characteristic energy barrier , not a time dependent distribution of barrier heights as typically found in many-body systems relaxing off-equilibrium. Finally, the quantum tunnelling description also tacitly assumes the existence of a static equilibrium state in which barriers are always positive. As we have clearly shown above, this is typically not the case and a dynamical approach is more appropriate.

The present scenario, where off-equilibrium phenomena dominate the anomalous low $`T`$ creep, could be experimentally verified by the discovery of “aging” phenomena  like those recently observed in .

Work supported by EPSRC, INFM-PRA(HOP)/PCI.


format: text

text:
# Adjusting the melting point of a model system via Gibbs-Duhem integration: application to a model of Aluminum

## I Introduction

The ability of a simulation to successfully predict the properties of real materials is primarily dependent upon the accuracy of the model interaction potential used. The construction of model interactions generally involves the optimization of the parameters of the potential with respect to mechanical properties of the material (crystal lattice constants, elastic constants, liquid density, etc.) as determined from experiment or ab initio calculations. Non-mechanical properties (i.e. those not obtainable as local averages over coordinates) such as phase transition temperatures are difficult to include in such optimization procedures and are generally calculated for the optimized model a posteriori, and the agreement of such quantities with experiment is not guaranteed. However, for some applications that deal directly with such properties, such as in the study of solid-liquid interfaces in which the melting temperature plays an obviously important role, it is desirable to develop efficient procedures for including such non-mechanical properties in the optimization. In this work, we outline a general procedure for adjusting the potential parameters for a system designed to model a real system to improve the agreement of the melting point of that system with the experimental value. As an example we present an application to an existing embedded atom model of aluminum .

Our re-parameterization scheme is based on the powerful Gibbs-Duhem integration method developed by Kofke. In this technique, the derivative along the coexistence curve of any coexistence condition (such as melting temperature or pressure) with respect to any parameter of the potential can be determined by an appropriate configurational average at a previously determined melting point. The method is generated by the integration of a generalized Gibbs-Duhem equation and the steps are analogous to the derivation of the familiar Clapeyron equation for the slope of the P-T coexistence curve. Gibbs-Duhem integration has been shown to be quite successful in efficiently determining the coexistence conditions for entire classes of potentials. For example, the phase diagram for the class of repulsive inverse power potentials, $`u(r)=ϵ\left(\frac{\sigma }{r}\right)^n`$, was determined by starting with the known hard-sphere ($`n=\mathrm{}`$) coexistence and integrating the derivative of coexistence curve with respect to the parameter $`s1/n`$. The method has also proved useful in a variety of other applications .

In the current application, one begins with a model potential, parameterized for a real system in the usual way with respect to mechanical properties of the real system. The melting temperature (or pressure) for the model system is then calculated by thermodynamic integration. Once this is done, the derivative of the melting temperature (or pressure) with respect to all parameters of the system can be determined via separate simulations on the coexisting fluid and solid using the Gibbs-Duhem procedure. The calculated derivatives allow an assessment of the effect of each individual parameter on the melting point. From this information an appropriate scheme to adjust the parameters to improve the melting point can be devised in such a way that the agreement with the other experimental properties is not unacceptably compromised. The Gibbs-Duhem integration and our re-parameterization scheme is outlined in more detail in the next section.

As a test application of this procedure, we examine an embedded atom model of aluminum developed by Mei and Davenport. This particular model was chosen for three reasons: First, the importance of aluminum as a material makes the development of an accurate model potential for simulation purposes desirable. Second, the large number of parameters and complicated nature of the embedded atom potential increases the need for a systematic, as opposed to ad hoc procedure for adjusting the melting point. In addition, the zero-pressure melting point for the Mei-Davenport potential has been previously determined to be 800$`\pm 9`$K somewhat lower than the experimental melting point of aluminum at 933K. (Note that, the melting point determined by Mei and Davenport was calculated for a 256 particle system - the actual value for this potential in the thermodynamic limit is slightly higher at 826.4$`\pm `$1.3K.) We find that, in this specific case, only one of the parameters of the potential has any significant effect on the potential and that changing this parameter according to the Gibbs-Duhem procedure yields a new model with the correct experimental melting point with no significant change in the quality of the agreement of the quantities with respect to which the model was originally optimized. Details of this calculation as well as a description of the model can be found in Section II below.

## II Gibbs-Duhem Integration and Melting Point Optimization

The technique of Gibbs-Duhem integration has been well described previously by Kofke, but in the interests of completeness and clarity we repeat the basic derivation here. Consider a single-component system with an arbitrary interaction potential, $`𝒰(\{𝐑_𝐢\},\{X_i\})`$, where the $`𝐑_𝐢`$ are the atomic coordinates and the $`X_i`$ are the parameters that define the potential - no restriction to pairwise additivity need be assumed. Assume there are two phases $`\alpha `$ and $`\beta `$ in coexistence at a temperature $`T`$ and pressure $`P`$. On the surface of coexistence, the chemical potentials (molar Gibbs free energies) of the two phases must be equal. To quantify how changes in $`P`$,$`T`$, and $`X_i`$ will affect the chemical potential one can define a generalized Gibbs-Duhem equation,

$$d\mu =sdT+vdP+\underset{i}{}\lambda _idX_i$$
(1)
where $`\mu `$ is the chemical potential, $`s`$ and $`v`$ are the entropy and volume per particle, respectively, and the $`\lambda _i`$ are generalized thermodynamic variables conjugate to the potential parameters, $`X_i`$, defined as

$$N\lambda _i\left(\frac{G}{X_i}\right)_{T,P,X_{j,ji}}$$
(2)
Now as one moves infinitesimally away from the original coexistence point $`(P,T,\{\lambda _i\})`$ to another point $`(P+dP,T+dT,\{\lambda _i+d\lambda _i\})`$ on the surface of coexistence, the change in $`\mu `$ must be identical in both phases. This condition together with Eq. 1 gives

$`\mu _\alpha \mu _\beta `$ $`=`$ $`(s_\alpha s_\beta )dT+(v_\alpha v_\beta )dP`$ (4)
$`+{\displaystyle \underset{i}{}}(\lambda _{\alpha ,i}\lambda _{\beta ,i})dX_i=0,`$
where $`\mu _\alpha `$ and $`\mu _\beta `$ are the chemical potentials for each of the respective phases. Assuming constant pressure ($`dP=0`$), since we are interested here in changes in the transition temperature, the previous equation can be rearranged to give

$$\left(\frac{T}{X_i}\right)_{P,X_{j;ji};\text{coex}}=\frac{T(\lambda _{\alpha i}\lambda _{\beta i})}{h_\alpha h_\beta }=\frac{T\mathrm{\Delta }\lambda _i}{\mathrm{\Delta }h}$$
(5)
where we have also assumed that at coexistence, $`\mathrm{\Delta }s=\mathrm{\Delta }h/T`$, where $`\mathrm{\Delta }h`$ is the latent heat per particle for the phase transition. (Note that, the corresponding equations for the $`\left(\frac{P}{\lambda _i}\right)_{T,\lambda _{j;ji};\text{coex}}`$ can be easily obtained by replacing $`\mathrm{\Delta }h/T`$ in Eq. 5 with $`\mathrm{\Delta }v`$.)

The $`\lambda _i`$ can be related to mechanical averages that can be easily calculated in a molecular dynamics or Monte Carlo simulation. First, the Gibbs free energy is related to the isothermal-isobaric distribution, $`\mathrm{\Delta }(N,P,T)`$ , as follows.

$$G=k_bT\mathrm{ln}\mathrm{\Delta }(N,P,T)$$
(6)
which for a classical system with interaction potential $`𝒰(\{𝐫_i\})`$ is given by

$$\mathrm{\Delta }(N,P,T)=\frac{1}{\mathrm{\Lambda }^{3N}N!}_0^{\mathrm{}}𝑑Vd^N𝐫\mathrm{exp}(\beta 𝒰\beta PV),$$
(7)
where $`V`$ is the volume. Taking the derivative of Eq. 6 with respect to the parameter $`X_i`$ gives

$`\left({\displaystyle \frac{G}{X_i}}\right)_{T,P,X_{j;ji}}`$ $`=`$ $`kT\left({\displaystyle \frac{\mathrm{ln}\mathrm{\Delta }}{X_i}}\right)_{T,P,X_{j;ji}}`$ (8)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }}}{\displaystyle _0^{\mathrm{}}}𝑑V{\displaystyle d^N𝐫\left(\frac{𝒰}{X_i}\right)_{X_j;ji}}`$ (10)
$`\times \mathrm{exp}(\beta 𝒰\beta PV).`$
Using Eq. 2 we have

$$\lambda _i=\left(\frac{dG}{dX_i}\right)_{T,P,X_{j;ji}}=\left(\frac{𝒰}{X_i}\right)_{X_{j;ji}}_{N,P,T}$$
(11)
Using Eq. 5 and Eq. 11 we can calculate the change necessary in each of the parameters $`\{X_i\}`$, to alter the melting temperature of our interaction potential by some arbitrary amount. As long as the changes are not too large, the calculation can be greatly simplified by the assumption that the deviation remains linear.

$$T_m(\{X_i\})=T_{m,0}+\underset{i}{}\left(\frac{T}{X_i}\right)_{X_{j,ji},\text{coex}}(X_iX_{i,0}),$$
(12)
where $`T_{m,0}`$ is the melting point of the original model with parameters $`\{X_{i,0}\}`$. (If the linear approximation does not hold - and it should always be checked - one could integrate the differential equation using an appropriate numerical technique.) From a single simulation, one could calculate the necessary change in all the different parameters of the potential corresponding to a particular change in the melting temperature. One could then, use this information to construct a cost function to re-optimize the potential that includes the previously used set $`\{A_j^{\text{exp.}}\}`$ of experimentally determined (mechanical) properties as well as the experimental melting temperature $`T_m^{\text{exp.}}`$ as optimization targets. For example, a linear least squares procedure could be utilized using appropriately chosen set of weight functions $`w_i`$:

$`f(\{X_i\})`$ $`=`$ $`w_T[T_m^{\text{exp.}}T_m(\{X_i\})]^2`$ (14)
$`+{\displaystyle \underset{j}{}}w_i[A_j^{\text{exp.}}A_j(\{X_i\})]^2.`$
(In the example application discussed in the next section, we found that only one of the several parameters of the potential had any significant effect on the melting point, which greatly simplified the re-parameterization by making it possible to modify the potential without having to directly use such a cost function. It is not expected, however, that this will be true, in general.)

## III Results for a model of Aluminum: determining the original melting point

As an application of the method we examine an embedded atom model  for aluminum developed by Mei and Davenport . The parameters in their potential were fit in order to optimize the potential with respect to a variety of experimental properties such as the cohesive energy ($`E_c`$), lattice constant ($`\sqrt{2}r_0`$), un-relaxed vacancy-formation energy, and elastic constants of the static fcc crystal at zero temperature. The total potential energy has the form

$$𝒰=\underset{i}{}F(\rho _i)+\frac{1}{2}\underset{i,ji}{}\varphi (r_{\mathrm{i}\mathrm{j}}),$$
(15)
where

$`F(\rho )`$ $`=`$ $`E_c\left[1{\displaystyle \frac{\alpha }{\beta }}\mathrm{ln}\left({\displaystyle \frac{\rho }{\rho _e}}\right)\right]\left({\displaystyle \frac{\rho }{\rho _e}}\right)^{\alpha /\beta }`$ (19)
$`+{\displaystyle \frac{1}{2}}\varphi _0{\displaystyle \underset{m=1}{\overset{3}{}}}s_mexp[(\sqrt{m}1)\gamma ]`$
$`\times \left[1+(\sqrt{m}1)\delta \sqrt{m}\left({\displaystyle \frac{\delta }{\beta }}\right)\mathrm{ln}\left({\displaystyle \frac{\rho }{\rho _e}}\right)\right]`$
$`\times \left({\displaystyle \frac{\rho }{\rho _e}}\right)^{\sqrt{m}\gamma /\beta },`$
and

$`\rho _i`$ $`=`$ $`{\displaystyle \underset{j(i)}{}}f(r_{ij}),`$ (20)
$`f(r)`$ $`=`$ $`\rho _e{\displaystyle \underset{l=0}{\overset{5}{}}}{\displaystyle \frac{c_l}{12}}\left({\displaystyle \frac{r_0}{r}}\right)^l,`$ (21)
$`\varphi (r)`$ $`=`$ $`\varphi _0\left[1+\delta \left({\displaystyle \frac{r}{r_0}}1\right)\right]exp\left[\gamma \left({\displaystyle \frac{r}{r_0}}1\right)\right],`$ (22)
and is used in conjunction with the following cutoff function:

$`q(r)=\{\begin{array}{cc}\hfill 1& rr_n,\hfill \\ \hfill (1x)^3(1+3x+6x^2)& r_n<r<r_c,\hfill \\ \hfill 0& rr_c\hfill \end{array}`$ (26)
$`x=(rr_n)/(r_cr_n)`$ (27)
with a cutoff distance $`(r_c)`$ between the third and fourth neighbor shells of a static fcc crystal. Using values of $`r_n=1.75r_0`$ and $`r_c=1.95r_0`$, the functions $`q(r)f(r)`$ and $`q(r)\varphi (r)`$ go smoothly to zero at $`r=r_c`$. The parameters of the potential  are given in Table I. For this embedded atom model, mass is measured in amu, distance in Å, and energy in eV. The natural simulation time unit is calculated to be 10.181fs.

To begin, the melting temperature of the original model (at $`P=0`$) must be determined. Mei and Davenport perform a calculation of the free energies of the aluminum melt and fcc crystal using thermodynamic integration of the Gibbs free energy. For a 256 particle system, they obtain a melting temperature of 800$`\pm `$9K.

In repeating their melting point determination, we found that the value they obtained is not quite correct, due primarily to the small system size studied and a problem in the choice of thermodynamic integration path for the liquid phase. Our melting point determination was performed using the same basic methodology as Mei and Davenport, described below.

The Gibbs free energies per particle of the liquid and crystal as a function of temperature at constant pressure can be obtained by thermodynamic integration using the integral form of the Gibbs-Helmholtz equation

$$\frac{g(T,P)}{T}=\frac{g(T_0,P)}{T_0}_{T_0}^T\frac{e(\tau ,P)}{\tau ^2}𝑑\tau ,$$
(28)
where $`T_0`$ is a predetermined reference temperature and $`e(T,P)`$ is the average total energy per particle. The Gibbs free energy at the reference temperature must be obtained separately by thermodynamic integration from a suitable ideal reference state. To do this, the interaction potential is parameterized along a linear path between that of the reference potential, $`𝒰_0`$ and the full potential $`𝒰`$

$$𝒰(\xi )=\xi 𝒰+(1\xi )𝒰_0,$$
(29)
The Gibbs free energy per particle relative to that of the reference system can then be obtained by thermodynamic integration along the path:

$`g(T,P)`$ $`g(T,P;\xi =1)`$ (31)
$`=g(T_0,P;\xi =0)+{\displaystyle _0^1}𝑑\xi {\displaystyle \frac{e(\xi )}{\xi }}_\xi .`$
For the crystal, the reference system chosen is that of an Einstein crystal

$$𝒰_E(\{r_i\})=\frac{1}{2}m\omega _D^2\underset{i}{}(𝐫_i𝐫_{i0})^2,$$
(32)
where the $`\{𝐫_{i0}\}`$ are the ideal crystal lattice positions, $`m`$ is the mass and $`\omega _D`$ is the Debye frequency, which to minimize the difference between the reference and full system, should be chosen to give a similar mean squared displacement for the atoms at the temperature of interest. For this system at 296K, the optimum Debye frequency corresponds to a Debye temperature ($`T_D=\mathrm{}\omega _D/k`$) of 207K. Note that this frequency is different than the one used by Mei and Davenport.

The reference system for the liquid phase is an ideal gas, but the transformation must performed as a two step process in order to avoid the liquid/gas phase transition. Mei and Davenport follow the reversible expansion method used by Broughton and Li. The two step process is accomplished by first turning off the attractive part of the potential followed a volume expansion to reach the ideal gas limit. In step one, it is extremely important to turn off the attractive part of the potential in a way that will not drastically alter the effective radius of the potential. If care is not taken, the system will freeze during the integration of this step. Here, we write the interatomic potential as a linear interpolation between the actual potential and the reference system

$`\varphi (r;\zeta )`$ $`=`$ $`\varphi _{rep}(r)+\zeta \varphi _{att}(r).`$ (33)
As the value of $`\zeta `$ is varied from 1 to 0, the system is transformed from the original Mei and Davenport potential to that of a purely repulsive potential. Step two of the integration is a volume expansion. The free energy change in this step is given by

$$\mathrm{\Delta }g_{\text{step 2}}=k_bT_0_0^\rho \frac{d\rho ^{}}{\rho ^{}}\left[\frac{\beta P}{\rho ^{}}1\right].$$
(34)
The potential splitting that Mei and Davenport use has the problem in that their repulsive potential $`\varphi _{rep}=\varphi _0\delta \mathrm{exp}\left[\gamma \left(\frac{r}{r_0}1\right)\right]`$ has an effective radius that is much larger than the actual potential (See Fig. 1a) so that as the rest of the potential is turned off the system freezes as the effective packing fraction increases. To remedy this we use a WCA splitting where the repulsive part of the potential is equal to the the potential energy for radii less than the radius at the minimum of the potential and zero for larger radii. This prescription (illustrated in Fig. 1) gives an effective radius more similar to the full system and avoids the freezing transition.

In order to obtain energy curves needed in Eq. 28 as a function of temperature, molecular dynamics (MD) simulations were conducted at several different temperatures using the Nosé-Poincaré-Anderson algorithm  for isothermal-isobaric molecular dynamics. The fcc crystal was simulated at 50K intervals from 296K to 946K, while the liquid was studied over a smaller temperature range from 762K to 1152K using 30K intervals. All of the isothermal-isobaric MD simulation runs were equilibrated for 100,000 steps and sampled for 300,000 steps. From these simulations we obtained the average energy and density for both the crystal and liquid as a function of temperature and system size. Both the energies and densities were fit to second order polynomials. The coefficients for these polynomials are shown in Table II. These polynomials were used in the construction of the free energy curves as described in Eq. 28.

Using the density obtained from the constant $`NPT`$ simulations, the Gibbs free energy of the fcc crystal at $`T_0=296K`$ was calculated by running several simulations at constant N,V,T using the Einstein crystal reference state to evaluate the integrand of Eq. 31 at values of $`\xi `$ corresponding to a ten-point Gauss-Legendre quadrature. The simulations were performed using the Nosé-Poincaré algorithm . Figure 2 shows a plot of the integrand vs. $`\xi `$ for all of the system sizes studied. At a value of $`\xi =1`$ the system is governed solely by the

embedded atom potential. As the value of $`\xi `$ goes to zero, the potential is gradually changed to that of an Einstein crystal. For each value of $`\xi `$, the system was equilibrated for 100,000 steps and sampled for 100,000 steps.

Simulations for step one of the liquid free-energy at a temperature of $`1092K`$ were performed in the same manner as for the crystal. Here the attractive part of the potential (Eq. 33) was slowly turned off as the value $`\xi `$ was changed from one to zero. This is shown graphically in Figure 3(a). As with the solid, ten-point Gauss-Legendre integration was used to numerically compute this integral. The second step of the liquid free-energy calculation included a series of constant $`NVT`$ simulations at decreasing densities starting with the repulsive potential system from the conclusion of step one. Average values for the step two integrand (Eq. 34) are shown graphically in Figure 3(b). The simulations for both steps of the liquid free energy integrations were equilibrated for 100,000 steps and sampled for 100,000 steps at each value of the integrand. The numerical integration for step two was performed using the ten-point Simpsons quadrature. This method was chosen over the Gauss-Legendre quadrature due to the inaccuracy of the sampling at low density. The value of the integrand at zero density was obtained by an analytical calculation of the second virial coefficient ($`B_2`$).

Free energies and the melting point are shown in Table III as a function of particle number. A graph of the

melting temperature vs. $`1/N`$ (Figure 4) shows that at infinite particle number the melting point for the embedded atom potential proposed by Mei and Davenport approaches $`826.4K\pm 1.3K`$. It should be noted that the major error in the melting point calculated by Mei and Davenport for their potential 800$`\pm 9K`$ is primarily due to the small size of their system. The problems with their potential splitting appears to have had little effect, probably due to cancellation of errors, as the melting point that we determine here for the 256 particle system agrees with theirs within the simulation error. Recently, Morris, Wang, Ho and Chan  argued on the basis of the stability of crystal-liquid interfaces that the melting point of the Mei-Davenport aluminum potential is actually significantly lower (around 725K). Our results do not support this conclusion and as a check we have carefully set up stable stress-free interfaces at our calculated melting point. The systems exhibit melting (freezing) as the temperature is raised (lowered) away from our calculated melting point. The lower temperature transitions of Morris, et al. were most likely due to significant un-relaxed stress in the crystal.

## IV Results for a model of Aluminum: Re-parameterizing the potential

To adjust the melting temperature of the embedded atom potential for aluminum, we first calculate $`\frac{T_m}{X_i}`$ for each of the parameters in the potential. (Note, for simplicity and consistency the expansion coefficients, $`c_l`$, in Eq. 21 were kept constant). The derivatives were calculated in a single (constant $`NPT`$) simulation at the melting temperature calculated in the previous section with $`P=0`$, for each of several system sizes ($`N=`$256, 500, 864, 2048 and 4000). The system was equilibrated for 100,000 steps followed by 100,000 steps for averageing. From this simulation it was determined that $`\frac{T_m}{X_i}`$ was significant only for the parameter $`\delta `$ \- the other parameters of the potential have little effect on the melting point. The complicated nature of the potential makes it difficult to assign any physical explanation to the sensitivity of the melting point to $`\delta `$ relative to the other parameters.

Next, a series of simulations were performed for each system size to integrate along the coexistence curve from the initial calculated melting point of the potential to the true experimental melting temperature. At each temperature along the coexistence curve, the system was equilibrated for 100,000 steps and sampled over 300,000 steps. In our experiments we use a fourth-order predictor corrector integrator to carry out the integration along the coexistence curve as a function of $`\delta `$. Fig 5 shows this integration graphically with final values of $`\delta `$ corresponding to a melting temperature of $`T_m=933K`$. These values are listed in Table III. (Note that, although we did an accurate numerical integration along the coexistence

curve, the results indicate that using the approximation that the deriviative is a constant in the region of interest would have been correct to within the simulation error.)

In order to confirm that the melting point did indeed change as expected, the melting point calculation was repeated using the newly calculated value of $`\delta `$. The new melting temperature for the embedded atom potential with a value of $`\delta =8.45`$ corresponding to a system size of N=4000, was calculated to be $`931.5K\pm 1.5K`$. Melting temperatures for the other system sizes are listed in Table III. The experimental value for the melting temperature of aluminum is 933.47K.

Mei and Davenport initially determined $`\delta `$ (and the other parameters in the potential) by fitting the potential to certain physical constants. For the new value of $`\delta `$ (8.45), we have re-calculated a variety of physical properties of aluminum. These quantities, for the original potential, the reparameterized potential and the experimental values are collected in Table IV. From this data,

we see that, in comparison to the original potential with $`\delta =7.35`$, the new potential more closely models the experimental values  of the ($`T=0K`$) elastic constants $`C_{11}`$ and $`C_{22}`$ while $`C_{44}`$ becomes slightly worse in comparison to its experimentally determined value. In addition, mostly as a consequence of the improved melting point, the latent heat is considerably improved.

## V Summary

We have outlined an application of the Gibbs-Duhem integration method of Kofke, with which a model interaction potential can be re-parameterized, including the experimental melting point in the optimization protocol. The melting temperature of a potential can then be adjusted similar to the tuning of other parameters in the potential. Since non-mechanical properties, such as the melting point, are not generally included in potential optimization and the agreement of such quantities with experiment is not guaranteed, such a procedure will be useful in situations, such as in the simulation of crystal-melt interfaces, where the having the correct melting point is highly desirable. The method is general and can easily be extended to a variety of systems. As an example of the utility of the method, we apply the procedure to reparameterizing a popular model of aluminum for which the melting point has been calculated to be over 100K below the accepted experimental value. We demonstrate that the reparamterized potential has a melting point that agrees within the statistical error with the experimental value of 933K and that reparameterization does not degrade the quality of the potential with respect a variety of properties, and that, in fact, quantities such as the elastic constants $`C_{11}`$ and $`C_{12}`$ and the latent heat, agreement is improved.

## VI Acknowledgements

We gratefully acknowledge J. W. Davenport for helpful conversations, as well as the Kansas Center for Advanced Scientific Computing for the use of their computer facilities. We also would like to thank the National Science Foundation (under grant CHE-9500211) as well as the University of Kansas General Research Fund for their generous support of this research.


text:
# 1 Introduction

## 1 Introduction

In studying the Mikheyev-Smirnov-Wolfenstein (MSW) effect  due to the coherent forward scattering of neutrinos by electrons in matter, it is often instructive to consider first special cases where the electron density is taken to be a simple function of distance. It is the purpose of the present paper to investigate perhaps the simplest case: the case where the electron density is a linear function of distance.

The problem of the linear electron density is formulated in Sec. 2. The case of two neutrino species has a long history , and the solution, as reviewed in Sec. 3, can be expressed in terms of parabolic cylinder functions—see, for example, Chapter VIII of , or equivalently confluent hypergeometric functions. However, the solution in this form is specific to the case of two neutrino species, and is not convenient for generalizations to more neutrino species. Physically, this generalization is essential because there are at least three types of neutrinos. Therefore, Sec. 4 is devoted to treating in a different way the MSW differential equations for linear electron density and two neutrino species. On the one hand, this alternative method must lead to the same solutions as those in Sec. 3; on the other hand, this new treatment can be readily generalized to any number of neutrino species. This case of linear electron density but any number of neutrino species forms the main content of the present paper, and various aspects of this case are treated in Sec. 5 and Sec. 6.

## 2 Formulation of the Problem

Let there be $`N`$ types of neutrinos, denoted by $`\nu _1`$, $`\nu _2`$, $`\mathrm{}`$ $`\nu _N`$, where $`\nu _1`$ is the neutrino of the first generation, i.e., the one that forms the SU(2) doublet with the electron. It is assumed that $`\nu _1`$ is the only neutrino which interacts differently with the electron because of the exchange of the intermediate boson $`W`$, while the others neutrinos $`\nu _2`$, $`\nu _3`$, $`\mathrm{}`$ $`\nu _N`$ all have the same interaction with the electron. Thus, the neutrino mass matrix $`M`$  is an $`N\times N`$ matrix. The eigenvalues of $`M`$ give the $`N`$ neutrino masses $`\mu `$.

In analyzing the MSW effect, the neutrino masses are usually taken to be much smaller than the momentum $`p`$ of the neutrino. Under this assumption, because

$$(p^2+\mu ^2)^{1/2}p+\frac{1}{2p}\mu ^2,$$
(2.1)
it is $`M^2`$ that enters in the differential equation for the MSW effect. Let $`\mathrm{\Psi }(x)`$ be the $`N`$-component neutrino wave function, then this differential equation is

$$i\frac{d}{dx}\mathrm{\Psi }(x)=\left[W(x)+\frac{1}{2p}M^2\right]\mathrm{\Psi }(x)$$
(2.2)
where $`W(x)`$ is an $`N\times N`$ matrix whose only non-zero element is

$$\left[W(x)\right]_{11}=\sqrt{2}G_\mathrm{F}N_e(x),$$
(2.3)
with $`G_\mathrm{F}`$ the Fermi weak-interaction constant and $`N_e(x)`$ the electron density at the point $`x`$.

The terminology “linear electron density” is used to mean that $`N_e(x)`$ is a linear function of $`x`$.

Since $`N_e(x)`$ is the density of electrons, it cannot be negative. Therefore, the MSW differential equation (2.2) is physically meaningful only for the half-line of $`x`$ where $`N_e(x)0`$. On the other hand, when the neutrino or the electron, but not both, is replaced by its antiparticle, the quantity $`\left[W(x)\right]_{11}`$ of Eq. (2.3) changes sign. Therefore, the complementary half-line of $`x`$ describes this slightly different physical situation. For this reason, Eq. (2.2) is to be studied for the entire range of $`x`$, from $`\mathrm{}`$ to $`+\mathrm{}`$.

For the present case of the linear electron density, Eq. (2.2) can be reduced, for a given value of $`p`$, to the dimensionless standard form

$$i\frac{d}{dt}\psi (t)=A(t)\psi (x),$$
(2.4)
where

$$\psi (t)=\left[\begin{array}{c}\psi _1(t)\\ \psi _2(t)\\ \psi _3(t)\\ \mathrm{}\\ \psi _N(t)\end{array}\right]$$
(2.5)
and

$$A(t)=\left[\begin{array}{ccccc}t& a_2& a_3& \mathrm{}& a_N\\ a_2& b_2& 0& \mathrm{}& 0\\ a_3& 0& b_3& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}\\ a_N& 0& 0& \mathrm{}& b_N\end{array}\right].$$
(2.6)
This is accomplished as follows.

(i) To change the independent variable from $`x`$ to $`t`$, there is a shift in origin and a rescaling with possibly a reversal of sign.

(ii) To change the dependent variable from $`\mathrm{\Psi }`$ to $`\psi `$, there is a rotation in the second to $`N`$th component and an introduction of exponential factors with possibly some minus signs.

Furthermore, from (i) and (ii), the elements of the matrix $`A(t)`$ can be chosen to satisfy the conditions

$$\underset{j=2}{\overset{N}{}}b_j=0,$$
(2.7)
$$b_2b_3b_4\mathrm{}b_{N1}b_N,$$
(2.8)
and

$$a_j0\text{for}j=2,3,4,\mathrm{}N.$$
(2.9)
Consider the following special cases.

(a) If, for some $`j`$, say $`j_0`$, $`a_{j_0}=0`$, then it is seen from Eqs. (2.5) and (2.6) that $`\psi _{j_0}`$ is decoupled from the other $`\psi _j`$’s. Thus, this special case of $`N`$ types of neutrinos is reduced to a problem with $`N1`$ types of neutrinos.

(b) If, again for some $`j`$, say $`j_0`$, $`b_{j_0}=b_{j_0+1}`$, then a rotation can be carried out between $`\psi _{j_0}`$ and $`\psi _{j_0+1}`$ such that, after this additional rotation, the new $`a_{j_0}`$ is zero. Thus, this special case of $`b_{j_0}=b_{j_0+1}`$ can be reduced to the above case of $`a_{j_0}=0`$, and hence again this second special case of $`N`$ types of neutrinos is reduced to a problem with $`N1`$ types of neutrinos.

It is therefore sufficient to study the ordinary differential equation (2.4) with Eqs. (2.5) and (2.6) under the condition (2.7) together with

$$b_2<b_3<b_4<\mathrm{}<b_{N1}<b_N$$
(2.10)
and

$$a_j>0\text{for}j=2,3,4,\mathrm{}N.$$
(2.11)
In view of the inequality (2.10), it turns out to be convenient to define symbolically

$$b_1=\mathrm{}\text{and}b_{N+1}=+\mathrm{}.$$
(2.12)
## 3 Case $`N=2`$

Let us review first the well-known case of the MSW effect for two types of neutrinos . By Eqs. (2.2)–(2.7), the MSW equations are

$$i\frac{d}{dt}\left[\begin{array}{c}\psi _1(t)\\ \psi _2(t)\end{array}\right]=\left[\begin{array}{cc}t& a_2\\ a_2& 0\end{array}\right]\left[\begin{array}{c}\psi _1(t)\\ \psi _2(t)\end{array}\right],$$
(3.1)
or more explicitly

$`i{\displaystyle \frac{d}{dt}}\psi _1(t)`$ $`=`$ $`t\psi _1(t)+a_2\psi _2(t),`$ (3.2)
$`i{\displaystyle \frac{d}{dt}}\psi _2(t)`$ $`=`$ $`a_2\psi _1(t),`$ (3.3)
with $`a_2>0`$. A second-order ordinary differential equation for $`\psi _1(t)`$ is obtained by applying $`d/dt`$ to Eq. (3.2) and using Eq. (3.3):

$$\frac{d^2\psi _1(t)}{dt^2}it\frac{d\psi _1(t)}{dt}+(a_2^2i)\psi _1(t)=0.$$
(3.4)
In order to remove the first-derivative term, let

$$\psi _1(t)=e^{it^2/4}\varphi _1(t)$$
(3.5)
Then the equation for $`\varphi _1(t)`$ is

$$\frac{d^2\varphi _1(t)}{dt^2}+(\frac{1}{4}t^2+a_2^2\frac{1}{2}i)\varphi _1(t)=0$$
(3.6)
Two linearly independent solutions of this Eq. (3.6) are the parabolic cylinder functions

$$D_\rho (\pm e^{i\pi /4}t),$$
(3.7)
where

$$\rho =ia_2^21.$$
(3.8)
Parabolic cylinder functions are special cases of the confluent hypergeometric function , the relation being

$$D_\rho (z)=2^{(\rho 1)/2}e^{z^2/4}z\mathrm{\Psi }(\frac{1}{2}\frac{1}{2}\rho ,\frac{3}{2};\frac{1}{2}z^2).$$
(3.9)
Since the confluent hypergeometric functions $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ satisfy the same second-order differential equation, the general solution of Eq. (3.6) is

$$\psi _1(t)=t\left[C\mathrm{\Phi }(1+\frac{1}{2}ia_2^2,\frac{3}{2};\frac{1}{2}it^2)+C^{}\mathrm{\Psi }(1+\frac{1}{2}ia_2^2,\frac{3}{2};\frac{1}{2}it^2)\right].$$
(3.10)
This is one convenient form for the solution for $`N=2`$.

## 4 Case $`N=2`$—an Alternative Approach

In the existing treatment in the literature for linear electron density and two types of neutrinos  as reviewed in Sec. 3, the crucial step is to recognize that the second-order differential equation (3.4) can be solved exactly in terms of known higher transcendental functions, either parabolic cylinder functions or confluent hypergeometric functions. More generally, for $`N`$ types of neutrinos, the corresponding differential equation is of $`N`$th order. Even for $`N=3`$, the third-order differential equation is not one for any well-known transcendental function. Therefore, in order to be able to generalize the treatment of $`N=2`$ to larger values of $`N`$, we must recast the solution of Sec. 3 so that parabolic cylinder functions and confluent hypergeometric functions do not play an essential role.

A useful question to ask is the following: In what way is the linear electron density especially simple? The answer must be sought in Eq. (2.6), from which it is seen that the independent variable $`t`$ appears only in one matrix element, and furthermore, it appears only linearly in that element. This implies that, if Fourier transform is applied to the differential equation (2.4), the differentiation with respect to the Fourier-transform variable appears only once. Hence it is expected that an explicit expression can be obtained for the Fourier transform of $`\psi `$.

Let

$$F(\zeta )=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑te^{i\zeta t}\psi _1(t),$$
(4.1)
then it follows from Eq. (3.4) that $`F(\zeta )`$ satisfies the first-order differential equation

$$\zeta ^2F(\zeta )\frac{d}{d\zeta }[i\zeta F(\zeta )]+(a_2^2i)F(\zeta )=0,$$

where we have omitted all terms from $`t=\pm \mathrm{}`$. This differential equation simplifies immediately to

$$i\zeta \frac{dF(\zeta )}{d\zeta }(\zeta ^2a_2^2)F(\zeta )=0,$$

or

$$\frac{1}{F(\zeta )}\frac{dF(\zeta )}{d\zeta }=\frac{i}{\zeta }(a_2^2\zeta ^2).$$
(4.2)
Integration over $`\zeta `$ gives

$$F(\zeta )=\text{const.}e^{i\zeta ^2/2}\zeta ^{ia_2^2}.$$
(4.3)
From the inequality (2.11), it is seen that the function on the right-hand side of Eq. (4.3) has a singularity at

$$\zeta =0=b_2.$$
(4.4)
Therefore the constant in (4.3) can take on different values for $`\zeta `$ positive and for $`\zeta `$ negative. In other words, the differential equation (4.2) is really two differential equations, one for $`\zeta >0`$ and the other for $`\zeta <0`$, consistent with the fact that the right-hand side of Eq. (4.2) has a singularity at $`\zeta =0`$. With this observation, it is natural to define

$$F_1(\zeta )=\{\begin{array}{cc}e^{i\zeta ^2/2}(\zeta )^{ia_2^2}\hfill & \text{for }\zeta <0,\hfill \\ 0\hfill & \text{for }\zeta >0,\hfill \end{array}$$
(4.5a)
and

$$F_2(\zeta )=\{\begin{array}{cc}0\hfill & \text{for }\zeta <0,\hfill \\ e^{i\zeta ^2/2}\zeta ^{ia_2^2}\hfill & \text{for }\zeta >0.\hfill \end{array}$$
(4.5b)
Inverting the Fourier transform (4.1), this choice leads to

$`\psi _1^{(1)}(t)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^0}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}(\zeta )^{ia_2^2},`$
$`\psi _1^{(2)}(t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}\zeta ^{ia_2^2}.`$ (4.6)
With the notation (2.12), these two formulas (4) can be written as

$$\psi _1^{(n)}(t)=_{b_n}^{b_{n+1}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}|\zeta |^{ia_2^2},$$
(4.7)
for $`n=1,2`$.

It remains to show that both $`\psi _1^{(1)}(t)`$ and $`\psi _1^{(2)}(t)`$ are confluent hypergeometric functions of the correct parameters and argument. For this purpose, it is convenient to define

$`\psi _c(t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\zeta \mathrm{cos}(\zeta t)e^{i\zeta ^2/2}\zeta ^{ia_2^2},`$
$`\psi _s(t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\zeta \mathrm{sin}(\zeta t)e^{i\zeta ^2/2}\zeta ^{ia_2^2},`$ (4.8)
so that it follows from Eqs. (4) that

$`\psi _1^{(1)}(t)`$ $`=`$ $`\psi _c(t)+i\psi _s(t),`$
$`\psi _1^{(2)}(t)`$ $`=`$ $`\psi _c(t)i\psi _s(t).`$ (4.9)
It is found that

$$\psi _c(t)=e^{i\pi /4}e^{\pi a_2^2/4}\mathrm{\hspace{0.17em}2}^{(1+ia_2^2)/2}\mathrm{\Gamma }(\frac{1}{2}+\frac{1}{2}ia_2^2)\mathrm{\Phi }(\frac{1}{2}+\frac{1}{2}ia_2^2,\frac{1}{2};\frac{1}{2}it^2)$$
(4.10)
and

$$\psi _s(t)=ie^{\pi a_2^2/4}\mathrm{\hspace{0.17em}2}^{ia_2^2/2}\mathrm{\Gamma }(1+\frac{1}{2}ia_2^2)t\mathrm{\Phi }(1+\frac{1}{2}ia_2^2,\frac{3}{2};\frac{1}{2}it^2).$$
(4.11)
There are various ways to verify Eqs. (4.10) and (4.11), including carrying out power series expansions in $`t`$ for the left-hand and right-hand sides.

Finally, we note from Eq. (7) on p. 257 of reference  that

$$\mathrm{\Psi }(a,c;x)=\frac{\mathrm{\Gamma }(1c)}{\mathrm{\Gamma }(ac+1)}\mathrm{\Phi }(a,c;x)+\frac{\mathrm{\Gamma }(c1)}{\mathrm{\Gamma }(a)}x^{1c}\mathrm{\Phi }(ac+1,2c;x).$$
(4.12)
Therefore, the results of Sec. 3 and this section are the same.

## 5 General values of $`N`$

The procedure of Sec. 4 for $`N=2`$ can be generalized in a straightforward way to larger values of $`N`$. Indeed, this is the major advantage over the previously known ones as reviewed in Sec. 3. This generalization to arbitrary values of $`N`$ is to be carried out in this section. Thus, the differential equations (2.4) need to be solved under the constraints (2.7), (2.10), and (2.11).

By Eqs. (2.5) and (2.6), the Eqs. (2.4) are more explicitly

$$i\frac{d\psi _1(t)}{dt}=t\psi _1(t)+\underset{j=2}{\overset{N}{}}a_j\psi _j(t)$$
(5.1)
and, for $`k=2,3,4\mathrm{}N`$,

$$\left(i\frac{d}{dt}b_k\right)\psi _k(t)=a_k\psi _1(t).$$
(5.2)
In order to get a differential equation for $`\psi _1(t)`$, apply the operator

$$\underset{k=2}{\overset{N}{}}\left(i\frac{d}{dt}b_k\right)$$

to Eq. (5.1). By Eq. (5.2), this gives

$$\left[\underset{k=2}{\overset{N}{}}\left(i\frac{d}{dt}b_k\right)\right]\left(i\frac{d}{dt}+t\right)\psi _1(t)=\underset{j=2}{\overset{N}{}}a_j^2\underset{\begin{array}{c}k=2\\ kj\end{array}}{\overset{N}{}}\left(i\frac{d}{dt}b_k\right)\psi _1(t).$$
(5.3)
Eq. (5.3) is an $`N`$th-order ordinary differential equation for $`\psi _1(t)`$; it reduces to Eq. (3.4) when $`N=2`$.

Following Sec. 4, define the Fourier transform $`F(\zeta )`$ of $`\psi _1(t)`$ by Eq. (4.1), then the first-order differential equation for $`F(\zeta )`$ is

$$\left[\underset{k=2}{\overset{N}{}}(\zeta b_k)\right]\left(\zeta i\frac{d}{d\zeta }\right)F(\zeta )=\underset{j=2}{\overset{N}{}}a_j^2\underset{\begin{array}{c}k=2\\ kj\end{array}}{\overset{N}{}}(\zeta b_k)F(\zeta ),$$
(5.4)
or

$$\left(\zeta i\frac{d}{d\zeta }\right)F(\zeta )=\underset{j=2}{\overset{N}{}}\frac{a_j^2}{\zeta b_j}F(\zeta ),$$
(5.5)
or

$$\frac{1}{F(\zeta )}\frac{dF(\zeta )}{d\zeta }=i\left(\zeta +\underset{j=2}{\overset{N}{}}\frac{a_j^2}{\zeta b_j}\right).$$
(5.6)
This Eq. (5.6) is the generalization of the previous Eq. (4.2) for $`N=2`$. Integration over $`\zeta `$ gives the generalization of Eq. (4.3):

$$F(\zeta )=\text{const.}e^{i\zeta ^2/2}\underset{j=2}{\overset{N}{}}(\zeta b_j)^{ia_j^2}.$$
(5.7)
From (2.11), the function on the right-hand side of this Eq. (5.7) has singularities at

$$\zeta =b_j$$
(5.8)
for $`j=2,3,\mathrm{}N`$. Therefore, define for $`n=1,2,3\mathrm{}N`$

$$F_n(\zeta )=\{\begin{array}{cc}e^{i\zeta ^2/2}_{j=2}^N|\zeta b_j|^{ia_j^2}\hfill & \text{for }b_n<\zeta <b_{n+1},\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$$
(5.9)
Inverting the Fourier transform (4.1) then gives the desired $`N`$ linearly independent solutions of the differential equation (5.3) as

$$\psi _1^{(n)}(t)=_{b_n}^{b_{n+1}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}\underset{j=2}{\overset{N}{}}|\zeta b_j|^{ia_j^2}$$
(5.10)
for $`n=1,2,3,\mathrm{}N`$. In both Eq. (5.9) and Eq. (5.10), the notation of (2.12) has been used. The general solution of (5.3) is of course

$$\psi _1(t)=\underset{n=1}{\overset{N}{}}C_n\psi _1^{(n)}(t)$$
(5.11)
with arbitrary constants $`C_n`$.

The other components of the $`\psi (t)`$ of (2.5) can be easily obtained also, and the result is

$$\psi (t)=\underset{n=1}{\overset{N}{}}C_n\psi ^{(n)}(t),$$
(5.12)
where

$$\psi ^{(n)}(t)=_{b_n}^{b_{n+1}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}\left(\underset{j=2}{\overset{N}{}}|\zeta b_j|^{ia_j^2}\right)\left[\begin{array}{c}1\\ \frac{a_2}{\zeta b_2}\\ \frac{a_3}{\zeta b_3}\\ \mathrm{}\\ \frac{a_{N1}}{\zeta b_{N1}}\\ \frac{a_N}{\zeta b_N}\end{array}\right].$$
(5.13)
## 6 Limiting Behaviors for Large Distances

The next task is to obtain the limiting behaviors of the various components of the wave function when the distance $`x`$ is large, either positive or negative. In other words, the problem is to find the limiting behaviors of the $`\psi _j^{(n)}(t)`$, as given explicitly by Eq. (5.13), both for $`t\mathrm{}`$ and $`t\mathrm{}`$, with all the $`a`$’s and $`b`$’s fixed. It is important to remember that these two limits correspond to different physical problems, as discussed after Eq. (2.3).

The consideration here will be limited to the part of the asymptotic behavior that does not vanish as $`t\pm \mathrm{}`$. This is the physically interesting part. There are two possible types of contributions, from points of stationary phase and from end points of integration.

### 6.1 Points of Stationary Phase

From Eq. (5.13), the points of stationary phase are determined by

$$\frac{}{\zeta }(\zeta t\frac{1}{2}\zeta ^2)=0$$
(6.1)
or

$$\zeta =t.$$
(6.2)
In Eq. (6.1), the additional phase due to the factor

$$\underset{j=2}{\overset{N}{}}|\zeta b_j|^{ia_j^2}$$

is not included because the $`a_j`$ and $`b_j`$ are all fixed while $`t\pm \mathrm{}`$. Eq. (6.2) implies that this point of stationary point is relevant only to:

* $`\psi _1^{(1)}(t)`$ as $`t\mathrm{}`$, and
* $`\psi _1^{(N)}(t)`$ as $`t\mathrm{}`$

in view of Eq. (5.13). For example, when $`j>1`$, $`\psi _j^{(1)}(t)`$ as $`t\mathrm{}`$ and $`\psi _j^{(N)}(t)`$ as $`t\mathrm{}`$ both behave as $`1/t`$ in absolute value so far as the contribution from this point of stationary phase (6.2) is concerned.

### 6.2 End Points of Integration

It is seen from Eq. (5.13) that, when $`k2`$, there is an extra factor of

$$\frac{a_k}{\zeta b_k}$$
(6.3)
associated with $`\psi _k^{(n)}(t)`$. But the range of integration for this $`\psi _k^{(n)}(t)`$ as given by Eq. (5.13) is from $`b_n`$ to $`b_{n+1}`$. Therefore, the contribution from these end points of integration can lead to a non-zero answer only when the index $`k`$ appearing in the expression (6.3) agrees with either $`n`$ or $`n+1`$. In other words, there are non-zero contributions as $`t\pm \mathrm{}`$ only to $`\psi _k^{(k1)}(t)`$ \[i.e., $`n=k1`$\] and $`\psi _k^{(k)}(t)`$ \[i.e., $`n=k`$\]. These particular components are given by

$`\psi _k^{(k1)}(t)={\displaystyle _{b_{k1}}^{b_k}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}\left[{\displaystyle \underset{j=2}{\overset{k1}{}}}(\zeta b_j)^{ia_j^2}\right]\left[{\displaystyle \underset{j=k}{\overset{N}{}}}(b_j\zeta )^{ia_j^2}\right]{\displaystyle \frac{a_k}{b_k\zeta }},`$
$`\psi _k^{(k)}(t)={\displaystyle _{b_k}^{b_{k+1}}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}\left[{\displaystyle \underset{j=2}{\overset{k}{}}}(\zeta b_j)^{ia_j^2}\right]\left[{\displaystyle \underset{j=k+1}{\overset{N}{}}}(b_j\zeta )^{ia_j^2}\right]{\displaystyle \frac{a_k}{\zeta b_k}}.`$
These Eq. (6.2) are exact.

Since the important contributions come from the vicinity of $`\zeta =b_k`$, all the $`\zeta `$’s in Eq. (6.2) can be replaced approximately by $`b_k`$ except in the factors $`e^{i\zeta t}`$, $`b_k\zeta `$, and $`\zeta b_k`$. Therefore

$`\psi _k^{(k1)}(t)`$ $``$ $`e^{ib_k^2/2}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right](a_k){\displaystyle ^{b_k}}𝑑\zeta e^{i\zeta t}(b_k\zeta )^{1+ia_k^2},`$
$`\psi _k^{(k)}(t)`$ $``$ $`e^{ib_k^2/2}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right]a_k{\displaystyle _{b_k}}𝑑\zeta e^{i\zeta t}(\zeta b_k)^{1+ia_k^2},`$ (6.5)
or

$`\psi _k^{(k1)}(t)`$ $``$ $`e^{ib_k^2/2}e^{ib_kt}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right](a_k){\displaystyle _0^{\mathrm{}}}𝑑xe^{ixt}x^{1+ia_k^2},`$
$`\psi _k^{(k)}(t)`$ $``$ $`e^{ib_k^2/2}e^{ib_kt}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right]a_k{\displaystyle _0^{\mathrm{}}}𝑑xe^{ixt}x^{1+ia_k^2}.`$ (6.6)
This integral can be evaluated exactly in terms of the gamma function.

### 6.3 Results

Fig. 1 shows which ones of the various $`\psi _k^{(n)}(t)`$ have non-vanishing behaviors for $`t\mathrm{}`$ and $`t\mathrm{}`$.

These non-vanishing behaviors are:

For $`t`$ positive and large,

$`\psi _1^{(1)}(t)`$ $``$ $`\sqrt{2\pi }e^{i\pi /4}t^{i\alpha }e^{it^2/2},`$ (6.7)
$`\psi _k^{(k1)}(t)`$ $``$ $`e^{ib_k^2/2}e^{ib_kt}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right](a_k)e^{\pi a_k^2/2}\mathrm{\Gamma }(ia_k^2)t^{ia_k^2},`$ (6.8)
$`\psi _k^{(k)}(t)`$ $``$ $`e^{ib_k^2/2}e^{ib_kt}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right]a_ke^{\pi a_k^2/2}\mathrm{\Gamma }(ia_k^2)t^{ia_k^2};`$ (6.9)
while, for $`t`$ negative and large,

$`\psi _1^{(N)}(t)`$ $``$ $`\sqrt{2\pi }e^{i\pi /4}|t|^{i\alpha }e^{it^2/2},`$ (6.10)
$`\psi _k^{(k1)}(t)`$ $``$ $`e^{ib_k^2/2}e^{ib_kt}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right](a_k)e^{\pi a_k^2/2}\mathrm{\Gamma }(ia_k^2)|t|^{ia_k^2},`$ (6.11)
$`\psi _k^{(k)}(t)`$ $``$ $`e^{ib_k^2/2}e^{ib_kt}\left[{\displaystyle \underset{\begin{array}{c}j=2\\ jk\end{array}}{\overset{N}{}}}|b_jb_k|^{ia_j^2}\right]a_ke^{\pi a_k^2/2}\mathrm{\Gamma }(ia_k^2)|t|^{ia_k^2}.`$ (6.12)
All the other components approach zero as $`t\mathrm{}`$ and as $`t\mathrm{}`$. In the asymptotic formulas (6.7) and (6.10), $`\alpha `$ is the quantity

$$\alpha =\underset{j=2}{\overset{N}{}}a_j^2.$$
(6.13)
In the formulas (6.8), (6.9), (6.11) and (6.12),

$$2kN,$$
(6.14)
where $`N`$ as always is the number of neutrino species.

## 7 Discussion

When we started to investigate the MSW differential equations for three neutrino species in the case of the linear electron density, we were mostly interested in various possibilities of finding approximate solutions. Therefore, it was quite a surprise to us that these coupled differential equations can be solved exactly not only for three, but also for any number of neutrino species.

In the work of Wolfenstein, Mikheyev, Smirnov  and others  on the sun taking into account two species of neutrinos, it has been found that most of the effect takes place in a fairly narrow region around a particular value of the electron density. Because of this, it is quite accurate to use a linear approximation to the electron density.

For more than two species of neutrinos, it is no longer true in general that there is a narrow region for most of the activity. Nevertheless, there are a number of circumstances where this is true. However, the conditions for this to hold has not yet been studied systematically. This is one direction for future work.

Under the assumption of the electron density being a linear function of distance, the exact, general solution of the MSW differential equation is given by Eqs. (5.12) and (5.13). This solution is in the form of a number of single integrals. When the number of neutrino species is more than 2, these integrals cannot be evaluated in terms of known functions, and therefore their properties need to be investigated. A small step in this direction has been taken in Sec. 6, where the asymptotic behaviors of these integrals have been evaluated for large distances but with all the $`a`$’s and $`b`$’s held fixed. It is believed that, in so far as this case of linear electron density is applicable to the physically interesting case of solar neutrinos, the asymptotic evaluation of Sec. 6 is far from being sufficient. It is more likely that not only the distance, but also some of the parameters, the $`a`$’s and $`b`$’s, are large. This is a second direction for future work.

Acknowledgments

We are greatly indebted to Dr. Conrad Newton for collaboration at the early stage of this work. One of us (TTW) thanks the Theory Division at CERN for its kind hospitality.


text:
# Inclusive Heavy Hadron Decays and Light-Cone Dynamics

## 1 Introduction

The beauty hadron is the heaviest hadron containing a single heavy quark but lives for a long time, while the top quark is too heavy to build hadrons and has a very short lifetime. This confers a special role to beauty hadrons in heavy hadron physics: The longevity of the beauty hadron makes it phenomenologically interesting, and the heaviness of the beauty hadron ensures theoretical reliability.

Because of the long lifetime of the beauty hadron, many interesting phenomena have chance to be observed. We discuss here semileptonic or radiative inclusive decays of $`B`$ mesons. The discussion can be easily extended to inclusive decays of other heavy hadrons. Semileptonic or radiative inclusive decays of $`B`$ mesons are interesting because these decays are our main tool to determine the fundamental Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, to probe hadron structure and strong interactions, to test the standard model, and to search for new physics beyond the standard model.

The principal topic in the theoretical description of inclusive decays of $`B`$ mesons is the calculation of nonperturbative QCD contributions. There are essentially two components of nonperturbative QCD effects: one on dynamics — it is responsible for the confinement of quarks inside a hadron, the other on kinematics — it transforms parton kinematics into hadron kinematics. The past decade saw great progress in the QCD treatment of inclusive decays of heavy hadrons, including the construction of the heavy quark effective theory  and the developments of the heavy quark expansion approach  and the light-cone approach . In this talk I describe the light-cone approach and present a brief review of some of its applications. A comparison of the light-cone approach and the heavy quark expansion approach can be found in Refs. .

The light-cone approach uses the methods of light-cone expansion and heavy quark effective theory to address both dynamic and kinematic effects of nonperturbative QCD. This QCD approach has led to an improvement over the naive parton model  for inclusive $`B`$ decays, so that model-independent predictions and a control over theoretical uncertainties become possible.

## 2 Light-Cone Expansion

All nonperturbative QCD physics for inclusive semileptonic $`B`$ decays are incorporated in the hadronic tensor:

$$W_{\mu \nu }=\frac{1}{2\pi }d^4ye^{iqy}B\left|[j_\mu (y),j_\nu ^{}(0)]\right|B,$$
(1)
where $`q`$ stands for the momentum of the virtual $`W`$ boson. In general, the hadronic tensor can be expressed in terms of five scalar functions $`W_{15}(\nu ,q^2)`$ of the two independent Lorentz invariants, $`\nu =qP/M_B`$ and $`q^2`$ ($`M_B`$ and $`P`$ denote the mass and momentum of the $`B`$ meson, respectively),

$$W_{\mu \nu }=g_{\mu \nu }W_1+\frac{P_\mu P_\nu }{M_B^2}W_2i\epsilon _{\mu \nu \alpha \beta }\frac{P^\alpha q^\beta }{M_B^2}W_3+\frac{q_\mu q_\nu }{M_B^2}W_4+\frac{P_\mu q_\nu +q_\mu P_\nu }{M_B^2}W_5.$$
(2)
The theoretical task is to calculate the five structure functions.

Let’s look into the expression (1) for the hadronic tensor. The current commutator in Eq. (1) is causal, $`[j_\mu (y),j_\nu ^{}(0)]=0`$ for $`y^2<0`$. Moreover, $`e^{iqy}`$ in Eq. (1) oscillates rapidly, averaging out contributions to the integral except in the domain $`y^21/q^2`$. For inclusive semileptonic $`B`$ decays, the momentum transfer squared lies in the range $`M_{\mathrm{}}^2q^2(M_BM_{X_{\mathrm{min}}})^2`$, where $`M_{\mathrm{}}`$ is the charged lepton mass and $`M_{X_{\mathrm{min}}}`$ is the minimal hadronic invariant mass in the final state, which is the $`D`$ meson (pion) mass for the $`bc`$ ($`bu`$) transition. Because of the heaviness of the $`B`$ meson, the decays occur mostly at large momentum transfer. Therefore, the combined consideration of causality, less rapid oscillations and kinematics leads us to realize that the hadronic tensor is dominated by the space-time near the light cone $`y^20`$.

Another argument for the light-cone dominance is as follows. One can express the commutator of two currents in terms of the bilocal operator:

$$[j_\mu (y),j_\nu ^{}(0)]=2(S_{\mu \alpha \nu \beta }i\epsilon _{\mu \alpha \nu \beta })\left[^\alpha \mathrm{\Delta }_q(y)\right]\overline{b}(0)\gamma ^\beta 𝒫exp[ig_s_y^0𝑑z^\mu A_\mu (z)]b(y),$$
(3)
where $`S_{\mu \alpha \nu \beta }=g_{\mu \alpha }g_{\nu \beta }+g_{\mu \beta }g_{\nu \alpha }g_{\mu \nu }g_{\alpha \beta }`$ and $`\mathrm{\Delta }_q(y)`$ is the Pauli-Jordan function. $`𝒫`$ denotes path ordering. The coefficient function of the bilocal operator in Eq. (3) is finite except on the light cone $`y^2=0`$ where it is singular. Therefore, the light-cone singularity also implies that the dominant contribution to the hadronic tensor comes from the space-time separations in the neighborhood of the light cone.

The light-cone dominance justifies the light-cone expansion in matrix elements of increasing twist. The light-cone expansion provides a formal and powerful way of organizing the nonperturbative QCD effects and singling out the leading term. The leading nonperturbative QCD contribution to inclusive semileptonic decays of $`B`$ mesons resides in the distribution function

$$f(\xi )=\frac{1}{4\pi }\frac{d(yP)}{yP}e^{i\xi yP}B|\overline{b}(0)y/𝒫exp[ig_s_y^0𝑑z^\mu A_\mu (z)]b(y)|B|_{y^2=0}.$$
(4)
The five structure functions, a priori independent, are then related to a single distribution function in leading twist approximation :

$`W_1`$ $`=`$ $`2[f(\xi _+)+f(\xi _{})],`$ (5)
$`W_2`$ $`=`$ $`{\displaystyle \frac{8}{\xi _+\xi _{}}}[\xi _+f(\xi _+)\xi _{}f(\xi _{})],`$ (6)
$`W_3`$ $`=`$ $`{\displaystyle \frac{4}{\xi _+\xi _{}}}[f(\xi _+)f(\xi _{})],`$ (7)
$`W_4`$ $`=`$ $`0,`$ (8)
$`W_5`$ $`=`$ $`W_3,`$ (9)
where $`\xi _\pm =(\nu \pm \sqrt{\nu ^2q^2+m_q^2})/M_B`$ .

It has been shown  that inclusive radiative decays of $`B`$ mesons $`BX_s\gamma `$ are also dominated by light-cone dynamics. It turns out that the distribution function is universal: The same distribution function encodes the leading nonperturbative QCD contribution to inclusive radiative decays of $`B`$ mesons. The universality originates from the fact that the primary object of analysis in long-distance effects is the same bilocal operator matrix element dictated by light-cone dynamics. In leading twist approximation, the photon energy spectrum is then given by

$$\frac{d\mathrm{\Gamma }(BX_s\gamma )}{dE_\gamma }=\frac{G_F^2\alpha m_b^2}{2\pi ^4M_B}|V_{tb}V_{ts}^{}|^2|C_7^{(0)}(M_W)|^2E_\gamma ^3f\left(\frac{2E_\gamma }{M_B}\right).$$
(10)
## 3 Properties of the Distribution Function

The $`b`$-quark distribution function $`f(\xi )`$ for the $`B`$ meson is a key object. Like the well-known parton distribution functions for the nucleon in deep inelastic lepton-nucleon scattering, the knowledge of the $`b`$-quark distribution function would help us greatly in understanding the nature of confinement and the structure of the $`B`$ meson. I will survey our present knowledge of the $`b`$-quark distribution function and look forward to likely advances in the future.

The distribution function is gauge invariant. It obeys positivity and the support of it is $`0\xi 1`$. The distribution function is exactly normalized to unity

$$_0^1𝑑\xi f(\xi )=1.$$
(11)
This normalization does not get renormalized as a consequence of $`b`$-flavored quantum number conservation. The distribution function contains the free quark decay as a limiting case with $`f(\xi )=\delta (\xi m_b/M_B)`$. The distribution function $`f(\xi )`$ has a simple physical interpretation: It is the probability of finding a $`b`$-quark with momentum $`\xi P`$ inside the $`B`$ meson with momentum $`P`$.

In addition, the mean $`\mu `$ and the variance $`\sigma ^2`$ of the distribution function were deduced  using the techniques of the operator product expansion and the heavy quark effective theory (HQET):

$$\mu _0^1𝑑\xi \xi f(\xi )=\frac{m_b}{M_B}\left(1+\frac{5E_b}{3}\right),$$
(12)
$$\sigma ^2_0^1𝑑\xi (\xi \mu )^2f(\xi )=\left(\frac{m_b}{M_B}\right)^2\left[\frac{2K_b}{3}\left(\frac{5E_b}{3}\right)^2\right],$$
(13)
where $`E_b=K_b+G_b`$ and $`K_b`$ and $`G_b`$ are the dimensionless HQET parameters of order $`(\mathrm{\Lambda }_{\mathrm{QCD}}/m_b)^2`$, which are often referred to by the alternate names $`\lambda _1=2m_b^2K_b`$ and $`\lambda _2=2m_b^2G_b/3`$. The parameter $`\lambda _2`$ can be extracted from the $`B^{}B`$ mass splitting: $`\lambda _2=(M_B^{}^2M_B^2)/4=0.12`$ GeV2. The parameter $`\lambda _1`$ suffers from large uncertainty.

The mean value and variance of the distribution function characterize the location of the “center of mass” of the distribution function and the square of its width, respectively. They specify the primary shape of the distribution function. From Eqs. (12) and (13) we know that the distribution function is sharply peaked around $`\xi =\mu m_b/M_B`$ close to 1 and its width of order $`\mathrm{\Lambda }_{\mathrm{QCD}}/M_B`$ is narrow, suggesting that the distribution function is close to the delta function form in the free quark limit. Note that the $`b`$-quark distribution function has different forms for different beauty hadrons.

Nonperturbative QCD methods such as lattice gauge theory, light-cone field theory, and QCD sum rules could help determine further the form of the distribution function.

The $`b`$-quark distribution function can also be extracted directly from experiment. The $`BX_s\gamma `$ photon energy spectrum  and the $`BX_q\mathrm{}\nu `$ spectra $`d\mathrm{\Gamma }/d\xi _+`$  share a common feature — namely they offer the intrinsically most sensitive probe of long-distance strong interactions because these spectra correspond to a discrete line solely on kinematic grounds in the absence of gluon bremsstrahlung and long-distance strong interactions. Indeed, our calculation based on the light-cone expansion shows that they are explicitly proportional to the nonperturbative distribution function. Therefore, the shapes of these spectra directly reflect the inner long-distance dynamics of the reactions and measurements of these spectra are idealy suited for direct extraction of the distribution function from experiment.

The universality of the distribution function implies great predictive power: Once the distribution function is measured from one process, it can be used to make predictions in all other processes in a model-independent manner.

## 4 Applications

I have described the QCD approach to inclusive decays of heavy hadrons from light-cone expansion and heavy quark effective theory. Because of the large mass of the decaying hadron, light-cone dynamics plays a governing role in inclusive decays of heavy hadrons. The light-cone expansion allows a rigorous and systematic ordering of nonperturbative QCD effects, and the identification of the leading term — the $`b`$-quark distribution function. Another large scale set by the $`b`$-quark mass allows the construction of the heavy quark effective theory, providing another complementary framework for organizing and parametrizing nonperturbative QCD effects. The additional properties of the distribution function have been learned exploiting the heavy quark effective theory.

This approach has been applied to calculate decay rates and distributions in semileptonic or radiative inclusive decays of $`B`$ mesons. Model-independent predictions have been made from the known normalization of the distribution function or the cancellation of the distribution function in the ratio of the decay rates. Some calculations of decay rates and distributions involve the modelling of the distribution function, since the form of it has as yet not been completely determined. However, for sufficiently inclusive quantities such as the total semileptonic decay rates, the results are nearly model-independent, since they are essentially only sensitive to the mean value and variance of the distribution function, which are known from the heavy quark effective theory. Direct extraction of the distribution function from precision measurements would eliminate the model dependence.

A crucial observation is that both dynamic and kinematic effects of nonperturbative QCD must be taken into account. The latter results in the extension of phase space from the quark level determined by the $`b`$-quark mass to the hadron level determined by the $`B`$ meson mass, thereby increasing decay rates. The calculations in the light-cone approach are able to include both dynamic and kinematic effects of nonperturbative QCD, and have shown that the net effect of nonperturbative QCD enhances the total semileptonic rates for $`BX_c\mathrm{}\nu `$  and $`BX_u\mathrm{}\nu `$ , contrary to the results  obtained in the heavy quark expansion approach. The heavy quark expansion approach  assumes quark-hadron duality, and cannot account for the rate due to the extension of phase space from the quark level to the hadron level. The light-cone approach provides a way around problems associated with the assumption of quark-hadron duality. The rate calculations in the light-cone approach have quantitatively shown  the importance of the inclusion of kinematic effects of nonperturbative QCD.

The interplay between nonperturbative and perturbative QCD effects has been accounted for in the light-cone approach, since confinement implies that free quarks are not asymptotic states of the theory and the separation of perturbative and nonperturbative effects cannot be done in a clear-cut way.

The following are some of the results.

(1) A new method for precise and model-independent determination of $`|V_{ub}|`$ has been proposed. For inclusive charmless semileptonic decays of $`B`$ mesons $`BX_u\mathrm{}\nu `$, the light-cone expansion and $`b`$-flavored quantum number conservation lead to the sum rule

$$S_0^1𝑑\xi _u\frac{1}{\xi _u^5}\frac{d\mathrm{\Gamma }}{d\xi _u}(BX_u\mathrm{}\nu )=|V_{ub}|^2\frac{G_F^2M_B^5}{192\pi ^3},$$
(14)
with the kinematic variable $`\xi _u=(q^0+|\stackrel{}{q}|)/M_B`$ in the $`B`$-meson rest frame. This sum rule receives no perturbative QCD correction, and avoids the dominant hadronic uncertainty. The sum rule (14) thus establishes a clean relationship between $`|V_{ub}|`$ and the observable $`S`$, allowing precise determination of $`|V_{ub}|`$. This determination of $`|V_{ub}|`$ is also model-independent in the sense that the sum rule is independent of phenomenological models. Moreover, this method is not only exceptionally clean theoretically, but also very efficient experimentally in background suppression. Applying the kinematic cut $`q^0>M_BM_D`$ or $`M_X<M_D`$, which leads to $`\xi _u>1M_D/M_B`$, one can discriminate between $`bu`$ signal and $`bc`$ background. The $`\xi _u`$ spectrum is unique. About $`80\%`$ of the spectrum satisfy $`\xi _u>1M_D/M_B`$.

(2) $`|V_{ub}|`$ has been determined  from the measured inclusive charmless semileptonic branching ratio of beauty hadrons. The determination of $`|V_{ub}|`$ from the measurement of the inclusive charmless semileptonic branching ratio in conjunction with the calculation of the total charmless semileptonic rate has larger theoretical uncertainty than the proposed determination of $`|V_{ub}|`$ from the measurement of the observable $`S`$ in conjunction with the sum rule (14), since the calculations of perturbative and nonperturbative QCD corrections to the total charmless semileptonic rate suffer from much larger theoretical uncertainty  than the sum rule (14), not to mention the fundamental uncertainty due to the assumption of quark-hadron duality if the total semileptonic rate is calculated in the heavy quark expansion approach.

(3) The $`BX_u\mathrm{}\nu `$ ($`\mathrm{}=e`$ or $`\mu `$) charged lepton energy spectrum , hadronic invariant mass spectrum , and lepton pair spectrum  have been analyzed in detail, resulting in the state-of-the-art descriptions of these spectra. It has been shown that the interplay between perturbative and nonperturbative QCD eliminates the singularities of the parton-level perturbative spectra.

(4) $`|V_{cb}|`$ has been determined  from the measured inclusive semileptonic branching ratio of beauty hadrons in conjunction with the calculation of the total semileptonic rate. The charged lepton energy spectrum in inclusive semileptonic decays of $`B`$ mesons $`BX_c\mathrm{}\nu `$ has been analyzed  in detail, resulting in the state-of-the-art description of the spectrum. The calculated spectrum was found  to be in agreement with the experimental data, providing an experimental test of the validity of the light-cone approach.

(5) The total semileptonic decay rate of the $`\mathrm{\Lambda }_b`$ baryon has been calculated. The calculated semileptonic branching ratio for $`\mathrm{\Lambda }_b`$ is consistent with the measurements.

(6) Nonperturbative QCD effects on the $`BX_s\gamma `$ photon energy spectrum have been calculated. The theoretically clean methods for the determinations of $`|V_{ts}|`$, $`|V_{ts}/V_{ub}|`$ and $`|V_{ts}/V_{cb}|`$ have been suggested.

## Acknowledgments

I am indebted to Emmanuel Paschos for collaboration. I would like to thank Stanley Brodsky, Chueng-Ryong Ji and Bruce McKellar for useful discussions. This work was supported by the Australian Research Council.