text:
# Observation of non-gaussian conductance fluctuations at low temperatures in Si:P(B) at the metal-insulator transition

## Abstract

We report investigations of conductance fluctuations (noise) in doped silicon at low temperatures (T$`<20`$K) as it is tuned through the metal-insulator transition (MIT). The scaled magnitude of noise, $`\gamma _H`$, increases with decrease in T following an approximate power law $`\gamma _HT^\beta `$. At low T, $`\gamma _H`$ diverges as $`n/n_c`$ crosses 1 from the metallic side. We find that the distribution function and second spectrum of the fluctuations show strong non-gaussian behavior below 20K as $`n/n_c`$ decreases through 1. In particular, the observed distribution function which is gaussian for $`n/n_c>>1`$, develops a log-normal tail as the transition is approached from the metallic side and eventually it dominates in the critical region.

The Anderson-Mott metal-insulator transition (MIT) has been one of the most studied areas in condensed matter . However, even after several decades of work , there remain many unanswered questions. A fundamental issue of considerable theoretical interest is that of conductance fluctuations. The two main aspects of this issue are the divergence of noise near the critical region and the onset of non-gaussian fluctuations. The origin of conductance fluctuation in the weak localization regime is dominated mainly by the Universal Conductance fluctuations(UCF), and is distinctly different from that in the strongly localized insulating regime ($`\sigma _{T0}=0`$) which show activated conductance. There are some experimental investigations of the magnitude of conductance fluctuations in thin films of $`In_xO_y`$  at room temperature and in Si inversion layers  at low temperatures which showed a divergence near the MIT. However, no experimental investigations have looked into the question of non-gaussian fluctuations at the MIT. In this report we present the results of the first experimental investigations of non-gaussian conductance fluctuations($`f<10Hz`$) near the MIT at low temperatures. We have investigated the above aspects of conductance fluctuations in bulk single crystals of doped Si as it is driven through the Anderson-Mott transition. Specifically, we find that the fluctuations, which have a gaussian distribution in the metallic regime becomes strongly non-gaussian with a log-normal distribution at the transition region. This also depends strongly on temperature. We enumerate below the important results:

1. In a 3-D system like doped Si ( Si:P(B)) the noise magnitude shows clear signatures of divergence as $`n/n_c`$ decreases through 1 and $`T0`$.
2. As $`n/n_c1`$ and the system crosses over to the insulating regime, the conductance fluctuations becomes strongly non-gaussian as evident both from the second spectra (defined later) and Probability Distribution Function (PDF).
3. The PDF as well as the second spectra also depend on the temperature and both of them becomes more non-gaussian as T is decreased.

The samples (P doped and in some cases B compensated single crystals of Si) were taken from the same batch that were used earlier for transport studies near the MIT  with a range of doping levels $`n`$ from 3.3$`\times 10^{18}/cc`$ to $`1\times 10^{19}/cc`$, covering both sides of the critical doping level ($`n_c`$) for MIT (see table, figure 1). They were initially characterised by conductivity ($`\sigma `$) measurements down to $`T=0.1K`$. Figure 1 shows the variation of $`\sigma `$ as a function of $`T^{1/2}`$. When fitted to the simple relation $`\sigma `$ = $`\sigma _0`$ \+ m$`T^{1/2}`$, samples with $`\sigma _0>0`$ are taken as metallic. In the samples with $`n/n_c1`$ the conductivity for $`T<10K`$ follows Mott’s variable range hopping (VRH) with the relation $`\sigma (T)=\sigma _cexp((T_0/T)^{1/4})`$. $`T_00`$ as $`nn_c`$ from the insulating side.

Conductance fluctuations (noise) was measured as a function of T for $`2K<T<300K`$ using a five-probe ac technique . The time series of voltage fluctuations were collected directly into a computer using a 16-bit A/D card and was digitally processed to obtain the power spectra $`S_V(f)`$. A typical power spectra is calculated from a time series of 1-5$`\times 10^5`$ data points at fixed temperatures which was controlled to within $`\pm `$1mK. Typically the power spectra $`S_V(f)`$ $``$$`1/f^\alpha `$ ($`\alpha 1`$) and it scales inversely as the sample volume. At low measurement bias $`S_V(f)`$ scales as $`V^2`$. However, beyond certain bias values, deviation from $`V^2`$ dependence was seen at lower temperatures, and will be discussed later on.

In all the samples, the noise magnitude increases as the temperature is decreased below 100K. For $`T>100K`$ the noise increases with T again and this arises from more classical mechanisms  which are outside the scope of this paper. Figure 2 shows a plot of $`\frac{\mathrm{\Delta }G^2}{G^2}`$ = $`𝑑f.S_V(f)/V^2`$ integrated over the experimental bandwidth for a representative insulating sample K139 ($`n/n_c`$=0.94) for T$`<`$10K. When measured with low enough biasing field, $`\frac{\mathrm{\Delta }G^2}{G^2}`$ shows an approximate inverse power-law dependence $`A/T^\beta `$ with temperature. The behavior presented in figure 2 is typical for all the samples on both sides of the MIT, right down to the lowest measured T.

At lower temperatures in all the samples, the measured noise is sensitive to the biasing field used for the measurement as can be seen in figure 2. It shows how the T dependence can change drastically when measured with different biases. As T decreases, the bias dependence of $`\frac{\mathrm{\Delta }G^2}{G^2}`$ becomes more severe. On increasing the measuring bias beyond a characteristic electric field $`E^{}`$(T) (see figure 2 inset), $`\frac{\mathrm{\Delta }G^2}{G^2}`$ decreases drastically as a function of the measuring bias following an inverse power-law. $`E^{}`$(T) in all samples decreased with T. We find that sample heating can explain only a small part of this noise suppression. The bias dependence noise has a strong impact on the measured temperature dependence of noise, and is clearly seen in all samples. Unless measured with a bias $`E<E^{}`$, the noise tends to ”flatten out” or saturate at lower temperatures, as seen in figure 2. The cause for this effect is currently under investigation and will not be elaborated any further in this paper. All the results presented henceforth have been taken with biasing fields $`E<E^{}`$ such that $`S_V(f)V^2`$ and has no bias dependence.

Figure 3 shows a plot of the Hooge parameter  $`\gamma _H`$ (=$`f.NS_V(f)/V^2`$) as a function of $`n/n_c`$ for different temperatures for $`f=1Hz`$, $`N`$ being the total number of carriers in the measured volume. For $`n/n_c>1`$ (metallic regime), $`\gamma _H`$ is comparatively low. As the doping level decreases and $`n/n_c1`$ from the metallic side, $`\gamma _H`$ increases slowly. For $`n/n_c`$ 1, $`\gamma _H`$ starts increasing rapidly and at the MIT, changes drastically over several orders of magnitude as it goes over to the insulating side. On the insulating side, the observed $`\gamma _H`$ depends strongly on the temperature (see inset 3a). At higher temperatures ($`T6K`$), the value of $`\gamma _H`$ passes through a sharp peak at a value of $`n/n_c`$ very close to but less than 1. It then drops sharply in the insulating side . Interestingly this is roughly the same temperature ($`T10K`$) at which the resistivity starts to become activated and follows the VRH law.

The inset (b) of figure 3 shows the variation of $`\beta `$ as a function of the doping level in the system. $`\beta `$ also has a strong n dependence close to the MIT. At the highest doping levels, $`\beta 0.5`$, slowly increasing as the concentration is decreased and shows a rather sharp increase while it passes through the MIT to reach a value of 6 in the sample with the lowest doping level measured. The varitaion of $`\beta `$ in the metallic regime as the MIT is approached has been discussed elsewhere .

On the weakly localized metallic side ($`n/n_c2`$), we have shown earlier  that the predominant contribution to conductance fluctuations arise from the mechanism of UCF . In this region, the measured noise was found to be mostly gaussian (both PDF and the second spectra). However, there is strong departure from such a gaussian behavior as $`n`$ is reduced through $`n_c`$. A build up of non-gaussian fluctuations as $`n/n_c`$ is reduced through 1 at low temperaure ($`T=4.2K`$) can be seen from the probability distribution function (PDF). This is obtained from the time series and is the probability of occurence $`P(\mathrm{\Delta }V)`$ of the voltage jump of magnitude $`\mathrm{\Delta }V`$ when the conductivity fluctuates in a current biased sample. In all the samples (spanning different regions of $`n/n_c`$) the PDF calculated from the time series of the conductance fluctuations can be expressed as a sum of a gaussian and a log-normal distribution: $`P(\mathrm{\Delta }V)=A_ge^{1/2(\frac{\mathrm{\Delta }V}{\sigma _g})^2}+\frac{A_{ln}}{\mathrm{\Delta }V}e^{1/2(\frac{ln(\mathrm{\Delta }V)M}{\sigma _{ln}})^2}`$. $`A_g`$ and $`A_{ln}`$ are the relative weights of the two contributions. Figure 4a and b show the PDF of 3 representative samples at 4.2K (in a liquid He bath). For the metallic sample (K348) the distrubution is completely gaussian with $`A_{ln}`$ =0. However, close to the the transition for the sample K242 one can see a gaussian distribution with a log-normal tail. For the insulating sample K240, there is no gaussian component ($`A_g=0`$) and the PDF is a log-normal distribution. The variation in the relative weights of the two contributions, as a function of $`n/n_c`$ , can be seen from the integrated area under the two distributions ($`\mathrm{\Delta }_G`$ for the gaussian part and $`\mathrm{\Delta }_{LN}`$ for the log-normal part) normalized by the total area under the PDF ($`\mathrm{\Delta }_{total}`$) in figures 4(d)and (e). In figure 4(d) it can be seen that at $`T=4.2K`$, $`\mathrm{\Delta }_G`$/$`\mathrm{\Delta }_{total}=1`$ away from the transition but close to the t ransition $`\mathrm{\Delta }_G`$/$`\mathrm{\Delta }_{total}=0`$. A complimentary behavior is seen for $`\mathrm{\Delta }_{LN}`$/$`\mathrm{\Delta }_{total}`$ . Figure 4(e) (done for sample K240) shows that as T is reduced from the high temperature the gaussian component is reduced while the log-normal component takes over.

Another important measure of non-gaussianity is the second spectra, often referred to as the ”noise in noise”. If the fluctuations are uncorrelated, then the second spectrum  defined as $`S^{(2)}(f)`$ = $`(\mathrm{\Delta }S^{(1)}(f))^2_{time}`$ should have a frequency independent ”white” frequency spectrum. $`S^{(2)}(f)`$ is a measure of higher order cumulants. A departure of $`S^{(2)}(f)`$ from a white spectrum will be an indication of non-gaussian spectra and finite higher order cumulants. Figure 5a shows the second spectrum of the samples at T=4.2K. The spectra ($`f<1Hz`$) fits the relation $`S^{(2)}(f)f^p`$. As a function of $`n/n_c`$, $`p`$ changes (see inset, figure 5a) from a small value of $`0.5`$ for $`n/n_c>1`$ (for white spectrum, $`p=0`$)to $`1.5`$ as the insulating state is reached, changing rapidly at $`n=n_c`$. Figure 5b shows how $`p`$ varies from a small value $`0.5`$ to a large value$`1.3`$ in the insulating sample K240 as the temperature is lowered. These observations clearly go to show that as the MIT is approached the dynamics that gives rise to the conductance fluctuation becomes increasingly non-gaussian and this is more so as the temperature is lowered.

