text:
# Double parton scatterings in b-quark pairs production at the LHC

## I Introduction

One of the main topics at the LHC is the production of $`b`$-quarks, both to search for CP violation, looking at $`b`$ decays, and to test QCD, by studying the production mechanism. Bottom quarks are also a large source of background to several processes of interest, as in various promising channels for Higgs detection. The production mechanism of heavy quarks in hadronic collisions pose, on the other hand, non trivial problems already at smaller energies, where also the simplest observable quantity, the integrated inclusive cross section, is not reproduced trivially.

The inclusive cross section of $`b`$-quarks production has been evaluated in pQCD at the next to leading order in $`\alpha _S`$. Unfortunately comparisons with the recent experimental data of the D0 Collaboration at TEVATRON have shown that the NLO pQCD calculations underestimate the cross section by a factor $`2,3`$, showing that NNLO corrections, whose explicit evaluation is still an open question, give a large contribution to the cross section.

A complementary approach to heavy quarks production, which keeps explicitly into account that transverse momenta and virtualities of the interacting partons become increasingly important in the kinematical regime of $`sm_b^2\widehat{s}\mathrm{\Lambda }^2`$, and includes terms at every order in $`\alpha _S`$ in the calculation of the cross section, is the $`k_t`$-factorization, where the interaction is factorized into un-integrated structure functions and off shell matrix elements . From the phenomenological point of view $`k_t`$-factorization is not inconsistent with HERA and TEVATRON data, allowing one to reproduce both the value of the integrated inclusive cross section and various differential distributions, including the correlation in the azimuthal angle between the produced $`b`$ quarks, where different approaches are less successfully compared with experiment (see  and references therein).

Interestingly, although the value of the integrated inclusive cross section cannot be obtained trivially, one may find several cases where the overall effect of higher order corrections amounts to a simple rescaling of the lowest order parton model result. There are in fact several distributions, derived either using the $`k_t`$-factorization approach or by working out the cross section at the NLO pQCD, which are rather similar (apart form the normalization factor) to those obtained with a simplest lowest order calculation. Hence, in a few cases, the whole effect of higher order corrections is (approximately) reduced to a single numerical value, the so-called $`K`$ factor:

$$K=\frac{\sigma (b\overline{b})}{\sigma _{LO}(b\overline{b})}.$$
(1)
where $`\sigma (b\overline{b})`$ is the inclusive cross section for $`b\overline{b}`$ production and $`\sigma _{LO}(b\overline{b})`$ the result of the lowest order calculation in pQCD.

When looking at extrapolations of the cross section at high energies, one finds that the result is affected by several uncertainties, as the knowledge of the parton structure functions at very small $`x`$ and the values of the heavy quark mass and of the running coupling constant. Although the expected inclusive cross section of $`b`$ production is hence still pretty uncertain at LHC energy, all estimates point in the direction of rather large values, as a consequence of the high parton luminosity. The fairly large flux of partons make it also plausible to expect a sizable rate of events, where two or more $`b\overline{b}`$ pairs are produced contemporarily by different partonic collisions in a given $`pp`$ interaction. Although at present stage all quantitative predictions for this much more structured interaction process are unavoidably pretty uncertain, the large cross sections foreseen at the LHC is, in our opinion, a strong motivation to make an attempt of giving a few quantitative indications on the production rate of multiple $`b\overline{b}`$ pairs trough multiparton interactions at the LHC, comparing with the rates to be expected by the more conventional single parton scattering mechanism.

Since the details of the elementary production of heavy quarks are still a matter of debate, we limit our considerations, for the production of multiple $`b\overline{b}`$ pairs, to the simplest cases, where the whole effect of higher order corrections is taken into account by the overall normalization factor. Given the lack of information on higher order corrections in the $`24`$ processes, we make moreover the assumption that the $`K`$ factors of the $`ggb\overline{b}b\overline{b}`$ and of the $`ggb\overline{b}`$ processes are equal. Hence we work out the $`ggb\overline{b}`$ process in the $`k_t`$-factorization approach, fixing the input parameters by comparing with the TEVATRON data, and extrapolate the cross section at LHC energies, identifying a few distributions where the effect of higher order corrections reduceds to a simple rescaling of the lowest order result. The value of the $`K`$-factor derived in this way is then used to renormalize the double $`(ggb\overline{b})^2`$ and the single $`ggb\overline{b}b\overline{b}`$ parton scattering cross sections, which we evaluate by working out all Feynman diagrams at order $`\alpha _S^4`$.

## II $`b\overline{b}`$ cross section at TEVATRON and LHC and $`K`$-factor

In the $`k_t`$-factorization approach the $`b\overline{b}`$ production cross section is expressed as

$$\sigma (ppb\overline{b})=\frac{d^2q_{t1}}{\pi }\frac{d^2q_{t2}}{\pi }𝑑x_1𝑑x_2f(x_1,q_{t1},\mu )f(x_2,q_{t2},\mu )\widehat{\sigma }(x_1,q_{t1};x_2,q_{t2};\mu )$$
(2)
where $`f(x,q_t,\mu )`$ is the unintegrated structure function, representing the probability to find a parton with momentum fraction $`x`$, transverse momentum $`q_t`$ at the factorization scale $`\mu `$, while $`\widehat{\sigma }`$ is the off-shell partonic cross section of the process $`g^{}g^{}Q\overline{Q}`$.

To work out the inclusive cross section we use two different prescriptions for constructing the $`k_t`$-distributions from the usual integrated parton densities. The first prescription is based on the conventional DGLAP evolutions equations, with virtual corrections re-summed in the survival probability factor $`T_a(k_t^2,\mu ^2)`$ . Hence the un-integrated structure function for the parton $`a`$ reads

$$f_a(x,k_t^2,\mu ^2)=T_a(k_t^2,\mu ^2)\left[\frac{\alpha _S(k_t^2)}{2\pi }_x^{1\delta }P_{aa^{}}(z)a^{}(\frac{x}{z},k_t^2)𝑑z\right]$$
(3)
where $`P_{aa^{}}(z)`$ is the splitting function, $`a^{}(x,k_t^2)`$ the integrated structure function and $`\delta `$ a cutoff parameter introduced to give sense to the integral. Although not written explicitly also $`T_a(k_t^2,\mu ^2)`$ depends on $`\delta `$, in such a way the $`f_a(x,k_t^2,\mu ^2)`$ is a smooth function when $`\delta `$ becomes small.

As for the second prescription we follow ref., where the un-integrated structure functions are obtained from the leading order BFKL equation and are expressed as the convolution of the usual collinear gluon densities $`G(x,\mu ^2)`$ with the universal function $`𝒢(x,k_t^2,\mu ^2)`$