The divergence of the conductance fluctuation at the MIT has been predicted in some models. Since the conductance a variable with large random fluctuations, it has an important impact on the complete PDF. Some models  propose that in the large conductance limit (metallic side), the relative fluctuations are smaller and hence the higher cumulants of the fluctuations are negligible. This leads to a roughly gaussian PDF. With decrease in conductance, the relative fluctuations are larger and the higher cumulants start contributing. On the insulating side, this shows up as a log-normal PDF. Since the conductance decreases both as $`n/n_c`$ decreases through 1 as well as T decreases on the insulating side, our results are in excellent agreement with this model. However, there is one more important aspect in the non-gaussian fluctuation. The departure of the second spectra from a frequency independent white spectrum at low temperatures in samples with $`n/n_c1`$ clearly indicates that the non-gaussian fluctuation has a correlated dynamics associated with it. We have no clear interpretation as to what gives rise to it but we suggest that this may be an indication of a glass like freezing similar to what has been seen in spin glasses . Recent theories  have suggested that in the MIT region there is a glass like freezing of the charge system. Analysis of conductivity data close to the MIT using various scaling forms has shown that a comprehensive scaling description could be obtained if a glass-like transition is invoked at the MIT. Recent noise experiments on 2-D electron system in Si inversion layers  has shown that close to the MIT there is an appearance of large low-frequency dynamics which makes the noise blow up near the transition. This has been interpreted as a signature of freezing into electron glass. The divergence of noise at the MIT is unlikely to arise from the usual dynamics of atomic defect relaxations  as the strongest coupl ing mechanism that can produce conductivity fluctuation from atomic dynamics is UCF which itself saturates at a variance of $`e^2/h`$ (Mechanisms like UCF have limits to the fluctuations they can produce, i.e. when UCF is saturated). The large $`\gamma _H`$ seen in the insulating samples need much stronger sources of fluctuations. As the system goes through the MIT the dopant concentration is reduced and there is no reason why the defect dynamics should become more correlated at the MIT (for a $`n`$ $`3\times 10^{18}cm^3`$ the average defect separation is $`10nm`$). In view of this it is most likely that the fluctuation that we see close to the MIT at low temperatures is electronic in origin and is a manifestation of the slowing down of the time scale of charge fluctuations as the MIT is approached. On approach to the MIT the single particle density of states can show both temporal and spatial fluctuations and if the dynamics is glassy then the relaxation time becomes an exponential function of the correlation length . We propose that the large rise in the low frequency noise at the approach to transition is a manifestation of low frequency dynamics of the fluctuation in local density of states which are strongly correlated by coulomb interaction in the absence of strong screening.

The authors wish to thank CSIR, Govt. of India for a fellowship (S.K) and a sponsored project (A.K.R).

Present address: Cavendish Laboratory, University of Cambridge, Madingley Road, Cambridge, UK.

Figure Captions

Figure 1: Conductivity($`\sigma `$) vs. temperature of the samples Si:P(B). $`\sigma `$ of the 3 most metallic samples PS24,PS41 and D150 have been scaled down by a factor of 5 for convenience. The table shows the relevent sample properties.

Figure 2: $`\frac{\mathrm{\Delta }G^2}{G^2}`$ vs. temperature in the insulating sample K139 measured at different biasing voltages. When measured with the lowest bias, the noise scales as $`T^\beta `$ for all the samples (for K139, $`\beta =6`$). The inset shows how $`\mathrm{\Delta }G^2`$ deviates from a $`G^2`$ dependence on application of larger measuring biases.

Figure 3: The scaled noise $`\gamma _H`$ for different values of $`n/n_c`$. Inset (a) shows the same graph, magnified around $`n/n_c=1`$. Inset (b) shows the $`n/n_c`$ dependence of $`\beta `$.

Figure 4: 4(a),(b) and (c) shows the PDFs at 4.2K for K348, K242 and K240. 4(d) shows the relative contribution of gaussian and lognormal parts (see text) of the PDF as a function of n. 4(e) shows the same as a function of T in sample K240.

Figure 5: 5(a) Second spectrum $`S^{(2)}(f)`$ at 4.2K (see text). The graphs have been shifted vertically for visual clarity. The inset shows the sharp growth of frequency dependence at the critical concentration $`n=n_c`$ at 4.2K. 5(b) T dependence of the exponent $`p`$ of the second spectra in K 240.


format: text

text:
# Tagging the 𝒑𝒏→𝒅ϕ reaction by backward protons in 𝒑𝒅→𝒅ϕ𝒑_{𝒔𝒑} processes

## I Introduction

Data of elementary nucleon reactions involving a neutron in the entrance channel are scarce. One can use, for instance, a tagged neutron ($`n`$) beam  impinging on a proton ($`p`$) target or a proton beam impinging on a nuclear target. In the latter case one needs models to deconvolute the $`pn`$ reaction. Here we are interested in threshold-near meson production. Over the last decade the data basis for the $`pp`$ reactions has been impressively improved . Such a detailed knowledge of elementary hadron reactions is necessary for analyzing heavy-ion collisions. Threshold-near meson production is also interesting for testing predictions of QCD in the non-perturbative region, provided e.g., by chiral perturbation theory. While it seems to be straightforward to translate cross sections for reactions of the type $`ppppM`$ ($`M`$ is a meson) into the cross section $`nnnnM`$, the channel $`pnpnM`$ is less simple accessible. Model calculations point to a non-trivial energy dependence of the ratio $`\sigma (ppppM)/\sigma (npnpM)`$, e.g. for $`\varphi `$ meson production . Experiments showed, furthermore, a surprisingly large ratio $`\sigma (pdpd\eta )/\sigma (pppp\eta )`$ . Therefore, it is desirable to test directly predictions for the $`pn`$ reactions.

The spectator technique  offers one possibility to isolate to a large extent the quasi-free reaction at a neutron when a proton beam impinges on a deuteron target. Since the final state interactions in the $`pn`$ system differs for singlet and triplet states and can be large at threshold, it is of great advantage to consider the particular reaction $`pndM`$, as the deuteron $`d`$ is in one well defined state. Guided by this, in  the reactions $`pndM`$ with $`M=\eta ,\eta ^{},\omega ,\varphi ,a_0^+`$ have been considered.

The aim of the present note is to analyze the complete reaction $`pddMp_{sp}`$ and to test to what extent the spectator technique  really selects a quasi-free reaction at the neutron in the deuteron target. In other words, we are going to derive a factorization theorem showing kinematical conditions for which the total cross section $`\sigma (pddMp_{sp})`$ factorizes in a part depending on the target deuteron characteristics and a part $`\sigma (pndM)`$. In doing so we select as a sufficiently transparent example the $`\varphi `$ production. The reasons are obvious: (i) the elementary $`\varphi `$ production amplitude is simple , (ii) the $`\varphi `$ meson is interesting with respect to OZI rule  and hidden strangeness in the nucleon , and (iii) the $`\varphi `$ meson is important with respect to both the inclusive $`K^{}`$ production  and the $`e^+e^{}`$ decay channel to be studied at HADES in near future . An extension to $`\rho `$ and $`\omega `$ production goes along the same lines, but due to the role of $`s`$, $`u`$ channel baryon resonances  the elementary amplitude of the subprocess $`pnpnM`$ is by far more complicate and needs separate consideration, as other mesons too.

Our paper is organized as follows. In section II we present the kinematics, cross section and amplitudes of the complete reaction $`pdd\varphi p_{sp}`$. Numerical results are exhibited in section III, where polarization observables are discussed and compared with the ones of the quasi-free sub-reaction $`pnd\varphi `$. We also derive selection rules. The factorization of the cross section is considered in section IV. The summary can be found in section V.

## II kinematics, cross section, amplitudes

Let us consider the process of $`\varphi `$ meson production in the exclusive reaction envisaged in

$`p_p+p_d=p_d^{}+p_\varphi +p_p^{},`$ (1)
at kinetic energies of the incoming proton corresponding to those achievable at the Cooler Synchotron COSY $`T_p2.7`$ GeV. The deuteron is here supposed to be detected in the forward direction, i.e., in direction of the incoming proton with relativistic energy, $`T_d^{}1`$ GeV, and, contrarily, the final (spectator) proton is slowly moving in the backward direction ($`\theta _{sp}>90^0`$, $`p_p^{}<100`$ MeV). The main contribution to the process (1) comes then from the spectator mechanism, where the incoming proton interacts with the internal neutron from the deuteron and produces a meson and a deuteron in the final state. The second proton in the deuteron acts merely as a spectator. Schematically the spectator mechanism is depicted in Fig. 1.

In the laboratory system, where the target deuteron is at rest, the relevant momenta are defined by $`p_p=(E_p,𝒑_p)`$ and $`p_d=(M_d,\mathrm{𝟎})`$ as initial momenta of the proton and deuteron, $`p_p^{}=(E_p^{},𝒑_p^{})`$, $`p_d^{}=(E_d^{},𝒑_d^{})`$, and $`p_\varphi =(E_\varphi ,𝒑_\varphi )`$ as momenta of final proton, deuteron and $`\varphi `$ meson, respectively. The fivefold differential cross section reads

$`2E_d^{}{\displaystyle \frac{d^5\sigma }{d^3p_d^{}d\mathrm{\Omega }_p^{}}}={\displaystyle \frac{1}{8(2\pi )^5\sqrt{\lambda }}}\left[{\displaystyle \frac{|𝒑_p^{}|^2}{|B|𝒑_p^{}|CE_p^{}|}}\right]{\displaystyle \frac{1}{6}}{\displaystyle \underset{\mathrm{𝑠𝑝𝑖𝑛𝑠}}{}}||^2,`$ (2)
where the Källen function $`\lambda =4M_d^2|𝒑_p|^2`$ and the kinematical coefficients $`B=M_d+E_pE_d^{}`$ and $`C=|𝒑_p|\text{cos}\theta _p^{}|𝒑_d^{}|\text{cos}\theta _{dp}^{}`$ reflect the energy conservation. The invariant amplitude $``$ depends, besides the spin projections $`s_p,_d;s_p^{},_d^{},_\varphi `$ on the quantization z-axis (oriented along $`𝐩_p`$), also upon the initial energy and five independent kinematical variables, say $`𝐩_d^{}`$ and $`\mathrm{\Omega }_p^{}`$.

To compute the amplitude $``$, let us define the scattering operator $`\widehat{𝒪}_{_\varphi ,_d^{}}`$ (a $`44`$ matrix in the spinor space) as the truncated amplitude for the subprocess $`pnd\varphi `$ (hatched blob of the diagram in Fig. 1), i.e., as an operator which, when sandwiched between two free spinors $`u_p`$ and $`u_n`$, results in the amplitude of the real process $`pnd\varphi `$ with two on-mass shell nucleons. In the general case this operator is assumed to describe the process of meson production from two off-mass shell nucleons. With this definition, the amplitude $``$ may be written as

$`=\overline{u}_\alpha (p_p^{},s_p^{})(\widehat{p}_p^{}m)_{\alpha \beta }\mathrm{\Phi }__d^{c\beta }(p_n,p_p^{})\widehat{𝒪}_{_\varphi ,_d^{}}^{\delta c}u_\delta (p_p,s_p),`$ (3)
where $`m`$ is the nucleon mass and $`\mathrm{\Phi }__d^{c\beta }(p_n,p_p^{})`$ denotes the Bethe-Salpeter (BS) amplitude of the deuteron in the momentum space. In Eq. (3), the spinor indices are exhibited explicitly, with Latin and Greek characters referring to neutrons and protons, respectively (summation over indices occurring pairwise is supposed). An expression for the amplitude (3) can be found by specifying a model for the operator $`\widehat{𝒪}`$. However, in order to keep the correspondence between the reaction (1) and the process of $`\varphi `$ meson production in $`NN`$ collisions as close as possible, we highlight explicitly its dependence upon the BS amplitude of the final deuteron by introducing another operator $`\widehat{O}`$ as

$`\widehat{𝒪}_{_\varphi ,_d^{}}^{\delta c}=i{\displaystyle \frac{d^4p}{(2\pi )^4}\overline{\mathrm{\Phi }}__d^{}^{\alpha b}(1^{},2^{})\widehat{O}_{\alpha \delta }^{bc}(12;1^{}2^{},_\varphi )}.`$ (4)
Then from Eq. (4) and from Fig. 1 it is easy to recognize the operator $`\widehat{O}(12;1^{}2^{},_\varphi )`$ as the one describing the $`\varphi `$ meson production in the sub-reaction $`pnd\varphi `$ with an off-mass shell initial neutron. Such an operator has been introduced in Ref.  where its properties have been discussed in detail. We briefly recall here the main steps in computing $`\widehat{O}(12;1^{}2^{},_\varphi )`$ and the amplitude $``$ in Eq. (3). Usually, within the BS formalism in order to exhibit formulae in a familiar matrix form, one introduces new, “charge-conjugated” quantities such as $`\mathrm{\Psi }\mathrm{\Phi }\gamma _c`$, $`\widehat{𝒬}\gamma _c\widehat{𝒪}`$, $`\gamma _c\gamma _\mu ^T=\gamma _\mu \gamma _c`$, $`v(p)=[\overline{u}(p)\gamma _c]^T`$, $`\overline{v}(p)=[\gamma _cu(p)]^T`$. Then the amplitude $``$ can be written as

$`=\overline{v}(p_p,s_p)\widehat{𝒬}_{_\varphi ,_d^{}}\mathrm{\Psi }__d(p_n,p_p^{})(\widehat{p}_p^{}+m)v(p_p^{},s_p^{}),`$ (5)
where in Eq. (5) and in subsequent expressions one has usual matrix multiplications and traces. For further simplification we insert between $`\widehat{𝒬}`$ and $`\mathrm{\Psi }`$ in (5) the complete set of the Dirac free spinors $`\widehat{I}=\left(1/2m\right)\underset{r=1}{\overset{4}{}}\epsilon _ru_r(p)\overline{u}_r(p)`$, where $`\epsilon _r=+1`$ for $`r=1,2`$ and $`\epsilon _r=1`$ for $`r=3,4`$, to obtain