$$(x,k_t^2,\mu ^2)=_x^1𝑑\xi 𝒢(\xi ,k_t^2,\mu ^2)G(\frac{\xi }{x},\mu ^2)$$
(4)
The weight factors $`𝒢(\xi ,k_t^2,\mu ^2)`$ have a known analytic expression, in double-logarithmic approximation, in terms of Bessel functions and depend on the quantity $`\overline{\alpha }_s=3\alpha _s/\pi `$, which in the BFKL formalism is a fixed parameter, related to the pomeron intercept $`\alpha (0)=1+\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ in leading log approximation is $`\mathrm{\Delta }=4\overline{\alpha }_s\mathrm{log}2`$. Following we take the value $`\mathrm{\Delta }=0.35`$.

To generate the un-integrated structure functions we have used the parton distributions set GRV94 with factorization scale $`\mu _F^2=\widehat{s}`$. Hence in evaluating the cross sections, in the $`k_t`$-factorization approach, we have set the renormalization scale equal to the gluon virtuality. To obtain the cross section at the lowest order in pQCD we have used the MRS99 parton distributions , with factorization and renormalization scale equal to the transverse mass of the $`b`$-quark. Comparing the total cross section values in the two approaches, we have obtained for the $`K`$-factor the value $`K5.5`$.

In Fig.1 we plot the integrated cross section of $`b\overline{b}`$ production at TEVATRON($`\sqrt{s}=1.8`$TeV) and at LHC($`\sqrt{s}=14`$TeV), as a function of the minimum value of the transverse momentum $`p_t^{min}`$ of the $`b`$-quark. The dotted curves represent the cross section derived using the unintegrated gluon structure function, according with the BFKL prescription of Eq(4), whereas the dashed lines have been obtained by using the prescription in Eq(3). The continuous lines represent the result of the lowest order calculation multiplied by the $`K`$ factor. At TEVATRON energy the $`b`$-quark distributions are within the rapidity interval $`|y|<1`$ and are compared with the D0 experimental data . The same distributions, extrapolated at LHC energy, are then plotted as a function of $`p_t^{min}`$ within the pseudorapidity interval $`|\eta |<0.9`$, corresponding to the acceptance of the ALICE detector.

In Fig.2 we show the rapidity$`(y)`$ and pseudorapidity$`(\eta )`$ distributions normalized to one and within $`|\eta |<0.9`$. Here the continuous histograms refer to the result of the lowest order calculation, rescaled by the $`K`$ factor, whereas the dashed histograms represent the distributions evaluated with the $`k_t`$-factorization approach. As one may see also in this case the whole effect of higher orders reduces to a simple rescaling.

## III $`b\overline{b}b\overline{b}`$ cross section

The leading order QCD process to produce two pairs of heavy quarks is given by the single parton scattering term at the fourth order in the coupling constant. A competing mechanism at the LHC energy is the double parton scattering. We compare the two mechanisms in proton-proton collisions in the kinematical range of the ALICE and of the LHCb detectors, namely at center-of-mass energies of $`5.5`$ and $`14`$ TeV, within the pseudorapidity regions $`|\eta |<0.9`$ and $`1.8<\eta <4.9`$, down to very low transverse momenta.

The single scattering pQCD sub-processes at the lowest order in $`\alpha _s`$, in $`ppb\overline{b}b\overline{b}`$, are the quarks initiated process, $`q\overline{q}b\overline{b}b\overline{b}`$, whose amplitude is given by the sum of $`14`$ Feynman diagrams for each flavor in the initial state, and gluon fusion, $`ggb\overline{b}b\overline{b}`$, represented by $`76`$ diagrams altogether, the latter amplitude giving the dominant contribution to the cross section at small $`x`$. To evaluate the cross section we have generated the matrix elements of the partonic amplitudes with MadGraph  and HELAS  and we have used the MRS99 parton distributions , with the factorization scale equal to the renormalization scale $`\mu _F=\mu _R`$, which we have kept fixed at the value of the transverse mass of the produced $`b`$ quark. For the mass of the bottom quark we have used the value $`m_b=4.6`$ GeV. The multi-dimensional integrations have been performed by VEGAS  and the resulting cross section has been finally multiplied by the $`K`$ factor obtained as described in the previous section.

The evaluation of the double parton scattering contribution to the cross section is considerably more uncertain because of the unknown non perturbative input to the process, given by the two-body parton distribution functions $`\mathrm{\Gamma }(x_1,x_2,\beta )`$, where $`x_{1,2}`$ are the fractional momenta of the two partons belonging to the same hadron and $`\beta `$ their distance in transverse space. Although not explicitly written, the distributions depend also on the scale factors characterizing each elementary interaction and on the different kinds of partons involved. Given the large parton population at low $`x`$, to proceed further we make the usual simplifying assumption of neglecting correlations in fractional momenta and we factorize the two-body parton distribution as $`\mathrm{\Gamma }(x_1,x_2,\beta )=G(x_1)G(x_2)F(\beta )`$, where $`G(x)`$ are the usual one-body parton distributions and $`F(\beta )`$ is a function normalized to 1 and representing the parton pair density in transverse space. With these assumptions the cross section acquires the simplified form

$$\sigma _D(b\overline{b};b\overline{b})=\frac{1}{2}\underset{ij}{}\mathrm{\Theta }^{ij}\sigma _i(b\overline{b})\sigma _j(b\overline{b})$$
(5)
where the indices $`i,j`$ label the different cases where each $`b\overline{b}`$ pair is originated either by a $`q\overline{q}`$ annihilation, discriminating the cases of sea and valence, or by two gluons and $`\sigma _i(b\overline{b})`$ represents the inclusive cross sections for $`b\overline{b}`$ production in a hadronic collision, with the index $`i`$ labelling a definite parton process. The weight factors $`\mathrm{\Theta }^{ij}`$ have dimension an inverse cross section and result from integrating the product of the two-body parton distributions in transverse space, while the factor $`1/2`$ is a consequence of the symmetry of the expression for exchanging $`i`$ and $`j`$. The dependence of $`\mathrm{\Theta }^{ij}`$ on the indices $`i,j`$ accounts for the possibility, for different pairs of partons in the hadron structure, to be characterized by different values of their relative average transverse distance  . Notice that by measuring the double parton collisions one has access to a new information on the hadron structure, summarized in these weight factors, which cannot be obtained in hard processes with a single parton interaction only.

The experimental information on the double parton scatterings is due to the four-jet production measurement in $`pp`$ collisions at $`\sqrt{s}=63`$ GeV, performed by the AFS Collaboration, and to the study of final states with three minijets and one photon in $`p\overline{p}`$ collisions at $`\sqrt{s}=1800`$ GeV, due to CDF . In both cases the cross section was expressed as

$$\sigma _D=\frac{m}{2}\frac{\sigma _S(A)\sigma _S(B)}{\sigma _{eff}}$$
(6)
where $`m=1`$ if the two parton processes $`A`$ and $`B`$ are identical, while $`m=2`$ if they are different and $`\sigma _S`$ is the single scattering inclusive cross section. The overall output of the experiment hence reduces to the value of a single parameter, the scale factor $`\sigma _{eff}`$, whose value is $`\sigma _{eff}=5`$ mb for AFS, while $`\sigma _{eff}=14.5`$ mb for CDF. The two experimental results are not inconsistent, given the different content of partons in the two cases, mainly valence quarks in the former case and mostly gluons and sea quarks in the latter; the experimental indication hence pointing in the direction of a sizable dependence of the factors $`\mathrm{\Theta }^{ij}`$ in (5) on the different elementary processes. Interestingly, the measurement of $`\sigma _{eff}`$ for different final states, as a function of the c.m. energy and cuts applied, allows one obtaining the values of the scale factors $`\mathrm{\Theta }^{ij}`$, letting in this way access to the three dimensional structure of the proton.

The dominant contribution to $`b\overline{b}b\overline{b}`$ production at the LHC is gluon fusion, so, for the present purposes, the sum in (5) may be approximated well by a single term, where the scale factor could not be strongly different with respect to the case of the CDF experiment. Hence, to evaluate the double scattering cross section, we have used the simplest expression

$$\sigma _D(b\overline{b}b\overline{b})=\frac{\sigma (b\overline{b})^2}{2\sigma _{eff}}.$$
(7)
where for $`\sigma _{eff}`$ we have taken the value reported by CDF. Notice that since $`\sigma _D`$ is proportional to $`\sigma _S^2`$, the effect of higher order corrections is enhanced on $`\sigma _D`$:

$`\sigma _S=K\sigma _S^{LO}`$ (8)
$`\sigma _D=K^2\sigma _D^{LO}`$ (9)
where $`\sigma _{S,D}^{LO}`$ refers to the lowest order expressions of the cross section.

## IV  Results

In Fig.3 we plot the expected rise of the total $`b\overline{b}b\overline{b}`$ production cross section as a function of the c.m. energy. The continuous curves refer to the double parton scattering contribution, while the dotted curves to single scattering. In each case the lower curve refers to the value $`K=2.5`$, while the higher curve to $`K=5.5`$, which are the typical estimate of the NLO-QCD and the result of our calculation within the $`k_t`$-factorization approach. Notice that at the LHC the double parton scattering gives a contribution to the integrated cross section about ten times larger than the single scattering.

To see how the cross section depends on the transverse momenta we have plotted in Fig.4 the two contributions to the integrated cross section, as a function of $`p_t^{min}`$, the minimum value of the transverse momenta of the $`b`$ quarks (which we require to be all inside the pseudorapidity interval $`|\eta |<0.9`$), for the center-of-mass energy values of 14 and 5.5 TeV. The continuous histograms refer to the double parton scattering contribution, while the dotted histograms to single scattering. The double parton cross section decreases faster with $`p_t^{min}`$ than the single parton cross section, the two contributions being of the same order at $`p_t^{min}=810`$ GeV. Pseudorapidity, and rapidity distributions at 14 and 5.5 TeV are plotted in Fig.5 ( always requiring for the two $`b`$-quarks $`|\eta |<0.9`$), where continuous and dashed histograms have the same meaning as in the previous cases.

To see how the results depend on rapidity, we have plotted in Fig6 the same rapidity ($`y`$) and pseudorapidity ($`\eta `$) distributions at $`\sqrt{s}=14`$ TeV, requiring both $`b`$-quarks to be in the pseudorapidity interval $`1.8<\eta <4.9`$, which corresponds to the acceptance of the LHCb experiment. In the same figure we also compare the two contributions of single and double parton scattering, integrated within the rapidity acceptance of the LHCb, as a function of $`p_t^{min}`$.

The overall indication which one obtains from the present study is that double parton scatterings dominate the $`b\overline{b}b\overline{b}`$ integrated cross section by a large factor, both in the central rapidity region and at the larger rapidity values of the LHCb experiment. In both cases the contribution of the single parton scattering term becomes important only after applying cuts to the transverse momenta of the order of 8-10 GeV. The large values expected for the $`b\overline{b}b\overline{b}`$ cross section, which at $`14`$ TeV are of the order of one $`\mu `$b, inside the ALICE and the LHCb detectors, and the localization of the double scattering contribution at relatively low $`p_t`$ values, suggest that heavy quark pairs production at the LHC might represent an efficient tool for studying the gluon initiated double parton scattering process.

Acknowledgment

This work was partially supported by the Italian Ministry of University and of Scientific and Technological Researches (MIUR) by the Grant COFIN2001.


format: text

text:
# Reaction Pathways Based on the Gradient of the Mean First-Passage Time

## Introduction.

In many physical, chemical, or biological reactions, initial (reactant) and final (product) states are known but reaction pathways are not. Examples range in complexity from single-particle Brownian motion to conformational changes of proteins, such as protein folding Pande et al. (1998); Eastman et al. (2001). Finding reaction pathways is one of the most fundamental challenges in chemistry and molecular biologyElber (1996); Straub (2001). It is important for understanding reaction mechanisms and for the calculation of reaction rates Hänggi et al. (1990).

In most cases of interest, reactions take place at finite temperature and therefore are stochastic: every reaction event follows a different path and takes a different amount of time. Among all possible paths from reactant to product, one seeks the path that best characterizes the reaction not (a). The selection criteria are twofold. First, one seeks minimal paths, free of unnecessary fluctuations. Second, and more important, one seeks paths that are representative of all reaction events so that individual reaction events can be considered to be following noisy paths around them. It is, however, challenging to formulate these criteria rigorously.

The most widely used formulation of reaction path is probably the steepest-descent path, which is constructed by identifying saddle points of the potential energy surface and then following the steepest descent from the saddle points such that energy barriers along the path are minimized Elber (1996). The steepest-descent path, however, does not involve temperature. Since reactions are driven by thermal fluctuations, reaction paths should depend on temperature. For example, if there is a direct path with high energy barriers and a roundabout path with low energy barriers, at a temperature higher than the barriers reaction will occur most likely along the direct path rather than the roundabout path while the steepest-descent path will be the roundabout path regardless of temperature.

Understanding this drawback of the steepest-descent path has led to alternative formulations of reaction path. One approach is to select the path of maximum flux Berkowitz et al. (1983); Huo and Straub (1997). Another formulation focuses on most probable paths Pratt (1986); Elber and Shalloway (2000). In this approach, an ensemble of reaction events of a fixed time interval is considered. A probability is then assigned to each event, and the path followed by the most probable event is taken as the reaction path not (b).

The above methods succeeded to a certain extent in elucidating reaction mechanisms, but they do not fully satisfy the criterion of representativity. The method of most probable path comes closest to satisfying the criterion, but it is not clear how to choose the time interval beforehand and whether an ensemble of reaction events of a single time interval suffices to represent the reaction.

In this paper we present a new formulation of reaction path. While previous approaches attempted to quantify paths, we use the concept of reaction coordinate which quantifies states. Reaction coordinate is a function that describes where in the progress of reaction a state is located. The most natural measure of the progress of reaction is provided by the mean first-passage time (MFPT), namely the average amount of time that the system starting from the state takes to reach the given product. The MFPT depends on the energy landscape, the temperature, as well as the boundary conditions, and most important, it is an average over all reaction events. Surprisingly, to our knowledge, the MFPT has never been used as a reaction coordinate before. Once the MFPT is determined, one can choose as reaction paths the paths along which the MFPT decreases most rapidly, which complies with the criterion of minimal path.

## Theory.

A practical scheme of determining reaction coordinates and reaction paths based on MFPTs can indeed be stated, as described below and illustrated in Fig. 1.

Consider a system undergoing a stochastic reaction and let $`𝒮`$ be the set of all states that the system can access. If $`𝒮`$ is continuous the reaction can be described by a Fokker-Planck-type equation, and if $`𝒮`$ is discrete it can be described by a master equation Gardiner (1985). The reactant and the product are specified by disjoint subsets, $``$ and $`𝒫`$ respectively, of $`𝒮`$. The reaction can be considered a first-passage process Redner (2001) because a reaction event ends as soon as the system reaches any state in $`𝒫`$. The MFPT $`\tau (𝐱)`$ from state $`𝐱`$ to the product $`𝒫`$ is then calculated for all states $`𝐱`$ in $`𝒮`$ (this involves solving an inhomogeneous partial differential equation when $`𝒮`$ is continuous and inverting a transition matrix when discrete MFP ; Bauer et al. (1988), as demonstrated later) and is used as a reaction coordinate. The location of the reactant set $``$ does not affect the calculation of the reaction coordinate $`\tau (𝐱)`$; it is involved only in determining reaction paths.

The scheme of constructing reaction paths from MFPTs depends on the character of the set $`𝒮`$. When $`𝒮`$ is continuous and described by a Cartesian coordinate system, reaction paths are constructed following the direction of $`\tau `$, along which $`\tau `$ decreases most rapidly. Thus, a reaction path $`𝐱(l)`$, parameterized through an arc length $`l`$, satisfies

$$\frac{dx_i}{dl}=v_i\left(_jv_jv_j\right)^{1/2},v_i=\frac{\tau }{x_i}.$$
(1)
Often, reactions are better described with non-Cartesian coordinates such as angles. In such cases reaction paths can be determined via a transformation to a Cartesian coordinate system, and the resulting equation of reaction path is

$$\frac{dx_i}{dl}=v_i\left(_jv_jv_j\right)^{1/2},v_i=_jg_{ij}^1\frac{\tau }{x_j}.$$
(2)
Here $`g_{ij}^1`$ is the inverse matrix of the metric tensor $`g_{ij}`$ and $`l`$ is now the arc length with respect to the non-Cartesian coordinate system ($`dl^2=_idx_idx_i`$).

When $`𝒮`$ is discrete and the reaction is governed by a master equation with transition rates $`k_{yx}`$ (from state $`x`$ to state $`y`$), the MFPT $`\tau _x`$ from state $`x`$ to the product $`𝒫`$ is again employed as a reaction coordinate. But in order to determine reaction paths an analogue of metric is required, as reaction paths for continuous systems are determined via gradients which involve metric. The most obvious choice for an analogue of metric is the transition rates $`k_{yx}`$ themselves, and we suggest the scheme that a reaction path going through state $`x`$ chooses the next state $`y`$ such that $`k_{yx}(\tau _x\tau _y)`$ is maximized. For the transition step from $`x`$ to $`y`$, the transition time $`1/k_{yx}`$ may be interpreted as a cost, and the decrease $`\tau _x\tau _y`$ in the MFPT as a gain. The scheme then amounts to maximizing the ratio between these two times, namely the gain-cost ratio, for each step.

According to the above scheme, a reaction path is constructed starting from each state in the reactant set $``$. In general, multiple reaction paths are obtained unless the reactant set is narrowed down to a single state. Some reaction paths may overlap if they go through a common state.

## Examples.

To demonstrate the above scheme and its outcomes, we consider first a single-particle Brownian motion on a two-dimensional potential surface described with a Cartesian coordinate system ($`x_1,x_2`$). We take the three-hole potential

$`U(x_1,x_2)`$ $`=`$ $`3e^{x_1^2(x_25/3)^2}+3e^{x_1^2(x_21/3)^2}`$ (3)
$`5e^{(x_11)^2x_2^2}5e^{(x_1+1)^2x_2^2}`$
which was also studied by others regarding the temperature dependence of reaction paths Huo and Straub (1997); Elber and Shalloway (2000).

The potential features two deep holes and one shallow hole (Fig. 2). The two deep holes are considered as the reactant and the product. Roughly two possible pathways can be seen; the upper path is longer than the lower one but has lower energy barriers. It is therefore expected that the upper path will be taken at low temperature and that the lower path will be taken at high temperature.

The Brownian motion can be described in terms of the probability distribution $`p(𝐱,t)`$ and the probability current $`𝐉(𝐱,t)`$. In the strong friction regime, they satisfy the Smoluchowski equation Gardiner (1985)

$`{\displaystyle \frac{}{t}}p(𝐱,t)=𝐉(𝐱,t)`$
$`𝐉(𝐱,t)=(\beta \gamma )^1e^{\beta U(𝐱)}[e^{\beta U(𝐱)}p(𝐱,t)],`$ (4)
where $`\gamma `$ is the friction coefficient. The MFPT $`\tau (𝐱)`$ is then obtained by solving the inhomogeneous partial differential equation

$$[e^{\beta U(𝐱)}\tau (𝐱)]=\beta \gamma e^{\beta U(𝐱)}$$
(5)
with appropriate boundary conditions MFP . We take the region $`(4x_14,3x_24)`$ as the whole set $`𝒮`$, and assume that its boundary is reflecting, namely the probability current $`𝐉`$ is tangential to the boundary. For the reactant $``$ we take the single point $`(1,0)`$, and for the product $`𝒫`$ we take the circular region of radius 0.1 centered at $`(1,0)`$. Because a reaction event ends when the particle reaches the product region $`𝒫`$, the boundary of $`𝒫`$ is absorbing, with the probability distribution $`p`$ vanishing at the boundary. Boundary conditions for $`p`$ and $`𝐉`$ lead to corresponding boundary conditions for $`\tau `$: $`\tau `$ vanishes at absorbing boundaries and $`\tau `$ is tangential to reflecting boundaries MFP .

The solutions of Eq. (5), numerically obtained with MATLAB Mat , for two different temperatures are shown in Fig. 3.

The differences between the two temperatures are dramatic. At the high temperature ($`\beta =1`$) the arrows denoting the directions of $`\tau `$ flow more or less directly towards the product. At the low temperature ($`\beta =8`$), on the other hand, the flow is significantly distorted so that energy barriers are avoided, with a singular point produced around $`(1.5,0.5)`$. Also, the reaction coordinate $`\tau `$ drops rapidly when barriers are crossed, as indicated by the contours of $`\tau `$ packed around saddle points.

Fig. 4 shows the reaction paths found at various temperatures.

As was expected, lower paths are taken at high temperature and upper paths are taken at low temperature. At intermediate temperature, such as $`\beta =4`$, reaction paths lie in between, indicating that the upper and lower paths are equally favorable.

In order to demonstrate the scheme for reaction paths of discrete systems, we consider light-harvesting complexes Blankenship (2002); LH ; Damjanović et al. (2002), choosing photosystem I of cyanobacterium Synechococcus elongatus as an example. Photosystem I is a protein-pigment complex, embedded in the cell membrane, that contains 96 chlorophylls. The aggregate of these chlorophylls is responsible for the first step in photosynthesis, namely absorption of light and subsequent migration of the resulting electronic excitation towards the special pair of chlorophylls, called P700 and located at the geometrical center of the aggregate, where the next step in photosynthesis, the charge separation, occurs. This excitation migration can be considered a reaction. Assuming that no more than a single chlorophyll is simultaneously excited in the system, the states are specified by the excited chlorophylls. The reactant is specified by the chlorophyll that initially absorbed a photon; the product is specified the P700 pair of chlorophylls. Reaction paths of this system provide representative and minimal pathways of the excitation migration from the initially excited chlorophyll to the P700 pair.

Since we are interested in first passage of excitation to P700, it is convenient to consider a subsystem of 94 chlorophylls, excluding the P700 pair of chlorophylls. The migration of excitation in this subsystem can be described by the master equation

$$\frac{d}{dt}p_x(t)=_yK_{xy}p_y(t),$$
(6)
where $`p_x(t)`$ is the probability that the excitation resides at chlorophyll $`x`$ at time $`t`$ and $`K_{xy}`$ is a $`94\times 94`$ transition matrix Ritz et al. (2001); Park et al. . We build the transition matrix $`K_{xy}`$ by using the inter-chlorophyll excitation transfer rates that were calculated in Ref. Şener et al. (2002) based on the theory developed in Damjanović et al. (1999) and the recently obtained structure of photosystem I Jordan et al. (2001). The MFPT $`\tau _x`$ from chlorophyll $`x`$ to the P700 pair is then Park et al.

$$\tau _x=_y\varphi _yK_{yx}^1/\varphi _x,\varphi _x=_y\xi _yK_{yx}^1.$$
(7)
Here $`K_{xy}^1`$ is the inverse matrix of $`K_{xy}`$ and $`\xi _x`$ is the excitation transfer rate from chlorophyll $`x`$ to P700. The reaction paths constructed from the MFPT are shown in Fig. 5.

A detailed discussion of this system is reported in Park et al. .

## Conclusion.

We have presented a new method to construct reaction paths based on the MFPT gradient that incorporates all reaction events, and illustrated how it captures important aspects of reactions, most notably temperature effects. Unlike previous attempts, the present method provides paths starting from all states of the system. We believe that the MFPT is the most natural choice for reaction coordinate, and expect that our approach will be used in many future studies of reactions. The generality of MFPT permits applications of the present method to a wide range of phenomena. The present method is particularly suitable for large reaction networks that are completely characterized by the method through a pathway graph that connect each state to the product.

###### Acknowledgements.

We thank Paul Grayson and Melih K. Şener for useful discussions. This work has been supported by National Institutes of Health grant PHS 5 P41 RR05969.


text:
# 1 Diagram leading to 𝜏→𝜇𝜂. The small circle is the lepton-flavor violating vertex of Babu and Kolda.

## Abstract

The existence of large $`\nu _\mu \nu _\tau `$ mixing suggests the likelihood of large smuon-stau mixing in supersymmetric models, leading to $`\mu `$ and $`\tau `$ number violation. In addition to interesting signatures in slepton and neutralino production and decay, this will lead to rare $`\tau `$ decays, such as $`\tau \mu \gamma `$. Recently, it has been pointed out that the $`\tau 3\mu `$ branching ratio could be substantial in the large $`\mathrm{tan}\beta `$ region of parameter space, due to an induced $`\mu \tau `$Higgs vertex. In this paper, another signature, $`\tau \mu \eta `$ is considered. In the large $`\mathrm{tan}\beta `$ region, it is shown that the branching ratio of $`\tau \mu \eta `$ is 8.4 times the branching ratio of $`\tau 3\mu `$, independent of any unknown parameters, and it will thus give the most stringent bound on Higgs-mediated lepton flavor violation, and may provide its first signature. In the other regions of parameter space, where $`\tau \mu \gamma `$ is the most prominent decay, the branching ratio for $`\tau \mu \eta `$ is always substantially lower.

$`\tau \mu \eta `$ in Supersymmetric Models

Marc Sher
Nuclear and Particle Theory Group,

College of William and Mary, Williamsburg, Virginia 23187

The flavor physics of quarks and leptons is one of the most prominent mysteries in particle physics. The most surprising development in flavor physics in the past decade has been the observation of very large mixing between muon and tau neutrinos. The mixing will, at some level, lead to mixing in the charged lepton sector, giving violations of muon and tau lepton number conservation.

In the Standard Model supplemented with right-handed neutrinos, such violation is generally small. However, a much bigger effect will occur in supersymmetric models, by inducing mixing between the scalar muon and the scalar tau. In the most general supersymmetric standard model, even if neutrino mixing were not present, one would have large mixing between all scalar leptons and scalar quarks. This would lead to very large flavor (quark and lepton) changing neutral currents, which are not observed. It is thus generally assumed that the squark and slepton masses are equal at the unification scale; an assumption that occurs naturally in many models, such as supergravity or gauge-mediated supersymmetry breaking.

Even if the masses of the smuon and stau are equal at the unification scale, without mixing, the presence of non-diagonal neutrino mass terms (either from different Dirac mass terms or non-diagonal right-handed neutrino mass terms) will affect the masses through renormalization group running, and will generate mixing terms. In this paper, we will consider only mixing between left-handed smuons and staus; in most models, such mixing is the largest. It should be kept in mind that solar neutrino oscillation experiments also indicate large mixing between muon and electron neutrinos, and in models with inverted hierarchies, there could be substantial mixing between left-handed smuons and selectrons, although the very strong bounds on $`\mu e\gamma `$ will constrain these effects.

The mixing is characterized by the parameter $`\delta _{23}M_{23}^2/\stackrel{~}{m}^2`$, where $`M_{23}`$ is the off-diagonal term in the slepton mass matrix and $`\stackrel{~}{m}`$ is the slepton mass scale. The value of $`\delta _{23}`$ is extremely model-dependent, of course, but in models in which the mixing arises entirely through renormalization group running, its value is typically between $`0.1`$ and $`1`$.

The most studied tau-number and muon-number violating process is $`\tau \mu \gamma `$. Many of these works consider various models for $`\delta _{23}`$. Normalizing the rate to the current bound

$$BR(\tau \mu \gamma )=1.1\times 10^6\left(\frac{\delta _{23}}{1.4}\right)^2\left(\frac{100\mathrm{GeV}}{\stackrel{~}{m}}\right)^2$$
(1)
In the constrained MSSM, where the parameter space is restricted to a manageable dimensionality, this process dominates in the low $`m_{1/2}`$ region.

In addition, one can consider tau-number and muon-number violation in production and decay of sleptons, neutralinos and charginos. As an example, the process $`\chi _2\chi _1\mu \tau `$ where $`\chi _{1,2}`$ are neutralinos, can be searched for at the LHC. In the constrained MSSM, this process dominates in the $`m_{1/2}>m_0`$ region of parameter-space, and is thus complementary to $`\tau \mu \gamma `$.

Recently, Babu and Kolda pointed out that $`\tau 3\mu `$ was a promising signature in models with a large value of $`\mathrm{tan}\beta `$. Earlier, they had noticed that squark mixing would induce a flavor non-diagonal quark-quark-Higgs Yukawa coupling, and had examined the consequences for rare B-decays. The same mechanism, however, will also induce a $`\mu \tau \mathrm{Higgs}`$ vertex, and thus directly to $`\tau 3\mu `$, through tree-level Higgs exchange (either the $`h`$, $`H`$, or $`A`$). The branching ratio, for a reasonable choice of mass parameters, is

$$\mathrm{BR}(\tau 3\mu )=(1\times 10^7)\times \left(\frac{\mathrm{tan}\beta }{60}\right)^6\times \left(\frac{100\mathrm{GeV}}{m_A}\right)^4$$

where $`m_A`$ is the pseudoscalar mass. This result is very insensitive to the SUSY spectrum, with the exception that it can increase by up to a factor of four for large $`\mu `$. This branching ratio will be accessible at B-factories in the near future.

In this paper, another signature of lepton-number violation is discussed: $`\tau \mu \eta `$. It will be shown that, in the large $`\mathrm{tan}\beta `$ region discussed in the previous paragraph, the branching ratio is much higher than $`\tau 3\mu `$, and may be a much more sensitive test of Higgs-induced lepton flavor violation.

In the Babu-Kolda model, the lepton flavor violating interaction is given by

$$_{LFV}=(2G_F^2)^{1/4}\frac{m_\tau \kappa _{32}}{\mathrm{cos}^2\beta }(\overline{\tau }_R\mu _L)[H^0+iA^0]+\mathrm{h}.\mathrm{c}.$$
(2)
where $`H^0`$ and $`A^0`$ are the heavier scalar and the pseudoscalar, respectively, and we have chosen the generally preferred region of parameter-space in which $`\mathrm{sin}(\alpha \beta )1`$. Here, $`\kappa _{32}`$ is given in Ref.  and depends on loop integrals (but is relatively insensitive to SUSY parameters). For SUSY parameters $`\mu =M_1=m_2=m_{\stackrel{~}{l}}=m_{\stackrel{~}{nu}}`$, $`M_R=10^{14}`$ GeV and the off-diagonal Dirac neutrino coupling equal to the top quark Yukawa coupling (as expected in $`SO(10)`$ models), one has $`\kappa _{32}=4\times 10^4`$.

With this interaction, one can have a $`\tau `$ convert into a $`\mu `$ and a virtual $`H^0`$ or $`A^0`$, which then converts into a $`\mu ^+\mu ^{}`$ pair. This gives the branching ratio mentioned above. However, one could equally well have an $`A^0`$ convert into a strange quark pair, which then becomes an $`\eta `$, as shown in Figure 1, giving $`\tau \mu \eta `$. This will have both advantages and disadvantages. The two body phase space is a major advantage, and the extra color factor and slightly bigger coupling (since $`m_s>m_\mu `$) are also advantages. The disadvantage is in converting the strange quarks into an $`\eta `$.

A general discussion of $`\tau \mu \eta `$ can be found in Ref. . The relevant matrix element is

$$<0|\overline{s}\gamma _5s|\eta >=\sqrt{6}F_\eta ^8\frac{m_\eta ^2}{m_u+m_d+4m_s}$$
(3)
Using this matrix element, it is straightforward to calculate the branching ratio. If one divides by the $`\tau 3\mu `$ branching ratio, the unknown parameters all cancel, and the result is (neglecting the muon mass)

$$\frac{\mathrm{\Gamma }(\tau \mu \eta )}{\mathrm{\Gamma }(\tau 3\mu )}=54\pi ^2\left(\frac{F_\eta ^8}{m_\mu }\right)^2\left(\frac{m_\eta }{m_\tau }\right)^4\left(1\frac{m_\eta ^2}{m_\tau ^2}\right)^2$$
(4)
With $`F_\eta ^8150`$ MeV, this ratio is 8.4, giving a branching ratio of

$$\mathrm{BR}(\tau \mu \eta )=(0.84\times 10^6)\times \left(\frac{\mathrm{tan}\beta }{60}\right)^6\times \left(\frac{100\mathrm{GeV}}{m_A}\right)^4$$

One can get this result approximately without doing a calculation. Imagine that final state interactions are turned off and that the strange quarks propagate as free particles. One expects the ratio of $`\tau \mu \overline{s}s`$ to $`\tau 3\mu `$ to have a factor of 3 for color and a factor of $`(m_s/m_\mu )^2`$ for the Yukawa coupling. The cross diagram in the muon case turns out to lower the rate by a factor of $`3/2`$, so the overall ratio is $`\frac{9m_s^2}{2m_\mu ^2}10`$. Since the only two body decays would be $`\mu \eta `$ and $`\mu \eta ^{}`$, and the latter is suppressed much more by phase space, the $`\mu \eta `$ rate should dominate this process.

Since the experimental bound on the branching ratio for $`\tau \mu \eta `$ is $`9.6\times 10^6`$, and that for $`\tau 3\mu `$ is $`1.9\times 10^6`$, it is clear that $`\tau \mu \eta `$ puts stronger constraints on the model. In order to reach the interesting region of parameter space ($`\mathrm{tan}\beta 60`$ and $`m_A100`$ GeV), the bound on $`\tau 3\mu `$ would need to be improved by a factor of $`20`$, whereas the bound on $`\tau \mu \eta `$ would need to be improved by a factor of $`10`$.

Could these improvements be made? Both the $`\tau 3\mu `$ and $`\tau \mu \eta `$ bounds are based on the CLEO-II sample of $`4.7\mathrm{fb}^1`$, and the CLEO experiment has now accumulated a total of 5 times that luminosity (which would give a total of 24 million tau pairs). So in the absence of backgrounds, the current bounds could improve by a factor of five. The efficiency for $`\tau 3\mu `$ is listed as $`15\%`$; the efficiency for $`\tau \mu \eta `$ is about $`3\%`$, when one includes the fact that they only search for the $`\gamma \gamma `$ channel for the $`\eta `$ (thus the factor of five difference in the current bounds). Including the three-pion decay, or increasing the fiducial area for finding photons, could improve that efficiency substantially. Over the next years, BABAR and BELLE will reach $`500\mathrm{fb}^1`$, which could easily reach the needed sensitivity, depending on the point at which they become background-limited. Even if $`\mathrm{tan}\beta `$ is somewhat smaller, or $`m_A`$ larger, the necessary sensitivity could possibly be reached at LHC, SuperKEKB or a tau-charm factory. Note that the $`\tau 3\mu `$ decay could still be dominant in the small region of parameter space in which $`M_H<<M_A`$.

Are there any other processes that can give $`\tau \mu \eta `$? One can have the box diagram of Figure 2, which will also yield $`\tau `$ decays into $`\mu `$ plus other mesons, including the $`\pi ,\rho `$ and $`\varphi `$. If we take the special case in which the neutralinos are pure photino, the rate for $`\tau \mu \pi `$ is given by

$$\frac{\mathrm{\Gamma }(\tau \mu \pi )}{\mathrm{\Gamma }(\tau \mu \gamma )}=\frac{32\pi \alpha F_\pi ^2m_{\stackrel{~}{l}}^8(I_1^2+I_2^2)}{81m_\tau ^2m_{\stackrel{~}{\gamma }}^8M_3^2(x)}$$
(5)
where $`xm_{\stackrel{~}{\gamma }}/m_{\stackrel{~}{l}}`$, $`m_{\stackrel{~}{\gamma }}`$ is the photino mass, $`m_{\stackrel{~}{l}}`$ is the average slepton mass, $`M_3(x)`$ is given in Ref. , and the integrals are

$$(I_1,I_2)=\frac{1}{16\pi ^2}_0^{\mathrm{}}𝑑x\frac{(x,\frac{1}{4}x^2)}{(x1)^2(xa)(xb)(xc)}$$
(6)
Here, $`a,b,c`$ are the ratios of the squark,smuon and stau masses to the photino mass. In giving this expression, we have used the fact that $`m_\pi ^2/(m_u+m_d)`$ is numerically close to $`m_\tau `$. In the case of $`\tau \mu \eta `$, there will be a suppression of a factor of $`4`$ from the $`s`$-quark charge, an increase of a factor of $`(F_\eta ^8/F_\pi )^23.0`$ from the decay contants, and the coefficient of the $`I_1^2`$ term will be decreased by a factor of $`(3m_\eta ^4/8m_s^2)`$ relative to $`m_\tau ^2`$, or a factor of about $`3`$. This ratio has been evaluated for the entire SUSY parameter space, assuming sparticle masses in the range of $`51000`$ GeV, and is always less than $`10^2`$. As a result, $`\tau \mu +`$meson arising from the box diagram is negligible.

The existence of large $`\nu _\mu \nu _\tau `$ mixing implies mixing between the left-handed smuon and stau. While one could look for this mixing directly in neutralino/slepton interactions, one can also look at $`\tau `$ decays. The decay $`\tau \mu \gamma `$ is one signature, however Babu and Kolda have noted the mixing will also lead, especially in the large $`\mathrm{tan}\beta `$ region, to $`\tau 3\mu `$. In this paper, it has been pointed out that $`\tau \mu \eta `$ will also occur in this Higgs-mediated model, with a branching ratio $`8.4`$ times bigger, and is thus more sensitive. In other models, where $`\tau \mu \gamma `$ is the main signature, the $`\tau \mu \eta `$ rate is substantially smaller.

I thank Carl Carlson, Chris Carone and Jon Urheim for useful discussions.


text:
# Nonlinear Spinor Field in Anisotropic Universes*footnote **footnote * Talk, given at the International Conference “Scientific Reading devoted to 90 years anniversary of Professor Yakov Petrovich Terletskii”, July 1-3, 2002, Russian Peoples’ Friendship University, Moscow, Russia.

## 1 Introduction

Several authors studied the nonlinear spinor fields (NLSF) since Ivanenko  showed that a relativistic theory imposes in some cases a fourth order self-coupling. Nonlinear spinor field in an anisotropic universe, namely in a Bianchi-type I universe (BI) is studied by us in a series of papers . In these papers we considered the nonlinear spinor field, as well as a system of interacting spinor and scalar fields. Beside the spinor fields, we also study the role of a $`\mathrm{\Lambda }`$ term in the evolution of the universe. For the details of a NLSF in a BI universe one can consult . In the light of results obtained in the previously mentioned papers the study of NLSF in other anisotropic universes presents great interest. In this report we consider the self-consistent system of nonlinear spinor and BVI gravitational fields. The results have been compared with those obtained for the BI universe.

## 2 Fundamental equations and general solutions

We choose the Lagrangian for the self-consistent system of spinor and gravitational fields in the form

$$L=\frac{R}{2\kappa }+\frac{i}{2}\left[\overline{\psi }\gamma ^\mu _\mu \psi _\mu \overline{\psi }\gamma ^\mu \psi \right]M\overline{\psi }\psi +L_N$$
(1)
with $`R`$ being the scalar curvature and $`\kappa `$ being the Einstein’s gravitational constant. The nonlinear term $`L_N`$ describes the self-interaction of spinor field and can be presented as some arbitrary functions of invariants generated from the real bilinear forms of spinor field having the form:

$`S`$ $`=`$ $`\overline{\psi }\psi (\mathrm{scalar}),`$ (3)
$`P`$ $`=`$ $`i\overline{\psi }\gamma ^5\psi (\mathrm{pseudoscalar}),`$ (4)
$`v^\mu `$ $`=`$ $`(\overline{\psi }\gamma ^\mu \psi )(\mathrm{vector}),`$ (5)
$`A^\mu `$ $`=`$ $`(\overline{\psi }\gamma ^5\gamma ^\mu \psi )(\mathrm{pseudovector}),`$ (6)
$`Q^{\mu \nu }`$ $`=`$ $`(\overline{\psi }\sigma ^{\mu \nu }\psi )(\mathrm{antisymmetric}\mathrm{tensor}),`$ (7)
where $`\sigma ^{\mu \nu }=(i/2)[\gamma ^\mu \gamma ^\nu \gamma ^\nu \gamma ^\mu ]`$. Invariants, corresponding to the bilinear forms, are

$`I`$ $`=`$ $`S^2,`$ (9)
$`J`$ $`=`$ $`P^2,`$ (10)
$`I_v`$ $`=`$ $`v_\mu v^\mu =(\overline{\psi }\gamma ^\mu \psi )g_{\mu \nu }(\overline{\psi }\gamma ^\nu \psi ),`$ (11)
$`I_A`$ $`=`$ $`A_\mu A^\mu =(\overline{\psi }\gamma ^5\gamma ^\mu \psi )g_{\mu \nu }(\overline{\psi }\gamma ^5\gamma ^\nu \psi ),`$ (12)
$`I_Q`$ $`=`$ $`Q_{\mu \nu }Q^{\mu \nu }=(\overline{\psi }\sigma ^{\mu \nu }\psi )g_{\mu \alpha }g_{\nu \beta }(\overline{\psi }\sigma ^{\alpha \beta }\psi ).`$ (13)
According to the Pauli-Fierz theorem  among the five invariants only $`I`$ and $`J`$ are independent as all other can be expressed by them: $`I_V=I_A=I+J`$ and $`I_Q=IJ.`$ Therefore, we choose the nonlinear term $`L_N=F(I,J)`$, thus claiming that it describes the nonlinearity in the most general of its form.

We choose the anisotropic inhomogeneous universe given by a Bianchi-type VI (BVI) metric:

$$ds^2=dt^2a^2e^{2mz}dx^2b^2e^{2nz}dy^2c^2dz^2,$$
(14)
with $`a,b,c`$ being the functions of time only. Here $`m,n`$ are some arbitrary constants and the velocity of light is taken to be unity. Note that the suitable choice of $`m,n`$ as well as the metric functions $`a,b,c`$ in the BVI given by (14) evokes the following Bianchi-type universes. Thus

* for $`m=n`$ the BVI metric transforms to a Bianchi-type V (BV) one, i.e.,
$`m=n`$, BVI $``$ BV $``$ open FRW;
* for $`n=0`$ the BVI metric transforms to a Bianchi-type III (BIII) one, i.e.,
$`n=0`$, BVI $``$ BIII;
* for $`m=n=0`$ the BVI metric transforms to a Bianchi-type I (BI) one, i.e.,
$`m=n=0`$, BVI $``$ BI;
* for $`m=n=0`$ and equal scale factor in all three direction the BVI metric transforms to a Friedmann-Robertson-Walker (FRW) universe, i.e.,
$`m=n=0`$ and $`a=b=c`$, BVI $``$ FRW.

Variation of the Lagrangian (1) with respect to field functions $`\psi (\overline{\psi })`$ gives the nonlinear spinor field equations:

$`i\gamma ^\mu _\mu \psi M\psi +𝒟\psi +i𝒢\gamma ^5\psi `$ $`=`$ $`0`$ (16)
$`i_\mu \overline{\psi }\gamma ^\mu +M\overline{\psi }𝒟\overline{\psi }i𝒢\overline{\psi }\gamma ^5`$ $`=`$ $`0`$ (17)
where $`𝒟=2SF_I`$ and $`𝒢=2PF_J`$.

Varying (1) with respect to metric function ($`g_{\mu \nu }`$) we find the Einstein equations for $`a,b,c`$ which for the metric (14) read

$`{\displaystyle \frac{\ddot{b}}{b}}+{\displaystyle \frac{\ddot{c}}{c}}+{\displaystyle \frac{\dot{b}}{b}}{\displaystyle \frac{\dot{c}}{c}}{\displaystyle \frac{n^2}{c^2}}`$ $`=`$ $`\kappa T_1^1,`$ (19)
$`{\displaystyle \frac{\ddot{c}}{c}}+{\displaystyle \frac{\ddot{a}}{a}}+{\displaystyle \frac{\dot{c}}{c}}{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{m^2}{c^2}}`$ $`=`$ $`\kappa T_2^2,`$ (20)
$`{\displaystyle \frac{\ddot{a}}{a}}+{\displaystyle \frac{\ddot{b}}{b}}+{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{\dot{b}}{b}}+{\displaystyle \frac{mn}{c^2}}`$ $`=`$ $`\kappa T_3^3,`$ (21)
$`{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{\dot{b}}{b}}+{\displaystyle \frac{\dot{b}}{b}}{\displaystyle \frac{\dot{c}}{c}}+{\displaystyle \frac{\dot{c}}{c}}{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{m^2mn+n^2}{c^2}}`$ $`=`$ $`\kappa T_0^0,`$ (22)
$`m{\displaystyle \frac{\dot{a}}{a}}n{\displaystyle \frac{\dot{b}}{b}}(mn){\displaystyle \frac{\dot{c}}{c}}`$ $`=`$ $`\kappa T_3^0.`$ (23)
Here overdots denote differentiation with respect to time ($`t`$) and $`T_\mu ^\nu `$ is the energy-momentum tensor of the material field and has the form

$$T_\mu ^\rho =\frac{i}{4}g^{\rho \nu }\left(\overline{\psi }\gamma _\mu _\nu \psi +\overline{\psi }\gamma _\nu _\mu \psi _\mu \overline{\psi }\gamma _\nu \psi _\nu \overline{\psi }\gamma _\mu \psi \right)\delta _\mu ^\rho L_{sp}.$$
(24)
Here $`L_{sp}`$ is the spinor field Lagrangian, which on account of spinor field equations (2) takes the form:

$$L_{sp}=𝒟S𝒢𝒫+F.$$
(25)
In the expressions above $`_\mu `$ denotes the covariant derivative of spinor, having the form :

$$_\mu \psi =_\mu \psi \mathrm{\Gamma }_\mu \psi $$
(26)
where $`\mathrm{\Gamma }_\mu (x)`$ are spinor affine connection matrices defined by the equality

$$\mathrm{\Gamma }_\mu (x)=(1/4)g_{\rho \sigma }(x)\left(_\mu e_\delta ^be_b^\rho \mathrm{\Gamma }_{\mu \delta }^\rho \right)\gamma ^\sigma \gamma ^\delta .$$

For the metric element (14) it gives

$`\mathrm{\Gamma }_0`$ $`=`$ $`0,`$ (27)
$`\mathrm{\Gamma }_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\dot{a}\overline{\gamma }^1\overline{\gamma }^0m{\displaystyle \frac{a}{c}}\overline{\gamma }^1\overline{\gamma }^3\right]e^{mz}`$ (28)
$`\mathrm{\Gamma }_2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\dot{b}\overline{\gamma }^2\overline{\gamma }^0+n{\displaystyle \frac{b}{c}}\overline{\gamma }^2\overline{\gamma }^3\right]e^{nz},`$ (29)
$`\mathrm{\Gamma }_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{c}\overline{\gamma }^3\overline{\gamma }^0`$ (30)
It is easy to show that

$$\gamma ^\mu \mathrm{\Gamma }_\mu =\frac{1}{2}\frac{\dot{\tau }}{\tau }\overline{\gamma }^0+\frac{mn}{2c}\overline{\gamma }^3,$$

where we define

$$\tau =abc.$$
(31)
The Dirac matrices $`\gamma ^\mu (x)`$ of curved space-time are connected with those of Minkowski one as follows:

$$\gamma ^0=\overline{\gamma }^0,\gamma ^1=\overline{\gamma }^1e^{mz}/a,\gamma ^2=\overline{\gamma }^2/be^{nz},\gamma ^3=\overline{\gamma }^3/c$$

with

$`\overline{\gamma }^0=\left(\begin{array}{cc}I& 0\\ 0& I\end{array}\right),\overline{\gamma }^i=\left(\begin{array}{cc}0& \sigma ^i\\ \sigma ^i& 0\end{array}\right),\gamma ^5=\overline{\gamma }^5`$ $`=`$ $`\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right),`$ (38)
where $`\sigma _i`$ are the Pauli matrices:

$`\sigma ^1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma ^2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma ^3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).`$ (45)
Note that the $`\overline{\gamma }`$ and the $`\sigma `$ matrices obey the following properties:

$`\overline{\gamma }^i\overline{\gamma }^j+\overline{\gamma }^j\overline{\gamma }^i=2\eta ^{ij},i,j=0,1,2,3`$ (46)
$`\overline{\gamma }^i\overline{\gamma }^5+\overline{\gamma }^5\overline{\gamma }^i=0,(\overline{\gamma }^5)^2=I,i=0,1,2,3`$ (47)
$`\sigma ^j\sigma ^k=\delta _{jk}+i\epsilon _{jkl}\sigma ^l,j,k,l=1,2,3`$ (48)
where $`\eta _{ij}=\{1,1,1,1\}`$ is the diagonal matrix, $`\delta _{jk}`$ is the Kronekar symbol and $`\epsilon _{jkl}`$ is the totally antisymmetric matrix with $`\epsilon _{123}=+1`$. Let us consider the spinors to be functions of $`t`$ and $`z`$ only, such that

$$\psi (t,z)=v(t)e^{ikz},\overline{\psi }(t,z)=\overline{v}(t)e^{ikz}$$
(49)
Inserting (49) into (2) for the nonlinear spinor field we find

$`\overline{\gamma }^0\left(\dot{v}+{\displaystyle \frac{\dot{\tau }}{2\tau }}v\right)\left({\displaystyle \frac{mn}{2c}}i{\displaystyle \frac{k}{c}}\right)\overline{\gamma }^3v+i\mathrm{\Phi }v+𝒢\overline{\gamma }^5v`$ $`=`$ $`0`$ (51)
$`\left(\dot{\overline{v}}+{\displaystyle \frac{\dot{\tau }}{2\tau }}\overline{v}\right)\overline{\gamma }^0\left({\displaystyle \frac{mn}{2c}}+i{\displaystyle \frac{k}{c}}\right)\overline{v}\overline{\gamma }^3i\mathrm{\Phi }\overline{v}𝒢\overline{v}\overline{\gamma }^5`$ $`=`$ $`0.`$ (52)
Here we define $`\mathrm{\Phi }=M𝒟`$. Introduce a new function $`u_j(t)=\sqrt{\tau }v_j(t)`$, for the components of the NLSF from (2) one obtains

$`\dot{u}_1+i\mathrm{\Phi }u_1\left[{\displaystyle \frac{mn}{2c}}i{\displaystyle \frac{k}{c}}+𝒢\right]u_3`$ $`=`$ $`0,`$ (54)
$`\dot{u}_2+i\mathrm{\Phi }u_2+\left[{\displaystyle \frac{mn}{2c}}i{\displaystyle \frac{k}{c}}𝒢\right]u_4`$ $`=`$ $`0,`$ (55)
$`\dot{u}_3i\mathrm{\Phi }u_3\left[{\displaystyle \frac{mn}{2c}}i{\displaystyle \frac{k}{c}}𝒢\right]u_1`$ $`=`$ $`0,`$ (56)
$`\dot{u}_4i\mathrm{\Phi }u_4+\left[{\displaystyle \frac{mn}{2c}}i{\displaystyle \frac{k}{c}}+𝒢\right]u_2`$ $`=`$ $`0.`$ (57)
Using the spinor field equations (2) and (2) it can be shown that the bilinear spinor forms, defined by (1), i.e.,

$`S`$ $`=`$ $`\overline{\psi }\psi =\overline{v}v,P=i\overline{\psi }\overline{\gamma }^5\psi =i\overline{v}\overline{\gamma }^5v,A^0=\overline{\psi }\overline{\gamma }^5\overline{\gamma }^0\psi =\overline{v}\overline{\gamma }^5\overline{\gamma }^0v,`$ (58)
$`A^3`$ $`=`$ $`\overline{\psi }\overline{\gamma }^5\overline{\gamma }^3\psi =\overline{v}\overline{\gamma }^5\overline{\gamma }^3v,V^0=\overline{\psi }\overline{\gamma }^0\psi =\overline{v}\overline{\gamma }^0v,V^3=\overline{\psi }\overline{\gamma }^3\psi =\overline{v}\overline{\gamma }^3v,`$ (59)
$`Q^{30}`$ $`=`$ $`i\overline{\psi }\overline{\gamma }^3\overline{\gamma }^0\psi =i\overline{v}\overline{\gamma }^3\overline{\gamma }^0v,Q^{21}=\overline{\psi }\overline{\gamma }^0\overline{\gamma }^3\overline{\gamma }^5\psi =i\overline{\psi }\overline{\gamma }^2\overline{\gamma }^1\psi =i\overline{v}\overline{\gamma }^2\overline{\gamma }^1v,`$ (60)
obeying the following system of equations:

$`\dot{S}_02{\displaystyle \frac{k}{c}}Q_0^{30}2𝒢A_0^0`$ $`=`$ $`0,`$ (62)
$`\dot{P}_02{\displaystyle \frac{k}{c}}Q_0^{21}2\mathrm{\Phi }A_0^0`$ $`=`$ $`0,`$ (63)
$`\dot{A}_0^0{\displaystyle \frac{mn}{c}}A_0^3+2\mathrm{\Phi }P_0+2𝒢S_0`$ $`=`$ $`0,`$ (64)
$`\dot{A}_0^3{\displaystyle \frac{mn}{c}}A_0^0`$ $`=`$ $`0,`$ (65)
$`\dot{V}_0^0{\displaystyle \frac{mn}{c}}V_0^3`$ $`=`$ $`0,`$ (66)
$`\dot{V}_0^3{\displaystyle \frac{mn}{c}}V_0^0+2\mathrm{\Phi }Q_0^{30}2𝒢Q_0^{21}`$ $`=`$ $`0,`$ (67)
$`\dot{Q}_0^{30}+2{\displaystyle \frac{k}{c}}S_02\mathrm{\Phi }V_0^3`$ $`=`$ $`0,`$ (68)
$`\dot{Q}_0^{21}+2{\displaystyle \frac{k}{c}}P_0+2𝒢V_0^3`$ $`=`$ $`0,`$ (69)
where we denote $`F_0=\tau F`$. Combining these equations together and taking the first integral one gets

$`(S_0)^2+(P_0)^2+(A_0^0)^2(A_0^3)^2+(V_0^0)^2+(V_0^3)^2+(Q_0^{30})^2+(Q_0^{21})^2=C=\mathrm{Const}`$ (71)
Before dealing with the Einstein equations (2) let us go back to (2). From the first and the third equations of the system (2) one finds

$$\dot{u}_{13}=(𝒢Q)u_{13}^22i\mathrm{\Phi }u_{13}+(𝒢+Q),$$
(72)
where, we denote $`u_{13}=u_1/u_3`$ and $`Q=[mn2ik]/2c`$. The equation (72) is a Riccati one  with variable coefficients. Further, setting $`u_{13}=v_{13}\mathrm{exp}[2i\mathrm{\Phi }(t)𝑑t]`$, from (72) one obtains

$$\dot{v}_{13}=(𝒢Q)v_{13}^2e^{2i{\scriptscriptstyle \mathrm{\Phi }(t)𝑑t}}+(𝒢+Q)e^{2i{\scriptscriptstyle \mathrm{\Phi }(t)𝑑t}}.$$
(73)
The general solution of (73) can be written as

$$v_{13}=\left[(𝒢Q)e^{2i{\scriptscriptstyle \mathrm{\Phi }(t)𝑑t}}𝑑t+C(t)\right]^1,$$
(74)
with the integration constant be defined from

$$\left[(𝒢Q)e^{2i{\scriptscriptstyle \mathrm{\Phi }(t)𝑑t}}𝑑t+C(t)\right]^2𝑑C=(𝒢+Q)e^{2i{\scriptscriptstyle \mathrm{\Phi }(t)𝑑t}}𝑑t.$$
(75)
Thus given the nonlinear term in the Lagrangian and solution of the Einstein equations, one finds the relation between $`u_1`$ and $`u_3`$ ($`u_2`$ and $`u_4`$ as well), hence the components of the spinor field.

Now we study the Einstein equation (2). In doing so, we write the components of the energy-momentum tensor, which in our case read

$`T_0^0`$ $`=`$ $`mSF+{\displaystyle \frac{k}{c}}V^3,`$ (77)
$`T_1^1`$ $`=`$ $`T_2^2=𝒟S+𝒢PF,`$ (78)
$`T_3^3`$ $`=`$ $`𝒟S+𝒢PF{\displaystyle \frac{k}{c}}V^3,`$ (79)
$`T_3^0`$ $`=`$ $`kV^0.`$ (80)
Let us demand the energy-momentum tensor to be conserved, i.e.,

$$T_{\nu ;\mu }^\mu =T_{\nu ,\mu }^\mu +\mathrm{\Gamma }_{\beta \mu }^\mu T_\nu ^\beta \mathrm{\Gamma }_{\nu \mu }^\beta T_\beta ^\mu =0$$
(81)
Taking into account that $`T_\mu ^\nu `$ is a function of $`t`$ only, from (81) we find

$`\mathrm{\Phi }\dot{S}_0𝒢\dot{P}_0+{\displaystyle \frac{k}{c}}\dot{V}_0^3{\displaystyle \frac{k}{c}}{\displaystyle \frac{mn}{c}}V_0^0`$ $`=`$ $`0,`$ (83)
$`\dot{V}_0^0{\displaystyle \frac{mn}{c}}V_0^3`$ $`=`$ $`0.`$ (84)
As one can easily verify, the equations (81) are consistent with those of (2).

Let us now deal with the Einstein equations (2). In view of (2), from (23) one obtains

$$a^mb^n/c^{mn}=𝒩\mathrm{exp}[\kappa kV^0𝑑t],𝒩=\mathrm{const}.,$$
(85)
whereas, subtraction of (20) and (19) leads to

$$\frac{d}{dt}\left[\tau \frac{d}{dt}\left\{\mathrm{ln}\left(\frac{a}{b}\right)\right\}\right]=\frac{m^2n^2}{c^2}\tau .$$
(86)
Note that, the two other equations, obtained by subtracting (21) from (19) and (21) from (20), respectively, are identical to (86), that can be easily verified inserting (23) and (2) into the equations in question. Finally, summation of (19), (20), (21) and (22), multiplied by 3, leads to the equation for $`\tau `$, which in view of (2) takes the form

$$\frac{\ddot{\tau }}{\tau }=2\frac{m^2mnn^2}{c^2}\frac{\kappa }{2}\left[3(MS+𝒟S+𝒢P2F)+2\frac{k}{c}V^3\right].$$
(87)
For the right-hand side of the equation (87) are some functions of $`\tau `$ only, the solution to this equation is well known . It can be shown that the quantities, related to the spinor field are indeed the functions of $`\tau `$. Now, if we assume the metric function $`c`$ also to be a function of $`\tau `$, both (86) and (87) can be solved explicitly. Thus both the nonlinear spinor and the gravitational field equations are solved in general. In what follows, we consider some special cases and compare the solutions obtained with those for a BI universe. In can be emphasized that the introduction of inhomogeneity both in gravitational (through $`m`$ and $`n`$) and spinor (through $`k`$) significantly complicates the whole picture.

To begin with we consider the linear case setting $`F(I,J)=0`$. As one sees, the structures of the spinor field equation (2), as well as the gravitational ones (86) and (87) in this case remain unaltered. Thus, the removal of the nonlinearity has little to offer in our cause. As it was mentioned earlier, the introduction of inhomogeneity is in the root of all these troubles. So, for some break-through we demand the spinor field completely space-independent setting $`k=0`$. In this case for the components of the energy-momentum tensor immediately we find

$`T_0^0`$ $`=`$ $`mSF,`$ (89)
$`T_1^1`$ $`=`$ $`T_2^2=T_3^3=𝒟S+𝒢PF,`$ (90)
$`T_3^0`$ $`=`$ $`0.`$ (91)
In view of (2) from (85) we obtain

$$a^mb^n/c^{mn}=𝒩,$$
(92)
but the equations (2), (86) and (87) are still unchanged. Thus we see that even the elimination of spinor field inhomogeneity has little to offer.

Finally, we consider the case, when the spinor field is independent of space i.e., $`k=0`$ and the gravitational field is given by a BV space-time with $`m=n`$ in (14). In this case for the spinor field we find

$`\dot{u}_1+i\mathrm{\Phi }u_1𝒢u_3`$ $`=`$ $`0,`$ (94)
$`\dot{u}_2+i\mathrm{\Phi }u_2𝒢u_4`$ $`=`$ $`0,`$ (95)
$`\dot{u}_3i\mathrm{\Phi }u_3+𝒢u_1`$ $`=`$ $`0,`$ (96)
$`\dot{u}_4i\mathrm{\Phi }u_4+𝒢u_2`$ $`=`$ $`0.`$ (97)
The bilinear spinor forms in this case satisfy

$`\dot{S}_02𝒢A_0^0`$ $`=`$ $`0,`$ (99)
$`\dot{P}_02\mathrm{\Phi }A_0^0`$ $`=`$ $`0,`$ (100)
$`\dot{A}_0^0+2\mathrm{\Phi }P_0+2𝒢S_0`$ $`=`$ $`0,`$ (101)
$`\dot{A}_0^3`$ $`=`$ $`0,`$ (102)
$`\dot{V}_0^0`$ $`=`$ $`0,`$ (103)
$`\dot{V}_0^3+2\mathrm{\Phi }Q_0^{30}2𝒢Q_0^{21}`$ $`=`$ $`0,`$ (104)
$`\dot{Q}_0^{30}2\mathrm{\Phi }V_0^3`$ $`=`$ $`0,`$ (105)
$`\dot{Q}_0^{21}+2𝒢V_0^3`$ $`=`$ $`0,`$ (106)
with the relations

$`(S_0)^2+(P_0)^2+(A_0^0)^2`$ $`=`$ $`B_1,`$ (108)
$`A_0^3`$ $`=`$ $`B_2,`$ (109)
$`V_0^0`$ $`=`$ $`B_3,`$ (110)
$`(V_0^3)^2+(Q_0^{30})^2+(Q_0^{21})^2`$ $`=`$ $`B_4,`$ (111)
with $`B_i`$ being the constant of integration.

For the gravitational field we find,

$$a=(𝒩)^{(1/m)}b,$$
(112)
with the equations

$$\frac{d}{dt}\left[\tau \frac{d}{dt}\left\{\mathrm{ln}\left(\frac{b}{c}\right)\right\}\right]=\frac{2m^2}{c^2}\tau ,$$
(113)
and

$$\frac{\ddot{\tau }}{\tau }=2\frac{m^2}{c^2}\frac{\kappa }{2}\left[3(MS+𝒟S+𝒢P2F)\right].$$
(114)
Note that the spinor field equation (2) completely coincides with those for a BI metric. If the nonlinear term in the Lagrangian is given as $`F=F(I)`$, then the components of the spinor field can be given as

$`\psi _1(t)`$ $`=`$ $`(C_1/\sqrt{\tau })\mathrm{exp}[i{\displaystyle (M𝒟)𝑑t}],`$ (116)
$`\psi _2(t)`$ $`=`$ $`(C_2/\sqrt{\tau })\mathrm{exp}[i{\displaystyle (M𝒟)𝑑t}],`$ (117)
$`\psi _3(t)`$ $`=`$ $`(C_3/\sqrt{\tau })\mathrm{exp}[i{\displaystyle (M𝒟)𝑑t}],`$ (118)
$`\psi _4(t)`$ $`=`$ $`(C_4/\sqrt{\tau })\mathrm{exp}[i{\displaystyle (M𝒟)𝑑t}],`$ (119)
with $`C_1,C_2,C_3,C_4`$ being the integration constants, such that

$$C_1^2+C_2^2C_3^2C_4^2=C_0,$$

with

$$S=C_0/\tau .$$

In case, the nonlinear term is given by $`F=F(J)`$, the components of the spinor field have the form

$`\psi _1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\tau }}}\left(D_1e^{i\sigma }+iD_3e^{i\sigma }\right),`$ (121)
$`\psi _2`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\tau }}}\left(D_2e^{i\sigma }+iD_4e^{i\sigma }\right),`$ (122)
$`\psi _3`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\tau }}}\left(iD_1e^{i\sigma }+D_3e^{i\sigma }\right),`$ (123)
$`\psi _4`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\tau }}}\left(iD_2e^{i\sigma }+D_4e^{i\sigma }\right).`$ (124)
Here $`\sigma =int𝒢dt`$, and the integration constants $`D_i`$ obey

$$2(D_1^2+D_2^2D_3^2D_4^2)=D_0,$$

with $`D_0`$ to be determined from

$$P=D_0/\tau .$$

Contrary to the BI metric, where the metric functions $`a,b,c`$ and $`\tau `$ are easily determined given the concrete form of nonlinearity, in the case considered, the process is much complicated with the space inhomogeneity being an active player. Thus we see that the introduction of inhomogeneity in the space-time has a far reaching effect in the evolution of both spinor and gravitational field.