$`={\displaystyle \frac{1}{2m}}{\displaystyle \underset{s_n}{}}T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)\overline{u}(p_n,s_n)\mathrm{\Psi }__d(p_n,p_p^{})(\widehat{p}_p^{}+m)v(p_p^{},s_p^{}).`$ (6)
The amplitude $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)=\overline{v}(p_p,s_p)\widehat{𝒬}_{_\varphi ,_d^{}}u(p_n,s_n)`$ coincides with the matrix element of the real process $`pnd\varphi `$, defined by the Mandelstam variables $`\nu _{pn}=(p_p+p_n)^2`$, $`t=(p_pp_d^{})^2`$. Here is worth emphasizing that the kinematics of the process (1) differs from the one of the real process $`pnd\varphi `$. In spite of commonly defined variable $`t`$, the numerical values of momenta of the final meson and deuteron do not coincide in the two processes. Rather, the amplitude $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ in the process (1) corresponds to an amplitude of meson production from an off-mass shell neutron. If the off-mass shell effects were negligibly small, then $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ in both processes would depend solely upon the momentum transfer $`t`$ and the effective energy $`\nu _{NN}=(E_p+E_p^{})^2(𝐩𝐩^{})^2`$. As shown in Ref. , $`T_{s_p,s_n}^{_\varphi ,_d^{}}`$ can be written in the form

$`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)={\displaystyle \frac{i}{(2m)^2}}{\displaystyle \underset{r,r^{}=1}{\overset{4}{}}}{\displaystyle \frac{d^4p}{(2\pi )^4}A_{s_p,s_n;r,r^{}}^_\varphi (12;1^{}2^{})\overline{v}_r^{}(2^{})\overline{\mathrm{\Psi }}__d^{}(1^{},2^{})u_r(1^{})},`$ (7)
where the partial spin amplitudes $`A_{s_p,s_n;r,r^{}}^_\varphi (12;1^{}2^{})`$ can be calculated within an effective one-boson exchange model by an explicit evaluation of the elementary diagrams depicted in Fig. 2 in . We use the parameter set “B” in Ref.  to calculate the amplitudes $`A_{s_p,s_n;r,r^{}}^_\varphi (12;1^{}2^{})`$. Now the second part of the matrix element in Eq. (6), which includes the BS amplitude $`\mathrm{\Psi }__d`$ for the deuteron at rest, can be straightforwardly evaluated by using the explicit expressions for the $`{}_{}{}^{3}S_{1}^{++}`$ and $`{}_{}{}^{3}D_{1}^{++}`$ components of the deuteron BS amplitude . The result is

$`=\sqrt{8\pi M_d}{\displaystyle \underset{s_n}{}}T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)\left[u_S(|𝐩^{}|)C_{\frac{1}{2}s_p^{}\frac{1}{2}s_n}^{1_d}u_D(|𝐩^{}|)\sqrt{4\pi }{\displaystyle \underset{m,\mu }{}}C_{\frac{1}{2}s_p^{}\frac{1}{2}s_n}^{1\mu }C_{2m1\mu }^{1_d}Y_{2m}(𝐧)\right],`$ (8)
where $`M_d`$ is the deuteron mass, $`u_{S,D}`$ denote $`S`$ and $`D`$ waves the BS wave function , $`𝐧=𝐩/|𝐩|`$ and $`Y_{2m}(𝐧)`$ the spherical harmonics of second order. In deriving Eq. (8) we used the identity $`\sqrt{2}\underset{\nu _1\nu _2}{}C_{\frac{1}{2}\nu _1\frac{1}{2}\nu _2}^{1_d}\chi _{\frac{1}{2}\nu _1}\chi _{\frac{1}{2}\nu _2}^T=i(𝝈𝝃__d)\sigma _y`$, where $`C_{\frac{1}{2}\nu _1\frac{1}{2}\nu _2}^{1_d}`$ are the Clebsch-Gordan coefficients, $`\chi _{\frac{1}{2}\nu }`$ and $`𝝃`$ stand for the spin-$`\frac{1}{2}`$ functions and the deuteron polarization vector, respectively. Equation (8) shows that an exact factorization of the cross section into two parts (one depending upon the target deuteron wave function, another one upon the kinematics and dynamics of the subprocess $`pnd\varphi `$) generally does not occur. This is because the summation over the spins of the internal neutron $`s_n`$ which leads to interference terms in the square of the matrix element. However, below we show that under special kinematical conditions an approximate factorization is possible.

## III Numerical results

From Eqs. (6 \- 8), the differential cross section (2) and any other polarization observables of the process $`pdd\varphi p_{sp}`$ are calculable. The BS wave functions $`u_{S,D}`$ basically coincide with the non-relativistic $`S`$ and $`D`$ waves of the deuteron for small values of the momentum of the spectator proton (see ). Therefore, the relativistic and non-relativistic expressions for the amplitude $``$ have the same form, and one may perform also calculations by using deuteron wave functions obtained within a non-relativistic approach. In our calculations we use the BS solution , obtained with a realistic one-boson exchange kernel. Additionally, Paris and Bonn wave functions are considered to investigate the sensitivity of observables on the deuteron model.

All our calculations have been performed within kinematical conditions corresponding to the COSY Proposal #75 , i.e., the detected deuteron is supposed to be in the forward direction and the spectator proton in the backward hemisphere with $`\theta _p^{}120^0140^0`$. The momenta of the deuteron are chosen in the range $`p_d^{}1.9`$ \- 3 GeV/c, which means that the spectator proton is detected with low velocity, i.e., $`p_p^{}0`$ \- 0.04 GeV/c. These kinematic conditions imply that both the reaction (1) and the subprocess $`pnd\varphi `$ occur essentially near the threshold.

In Fig. 2, the differential cross section is presented as a function of the deuteron momentum for the specific choice of variables $`\theta _d^{}=0^0`$, and $`\theta _p^{}=130^o`$. (We changed the kinematics with $`\theta _d^{}`$ up to $`5^o`$ and find a very smooth behavior of all observable.) The kinetic energy of the initial proton is $`T_p=2.69`$ GeV. The maximum of the cross section at $`p_d^{}2.35`$ GeV/c corresponds to that kinematical condition where the momentum of the spectator proton is maximal, $`p_p^{}0.035`$ GeV/c. Together with the large spectator angles this implies that at these kinematical conditions the excess energy in the $`pnd\varphi `$ subprocess is minimal ($`\mathrm{\Delta }s^{1/2}4`$ MeV) and very near the threshold. Obviously, since at such kinematics the production amplitude is almost constant, the maximum in the differential cross section has a pure kinematical origin being entirely governed by the kinematical factor within the square brackets in Eq. (2). Fig. 2 highlights that different models for the deuteron wave function result in quantitatively different cross sections. This is because the integration over the internal momentum in the final deuteron covers the region of the minimum of the $`S`$ wave where different models provide essentially different wave functions (see discussion in ).

Let us now discuss polarization observables. Since the spin structure of the amplitude $``$ in Eq. (8) is basically determined by the spin structure of the amplitude $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ of the subprocess $`pnd\varphi `$, we proceed with an analysis of this amplitude. Bearing in mind that $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ is a Lorentz invariant quantity one can, for the sake of simplicity, investigate its structure in the center of mass of the detected deuteron (near the threshold this system coincides with the center of mass system of two nucleons), provided the Lorentz boost is performed along the $`z`$-axis. In such a case the spin structure of $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ is the same as in the laboratory system (no additional Wigner rotations are needed). Near the threshold, the final $`\varphi d`$ system is in a $`S`$ orbital state corresponding to $`I_f=0`$, $`J_D^\pi =1^+`$ ($`L_D=0,2`$), $`J_f^\pi =J_\varphi ^\pi +J_D^\pi =0^{},1^{},2^{}`$, where $`I`$, $`L`$ and $`J^\pi `$ are the total isospin, total radial angular momentum, angular momentum and parity, respectively. Thus from symmetry constraints the allowed configurations in the initial state are $`I_i=0`$, $`L_i=1,3,5\mathrm{}`$ and total spin $`S_i=0`$. The conservation law $`J_i=J_f`$ implies $`L_i=1`$, so that $`J_i^\pi =1^{}`$. Hence, near the threshold the main contribution into the amplitude $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ comes from initial spin configurations with $`s_p=s_n`$ and, since the two nucleons in the final state form a deuteron (with the total spin $`S=1`$), the most probable transitions near the threshold are expected to be those with nucleon spin-flips in the corresponding vertices and with the deuteron projections $`_d^{}=\pm 1`$. Consequently, the projections of the $`\varphi `$ meson are to be $`_\varphi =1`$.

Coming back to the complete reaction one can expect that the prevailing contribution into the full amplitude $``$ comes from configurations with $`_ds_p^{}=s_p`$, if (i) the effects of $`D`$ waves in the target deuteron may be disregarded near the threshold, (ii) the spectator mechanism describes the reaction appropriately, and (iii) the factorization (see below) holds. Hence, near the threshold the polarization observables of the reaction (1) are predicted to behave in the same manner as the ones of the subprocess $`pnd\varphi `$. Therefore, the study of various polarization observables allows to check the factorization.

Let us define the asymmetry of the reaction as

$`𝒜={\displaystyle \frac{d\sigma (s_p+s_n=1)+d\sigma (s_p+s_n=1)d\sigma (s_p+s_n=0)}{d\sigma (s_p+s_n=1)+d\sigma (s_p+s_n=1)+d\sigma (s_p+s_n=0)}},`$ (9)
where the spin of the intrinsic neutron is defined as $`s_n=_ds_p^{}`$. Then, from the above selection rules near the threshold, the asymmetry $`A`$ is predicted to be negative and close to -1. In Fig. 3, results of a complete calculation of the asymmetry $`A`$ are presented. The kinematical conditions are as in Fig. 2. It is seen that in the whole kinematical range of $`p_d^{}`$ the asymmetry is basically the same as in the subprocess $`pnd\varphi `$. Moreover, as expected, at $`p_d^{}`$ near the maximum the asymmetry approaches -1.

Similar results are obtained for other polarization observables. For instance, the tensor analyzing power defined as

$`T_{20}={\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle \frac{d\sigma (_d^{}=1)+d\sigma (_d^{}=1)2d\sigma (_d^{}=0)}{d\sigma (_d^{}=1)+d\sigma (_d^{}=1)+d\sigma (_d^{}=0)}},`$ (10)
is predicted to be $`1/\sqrt{2}`$ from the selection rules. Indeed, the full evaluation of Eqs. (6 -8) yields a value close to that (see Fig. 4).

Apart from the polarization observables which near the threshold are determined directly by the corresponding quantities for the subprocess $`pnd\varphi `$, one can define another set of observables with completely different behavior in the two reactions, i.e., in the subprocess $`pnd\varphi `$ and in the reaction (1). For this purpose let us define observables at fixed initial deuteron and proton spin projections, e.g., at $`_d=0`$ and $`s_p=\frac{1}{2}`$, respectively. Near the threshold this means that the spin projection of the neutron will be $`s_n=s_p^{}`$ and by measuring the spectator polarization one can form different combinations which determine (allowed or forbidden) transitions for the subprocess $`pnd\varphi `$. Together with conservation rules for the complete reaction this admits a study of different observables which vanish in the sub-reaction $`pnd\varphi `$. Correspondingly, let us define the quantity

$`W_{zz}={\displaystyle \frac{d\sigma (_\varphi =1)d\sigma (_\varphi =1)}{d\sigma (_\varphi =1)+d\sigma (_\varphi =1)}}|_{s_p=\frac{1}{2};_d=0}.`$ (11)
If one measures this quantity for the positive spin projection of the spectator, then in the subprocess $`pnd\varphi `$ one has super allowed transitions ($`s_p+s_n=0,_\varphi =\pm 1,_d^{}=1`$) and, since in (11) a summation over the final deuteron projections is assumed, $`W_{zz}`$ is predicted to approach zero at the threshold (see ).

A different situation occurs if one considers $`W_{zz}`$ at $`s_p^{}=1/2`$. In this case, at least for $`t=t_{min(max)}`$ (i.e., the forward/backward direction for the subprocess $`pnd\varphi `$), the conservation of the spin projection in the complete reaction (1) strongly suppresses the final states with $`_\varphi =_d^{}=\pm 1`$ and $`_\varphi =1,_d^{}=0`$, so that the only contributing final states remain these with $`_\varphi =1,_d^{}=0,1`$ and $`_\varphi =1,_d^{}=1`$, which imply $`W_{zz}1/3`$, in contrast to the vanishing $`W_{zz}`$ in $`pnd\varphi `$ . Figure 5 illustrates our full numerical evaluation for $`W_{zz}`$ at $`s_p^{}=1/2`$.

## IV Factorization

Let us discuss in more detail the behavior of the cross section Eq. (2) with the matrix element (8) near the threshold. As mentioned above, the summation over the spin projection of the intrinsic neutron $`s_n`$ in Eq. (8) prevents an exact factorization of the cross section into the cross section of the subprocess $`pnd\varphi `$ and a “target factor”. However, from Figs. 3 \- 5 we conclude that, since the deviations of the exactly calculated observables from those determined by the selection rules of the subprocess $`pnd\varphi `$ are small, a factorization should hold indeed with a good accuracy. For instance, the asymmetry in Fig. 3 demonstrates the validity of the selection rule $`s_p=s_n`$, i.e., the amplitude $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ is proportional to $`\delta _{s_p,s_n}`$. Similarly, Fig. 5 shows that the main contribution comes from transitions with $`|_d^{}|=1`$, and from the above selection rules the amplitude $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)`$ is also proportional to $`\delta _{_d^{},_\varphi }`$. Hence, near the threshold one may represent the amplitude as $`T_{s_p,s_n}^{_\varphi ,_d^{}}(\nu _{pn},t)=\delta _{_d^{},_\varphi }\delta _{s_p,s_n}\stackrel{~}{T}_{s_p,_d^{}}(\nu _{pn},t),`$ where $`\stackrel{~}{T}_{s_p,_d^{}}(\nu _{pn},t)`$ is the threshold-near amplitude of the process $`pnd\varphi `$. Then in Eq. (8) the summation over $`s_n`$ disappears and, as conjectured above, the unpolarized cross section Eq. (2) factorizes. Substituting this amplitude into Eq. (2), introducing the corresponding flux factors for the sub-reaction $`pnd\varphi `$ and performing the summation over all spin projections, the unpolarized cross section can be cast in the form

$`2E_d^{}{\displaystyle \frac{d^5\sigma }{d^3p_d^{}d\mathrm{\Omega }_p^{}}}={\displaystyle \frac{2}{(2\pi )^3}}{\displaystyle \frac{\nu _{pn}(\nu _{pn}4m^2)|𝒑_p^{}|^2}{|𝒑||B|𝒑_p^{}|CE_p^{}|}}\left(u_S^2(|𝐩_p^{}|)+u_D^2(|𝐩_p^{}|)\right){\displaystyle \frac{d\sigma (pnd\varphi )}{dt}},`$ (12)
where the differential cross section $`d\sigma (pnd\varphi )/dt`$ depends only upon the transferred momentum $`t`$ (determined by $`p_p`$ and $`p_d^{}`$) and the invariant mass $`\nu _{pn}`$ (determined by $`p_p`$ and $`p_p^{}`$ due to the balance $`E_n=M_dE_p^{}`$, $`𝐩_𝐧=𝐩_𝐩^{}`$). From Eq. (12) it is seen that by studying experimentally the amplitude (8), near the threshold, one can obtain directly information about the subprocess $`pnd\varphi `$. Here it is worth reminding that the factorization of the cross section holds only if the contribution of the $`D`$ wave is negligible, i.e., if in Eq. (8) one can put $`s_n=s_p`$ and arrives at

$`\sqrt{8\pi M_d}T_{s_p,s_p}^{_\varphi ,_d^{}}(\nu _{pn},t)\left[u_S(|𝐩_p^{}|)C_{\frac{1}{2}s_p^{}\frac{1}{2}s_p}^{1_d}u_D(|𝐩_p^{}|)\sqrt{4\pi }{\displaystyle \underset{m\mu }{}}C_{\frac{1}{2}s_p^{}\frac{1}{2}s_p}^{1\mu }C_{2m1\mu }^{1_d}Y_{2m}(𝐧)\right].`$ (13)
Eqs. (8, 12) and (13) show that the smaller the values of the spectator momentum are the higher the accuracy of the factorization is. This has been confirmed numerically by comparing exact results based on Eqs. (2, 8) and approximate results based on Eqs. (12, 13). In the whole kinematical region of interests the difference is less than 0.5%.

One may generally conjecture that a structure analog to Eq. (8) gives rise to factorization under suitable kinematical conditions, hence providing a theoretical foundation of quasi-free reactions at the neutron in the deuteron.

## V Summary

In summary, within a covariant Bethe-Salpeter approach we study the relation of the tagged quasi-free reaction $`pnd\varphi `$ and the complete process $`pdd\varphi p_{sp}`$ near the threshold. Guided by the very proximity of the numerical results for cross section and polarization observables based on the exact formulae and on the factorized ones, we find that, for special kinematical conditions, selection rules are operative immediately causing a factorization $`\sigma (pdd\varphi p_{sp})=F\sigma (pnd\varphi )`$ where $`F`$ depends on deuteron and kinematical quantities. This is equivalent to an omission of the contribution of the deuteron $`D`$ wave at small values of the spectator momentum.

As an extension of the present work, an investigation of the twice tagged reaction $`ddnnMp_{sp1}p_{sp2}`$ can be envisaged to study the quasi-free meson production by the $`nn`$ reaction proposed in .

###### Acknowledgements.

We are grateful to A.I. Titov for many valuable discussions. L.P.K. and S.S.S. acknowledge the warm hospitality of the nuclear theory group in the Research Center Rossendorf. The work has been supported by Heisenberg - Landau program, BMBF grant 06DR921, and RFBR 00-15-96737.


text:
# Shock breakouts in SNe Ib/c

## Shock breakouts in SNe Ib/c

S.I. Blinnikov1, D.K. Nadyozhin1, S.E.Woosley2 and E.I.Sorokina3
1 Institute for Theoretical and Experimental Physics, Moscow
2 UCO/Lick Observatory, University of California, Santa Cruz
3 Sternberg Astronomical Institute, Moscow

It was realized long ago that a shock propagating along the profile of decreasing density should accelerate . In compact presupernovae, like SNe Ib/c, the shock can become relativistic  and is able to produce a burst of X-ray and even $`\gamma `$-ray radiation . The case of a peculiar type Ic SN1998bw, probably related to GRB980425, has shown that SNe Ib/c may have the highest explosion energy and highest production of 56Ni among core-collapsing supernovae . Supernovae of type Ib/c are also interesting for theory due to problems with their light curve modeling. Their understanding may serve as a diagnostics of mass loss from massive stars at the latest phases of evolution.

Numerical modeling of shock breakout in SNe Ib/c was done previously using some simplifying approximations . Our method, realized in the code stella, allows us to get more reliable predictions for the outburst. Improvements in the theory done in the current work:

* multi-energy-group time-dependent radiation transfer,
* taking into account of (some) relativistic effects.

A representative presupernova in our runs was a WR star built by the code kepler  (model 7A). Late light curves and spectra for this model were studied in . We present UBV light curves found for the same model by stella in Fig. 1.

The Model 7A has mass 3.199 $`M_{}`$(including the mass of the collapsed core). Its radius prior to explosion is not strictly fixed because the outer mesh zones in kepler output actually model the strong stellar wind and are not in hydrostatic equilibrium. So, we fixed the radius by hand and we got four models with radii from 0.76 up to 2 $`R_{}`$ which are in hydrostatic equilibrium. Our results are summarized in the Table.

Explosions with kepler  gave maximum temperature of photons $`T5\times 10^5`$ K at shock breakout. We have much finer mesh zoning at the edge of the star (down to 1$`\times 10^{12}`$ $`M_{}`$) and better physics and we find higher values of $`T`$. The Fig. 2 shows the difference between effective and color temperatures (labels ‘e’ and ‘c’, respectively, in the Table). One should note that the peak values of luminosity and temperature given in the table and plots in Fig. 2 do not contain the light travel time correction.

Models, labeled as 3.5N were computed by D.K.Nadyozhin in 1992 with the equilibrium radiation diffusion hydrocode snv. Left plot in Fig. 3 shows one of his results. Dashed lines in Fig. 3 demonstrate the effect of averaging the light curve due the light travel time correction.

| Parameters of shock breakouts | | | | | |
| --- | --- | --- | --- | --- | --- |
| $`M`$ | $`R_0`$ | $`E_{\mathrm{kin}}`$ | $`L_\mathrm{p}`$ | $`T_\mathrm{p}`$ | $`\mathrm{\Delta }t`$ |
| $`M_{}`$ | $`R_{}`$ | foe | erg/s | $`10^6`$K | s |
| 3.2 | 0.76 | 1.24 | 4.2$`\times 10^{44}`$ | 4.2c | 0.021 |
| 3.2 | 1.00 | 1.32 | 5.8$`\times 10^{44}`$ | 4.3c | 0.026 |
| 3.2 | 1.23 | 1.30 | 6.8$`\times 10^{44}`$ | 4.3c | 0.043 |
| 3.2 | 2. | 1.39 | 9.4$`\times 10^{44}`$ | 4.3c | 0.12 |
| 3.2 | 2. | 4.36 | 3.6$`\times 10^{45}`$ | 5.3c | 0.028 |
| 3.2 | 2. | 8.86 | 4.8$`\times 10^{45}`$ | 7.2c | 0.020 |
| 3.5N | 0.76 | 1.30 | 1.4$`\times 10^{45}`$ | 5.1e | 0.028 |
| 3.5N | 1.23 | 1.30 | 8.1$`\times 10^{44}`$ | 3.5e | 0.067 |

$`M`$ and $`R_0`$ — presupernova mass and radius, respectively, in solar units
$`E_{\mathrm{kin}}`$ — the kinetic energy at infinity in $`10^{51}`$ ergs
$`L_p`$ and $`T_p`$ — the peak luminosity and temperature
$`\mathrm{\Delta }t`$ is the width of the light curves at 1 stellar magnitude below $`L_p`$

We conclude that photon spectra at shock breakouts in SNe Ib/c peak near $`3kT1`$ keV and the X-ray fluences predicted can be observed by future missions up to tens of Mpc.

We are grateful to MPA staff for permanent and generous support. Our work in US was supported by grants NSF AST-97 31569 and NASA - NAG5-8128, in Russia RBRF 00-02-17230 and RBRF 02-02-16500.


text:
# Observation of the spin-charge thermal isolation of ferromagnetic 𝐆𝐚_0.94𝐌𝐧_0.06𝐀𝐬 by time-resolved magneto-optical measurement

## Abstract

The dynamics of magnetization under femtosecond optical excitation is studied in a ferromagnetic semiconductor $`\mathrm{Ga}_{0.94}\mathrm{Mn}_{0.06}\mathrm{As}`$ with a time-resolved magneto-optical Kerr effect measurement with two color probe beams. The transient reflectivity change indicates the rapid rise of the carrier temperature and relaxation to a quasi-thermal equilibrium within a few ps, while a very slow rise of the spin temperature of the order of 500ps is observed. This anomalous behavior originates from the thermal isolation between the charge and spin systems due to the spin polarization of carriers (holes) contributing to ferromagnetism. This constitutes experimental proof of the half-metallic nature of ferromagnetic $`\mathrm{Ga}_{0.94}\mathrm{Mn}_{0.06}\mathrm{As}`$ arising from the large exchange energy of carriers compared to Fermi energy.

Carrier-induced ferromagnetism is one of the most visited spontaneous symmetry-breaking phenomena in the solid-state physics. The underlying physics of this phenomenon deeply involves fundamental aspects of quantum physics including interactions between localized and delocalized carriers, where many-body quantum correlation plays an essential role.

The discovery of diluted magnetic semiconductors (DMS’s) with high Curie temperature that have become available due to recent advances in molecular beam epitaxy and transition metal ion doping techniques opens new horizons in exploring of the carrier mediated ferromagnetismMunekata ; HOhno1 ; HOhno2 . In particular, the compatibility of DMS’s with conventional semiconductor technology is very promising to employ the spin degree of freedom in practical devices. Correspondingly, along with extensive efforts to clarify the microscopic mechanism of the ferromagnetism, much work has been devoted to the search for systems with higher $`T_c`$Matsukura ; Akai ; Dietl ; Inoue ; HOhno3 .

In the initial stage of the research for the origin of the ferromagnetism, two mechanisms of the interaction between the local moments were considered. The first was based on the Ruderman-Kittel-Kasuya-Yoshida (RKKY) mechanism, in which the ferromagnetic interaction of localized spins is mediated by nearly free carriers Matsukura . This model reproduces the observed concentration dependence of Curie temperatures in II-VI DMSHaury .

However, recent photoemission experiments in Mn doped GaAs have revealed a rather low density of state, which indicates localization of carriersOkabayashi  and thus we need to consider the mechanism beyond simple RKKY model. The second one was based on the double exchange mechanism which is suggested by a band calculation of In1-xMnxAs Akai . In this approach, the ferromagnetic coupling between localized spins is induced by the presence of carriers with strong d-band like character. These carriers originate from transition metal impurities and essentially contribute to the conductivity of the material. This model was also confronted with the difficulty from both the photoemission and soft X-ray absorption dataIshiwata , which indicate that the holes have very little d-component. Hence at present the researchers are looking for the third way which compromises the above two approach and many hypotheses have been proposed.

The pronounced spin polarization of conducting carriers in several ferromagnetic materials, which are often referred to as half-metals in which the density of state is nonzero for only one spin direction at the Fermi surface, makes these materials very attractive for spintronics applications Wolf . Among high-$`T_c`$ ferromagnetic semiconductors, half-metallic nature have been suggested for the Mn doped GaAs. Such a suggestion is supported by results of the polarized electro-luminescence YOhno  and tunneling magneto-resistance measurementsTanaka  along with ab initio calculation of the band structure of this material Jain . The half-metallicity is very important not only for the application but also for putting a strong constraint for the theories to explain the ferromagnetism. Nevertheless, since the exchange interaction energy is of the same order as Fermi and potential fluctuation energy, the existing experimental data do not allow us to arrive at decisive conclusion of the electronic structure of Mn doped GaAs and their half-metallicity.

Simultaneous probing of the spin and carrier dynamics by time-resolved magneto-optical Kerr effect (TR-MOKE) measurements gives us a chance to reveal the half-metallic nature of the material. The TR-MOKE technique has been first applied to the investigation of the picosecond photo-induced demagnetization of NiBeaurepaire . The phenomenon is analyzed by “three-temperature model ”which describes the dynamical process of thermal interaction among sub-systems(charge, spin, and lattice systems). We have employed the TR-MOKE technique to investigate the photo-induced spin and charge dynamics in ordered double perovskite Sr2FeMoO6Kise , which is known to show very large magneto-resistance at room temperature. In particular, we have discovered that the response of the spin system is much slower than the charge system and shown that such a spin-charge thermal isolation originates from the half-metallic nature of Sr2FeMoO6. In turn, the observation of the spin-charge thermal isolation may be suggested as an experimental criterion of the half-metallic nature of a material.

Although TR-MOKE is an unique method to probe ultrafast dynamics of magnetization, the extraction of the pure magnetization component from the MOKE signal is an important experimental problem because a photo-induced changes of refractive index also contributes to the pump-induced TR-MOKE signal. In the experiment with ferromagnetic metal films, the magnetization component of the signal has been extracted by comparing the results of the linear and nonlinear magneto-optical measurements in both reflection and transmission geometries Guidoni .

In this paper, we propose and develop a two color probe method of TR-MOKE applicable to opaque DMS’s, which allows us to extract magnetization component in reflection geometry. With this method, we investigate magnetization dynamics in $`\mathrm{Ga}_{0.94}\mathrm{Mn}_{0.06}\mathrm{As}`$ which has the highest Curie temperature, $`T_c`$=110K among ferromagnetic semiconductorsHOhno3 . The extracted magnetization dynamics reveals that spins thermalize much slowly in comparison with the holes indicating the half-metallicity of this material, supporting some band calculations Jain  based on a large exchange energy. A sample of $`\mathrm{Ga}_{0.94}\mathrm{Mn}_{0.06}\mathrm{As}`$ with 1.05$`\mu \mathrm{m}`$ thickness was grown by molecular beam epitaxy on a GaAs substrate and annealed at 280C.

The measurements of linear spectra and TR-MOKE were performed in the polar Kerr configuration under a magnetic field of 0.2T by alternating the field direction in a liquid N2 flow cryostat. Detailed description of the experiment was shown in the previous paperKise . In the TR-MOKE measurements, a Ti-sapphire regenerative amplifier system with a photon energy of 1.55eV, a pulse duration of 150 fs, and a repetition rate of 1 KHz was used as a light source. The pulses at the fundamental and doubled frequencies were used as a probe and a pump, respectively. The intensity of the pump pulse was 45$`\mu \mathrm{J}/\mathrm{cm}^2`$, which corresponds to the carrier density of 6.2$`\times `$10$`{}_{}{}^{19}\mathrm{cm}_{}^{3}`$. We also performed a TR-MOKE measurement at different probe photon energy under the same pump condition. Part of the white light continuum obtained with sapphire crystal was spectrally filtered and used for the probe pulse with 1.77eV photon energy. The stability of the light source limits the sensitivity of our measurements, typically 1 mdeg at 1.77eV, and 0.2mdeg at 1.55eV. The temporal resolution of 0.3ps was limited by the pump pulse duration.

Figure 1 shows the Kerr rotation and ellipticity spectrum at 77K. The arrows in Fig. 1 indicate the photon energies of the probe pulses (1.55eV and 1.77eV) used in the TR-MOKE measurement. Figure 2(a) shows the dependence of the differential reflectivity at 1.55eV under a 0.2T magnetic field at 95K on the time delay between the pump and probe pulses. One can observe a resolution-limited rapid reflectivity change, while a very slow relaxation is observed after several 100 ps. The temporal evolution of the Kerr rotation $`\mathrm{\Delta }\theta `$ and ellipticity $`\mathrm{\Delta }\eta `$ at 1.55eV and 1.77eV (Figures 2(b) and (c), respectively) are more complicated in comparison with that of Sr2FeMoO6Kise . In particular, $`\mathrm{\Delta }\theta `$ at 1.55eV is initially negative, however it changes its sign at later time. Furthermore, the strong dependence of the temporal profiles on the probe frequency clearly shows that TR-MOKE signals do not directly give information of magnetization. In order to analyze these signals, let us present the Kerr rotation and ellipticity in the following form:

$$\theta =f_\theta M,\eta =f_\eta M$$
(1)
where $`M`$ is the magnetization, while $`f_\theta `$ and $`f_\eta `$ depend on the electronic properties of the material and can be presented in terms of the refractive index and absorption coefficient. Correspondingly, the photo-induced change in $`\theta `$ and $`\eta `$ consists of two components:

$$\begin{array}{c}\mathrm{\Delta }\theta (t)f_\theta \mathrm{\Delta }M(t)+\mathrm{\Delta }f_\theta (t)M\\ \mathrm{\Delta }\eta (t)f_\eta \mathrm{\Delta }M(t)+\mathrm{\Delta }f_\eta (t)M\end{array}$$
(2)
At the probe frequencies of the present experiment, $`\theta <0`$ and $`\eta <0`$(see Fig. 1), i.e. both coefficients, $`f_\theta <0`$ and $`f_\eta <0`$. Thus, the reduction of magnetization leads the positive signal for both $`\mathrm{\Delta }\theta `$ and $`\mathrm{\Delta }\eta `$ from the first terms of Eq.(2). The second terms, $`\mathrm{\Delta }f_\theta `$ and $`\mathrm{\Delta }f_\eta `$ are the photo-induced change of Kerr signals associated with the complex refractive index change($`\mathrm{\Delta }\stackrel{~}{n}`$) caused by photo-excitation. At a probe energy of 1.55 eV (Fig. 2 (b)), we observe $`\mathrm{\Delta }\theta <0`$ and $`\mathrm{\Delta }\eta >0`$ right after the excitation. This indicates that the rapid changes in $`\theta `$ and $`\eta `$ are governed by the second term of Eq.(2).

Fig. 2(d) presents temporal profiles of $`\mathrm{\Delta }\theta `$ and $`\mathrm{\Delta }\eta `$ at probe 1.55eV (Fig. 2(b)) within the 2 ps time interval. By multiplying $`\mathrm{\Delta }\theta `$ by -4.6, both curves fit well. This indicates that in this time interval, both $`\mathrm{\Delta }\theta `$ and $`\mathrm{\Delta }\eta `$ are governed by a single dynamical parameter. Therefore, we can take $`\mathrm{\Delta }f_\eta (t)=C\mathrm{\Delta }f_\theta (t)`$ where $`C`$ is a constant in the whole temporal region. By substituting this equation into Eqs. (2) we readily arrive at the conclusion that the value of $`C\mathrm{\Delta }\theta \mathrm{\Delta }\eta `$ is proportional to the magnetization change $`\mathrm{\Delta }M(t)`$:

$$(Cf_\theta f_\eta )\mathrm{\Delta }M(t)=C\mathrm{\Delta }\theta (t)\mathrm{\Delta }\eta (t)$$
(3)
Eq.(3) shows that we can extract $`\mathrm{\Delta }M`$ from the measured values of MOKE and TR-MOKE signals and parameter $`C`$, which depend on the probe frequencyaboutC . To examine the validity of this procedure, we analyzed the data for the probe at 1.77eV to obtain parameter $`C`$ and extracted $`\mathrm{\Delta }M`$. Figure 3 shows, that the normalized $`C\mathrm{\Delta }\theta \mathrm{\Delta }\eta `$ at 1.55eV (open circles) and 1.77eV (open squares) show the same temporal behavior. The chosen values of $`C`$ are -4.6 and 1.2 for 1.55eV and 1.77eV, respectively. Our results clearly show that we can extract the magnetization component from TR-MOKE signals.

Now we discuss the temporal responses of $`\mathrm{\Delta }R/R`$ shown in Fig. 2(a) and $`\mathrm{\Delta }f_\theta `$ and $`\mathrm{\Delta }f_\eta `$ of eq.(2). In semiconductors, these changes are driven by the photo-induced change in the complex refractive index caused by the band filling effects and successive carrier heating. In $`\mathrm{Ga}_{1\mathrm{x}}\mathrm{Mn}_\mathrm{x}\mathrm{As}`$ as well as other low-temperature grown GaAs, the excess arsenic atoms form high density As antisites ($`10^{19}`$cm-3), which act as deep donor levels. The photo-generated electrons are immediately trapped in such levels and non-radiatively recombined within 1psBenjamin . It is known that the additional acceptor doping lowers the Fermi level and promotes the ionization of deep donor levels, which enhances electron trapping. In the present experiment, the high density Mn doping brings acceptor levels with density as high as 1019cm-3. Thus the trapping and non-radiative recombination of photo-excited carriers should be very efficient and the carrier lifetime is less than 1psStellmacher ; Haiml . The resolution-limited fast decay of $`\mathrm{\Delta }R/R`$ within 1ps corresponds to such a rapid non-radiative recombination process. The non-radiative process heats up the lattice and the preexisting holes. Comparing with the temperature dependence of the linear reflectivity spectrum, we estimate the temperature increase of the holes after the fast decay as about 3-5 K. Considering the small heat capacity of holes, which is four order smaller than that of the lattice, the estimated small temperature increase of holes indicates that holes and lattice should have reached a thermal equilibrium during the non-radiative recombination processes. The slow decay observed in $`\mathrm{\Delta }R/R`$ over several hundred ps is due to thermal diffusion.

After the non-radiative recombination, $`\mathrm{\Delta }\stackrel{~}{n}`$ induced by the temperature increase of holes governs the temporal responses of $`\mathrm{\Delta }R/R`$ and $`\mathrm{\Delta }f_\theta `$ and $`\mathrm{\Delta }f_\eta `$ as we have noted that all the signals are described with a single dynamical parameter shown in Fig. 2(d).

Now we discuss the origin of the slow thermalization of the extracted magnetization in Fig. 3. The observed slow thermalization of spin systems is rather unusual for GaAs based semiconductors because the valence band is strongly affected by spin-orbit interaction and the heavy hole and light hole bands are degenerated at the $`\mathrm{\Gamma }`$ point. The spin relaxation time for the hole in pure GaAs is known to be about 110fs Hilton . Moreover, with doping the magnetic impurityKikkawa , the spin relaxation time decreases due to the spin flip scattering by magnetic impurities.

However, our experimental findings indicate that the spin temperature rises slowly towards the hole and lattice temperature. According to the three temperature model Beaurepaire , the rising speed of spin temperature is determined by the heat capacities and the heat flow rates. The heat capacity of the spin system is 18$`\mathrm{J}/\mathrm{K}\mathrm{m}^3`$, which is much smaller than those of lattice and hole systems, $`1.0\times 10^6\mathrm{J}/\mathrm{K}\mathrm{m}^3`$ and $`1.8\times 10^2\mathrm{J}/\mathrm{K}\mathrm{m}^3`$, respectivelyBlakemore . Therefore, we need to consider the anomalous thermal isolation of spin systems. This resembles what we have observed in a half-metallic ferromagnet Sr2FeMoO6Kise .

After the rapid non-radiative relaxation, the subsequent carrier scattering mainly occurs near the Fermi surface. Therefore, the demagnetization rate of carriers strongly depends on their spin polarization near the Fermi level. In particular, if the density of state for one spin direction is much smaller than that of the other direction, the spin-flipping is inefficient. This situation occurs for a system where the exchange splitting of carriers’ band is larger than the Fermi energy, leading to the half-metallic electronic structure. Therefore, the observed slow demagnetization invokes this condition which holds for $`\mathrm{Ga}_{0.94}\mathrm{Mn}_{0.06}\mathrm{As}`$ as well as for Sr2FeMoO6.

The observed temporal responses can be numerically simulated by using the three-temperature modelBeaurepaire . The solid curve of Fig. 3 shows that this model can reproduce the temporal evolution of spin systems within a wide time interval. The rising time of the spin temperature is found to be 500ps. It is necessary to emphasize that the obtained coupling constant between the charge and spin systems is much smaller than that in Ni. The observed spin relaxation at longer time delay might be due to spin-lattice coupling, which stems from the magnetic anisotropyHubner .

In conclusion, we observe the magnetization dynamics in $`\mathrm{Ga}_{0.94}\mathrm{Mn}_{0.06}\mathrm{As}`$ over a wide temporal range from ps to ns using two-color TR-MOKE measurements. We discover that the spin system thermalizes much slower than the hole and lattice systems, indicating spin-charge thermal isolation. This experimental finding confirms that $`\mathrm{Ga}_{0.94}\mathrm{Mn}_{0.06}\mathrm{As}`$ has a half-metallic electronic structure, making it a prospective material for spintronic applications. The half metalicity indicates a large exchange energy of the carriers compared to the Fermi energy, which plays an essential role for the magnetism and transport properties of this material.

This work is supported by a Grant-in-Aid for Scientific for COE research (Phase Control of Spin-Charge-Photon Coupled Systems) and a Grant-in-Aid for COE Research (Quantum Dot and Its Application) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan.


text:
# Implementing the Keldysh formalism into the ab initio Gaussian Embedded Cluster Method for the calculation of quantum transport

## Abstract

We discuss the key steps that have to be followed to calculate quantum transport out of equilibrium by means of the ab initio Gaussian Embedded Cluster Method recently developed by the authors. Our main aim is to emphasize through several examples that, if a sufficiently large portion of the electrodes is included in the ab initio calculation, which does also incorporate an electrochemical potential difference $`\mu _L\mu _R=eV`$, there is no need to impose an electrostatic potential $`V`$ drop accross the system.

The effort devoted to investigate electronic transport through molecular bridges and metallic nanocontacts has sharply increased in recent yearsJoachim et al. (2000); Nitzan (2001); Agrait et al. . In particular, a variety of ab initio methods are being developed aiming to catch the essentials of these systemsLang (1995); Mujica et al. (2000); Taylor et al. (2001); Xue et al. ; Damle et al. (2001); Brandbyge et al. (2002); Damle et al. . These methods share most of the ingredients \[Landauer-Keldysh formalism for transport and density functional (DF) theory\], but differ in their numerical implementation, which translates into their computational efficiency and predictibility power. A critical problem not yet fully solved is the way these methods actually deal with out-of-equilibrium transportKeldysh (1965); Meir and Wingreen (1992); Hershfield et al. (1991). Most of the approaches that are being proposed impose from the outset an external and uniform electrostatic potential drop across the molecule or nanocontact Taylor et al. (2001); Damle et al. (2001); Brandbyge et al. (2002). The latter is justified when considering planar electrodes, which in most cases do not correspond to experimental electrode geometries. This choice simplifies the numerical implementation, but can condition the outcome of the calculation. In this work we discuss how the Keldysh formalism can be implemented into ab initio transport schemes and, in particular, into the Gaussian Embedded Cluster Method (GECM) recently developed by the authorsPalacios et al. (2001, 2002). Conclusive numerical evidence is presented which shows that, if a sufficiently large portion of the electrodes is incorporated into the ab initio calculations, there is no need to add an external electrostatic potential to the self-consistent potential. This provides a natural, consistent, way of solving the electrostatic problem for generic electrode geometries.

In a single-particle description, inherent to DF calculations, one obtains the current $`I`$ through the molecule or nanocontact at non-zero temperature and finite bias voltage $`V`$ by integrating the transmission probability:

$$I=\frac{e}{h}_{\mathrm{}}^{\mathrm{}}T(E,V)[f(E\mu _L)f(E+\mu _R)]dE,$$
(1)
where $`f(E\mu _{\mathrm{L},\mathrm{R}})`$ is the Fermi distribution for the left and right electrodes whose electrochemical potentials are denoted by $`\mu _\mathrm{L}`$ and $`\mu _\mathrm{R}`$, respectively, and $`\mu _\mathrm{L}\mu _\mathrm{R}=eV`$. The transmission coefficient $`T(E,V)`$ depends both on temperature and the applied bias, and is given by

$$T(E,V)=\mathrm{Tr}[\widehat{\mathrm{\Gamma }}_\mathrm{L}\widehat{G}^{(+)}\widehat{\mathrm{\Gamma }}_\mathrm{R}\widehat{G}^{()}],$$
(2)
where hats denote operators. The Green function operators and the gamma operators depend on energy and the applied potential (both arguments dropped for convenience). The latter are defined as $`\widehat{\mathrm{\Gamma }}_{\mathrm{L},\mathrm{R}}=\widehat{\mathrm{\Sigma }}_{\mathrm{L},\mathrm{R}}^{(+)}\widehat{\mathrm{\Sigma }}_{\mathrm{L},\mathrm{R}}^{()}`$, where $`\widehat{\mathrm{\Sigma }}_{\mathrm{L},\mathrm{R}}`$ denotes the self-energy operator for the part of the right and left semi–infinite electrodes which is not included in the ab initio calculation. The Green function operator is in its turn defined as

$$\left[(E\pm i\delta )\widehat{I}\widehat{F}\widehat{\mathrm{\Sigma }}_L^{(\pm )}\widehat{\mathrm{\Sigma }}_R^{(\pm )}\right]\widehat{G}^{(\pm )}=\widehat{I},$$
(3)
where $`\widehat{I}`$ is the unity operator. On the other hand, the density matrix out of equilibrium is obtained fromerr

$$P_{\alpha \beta }=\frac{i}{2\pi }_{\mathrm{}}^{\mathrm{}}dE\underset{\gamma \delta }{}\left[S_{\alpha \gamma }^1G_{\gamma \delta }^<S_{\delta \beta }^1\right]$$
(4)
where $`G_{\gamma \delta }^<`$ are the matrix elements in a non-orthogonal basis of the lesser Green function operator given by

$$\widehat{G}^<=i\widehat{G}^{(+)}\left[f(E\mu _L)\widehat{\mathrm{\Gamma }}_\mathrm{L}+f(E\mu _R)\widehat{\mathrm{\Gamma }}_R\right]\widehat{G}^{()}.$$
(5)
Most of the technical difficulties implicit in the evaluation of the above expressions have been discussed in previous worksTaylor et al. (2001); Damle et al. (2001); Palacios et al. (2002); Brandbyge et al. (2002). There is, however, an important conceptual issue that, in our opinion, remains unsolved and worth clarifying. In the presence of a voltage difference between electrodes one might be tempted to impose from the outset an external electrostatic potential drop across the molecule or nanocontact. As the potential profile depends strongly on the electrode geometry, one would have to solve Poisson equation with the boundary conditions appropriate for such geometry. A way around this problem is to assume a simple form for the electrodesTaylor et al. (2001); Damle et al. (2001); Brandbyge et al. (2002). However, adding an external potential through the homogeneous solution of the Poisson equation is unnecessary if a significant part of the metallic electrodes has already been included in the initial cluster and an electrochemical potential difference is maintained. The difference between left and right electrochemical potentials charges one electrode and discharges the other one, effectively creating an electrostatic potential drop across the constrictionLandauer (1957). In general, once self-consistency is attained, an electrostatic potential difference $`V`$ between atoms one or two layers inside opposite electrodes develops while they remain neutral. On the other hand, the potential difference between atoms on the surface of opposite electrodes is smaller than $`V`$ since they carry the charges. In summary, we assume all the charges responsible for the electrostatic potential across the constriction to be part of our cluster instead of imposing a usually unknown external electrostatic potential. For a large enough number of atoms in the electrodes this should be essentially correct.

The way we have actually implemented this self-contained, out-of-equlibrium scheme into the GECM initially proposed in Ref. Palacios et al., 2001 and later fully developed in Ref. Palacios et al., 2002 is hereafter described. A first step consists in carrying out a DF calculation by means of the GAUSSIAN98 codeFrisch et al.  of the region that contains the molecule or the set of atoms that form the contact between the electrodes plus a significant part of the electrodes. Once a reasonable accuracy has been achieved, the Green function for the infinite system is calculated and the density matrix given by Eq. 4 used in the subsequent DF process. The infinite electrode selfenergies are calculated by means of the Bethe lattice approximation and a minimal basis set with pseudopotentials was used in the calculations \[see Ref. Palacios et al., 2002 for details\]. In previous works, the DF calculations were carried out by means of the Becke three-parameter hybrid functional plus the correlation functional of Lee, Yang and Parr (B3LYP)Becke (1993). A difficulty specific to the out-of-equilibrium case comes from the fact that the density matrix in Eq. 4 is complex. This can be handled in two ways: i) treating exchange within DF, instead of carrying out a full Hartree-Fock calculation (this can be easily done by using BLYP instead of B3LYP), and ii) incorporating the complex density matrix into the GAUSSIAN98 code. In the present case we have chosen the former as the second procedure would have required major changes in the standard codelin . In order to get reasonably reliable results for finite bias, a root mean square error in the density matrix better than $`10^4`$ was mandatory.

Calculations have been carried out for gold and aluminum nanocontacts similar to those of the insets in Figures 1 and 3. Two types of fcc clusters have been considered: (a) two (111) planes containing 19 atoms each plus a four atom chain labeled $`M19M4M19`$ (a capacitor was described by taking out the four atom chain leaving a wide gap between the planes); (b) two (100)-oriented pyramidal clusters plus a 4 atom chain labeled $`M14M4M14`$ \[(100)-oriented pyramidal nanocontacts as $`M5M5`$, $`M14M14`$ and $`M30M30`$ have also been considered\]. Interatomic bulk distances have been taken for the whole cluster (2.86 $`\AA `$ for Al and 2.88 $`\AA `$ for Au). All results correspond to zero temperature and an external bias voltage in the range 0-5 V. The structural stability of the system models is an important issue when one is interested in the interpretation of experimental results. This is not the case here.

Figure 1 shows the average on-site energies ($`5d6s6p`$ for Au and $`3s3p`$ orbitals for Al) on all atoms of the nanocontacts $`M19M4M19`$ for zero and 5 V bias. This magnitude reflects only the electrostatic potential on each atom, not telling us anything about the chemical potential (charge) on them. The on-site energies in the planes are less dispersed in the case of gold than in aluminum. This is probably due to the simple $`s`$-character of the wavefunctions at the Fermi level that gold has (as opposed to the $`sp`$ character in aluminum). At 5 V bias, the major drop in the potential occurs at the chain/plane contacts. While the total potential drop in gold (aluminum) is 4.44 eV (4.6 eV) the drop between the first and the last chain atoms is only 1.67 eV (1.43 eV). In the case of zero bias, the results are symmetric with respect to the geometric center of symmetry, as expected. Instead, a similar symmetry is absent for 5 V. Namely, whereas the potential drop between the left electrode and the first atom in the chain is 1.92 eV (1.69 eV) for gold (aluminum) nanocontacts, it is only of 0.85 eV (1.48 eV) between the chain end and the right electrode. This feature seems common to all previously reported results and reflects not only bulk band structure features, but also details of the nanocontact geometry (see below). The transmission for $`V=0`$ and 5 V is shown in Figure 2. There are noticeable differences both in gold and aluminum. In the latter the peak structure is significantly changed, whereas in gold the gap below -0.5 eV is partially filled. It is worth noting that although the total current calculated by integrating in a window \[-$`V`$/2, $`V`$/2\] the transmission for $`V`$=0 does not differ much from that obtained with the full non-equilibrium approach (more important changes are expected in complex molecules), the differential conductance is significantly different.

The average on-site energies on all atoms of the nanocontacts $`M14M4M14`$ for zero and 5 V bias, are depicted in Figure 3. As above, the on-site energies in the outer planes are less dispersed in the case of gold than in aluminum. Now the potential drop along the chain is significantly different for the two metals: while in the case of gold the total potential drop is 4.29 eV and the drop between the first and the last chain atoms is only 0.162 eV, in aluminum these numbers change to 4.56 and 2.043, respectively. These results support our claim in the sense that what matters are not only the bulk features of the metals but also the nanocontact geometry. The asymmetry in the potential drop remarked above is also noted. Again, the differences in the transmission between zero and 5 V bias are significant. The peak structure in aluminum is noticeably reduced and the gap in gold filled over a range larger than 1 eV.

In Table I we report the fittings of the numerical results for the difference between the average energies of all atoms in the left and right planes in $`M19M4M19`$ nanocontacts, and the two outer planes in the pyramidal nanocontacts, versus the applied external potential. The fittings were done with $`aV^b`$, instead of assuming a linear relationship from the start. All regression coefficients were higher than 0.999. Although $`a`$ and $`b`$ are always rather close to 1, some significant features are worth of comment. There is a steady increase of $`a`$ and $`b`$ as the size of the outer planes increases. The values closest to one are found for the capacitors (the nanocontacts with a $`14\AA `$ gap). In the latter case the coefficient $`a`$ is only 7% (5%) smaller than 1 for gold (aluminum) nanocontacts, while the exponent cannot in practice be differentiated from one. The small deviations with respect to the ”ideal” behavior are due to the finite charge that those planes carry and/or to a yet insufficiently large electrode. As expected, in pyramidal nanocontacts the results steadily approach the ideal behavior as the size of the cluster increases. These results constitute the main message of the present work. The results for the current versus voltage were also fitted with $`I/G_0=bV^d`$ (see Table I). In this case the lowest regression coefficient is 0.99 (in some cases the results are better fitted with second order polynomia). The significant deviations from a linear behavior (Ohm’s law) are not surprising in systems as small as those investigated here and have been reported previously. The current is in general higher in aluminum than in gold, which is consistent with the larger conductance found for aluminum at the Fermi levelPalacios et al. (2002). Deviations from linearity are also more important in the case of aluminum, with the exponent increasing with the size of the pyramidal cluster. This is a consequence of a number of facts, including the increasing contributions from the $`\pi `$ channel discussed in Ref. Palacios et al., 2002.

Summarizing, we have discussed the practical procedures one should follow to implement the non-equilibrium Keldysh formalism into the Gaussian Embedded Cluster Method previously developed by the authors. In our view the most important outcome of this work is the conclusive evidence concerning the need of incorporating a sufficiently large portion of the electrodes in the ab initio calculation. If this is done, the artificial addition of an applied external potential to the potential in the Schrödinger equation is unnecessary.

Partial financial support by the spanish MCYT (grants BQU2001-0883, PB96-0085, and MAT03-04450-C03) and the Universidad de Alicante is gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

One of the key issues of present high energy physics is to check whether the Higgs field(s) account for the mechanism of the mass generation. In the standard model (SM) of the high energy physics, only one physical Higgs boson exists after symmetry spontaneously broken, while five Higgs bosons in the minimal supersymmetric standard model (MSSM).

Until now, the precise high energy experiments have tested the gauge bosons self interactions and the gauge bosons interactions with fermions (except top sector) of the SM with very good accuracy, while the Higgs sector is totally untouched . In the SM and MSSM, the couplings among Higgs bosons and top quarks play a very unique role because which are not suppressed by $`\frac{m_q}{m_W}`$. Very likely, these couplings are firstly detected at experiments. Besides the aspect of the study on the properties of the Higgs bosons, there is another strong motivation to examine the $`Ht\overline{t}`$ coupling: searching new physics beyond the SM. It is commonly thought that new physics might appear in the top sector for its huge mass.

Supersymmetric (SUSY) QCD interactions by exchanging gluinos and squarks might impact greatly on the Higgs-top quarks vertexes. If the Higgs bosons are heavy enough, the SUSY QCD effects can be studied in the decay channel $`Ht\overline{t}`$ (H is general notation for Higgs boson) . Otherwise, production channels offer alternative methods and should be carefully examined. At hadron colliders, Higgs boson can be associated production with top quark pair through $`gg`$ and/or $`q\overline{q}`$ sub-processes . Besides hadron colliders, as a supplementary machine, the linear colliders (LC) will provide ideal place to study such effects. Even in the seventies, the process $`e^+e^{}H\overline{t}t`$ in the SM at LC has been proposed, which can be used as finding Higgs boson and measuring the couplings of $`g_{H\overline{t}t}`$ . The calculations of the important QCD effects for the process are not accomplished until the end of the last century . In the MSSM, the neutral Higgs boson associated production with top quark pair at LC including the QCD effect are available soon after . However, the SUSY QCD effects on the process $`e^+e^{}t\overline{t}H`$ are not touched, which is the subject of this paper. In this paper, we focus on the lightest neutral Higgs boson $`h^0`$ associated production with top quark pair in the MSSM.

## 2 Basics and technical set-ups

The diagrams including the SUSY QCD corrections are shown in Fig.1-3. These diagrams and Feynman amplitudes are generated by FeynArts . Fig. 1 and 2 are the Feynman diagrams at tree and one-loop level, respectively. In Fig. 3, we draw the diagrams containing counter-terms. In this paper, we renormalize the quark wave function and mass in the on-shell (OS) renormaliztion scheme. Therefore in Fig. 3, we do not show the diagrams with counter-terms at external legs because in the OS renormaliztion scheme, such kind of diagrams are completely eliminated by the corresponding virtual diagrams.

The analytical results for Feynman amplitudes and phase space integrations are tedious but easily obtained through standard procedure, which are useless to almost all readers. Therefore, we decide to omit the explicit analytical results here 222Programs are available on request.. For the calculations of the amplitudes and the Monte Carlo phase space integrations, we rely heavily on the popular packages FormCalc and LoopTools  333A technical detail, LoopTools is not applicable for some $`D_0`$ with certain parameters that all $`(m_i^2+m_j^2p_{ij})^2<4m_i^2m_j^2`$, details of the definition and convention for $`D_0`$ can be found in manual of LoopTools or in Ref . For this exceptional case, the new program has been added according to Ref .. Before we shift to numerical discussions, it should be emphasized that the entire calculations and the cancellation of ultraviolet divergences have been checked carefully.

## 3 Numerical discussions and summary

In the MSSM, there are dozens of parameters and experiments constrains, for examples $`\mathrm{\Delta }\rho `$ and direct mass limits for gluino, squarks and Higgs bosons. We rely on the built-in programs of FormCalc to calculate the SUSY particles mass spectrum (we use two-loop results for Higgs mass calculations) and SUSY parameters as well as set constraints on known limits from experiments. Specifically, we set

$`\mathrm{\Delta }\rho <3\times 10^3`$ (1)
$`m_{\stackrel{~}{g}}>175GeV,`$ $`m_{\stackrel{~}{t}}>80GeV,`$ $`m_{\stackrel{~}{b}}>70GeV`$
$`m_{\stackrel{~}{q}}>150GeV,`$ $`m_{h^0}>91GeV.`$ (2)
For simplicity we have assumed the SUSY parameters $`M_Q=M_U=M_D=M_S`$ and $`A_t=A_b`$.

As for the first attempt, we would like to see the SUSY QCD effects for the non-mixing case, i.e. the non-diagonal matrix elements $`M_{LR}A_t\mu /\mathrm{tan}\beta `$ for stop mass matrix are zero. Our numerical results show that the SUSY QCD effects are small for such case, and the relative corrections $`\delta =(\sigma _{NLO}\sigma _{LO})/\sigma _{LO}`$ are at most about two percents, which are hard to detect at experiments. In the following we will focus on the mixing case.

After scanning the whole parameters space, in order to demonstrate the numerical results, in Fig. 4-6 we show the cross sections and relative corrections $`\delta `$ for a specific parameter set as a function of $`\mathrm{tan}\beta `$, center-of-mass energy $`\sqrt{s}`$ of $`e^+e^{}`$, as well as the mass of gluino. From the figures, we can see that the next-to-leading order (NLO) SUSY QCD effects decrease the LO cross sections and the magnitude can reach $`50\%`$. For such kind of large corrections, it should be noted that main contributions come from diagrams (5)-(8) of Fig. 2, and the lighter stop mass is about 100 GeV in this case. Similar to the SM QCD effects, the large corrections occur in the vicinity of the production threshold (Fig. 5), where the cross sections become small. Furthermore, from Fig. 6, it is obviously that the corrections are at the peak around $`m_{\stackrel{~}{g}}=190200`$ GeV, where the re-scattering enhancement $`t^{}\stackrel{~}{g}+\stackrel{~}{t}_1t+h^0`$ appears at $`m_{\stackrel{~}{g}}+m_{\stackrel{~}{t}_1}m_t+m_{h^0}`$.

In table 1, we show that cross sections at LO and NLO as well as the relative corrections as a function of $`A_t`$. The large SUSY QCD effects emerge only with positive large $`A_t`$ \[or $`M_{LR}`$\], which ensures that the large contributions from lighter stops don’t vanish.

From our parameters space scanning, in order to have the large SUSY QCD effects, the following requirements should be satisfied

1. there must be large mass splitting between stops, moreover $`M_{LR}`$ should be large and positive,
2. $`M_S`$ should not be large,
3. gluino mass should not be heavy.

As mentioned above, condition (1) prevents the cancellation between stops, and makes one stop mass light which is essential for large SUSY QCD effects, especially in the threshold regions. Conditions (2) and (3) ensure that SUSY effects do not decouple from the process.

To summarize, the SUSY QCD corrections to the lightest neutral Higgs boson associated production with top quark pair are studied in the MSSM at LC. Our studies show that generally the corrections are very moderate (say 10%) and under control, except for some rescattering effects which lead to a breakdown of perturbation theory and require a more detailed study, resummation or something of this kind. Contrary to the QCD corrections arising from gluons which are about 50% for $`m_{h^0}=120`$ GeV and $`\sqrt{s}=500`$ GeV \[near the production threshold\], the SUSY QCD effects can be also as large as about -50% in this kinemical region for the favorable model parameter space under OS renormalization scheme. Such kind of large SUSY QCD effects might be the probe to determine the sign of $`M_{LR}`$. We note that if nature does choose lighter stops, which can be copiously produced at hadron colliders and linear colliders. Comparison between direct productions and indirect virtual effects might give more information on the MSSM.

## Acknowledgements

The author would like to thank Dr. Thomas Hahn, Dr. Christian Schappacher and Dr. G.J. van Oldenborgh for the help on the packages used in this paper. This work was supported in part by the Nature Sciences and Engineering Research Council of Canada, the Alexander von Humboldt Foundation, National Nature Science Foundation of China.