format: text

text:
# Determination of the Fermi Surface of MgB2 by the de Haas-van Alphen effect

## Abstract

We report measurements of the de Haas-van Alphen effect for single crystals of MgB2, in magnetic fields up to 32 Tesla. In contrast to our earlier work, dHvA orbits from all four sheets of the Fermi surface were detected. Our results are in good overall agreement with calculations of the electronic structure and the electron-phonon mass enhancements of the various orbits, but there are some small quantitative discrepancies. In particular, systematic differences in the relative volumes of the Fermi surface sheets and the magnitudes of the electron-phonon coupling constants could be large enough to affect detailed calculations of Tc and other superconducting properties.

There has been rapid progress in understanding the physical properties of magnesium diboride since the discovery of superconductivity at 39 K just over two years ago. The accepted consensus is that it is an $`s`$-wave, phonon-mediated superconductor but with some highly unusual properties. The most important of these are the anomalously high $`T_c`$ and the existence of two, almost distinct, superconducting gaps. This understanding is based on many experiments and on theoretical calculations of the unusual electronic structure of MgB2 kortus01 ; liu01 ; choi02 ; mazin02 .

There have been two direct experimental probes of the Fermi surface (FS) structure of MgB2; angle-resolved photoemission spectroscopy uchiyama  and the de Haas-van Alphen effect (dHvA)yelland02 . Our previous dHvA study yelland02  showed good agreement between theory and experiment regarding the areas of the orbits that were observed and their electron-phonon mass enhancements. However, only three out of nine predicted dHvA orbits were observed (labelled 1,2 and 3 in Fig. 1) and so information was only obtained about two of the four FS sheets predicted by band calculations. It was not entirely clear whether the unobserved orbits were missing simply because of their relatively short mean-free-paths, or whether the topology of the other sheets was substantially different. This issue clearly affects calculations of many physical properties (especially $`T_c`$). In this Letter, we report new measurements of the dHvA effect in magnetic fields up to 32 Tesla for two single crystals. These new results give evidence for orbits on all four sheets of the Fermi surface.

Experiments were conducted on two different crystals grown by different groups. Sample B is one of the crystals studied previously yelland02  and was grown in Tokyo by high pressure synthesis using natural boron. Sample K was grown in Zürich by a similar technique but isotopically pure boron-10 was used as a starting ingredient b10note . Other crystals from the same batches were found to have $`T_c`$ (onset) of 38.0 K and 37.7 K for the Tokyo and Zürich crystals respectively. dHvA oscillations were observed by measuring the torque ($`\mathrm{\Gamma }`$) with a sensitive piezo-resistive cantilever technique bergemannphd . Both samples were studied extensively at Bristol in fields up to 20.5 T and temperatures down to 0.3 K, and in fields up to 32 T at the NHMFL.

The first harmonic of the oscillatory part of the torque for a 3D Fermi liquid is given by shoenberg ; wasserman

$$\mathrm{\Gamma }_{osc}\frac{1}{C}\frac{dF}{d\theta }B^{\frac{1}{2}}R_DR_TR_S\mathrm{sin}\left(\frac{2\pi F}{B}+\gamma \right)$$
(1)
where $`F`$ is the dHvA frequency \[$`F=(\mathrm{}/2\pi e)A`$, $`A`$ is the extremal orbit area in $`𝐤`$-space\]; $`C`$ is the curvature factor, $`\gamma `$ is the phase; $`R_D`$, $`R_T`$ and $`R_S`$ are the damping factors from impurity scattering, temperature and spin splitting respectively. The Dingle factor, $`R_D=\mathrm{exp}(\frac{\pi m_B}{eB\tau })`$, where $`m_B`$ is the unenhanced band mass shoenberg ; wasserman  and $`\tau `$ is the scattering time. An equivalent expression for $`R_D`$ is $`R_D=\mathrm{exp}(\frac{\pi \mathrm{}k_F}{eB\mathrm{}})`$, which shows clearly the increased damping with shorter mean-free-path $`\mathrm{}`$ and increased average Fermi wavevector ($`k_F`$) of the orbit. $`R_T=X/(\mathrm{sinh}X)`$ where $`X=\frac{2\pi ^2k__B}{\mathrm{}e}\frac{m^{}T}{B}`$, $`m^{}`$ is the quasi-particle effective mass. Finally, the spin splitting factor is given by $`R_S=\mathrm{cos}(\frac{\pi ngm_B(1+S)}{2m_e})`$ where $`(1+S)`$ is the orbitally averaged exchange-correlation (Stoner) enhancement factor, $`g`$ is the electron $`g`$-factor, $`m_e`$ is the free-electron mass and $`n`$ is an integer.

The electronic structure and dHvA orbits of MgB2 have been calculated by three different groups harima ; mazin02prb ; rosner02 . The calculated Fermi surface is shown in Fig. 1, together with the expected dHvA extremal orbits harimanote . The calculations kortus01  show that the electronic states near the Fermi level arise primarily from the boron atomic orbitals. The calculated Fermi surface is composed of four distinct sheets. Two of these arise from boron $`\sigma `$ orbitals and are quasi-two dimensional warped cylinders running along the $`c`$-direction, whereas the other two are tubular networks with larger $`c`$-axis dispersion, that are mainly formed from the boron $`\pi `$ orbitals. In total 9 primary extremal orbits have been predicted, and 7 of these are labelled in the Figure (two orbits with frequencies $`>`$30kT have been omitted). Calculations of the dHvA frequencies and masses of the various orbits by the three groups are all in good agreement (the differences are typically 100-200 T, and $`5`$% respectively).

The fast Fourier transforms (FFTs) of the raw torque data between 20 and 32 T at 1.4 K for crystal K are shown in Fig. 2, as a function of angle $`\theta `$ as the sample was rotated from $`Hc`$ to $`Ha`$ ($`Hc0^{}`$). In addition to the frequencies observed in our previous study yelland02  ($`F_1`$, $`F_2`$ and $`F_3`$) several additional peaks are visible. Some care is needed in interpreting these new frequencies. Because the cantilever is deflected slightly by the torque, the torque measurements are not made at constant angle. This generates spurious harmonics and combinations of the main dHvA oscillations. For crystal K only weak harmonics of F2 are observed (just above $`F_5`$– see Fig. 2), whilst for the larger crystal B, several harmonics of $`F_3`$ and a frequency corresponding to $`F_1+F_2`$ were observed.

The observed frequencies (omitting those assigned as harmonics or combinations) are shown in Fig. 3 as a function of $`\theta `$ for both crystals. The solid lines in Fig. 3 are fits of the observed $`F(\theta )`$ values to $`F(\theta )=_{i=1}^3a_i/\mathrm{cos}^i(\theta \theta _0)`$ ($`\theta _0=0`$ or $`90^{}`$), which we use to extrapolate the observed frequencies to the symmetry points (we denote these frequencies as $`F_n^0`$). The assignment of the FFT peaks to the orbits shown in Fig. 1 was achieved by comparing the values of the frequencies obtained by extrapolation, and their angular dependencies, with the calculations.

For crystal K, signals from 6 orbits associated with all 4 sheets of Fermi surface are observed. We are therefore able to verify experimentally the Fermi surface topology predicted by Kortus et al. kortus01  and shown in Fig. 1. For crystal B, in addition to the 3 orbits observed previously, 3 further frequencies are seen. One of these can be assigned to $`F_4`$. The frequency of $`F_3^{}`$ (F=4600 T) is close to that predicted for orbit 7 (F=4294 T)harima  (and has a similar effective mass), however, the frequency is also close to that expected for an orbit equivalent to $`F_3`$ but from the tube oriented along $`a^{}`$ (i.e., at $`60^{}`$ to $`a`$). A subsequent in-plane rotation study showed the latter to be the most likely origin. The origin of $`F_{3S}`$ is less clear, but it could arise from a slight warping of the same in-plane tube responsible for $`F_3`$.

Table 1 shows that the values of $`F^0`$ found for the two crystals are in good agreement. In total, we have studied crystals from six batches (5 from Tokyo and one from Zürich) and so far, the dHvA frequencies agree to within 30 T or 0.06 % of the basal area of the first Brillouin zone (50.2 kT). This suggests that although $`T_c`$ is slightly reduced compared with the best polycrystalline samples, any possible Mg deficiency is very reproducible, and probably small cooper03 .

The differences between the measured $`F^0`$ values and those predicted by theory is a significant fraction of $`F^0`$ in some cases, but we note that 100 T discrepancy only amounts to $``$ 0.3 % of the basal area of the Brillouin zone. The volumes of the tubes are proportional to the average of the two extremal areas rosner02  and these are both $`16`$% smaller than the calculations mazin02prb , implying a corresponding reduction in the number of holes in these two tubes rosner02 . This may have a significant effect on calculations of physical properties such as $`T_c`$, its dependence on alloying or pressure, and the London penetration depth.

It is instructive mazin02prb ; rosner02  to calculate the Fermi energy shift $`\mathrm{\Delta }E`$ which would bring the theoretical frequencies in line with experiment. As the dHvA band mass is defined as $`m_B=\frac{\mathrm{}^2}{2\pi }\frac{A}{E}`$, the necessary band shifts are given by $`\mathrm{\Delta }E=\frac{\mathrm{}e}{m_B}\mathrm{\Delta }F`$ (where $`\mathrm{\Delta }F=F_{\mathrm{Th}}F_{\mathrm{exp}}`$). There is remarkable consistency between the $`\mathrm{\Delta }E`$ values for the orbits, with the values roughly falling into two groups (Table 1). For the $`\sigma `$ sheet orbits (1,2,5,6) the average shift is $`83\pm 4`$ meV, whereas for the $`\pi `$ sheet orbits (3,4) it is $`61\pm 5`$ meV. Because of the high degree of reproducibility of the frequencies between samples, it is unlikely that this discrepancy is caused by sample impurities or non-stoichiometry. Instead, it seems to imply a shortcoming of the LDA calculationsmazin02prb ; rosner02 .

Quasiparticle effective masses were determined by performing field sweeps at different temperatures and fitting the amplitudes to Eq. 1. Our experimental values of $`m^{}`$ are compared with the calculated band masses in Table 1. In MgB2 the dominant source of mass enhancement is the electron phonon interaction. If we assume this is the only source of enhancement we can calculate an upper bound for the electron-phonon coupling constants $`\lambda ^{\mathrm{ep}}`$, from $`\lambda ^{\mathrm{ep}}=m^{}/m_B1`$. The results (Table 1) show that the values of $`\lambda ^{\mathrm{ep}}`$ on both the $`\sigma `$ sheets are approximately a factor three larger than those on the $`\pi `$ sheets.

A detailed comparison with the orbit-resolved theoretical values mazin02prb  of $`\lambda ^{\mathrm{ep}}`$ is also shown in the table. Generally our values are slightly smaller than the theoretical ones. The most significant differences are for the larger orbits on the $`\sigma `$ tubes for which $`\lambda ^{\mathrm{ep}}`$ values which are $`20\%`$ smaller than theory. The small shifts in Fermi level described above do make small differences to the calculated band masses mazin02prb , increasing the $`\lambda _{\mathrm{exp}}^{\mathrm{ep}}`$ values by $``$ 6%, but the differences remain significant. These differences could have a relatively large effect on the detailed calculations of $`T_c`$. For example, using the isotropic McMillan equation, with parameters appropriate to MgB2 choi02  we estimate a 20% reduction in $`\lambda ^{\mathrm{ep}}`$ would imply a 8 K reduction in $`T_c`$. It is likely that the reduced $`\lambda `$ values we observe are caused by phonon anharmonicity, which has been shown to reduce the average value of $`\lambda `$ by around 20 %choi02 .

Scattering rates ($`\tau ^1`$) were calculated by fitting the raw torque versus field data to Eq. 1. The results are also shown in Table 1. It can be seen that the mean free path ($`\mathrm{}`$) for sample K is significantly longer for orbits $`F_2`$, $`F_1`$ on the $`c`$-axis tube but not for $`F_3`$ on the in-plane tubular network. The reason for this is not clear, but it does not necessarily indicate a higher crystal quality. However, $`F_5`$ and $`F_6`$ are almost certainly observable in crystal K because of the larger values of $`\mathrm{}`$ for the two $`c`$-axis $`\sigma `$ tubes. These differences in $`\mathrm{}`$ may be at least partially responsible for a significant difference in $`B_{c2}`$ values for the two crystals. For crystal B, $`B_{c2}^c`$ = 3.3 T, $`B_{c2}^c`$ = 17 T, and for crystal K, $`B_{c2}^c`$ = 2.5 T, $`B_{c2}^c`$ = 12 T (all values quoted at $`T`$=0.3 K).

The longer mean free path of the smaller $`\sigma `$ tube ($`F_1`$ and $`F_2`$) in crystal K, has allowed us to track the amplitude of the signal over a very wide range of angle (up to 67 and 81 respectively). The data shown in Fig. 4 show pronounced minima at $`\theta =51\pm 1^{}`$, $`\theta =53\pm 1^{}`$, $`\theta =75.7\pm 0.5^{}`$ for $`F_1`$, $`F_2`$, and $`F_3`$ respectively, which we attribute to the spin-zero effect (the $`R_s`$ damping term in Eq. 1). Using the calculations of $`m_B(\theta )`$ by Harima harima  we can deduce the enhancement of the spin susceptibility ($`1+S`$) for these three orbits (at the spin zero angle). Taking $`g`$ = 2, we find that $`S`$=0.07, 0.12 and 0.45 for the three orbits respectively. The reason for the significantly larger enhancement on the $`\pi `$ band is not clear. Band structure calculations mazin02prb  predict $`S=0.31`$ and 0.26 on the $`F_2`$ and $`F_3`$ orbits respectively. The calculation therefore overestimates $`S`$ on the $`\sigma `$ sheet by a factor 2.5 and underestimates it on the $`\pi `$ sheet by a factor 1.7.

The solid lines in Fig. 4 show a fit to the data for $`F_1`$ and $`F_2`$ with Eq. 1, using the $`m_B`$ values of of Harima harima . There are 3 free parameters: the overall amplitude, $`\tau ^1`$ (assumed constant as a function of angle), and $`S`$. The fit is remarkably good and the values of $`\mathrm{}`$ were found to be 410 Å and 900 Å for $`F_1`$ and $`F_2`$ respectively which are close to the values obtained from a fit to the field dependence at a single angle given in Table 1. However, it was not possible to fit the data for $`F_3`$ because the abrupt fall for $`\theta 68^{}`$ cannot be explained simply in terms of the angular dependence of mB and a constant scattering rate (Dingle plots show that $`\tau `$ is approximately constant with angle).

In conclusion, our data strongly support the overall topology of the predicted electronic structure of MgB2 and the calculations of the electron-phonon coupling constants for the different orbits. Our data give direct evidence that the electron-phonon interaction is large on both $`c`$-axis $`\sigma `$ sheets and much smaller on both $`\pi `$ sheets. We have therefore obtained conclusive evidence in favor of the two key ingredients in the two-gap model of superconductivity in this compound, namely the Fermi surface topology and the disparity in the electron-phonon coupling for the $`\sigma `$ and $`\pi `$ bands.

We thank a number of people for useful conversations and insights: N.W. Ashcroft, S. Drechsler, H. Harima, S.M. Hayden, and I. Mazin. This work was supported by the EPSRC (U.K.), NEDO (Japan) and the NSF through grant number NSF-DMR-0084173. PJM gratefully acknowledges the support of the Royal Society (London).


text:
# hep-ph/0304231 LA-UR-03-2566 Using 𝜋₂(1670)→𝑏₁(1235)𝜋 to constrain hadronic models

## 1 Analyses of experimental data on $`\pi _2(1670)b_1(1235)\pi `$

Recently, the VES Collaboration published for the first time an upper bound of $`0.0019`$ on the branching fraction for Br$`[\pi _2b_1\pi ]`$, at the $`97.7`$% confidence level. This branching fraction is measured in $`37`$ GeV $`\pi ^{}`$ collisions on a nucleus, in the reaction $`\pi ^{}A\omega \pi ^{}\pi ^0A^{}`$ . This small branching fraction is consistent with a preliminary analysis performed by the E852 Collaboration  of data on the reaction $`\pi ^{}p\omega \pi ^{}\pi ^0p`$, in collisions of an $`18`$ GeV $`\pi ^{}`$ beam with a proton target.

The decay $`\pi _2b_1\pi `$ is allowed by conservation of parity, angular momentum, isospin and G-parity, and so its strength should be comparable with that of other decays which are allowed by the same quantum numbers, which are conserved to an extraordinary degree by the strong interactions. In order to show that the branching ratio is small for dynamical reasons, independent of any model, factors due to phase space and flavor should be removed. The standard expression for the partial width is

$$\mathrm{\Gamma }=\frac{p}{8\pi (2J_{\pi _2}+1)m_{\pi _2}^2}|p^Lf|^2$$
(1)
where $`m_{\pi _2}`$ and $`J_{\pi _2}`$ are the mass and total angular momentum of the decaying $`\pi _2`$, the decay momentum $`p`$ is measured in the rest frame of the $`\pi _2`$, the relative orbital angular momentum of the decay products is $`L`$, and $`p^Lf`$ is the decay amplitude. The amplitude with the phase space ($`p^L`$) and flavor ($`f`$) factors removed is $``$. In Table 1 we show the ratios of $`||^2`$ for the observed decay modes of the $`\pi _2`$ to that of the dominant decay mode ($`f_2\pi `$). A further refinement is to remove the dependence on the kinematics of the decays from the form factors of the initial and final mesons. With universal Gaussian wave functions for the mesons, this can be accomplished by defining $`=\mathrm{exp}(p^2/[12\beta ^2])\stackrel{~}{}`$, where $`\beta =0.4`$ GeV .

The ratios of the squares of these amplitudes with the flavor, phase space, and kinematic factors removed is also shown in Table 1. It is evident that the $`b_1\pi `$ decay is a factor of between 3 and 11 weaker than the other decay modes for dynamical reasons, making it anomalously small. This is emphasized by Fig. 1, which shows the $`||^2`$ ratios plotted logarithmically. Since there is only an experimental upper bound on the $`b_1\pi `$ mode, this suppression factor could be even larger. There is also evidence from recent analyses of E852 data  of a $`\pi _2(1670)`$ signal in the $`f_1\pi `$ and $`a_2\eta `$ final states. The discovery of additional final states will have the effect of further reducing the $`b_1\pi `$ branching fraction. We urge future experiments to put more restrictive bounds on the $`\pi _2b_1\pi `$ decay mode.

## 2 Models that can accommodate $`\pi _2(1670)b_1(1235)\pi `$

The decay $`\pi _2b_1\pi `$ is particularly clean in the sense that it is only sensitive to OZI allowed decays. This is because OZI-forbidden decay processes, which allow the creation of either the isovector $`\pi _2`$, $`b_1`$ or $`\pi `$ out of isoscalar gluons, are forbidden by isospin symmetry (see Fig. 2). The suppression of isospin symmetry breaking amplitudes is much greater than that of OZI forbidden amplitudes, the latter being about a factor of 10.

In non-relativistic quark-pair-creation models, where OZI-allowed meson decay processes are modeled by an initial $`q\overline{q}^{}`$ pair decaying to the two pairs $`q\overline{q}^{\prime \prime }`$ and $`q^{\prime \prime }\overline{q}^{}`$ (see Fig. 3), a simple selection rule arises when all the mesons have quark-spin $`S=0`$. If the $`q^{\prime \prime }\overline{q}^{\prime \prime }`$ pair is created with quark-spin $`S_{\text{pair}}=1`$, then conservation of quark-spin implies that the amplitude is zero . In the quark model, conventional mesons with $`S=0`$ have $`J^{PC}=0^+,\mathrm{\hspace{0.33em}1}^+,\mathrm{\hspace{0.33em}2}^+,\mathrm{\hspace{0.33em}3}^+,\mathrm{\hspace{0.33em}4}^+,\mathrm{\hspace{0.33em}5}^+,\mathrm{}`$, of which only states corresponding to the first three $`J^{PC}`$ have been established experimentally . The isovector resonances with these three $`J^{PC}`$ and in their radial ground states are $`\pi ,b_1`$ and $`\pi _2`$, respectively. The only kinematically allowed decay involving these three $`S=0`$ resonances is $`\pi _2b_1\pi `$. Moreover, all other kinematically allowed decays involving $`\pi ,b_1,\pi _2`$, and their isoscalar partners, are forbidden by the quantum numbers conserved by the strong interaction. The first explicit mention of the quark-spin selection rule or its application to $`\pi _2b_1\pi `$ was in Refs. , although it is implicit in Ref. .

No other strong decay involving conventional mesons composed of quarks other than $`u,d`$ quarks currently appears to be able to test the selection rule. Decays $`q\overline{q}q\overline{q}+q\overline{q}`$ with $`q\{s,c,b\}`$, where each meson is in its radial ground state with the $`S=0`$ quantum numbers $`J^{PC}=0^+,\mathrm{\hspace{0.33em}1}^+`$ or $`2^+`$, are forbidden for the same reasons as decays between the isoscalar resonances above. With the exception of the pseudoscalars, quark-model mesons with the open flavor structure $`K,D,D_s,B,B_s`$ or $`B_c`$, and lying on the un-natural parity sequence $`J^P=0^{},\mathrm{\hspace{0.33em}1}^+,\mathrm{\hspace{0.33em}2}^{},\mathrm{\hspace{0.33em}3}^+,\mathrm{\hspace{0.33em}4}^{},\mathrm{\hspace{0.33em}5}^+,\mathrm{}`$ are mixtures of $`S=0`$ and $`S=1`$ states, since $`S=1`$ components are no longer excluded by charge conjugation symmetry. In QCD, if one of the initial or final mesons in the decay has this open flavor structure, a second meson must also. This implies that the selection rule can only be tested in decays involving open flavor mesons if there are two open flavor pseudoscalar mesons involved. Since two pseudoscalars with an arbitrary relative angular momentum couple to the natural-parity sequence $`J^P=0^+,\mathrm{\hspace{0.33em}1}^{},\mathrm{\hspace{0.33em}2}^+,\mathrm{\hspace{0.33em}3}^{},\mathrm{\hspace{0.33em}4}^+,\mathrm{\hspace{0.33em}5}^{},\mathrm{}`$, the $`S=0`$ selection rule cannot be tested with decays involving open flavor mesons. It is, therefore, evident how central and unique the decay mode $`\pi _2b_1\pi `$ is for testing this selection rule.

The selection rule obtains when the $`\pi _2,b_1`$ and $`\pi `$ are treated non-relativistically as $`S=0`$ mesons in the quark model. Remarkably, relativistic interactions cannot introduce $`S=1`$ components in the $`\pi _2`$, $`b_1`$ and $`\pi `$ wave functions, so that the selection rule remains valid after relativistic interaction corrections to the quark model. This is because the $`q\overline{q}^{}`$ Fock state wave function of the $`\pi _2`$ can only have $`{}_{}{}^{1}D_{2}^{}`$ quantum numbers before relativistic interactions, and the interactions cannot change that. The analogous argument holds for $`b_1`$ and $`\pi `$. Even in the fully relativistic equal-time Bethe-Salpeter equation the selection rule is exact . It remains an open question whether a selection rule would be found in field theoretic calculations of $`\pi _2b_1\pi `$, e.g. in the lattice QCD, QCD sum rule, and Dyson-Schwinger equation approaches.

It has been pointed out that a success of the non-relativistic $`{}_{}{}^{3}P_{0}^{}`$ pair-creation model (Fig. 3), where $`S_{\text{pair}}=1`$, is the fact that the decay $`\pi _2b_1\pi `$ is predicted to vanish . Other decay models where $`S_{\text{pair}}=1`$, such as the non-relativistic chromo-electric string-breaking model where the pair has $`{}_{}{}^{3}S_{1}^{}`$ or $`{}_{}{}^{3}D_{1}^{}`$ quantum numbers  (Fig. 3), will also have this suppression. Both the $`{}_{}{}^{3}P_{0}^{}`$ and $`{}_{}{}^{3}S_{1}^{}`$ models involve a decay operator proportional to $`𝝈𝐩`$, where the $`𝝈`$ is the spin of the created quark anti-quark pair, and $`𝐩`$ is a momentum operator. It is not surprising that the $`{}_{}{}^{3}P_{0}^{}`$, $`{}_{}{}^{3}S_{1}^{}`$ and $`{}_{}{}^{3}D_{1}^{}`$ models obey the selection rule, since these all treat the quarks non-relativistically, as though they are heavy. This is a special case of a result that is shown in Appendix A: when each of the mesons participating in the decay is composed of a very heavy quark and anti-quark, the selection rule is exact to all orders of QCD perturbation theory.

Since ’t Hooft’s instanton-induced six-quark vertices only affects strong decays where all participating mesons have J=0, and their singlet flavor structure requires the presence of a strange quark (and anti-quark), decay models based on these vertices also predict vanishing $`\pi _2b_1\pi `$ decay .

## 3 One-boson emission models

The one-boson exchange (OBE) model describes the coarse features of the baryon spectrum as being due to confinement and the exchange of pseudoscalar  and scalar and vector  bosons between the quarks. For light-quark baryons an important pseudoscalar exchange potential comes from pion exchange. This model is not applied to meson spectroscopy. Two reasons are often given for this. The first is that if the light pseudoscalar bosons are the pseudo-Goldstone bosons of spontaneously-broken chiral symmetry, then it is inconsistent to also treat them as quark-anti-quark bound states and allow OBE to act between the quark and anti-quark. This argument would not appear to be applicable to heavier quark-anti-quark bound states such as the $`\pi _2(1670)`$ and $`b_1(1235)`$.

A second reason for not applying this model to the meson spectrum is that if one-boson-exchange in baryons is viewed microscopically, with the pion treated as a $`q\overline{q}`$ pair, an exchange of quarks in the process $`qq^{}q^{}q`$ can be viewed in one time ordering as an exchange of $`q\overline{q}^{}`$, which can be identified with a meson. In a meson the exchange of quarks occurs in the process $`q\overline{q}^{}\overline{q}^{}q`$, which in one time ordering is the exchange of a di-quark $`qq^{}`$ and not a meson. Exchange of mesons like the pion between quarks is, therefore, not expected to be important to the structure of mesons, even if it is important for baryons.

Once one admits a quark-pseudoscalar meson vertex as employed in baryon spectroscopy, this vertex naturally leads baryons to decay to a baryon plus a pseudoscalar meson, and mesons to decay to a meson and a pseudoscalar meson. For this reason the OBE model of baryon spectroscopy implies a one-boson emission decay model in baryons and in spatially excited mesons. This model should, therefore, be confronted with $`\pi _2b_1\pi `$.

In the $`{}_{}{}^{3}P_{0}^{}`$, $`{}_{}{}^{3}S_{1}^{}`$ and $`{}_{}{}^{3}D_{1}^{}`$ models, pionic decay of mesons proceeds via $`q\overline{q}^{}`$ pair decaying to the two final meson pairs $`q\overline{q}^{\prime \prime }`$ and $`q^{\prime \prime }\overline{q}^{}`$, one of which is identified with the pseudoscalar boson. As shown in Fig. 3, the one-pion emission model has either $`qq^{\prime \prime }\pi `$, or $`\overline{q}^{}\overline{q}^{\prime \prime }\pi `$. The lowest order one-pion coupling to the quark or anti-quark is given by the Lagrangian density

$$_\pi =i\frac{g_A^q}{2f_\pi }\overline{\psi }(x)\gamma _5\gamma _\mu ^\mu \stackrel{}{\pi }(x)\stackrel{}{\tau }\psi (x)+\text{H.c.}$$
(2)
An expansion of this axial current gives a decay operator of the form $`𝝈_q𝐤`$ (Eqs. 2 and 28 of Ref. ), where $`𝝈_q`$ is the spin of the quark emitting the pion, and $`𝐤`$ is the pion momentum. This means that the operator creating the boson is a vector operator in the space of the spin of the decaying meson, and so cannot link an initial $`S=0`$ meson to a final $`S=0`$ meson, so the selection rule is also valid for lowest order one-boson emission.

We conclude that the phenomenologically successful pair-creation model for light-light mesons (the $`{}_{}{}^{3}P_{0}^{}`$ model) , the chromo-electric string-breaking model ($`{}_{}{}^{3}S_{1}^{}`$ or $`{}_{}{}^{3}D_{1}^{}`$ model), instantons , and the lowest order one-boson emission model, which has successfully been applied to the decay of heavy-light mesons , are consistent with the experimental decay width of $`\pi _2b_1\pi `$.

## 4 Models possibly constrained by $`\pi _2(1670)b_1(1235)\pi `$

Higher order contributions in one-boson emission models contain terms that are not of the form $`𝝈_q𝐩`$, which violate the selection rule. An example is interactions where both a pseudoscalar boson is emitted, and a particle is exchanged between the quark and anti-quark in the initial meson (Eqs. 13, 38 and 39 of Ref. ). The amplitudes corresponding to the higher order contributions can be similar in size to those corresponding to the lowest order contribution111See Table 4 of Ref. . Note that the size of the part of the higher-order interaction that is not of the form $`𝝈𝐩`$ is not evaluated in Ref. .. This suggests that consistency with the small decay branch for $`\pi _2b_1\pi `$ can constrain models which do not obey the selection rule, such as the higher order contributions introduced in one-pseudoscalar-boson emission models  to cure problems with the lowest order contribution . It can also provide a viability check on proposed decay mechanisms. An example, depicted in Fig. 3, is where there is a single gluon exchanged between a quark in the decaying hadron and the vertex at which the quark pair is created. Although this one-gluon exchange quark pair creation decay mechanism violates the selection rule222One-gluon exchange involves both Coulomb and transverse interactions. The former has a simple $`𝝈𝐩`$ pair creation operator, but the latter involves both spin vector pair creation, and an additional term at the vertex where the quark or anti-quark emits a gluon (See Eqs. B5-B7 of Ref. ). This additional term includes a $`𝝈𝐩/m`$ contribution , so that the overall transverse gluon interaction has spin vector operators at both interaction vertices of the gluon, giving rise to a violation of the selection rule., it is found to be sub-dominant relative to the $`{}_{}{}^{3}P_{0}^{}`$ model , so that it is not expected to be constrained by $`\pi _2b_1\pi `$. If appreciable strength for $`\pi _2b_1\pi `$, inconsistent with experiment, is predicted by either higher order terms present in the one-boson emission decay mechanism, or by the one-gluon exchange pair creation decay mechanism, one of these decay models could be ruled out. This could distinguish between the OBE and one-gluon exchange models of the coarse features of the light baryon spectrum.

Even though the main models commonly applied to strong decays have been discussed, a comprehensive discussion of all proposed decay mechanisms has not been given. Such mechanisms should be confronted with the experimental data on $`\pi _2b_1\pi `$.

## 5 Further constraints due to $`\pi _2(1670)b_1(1235)\pi `$

In addition to aspects of the decay models discussed in the previous section, further breaking of the selection rule can arise from mixing with other Fock states. The mixing of mesons participating in the decay with non-$`q\overline{q}^{}`$ Fock states is constrained by the experimentally measured $`\pi _2b_1\pi `$ width. Examples of such mixing are mixing between the $`S=0`$ meson $`\pi _2`$ and the $`S=1`$ hybrid $`\pi _2`$ meson expected nearby in mass, and non-mesonic Fock states in the pseudo-Goldstone boson $`\pi `$.

## Acknowledgments

Helpful discussions with P. Eugenio and V.A. Dorofeev are gratefully acknowledged. This research is supported by the U.S. Department of Energy under contracts DE-FG02-86ER40273 (SC), and W-7405-ENG-36 (PRP).

## Appendix A Appendix: The quark-spin selection rule is exact for heavy quarks

The quark-gluon interaction in the QCD Lagrangian density (suppressing flavor and color) is

$$=g\overline{\psi }(x)\gamma _\mu A^\mu (x)\psi (x)+\text{H.c.}$$
(3)
Second quantize the free quark fields in the usual way,

$$\psi (x)=\frac{d^3p}{\sqrt{(2\pi )^32E(p)}}\underset{\nu }{}\left[a^\nu (p)u^\nu (p)e^{ipx}+b^\nu (p)v^\nu (p)e^{ipx}\right],$$
(4)
where $`a^\nu (p)`$ and $`b^\nu (p)`$ are the quark and anti-quark annihilation operators. Substituting Eq. 4 into Eq. 3 yields

$``$ $`=`$ $`g{\displaystyle }{\displaystyle \frac{d^3pd^3p^{}}{(2\pi )^3\sqrt{(2E(p))(2E(p^{}))}}}{\displaystyle \underset{\nu \nu ^{}}{}}[u^\nu (p)\gamma _0\gamma _\mu u^\nu ^{}(p^{})A^\mu (x)e^{i(p^{}p)x}a^\nu (p)a^\nu ^{}(p^{})`$ (5)
$`+`$ $`v^\nu (p)\gamma _0\gamma _\mu v^\nu ^{}(p^{})A^\mu (x)e^{i(pp^{})x}b^\nu (p)b^\nu ^{}(p^{})`$
$`+`$ $`u^\nu (p)\gamma _0\gamma _\mu v^\nu ^{}(p^{})A^\mu (x)e^{i(p+p^{})x}a^\nu (p)b^\nu ^{}(p^{})`$
$`+`$ $`v^\nu (p)\gamma _0\gamma _\mu u^\nu ^{}(p^{})A^\mu (x)e^{i(p+p^{})x}b^\nu (p)a^\nu ^{}(p^{})]+\text{H.c.}`$
The first and second terms describe the quark and anti-quark interactions with the gluon field, respectively, the third term describes creation of a quark-anti-quark pair, and the fourth term annihilation of a quark-anti-quark pair.

In the limit of very heavy quarks

$$u^\nu (p)=\sqrt{2m_Q}\left(\begin{array}{c}\chi ^\nu \\ 0\end{array}\right)v^\nu (p)=\sqrt{2m_Q}\left(\begin{array}{c}0\\ \chi ^\nu \end{array}\right),$$
(6)
where the $`\chi ^\nu `$ are the usual Pauli spinors. Then the first and second terms in Eq. 5 contain

$$u^\nu (p)\gamma _0\gamma _\mu u^\nu ^{}(p^{})=v^\nu (p)\gamma _0\gamma _\mu v^\nu ^{}(p^{})=2m_Q\chi ^\nu \chi ^\nu ^{}\delta _{\mu 0}=2m_Q\delta _{\nu \nu ^{}}\delta _{\mu 0},$$
(7)
so quark-gluon and anti-quark-gluon interactions do not change the spin of heavy quarks or anti-quarks. The third and fourth terms in Eq. 5 contain

$$u^\nu (p)\gamma _0\gamma _\mu v^\nu ^{}(p^{})=v^\nu (p)\gamma _0\gamma _\mu u^\nu ^{}(p^{})=2m_Q\chi ^\nu \sigma _i\chi ^\nu ^{}\delta _{\mu i}$$
(8)
where $`i\{1,2,3\}`$. Hence quark-anti-quark pair creation and annihilation involve a spin change described by the Pauli matrices $`\sigma _i`$.

The spin of a propagating heavy quark remains unchanged by quark-gluon interactions, according to the first and second terms of the interaction in Eq. 5, and Eq. 7. The exception to this is when the the quark travels in a Z-graph, which corresponds to quark-anti-quark pair creation and then annihilation via the third and fourth terms of the interaction in Eq. 5. However, these Z-graphs are suppressed by powers of $`1/m_Q`$, so that for very heavy quarks they do not contribute. The spin of a propagating heavy quark remains unchanged to all orders in QCD perturbation theory.

This implies that the spin of a quark or anti-quark is changed only when a quark-anti-quark pair is created or annihilated, through an operator of the form $`𝝈𝐀`$ (Eqs. 5 and 8). When an initial heavy-quark meson $`Q\overline{Q}^{}`$ pair undergoes an OZI allowed decay to the two final heavy-quark meson pairs $`Q\overline{Q}^{\prime \prime }`$ and $`Q^{\prime \prime }\overline{Q}^{}`$, the spin is only changed when the $`Q^{\prime \prime }\overline{Q}^{\prime \prime }`$ pair is created333 According to Eqs. 57 and 8, only the time-like component of the gluon field couples to the propagating quark or anti-quark, while only the space-like component couples in the case of quark-anti-quark pair creation or annihilation. The time-like component can change into the space-like component via gluon-loop diagrams in both covariant and Coulomb gauge, and also via ghost-loop diagrams in covariant gauge.. Also, since the individual mesons are composed of very heavy quarks, moving non-relativistically, they have a specific quark-spin (assuming no accidental mixing with states nearby in mass). It follows that the spin selection rule is exact to all orders in QCD perturbation theory when the mesons participating in the decay are built from very heavy quarks and anti-quarks. Light quark loops do not change these conclusions. For very heavy quarks, $`1/m_Q`$ corrections are negligible compared to higher order corrections in $`\alpha _s`$, because $`\alpha _s(m_Q)`$ depends only logarithmically on $`m_Q`$.


text: # 1 40<𝑃_𝑡^𝛾<50𝐺𝑒𝑉/𝑐.


text:
# A terahertz molecular switch

## Abstract

We present time-dependent results describing the current through a molecular device, modeled as a complex with two active centers connected to leads under bias. We show that, at a properly adjusted external voltage, a passing terahertz electromagnetic pulse may cause a transition between states of finite and negligible current, suggesting that the system might be useful as a nanoscopic switch in the terahertz range. A phase diagram defining the bias region in which the transition takes place within a short time is given. As described, the physical processes involved are of an entirely different nature than those in ordinary photodetectors.

Individual molecules as well as molecular complexes are increasingly being perceived as possible electronic device elements, constituting great promise in the process of miniaturization. Some successes have already been reported along this linedekk ; joac ; park ; lian . Among recent examples are molecules that exhibit big on-off current ratios and a large negative differential resistance in two-terminal transportchen ; chen1 , behaving much in the way some mesoscopic semiconductor tunneling heterostructures doesak ; tsu . Such performance may be understood decomposing the system into coupled molecular subunits that transit collectively through a conducting quantum resonance as the bias along the device is modifiedhan ; datta ; hettler . The current-voltage profile is then determined by the properties of this resonance, which in turn may depend strongly on the dynamic accumulation of electronic charge in the various subunits, thus affecting the tunneling rate for transport.

When a pair of neighboring subunits of a different kind are involved, an energy level associated with one site must reach alignment with one in the other site in order for a device resonance to be established. The narrow voltage range in which significant current flows in the experimental data of Chen et alchen  suggests that such mechanism may be at work in some molecules. Using this two site model, Han et al found qualitative agreement between the model predictions and the experimental datahan . A time dependent approach for the analysis of the same model allowed us to fit the data to an impressive degree of accuracy, and showed in addition that switching between distinct current states is possible, as reported below. Switching is extremely fast compared to that in organic complexes, where the effect has been assigned to local structural disorder induced by the external electric fieldgao . While in the latter switching times are of the order of 0.1 $`\mu `$s, we here report characteristic times in the order of 0.1 ps, six orders of magnitude shorter. The device is then sensitive to terahertz pulses, suggesting that the system is potentially useful as a nanoscopic detector of electromagnetic radiation in the THz range.

We consider a molecular complex connected to leads subject to an external dc bias potential. As shown schematically in Fig. 1, the complex is modeled by two active sites, connected to left and right particle reservoirs at chemical potentials $`\mu _L`$ and $`\mu _R`$, respectively. The device is described by the Hamiltonian,

$`H`$ $`=`$ $`{\displaystyle \underset{i0,1;\sigma }{}}\epsilon _in_{i\sigma }+t{\displaystyle \underset{i(0,1);\sigma }{}}(c_{i,\sigma }^{}c_{i+1,\sigma }+c_{i+1,\sigma }c_{i,\sigma }^{})+{\displaystyle \underset{\alpha =0,1;\sigma }{}}\epsilon _\alpha n_{\alpha ,\sigma }+v_L{\displaystyle \underset{\sigma }{}}(c_{1,\sigma }^{}c_{0,\sigma }+c_{0,\sigma }^{}c_{1,\sigma })`$ (1)
$`+v_R{\displaystyle \underset{\sigma }{}}(c_{1,\sigma }^{}c_{2,\sigma }+c_{2,\sigma }^{}c_{1,\sigma })+g{\displaystyle \underset{\sigma }{}}(c_{0,\sigma }^{}c_{1,\sigma }+c_{1,\sigma }^{}c_{0,\sigma })+U{\displaystyle \underset{\alpha =0,1}{}}n_\alpha _{}n_\alpha _{}.`$
where $`c_{i\sigma }^{}`$ ($`c_{i\sigma }`$ ) creates (destroys) an electron with spin $`\sigma `$ at site i and $`n_{i,\sigma }=`$ $`c_{i\sigma }^{}`$ $`c_{i\sigma }`$ represents the occupation of site i. The molecular sites are labeled 0 and 1 and have intrinsic energies $`\epsilon _0^o`$and $`\epsilon _1^o`$, respectively, and an intra-molecular coupling $`g`$. The coupling strength of the molecule to the left (right) lead is $`v_{L(R)},`$ while t represents the hoping parameter within the leads. The total site energy $`\epsilon _i`$ includes the fixed intrinsic energy at the site, the external radiation-induced voltage and the applied dc bias, the latter represented by a classical term linear in the spatial coordinate $`i`$. Finally, the effect of charge build-up at the molecular sites is accounted for by a non-linear term with coupling constant U, the last one in the above expression, that leaves out occupation with the same quantum numbers. The time-dependent wave function that evolves according to this Hamiltonian can be expanded in a tight-binding basis as

$$|\psi _{k\sigma }(\tau )=\underset{i}{}a_{i\sigma }^k(\tau )|\varphi _{i\sigma },$$
(2)
where $`|\varphi _{i\sigma }`$ is a Wannier state of spin $`\sigma `$ localized at site $`i,`$ and the coefficients $`a_{i\sigma }^k(\tau )`$ obey the non-linear equations

$`i\mathrm{}{\displaystyle \frac{da_{j,\sigma }^k}{d\tau }}`$ $`=`$ $`\epsilon _j(\tau )a_{j,\sigma }^k+t(a_{j1,\sigma }^k+a_{j+1,\sigma }^k)(j1,0,1,2),`$
$`i\mathrm{}{\displaystyle \frac{da_{1(2),\sigma }^k}{d\tau }}`$ $`=`$ $`\epsilon _{1(2)}(\tau )a_{1(2),\sigma }^k+v_{L(R)}a_{0(1),\sigma }^k+ta_{2(3),\sigma }^k,`$
$`i\mathrm{}{\displaystyle \frac{da_{0,\sigma }^k}{d\tau }}`$ $`=`$ $`(\epsilon _0(\tau )+Un_0)a_{0,\sigma }^k+v_La_{1(1),\sigma }^k+ga_{1,\sigma }^k,`$
$`i\mathrm{}{\displaystyle \frac{da_{1,\sigma }^k}{d\tau }}`$ $`=`$ $`(\epsilon _1(\tau )+Un_1)a_{1,\sigma }^k+v_Ra_{2,\sigma }^k+ga_{0,\sigma }^k.`$ (3)
In this expression the average site densities equal $`n_0=_{k,\sigma }\left|a_{0,\sigma }^k\right|^2,n_1=_{k,\sigma }\left|a_{1,\sigma }^k\right|^2`$, where the sum over $`k,\sigma `$ covers all occupied electron states. In writing Eqs. (3) we have adopted a Hartree model for the electron-electron interaction. As we will show in what follows, the terms proportional to the site densities generated by this interaction, nonlinear in nature, are of key importance in the behavior of the system.

Solutions are found discretizing equations (3) and using a half-implicit numerical method, which is second-order accurate and unitarymains ; orellana . It is assumed that the coefficients outside the structure are given at incident wavenumber $`k`$ by the expressions

$$a_{j,\sigma }^k=(Ie^{ikja}+R_je^{ikja})e^{i\epsilon \tau /\mathrm{}}(jaL)$$
(4)
$$a_{j,\sigma }^k=T_je^{ik^{}ja}e^{i\epsilon \tau /\mathrm{}}(jaL),$$
(5)
where the device is explicitly defined to lie between sites $`L`$ and $`L,`$ $`a`$ is the lattice constant, and $`k^{}`$ $`=[2m(\epsilon V)/\mathrm{}^2]^{1/2}`$ is the wavenumber of the transmitted particle, with $`V=\mu _L\mu _R`$. The incident amplitude $`I`$ is supposed to be spatially constant. Also, far from the molecule the amplitude of the reflected and transmitted waves $`R_j`$ and $`T_j`$ are supposed to be weakly dependent on site $`j`$. This permits to restrict the dependence to the linear term, which is found to be an adequate approximation provided the time step taken to discretize equations (3) is less than certain limit value that depends upon the parametrization of the system. In the numerical procedure the Wannier amplitudes obtained for one bias are used as the starting point for the next bias step. Once these coefficients are known the current is calculated from orellana ,

$$J_j=\frac{e}{h}t\underset{\sigma }{}_o^{k_f}Im(a_{j,\sigma }^ka_{j+1,\sigma }^k)𝑑k$$
(6)
We apply the above formalism to a specific situation. An excellent fit to the experimental I-V curve at T=60K of Chen et al.chen  is obtained in the absence of radiation within this model, by choosing the parameters $`\epsilon _0^o=0.52,\epsilon _1^o=0.92`$, $`g=0.008,t=2.0,U=0.8`$ and intrinsic level broadenings $`\mathrm{\Gamma }_L=\mathrm{\Gamma }_R=0.16`$, all in units of eV. Here $`\mathrm{\Gamma }_{L(R)}=\pi v_{L(R)}^2\rho (0)`$, with $`\rho (0)`$ the density of states at the Fermi energy. The energies at the active sites are represented by $`\epsilon _i`$ $`=\epsilon _i^o[V+\delta E\mathrm{sin}(2\pi \nu \tau )]/2,`$ $`(`$ $`i=0,1),`$ where $`\delta E`$ and $`\nu `$ are the amplitude and frequency of the incoming THz radiation, respectively. Energies are measured with respect to the Fermi energy at the emitter contact ($`\mu _L`$). Figure 2 shows our results, together with experimental data points from Ref. 5 shown as circles. Notice that the current peak is quite narrow, and there is a steep fall in the upper critical voltage edge $`V_c`$=2.16 V. As the bias is increased and the first site energy approaches the Fermi energy charge begins to flow, first into site 0, then to site 1 as well as shown in the inset and, schematically, in Fig. 1b. The effect of this charge is to raise the local potential at each subunit, thus affecting the transport coefficient of the overall quantum mechanical structure. The abrupt fall at $`V_c`$ is caused by a sudden loss of charge at site 1, and subsequent break of the energy level alignment that allows for resonant tunneling through the molecule.

Figure 3 shows the time evolution of the current flowing through the system (Fig. 3a) and the charge density at sites 0 (full line) and 1 (broken line) (Fig. 3b) when electromagnetic radiation of frequency $`\nu =10`$ THz and strength $`\delta E=1.8`$ mV is turned on at time =0. The offset bias $`\delta V=V_cV`$ has the value 0.2mV. Before the radiation enters the sample the system is in resonance, current flows and both sites have significant charge. At $`\tau `$=0 the external field is turned on and after a time of the order of a picosecond, site 1 looses its charge and the current drops to a low value. This behavior illustrates the fast switching response of the device as a THz pulse passes by. The actual transition time is controlled by the internal coupling constant, so that in our case $`T\mathrm{}/g0.08ps`$. There are small oscillations driven by the external field, which we have found to become larger in amplitude when the signal strength is increased. The fact that the system undergoes a few coherent oscillations before falling suggests that the random fluctuations caused by a finite temperature are not capable of overcoming the inertia of the device. Notice that in the low current state site 1 has lost its charge. This feature reveals that the system is bi-stable very close to $`V_c`$. When site 1 is charged, the local potential is dynamically adjusted to keep the system in resonance. A disturbance may cause the flow of charge out of the site, however, with the result that the resonance condition is lost and current stops.

Suppose now that the system is in a given state of high current, biased slightly below $`V_c`$. The transition to a low current state is then possible provided the amplitude of the external field is large enough to bring the system to criticality. Figure 4 shows the region in which the transition happens within 2 picoseconds time (labeled YES), separated from that region in which the radiation is too weak to cause the current to fall that fast, if at all (labeled NO). The frequency is $`\nu =10`$ THz and the bias is shown in terms of the offset from criticality $`\delta V=V_cV`$. A rough fit to the boundary is given by $`\delta E=0.014\times \delta V^{1/2}`$, all quantities in units of V. Our numerical methods allowed testing frequencies up to order 1/T, and within that range we found the boundary to shift to higher values as the frequency increases. In the figure, the filled circle is a point in the boundary for $`\nu =20`$THz while the triangle is for $`\nu =5`$ THz. This is to be expected since in the low frequency limit one should have $`\delta E\delta V.`$

In summary, we have shown that switching between different current states in transport through molecules is possible and may be triggered by a passing electromagnetic pulse with characteristic times in the terahertz range. Depending on the applied external bias, the passage of current is turned off by the radiation pulse, thus allowing the device to act as a tunable detector of terahertz radiation. Switching occurs because of the existence of bi-stable states near the I-V characteristic high-voltage edge. A similar bi-stable region is known to occur in transport through GaAlAs double barrier resonant tunneling devicesorellana ; gold ; mart ; inko . A single active energy level localized in the quantum well formed between the barriers is responsible for the bi-stability in such case. Near the bottom of the emitter conduction band a current carrying state through this level, lifted in energy by the charge trapped in the well, may coexist with a state of no current at the same bias with the well uncharged. In a two-site molecular device as the one treated here the effect derives from the dynamic coupling between the sites within the molecule, and thus has a different physical origin.

We thank Dr. J.E. Han for valuable discussions and bringing the work of Chen et al to our attention, and Dr. M. Orszag for comments on the manuscript. This work was supported in part by FONDECYT grants No. 1020269 and No. 1020829, Programa Milenio ICM P 99-135-F(PO) and Cátedra Presidencial en Ciencias (FC).


text:
# Supersymmetric quenched complexity in the Sherrington-Kirkpatrick model

## Abstract

By using the BRST supersymmetry we compute the quenched complexity of the TAP states in the SK model. We prove that the BRST complexity is equal to the Legendre transform of the static free energy with respect to the largest replica symmetry breaking point of its overlap matrix.

A key issue in the physics of complex systems is the computation of the entropy of the metastable states, normally called complexity in spin-glasses, and configurational entropy in structural glasses and supercooled liquids. A knowledge of the complexity is crucial for understanding the dynamics of a system when this is heavily influenced by strong metastability effects. Moreover, in some theoretical frameworks, the drop in the number of accessible states leads to an ergodicity breaking transition. In this context the complexity is essential also from the thermodynamic point of view, as in the Adam-Gibbs theory of the thermodynamic glass transition adam .

Despite its enormous theoretical relevance, there are few analytic calculations of the complexity in glassy systems, and this for a very good reason. In a nutshell, to find the complexity we have to compute the number of local minima (metastable states) of some state function, which is typically highly nontrivial. Just to fix ideas, we may think that this function is the Hamiltonian $`H`$. To compute the complexity, we must impose that the gradient of $`H`$ (the force) vanishes in the local minima, and we have to include as a normalization factor the second derivative of $`H`$ (the Hessian). Moreover, we may want to classify the metastable states according to the value of $`H`$, that is to their height in the landscape. Therefore, beside the force and the Hessian, we must include the state function itself in the calculation. Computing the complexity is thus a formidable technical task, since we have to deal with three very complicated functions: $`H`$, $`H`$, and $`^2H`$. In comparison, the calculation of the partition function, which just involves $`H`$, is an easy business.

This apparent difficulty in the calculation of the complexity stems from the fact that most methods treat $`H,H`$ and $`^2H`$ as three independent objects, when of course they are not. Every calculation which fails to capture the fact that it is essentially just one function, $`H`$, that we are dealing with, effectively wastes a crucial information. It would be therefore important to find a tool which exploits this information to simplify the calculation of the complexity. The Becchi-Rouet-Stora-Tyutin (BRST) supersymmetry brs ; tito  seems to be such a tool. As it was first noted in juanpe  for a particular model, a BRST calculation of the complexity is in fact equivalent to the one of the partition function. This is indeed what we expect from a method which does not treat $`H,H`$ and $`^2H`$ as independent functions.

The formal equivalence between complexity and standard thermodynamics found in juanpe  by means of the BRST supersymmetry is a very important theoretical issue. In the context of spin-glasses the existence of such a connection has been much investigated in the past bm1 ; bm2 ; bm3 ; ddy ; bmy ; monasson ; franzparisi ; potters ; crisatap . In a classic paper bm1 , Bray and Moore first calculated the complexity of the Sherrington-Kirkpatrick (SK) model sk  by counting the number of local minima of the Thouless-Anderson-Palmer (TAP) free energy tap , which in mean-field spin-glasses is the state function discussed above. The same authors also noted in bm2  some deep formal connections between TAP complexity and standard thermodynamics, while De Dominicis and Young showed in ddy  that TAP and static approaches were in fact equivalent, once some key hypothesis were made. These studies culminated in a remarkable work bmy , where Bray, Moore and Young uncovered a sort of Legendre transform relationship between TAP complexity and static free energy.

A method to compute the complexity which does not rely on the existence of a TAP free energy, was introduced by Monasson monasson , and by Franz and Parisi franzparisi . The basic idea is to introduce a coupling between different systems forcing them to live in the same metastable state. The free energy cost of such a constrained super-system is equal to the entropic contribution of the metastable states, that is the complexity. Within this approach, close connections between complexity and thermodynamics, similar to those found in the TAP context in bmy , were found. In particular, in spin-glass models with one step of replica symmetry breaking 1rsb ; mpv , the formulation of Monasson shows that the complexity is equal to the Legendre transform of the static free energy with respect to the breaking point $`x`$ of the overlap matrix monasson ; crisatap .

Despite all these investigations, it is fair to say that a general formal connection between complexity of the metastable states and static free energy, has not been proved yet. In particular, it is unclear how the Legendre transform method of monasson  should be used in systems with more than one step of replica symmetry breaking, as the SK model. In fact, none of the previous SK investigations bm2 ; bmy  succeeded in proving the existence of a sharp Legendre transform relationship as in one-RSB systems.

In this Letter we find for the first time an exact connection between complexity of the metastable states and static free energy in the SK model: we prove that the quenched TAP complexity obtained by means of the BRST supersymmetry is the Legendre transform of the static free energy with respect to the largest breaking point of its overlap matrix. Our result confirms the validity of the Legendre transform method of monasson ; franzparisi ; potters , and its consistency with the investigations of bm2 ; bm3 ; ddy ; bmy . Moreover, our findings strongly suggest that the BRST supersymmetry should be considered as an essential tool for the computation of the complexity in more general glassy systems.

The complexity of the TAP states with free energy density $`f`$, at inverse temperature $`\beta `$, is defined as bm1 ,

$$\mathrm{\Sigma }(\beta ,f)=\frac{1}{N}\mathrm{log}\underset{\alpha =1}{\overset{𝒩}{}}\delta \left[N\beta f\beta F_{TAP}(m^\alpha )\right]=\frac{1}{N}\mathrm{log}𝑑re^{Nr\beta f}\underset{\alpha =1}{\overset{𝒩}{}}e^{r\beta F_{TAP}(m^\alpha )},$$
(1)
where $`m^\alpha \{m_i^\alpha \}`$, are the local magnetizations at site $`i=1\mathrm{}N`$, in state $`\alpha =1\mathrm{}𝒩`$. A state $`m^\alpha `$ is defined as a local minimum of the TAP free energy $`F_{TAP}(m)`$ tap . If we define the thermodynamic potential $`\mathrm{\Psi }(\beta ,r)`$,

$$\mathrm{exp}\left(\beta Nr\mathrm{\Psi }\right)\underset{\alpha =1}{\overset{𝒩}{}}e^{r\beta F_{TAP}(m^\alpha )},$$
(2)
we can use the steepest descent method in (1), and obtain the complexity as the Legendre transform of $`\mathrm{\Psi }(\beta ,r)`$,

$$\mathrm{\Sigma }(\beta ,f)=\beta rf\beta r\mathrm{\Psi }(\beta ,r),$$
(3)
where the parameter $`r=r(\beta ,f)`$ is fixed by the equation,

$$\mathrm{\Psi }(\beta ,r)+r\frac{\mathrm{\Psi }(\beta ,r)}{r}=f.$$
(4)
From (2) we see that for $`r=1`$ the potential $`\mathrm{\Psi }`$ must be equal to the standard static free energy of the system $`F(\beta )`$, calculated in the TAP context. This calculation was first performed in ddy , where it was shown that the relation $`\mathrm{\Psi }(\beta ,r=1)=F(\beta )`$ only held if some suitable assumptions were made. In brst1  it was proved that the assumptions used in ddy  were in fact a general consequence of the BRST supersymmetry. However, the situation was less clear for generic values of $`r`$, since the calculations of bm3 ; bmy  for $`r1`$ seemed to explicitly break the BRST invariance nota . In what follows we perform a supersymmetric quenched calculation of $`\mathrm{\Psi }(\beta ,r)`$, and prove that this potential is intimately related to the static free energy $`F(\beta )`$ even for $`r1`$. The TAP free energy for the SK model is given by tap ,

$$\beta F_{TAP}(m)=\frac{\beta }{2}\underset{ij}{}J_{ij}m_im_j+\frac{1}{\beta }\underset{i}{}\varphi _0(m_i),$$
(5)
with,

$$\varphi _0(m)=\frac{1}{2}\mathrm{log}(1m^2)+m\mathrm{tanh}^1(m)\mathrm{log}2\frac{\beta ^2}{4}(1q)^2.$$
(6)
The variable $`q`$ is the self-overlap of the TAP states, $`q=\frac{1}{N}_im_i^2`$. The quenched couplings $`J`$ are random variables with Gaussian distribution and variance $`N`$. From (2) we have that the quenched potential $`\mathrm{\Psi }(\beta ,r)`$ is,

$$\beta r\mathrm{\Psi }(\beta ,r)=\frac{1}{nN}\mathrm{log}\overline{\rho (\beta ,r|J)^n},$$
(7)
with,

$$\rho (\beta ,r|J)=\underset{\alpha =1}{\overset{𝒩}{}}e^{r\beta F_{TAP}(m^\alpha )}=\underset{i}{}dm_i\delta (_iF_{TAP}(m))|det(_i_jF_{TAP}(m))|e^{\beta rF_{TAP}(m)}.$$
(8)
In (7) we have $`N\mathrm{}`$ and $`n0`$, and the over-bar indicates an average over the disorder. As usual, the modulus of the determinant will be dropped. This amounts to assume that at sufficiently low temperatures the largest part of TAP solutions are minima brst1 . Of course, any method which drops the modulus is doomed to fail if stable minima are subdominant with respect to unstable saddles. After introducing the commuting fields $`x_i`$ to implement the delta functions, and the anti-commuting (Grassmann) fields $`\overline{\psi }_i,\psi _i`$ for the determinant, we find,

$$\rho (\beta ,r|J)=𝒟m𝒟x𝒟\overline{\psi }𝒟\psi e^{\beta S(m,x,\overline{\psi },\psi )},$$
(9)
where the action $`S`$ is given by,

$$S(m,x,\overline{\psi },\psi )=\underset{i}{}x_i_iF_{TAP}(m)+\underset{ij}{}\overline{\psi }_i\psi _j_i_jF_{TAP}(m)rF_{TAP}(m).$$
(10)
By averaging $`\rho (\beta ,r|J)^n`$ over the disorder we obtain the following effective action,

$`\beta S`$ $`=`$ $`{\displaystyle \frac{\beta ^2}{2N}}\left[{\displaystyle \underset{ab}{\overset{n}{}}}\left({\displaystyle \underset{i}{\overset{N}{}}}x_i^ax_i^b\right)\left({\displaystyle \underset{j}{\overset{N}{}}}m_j^am_j^b\right)+{\displaystyle \underset{ab}{\overset{n}{}}}\left({\displaystyle \underset{i}{\overset{N}{}}}x_i^am_i^b\right)^2{\displaystyle \underset{ab}{\overset{n}{}}}\left({\displaystyle \underset{i}{\overset{N}{}}}\overline{\psi }_i^a\psi _i^b\right)^2\right]`$ (11)
$`+`$ $`{\displaystyle \frac{\beta ^2}{2N}}\left[{\displaystyle \frac{r^2}{2}}{\displaystyle \underset{ab}{\overset{n}{}}}\left({\displaystyle \underset{i}{\overset{N}{}}}m_i^am_i^b\right)^22r{\displaystyle \underset{ab}{\overset{n}{}}}\left({\displaystyle \underset{i}{\overset{N}{}}}m_i^ax_i^b\right)\left({\displaystyle \underset{j}{\overset{N}{}}}m_j^am_j^b\right)\right]`$
$`+`$ $`{\displaystyle \underset{a}{\overset{n}{}}}{\displaystyle \underset{i}{\overset{N}{}}}\left[x_i^a\varphi _1(m_i^a)+\overline{\psi }_i^a\psi _i^b\varphi _2(m_i^a)r\varphi _0(m_i^a)\right],`$
where $`\varphi _1=_m\varphi _0`$ and $`\varphi _2=_m^2\varphi _0`$. The scalar overlap $`q`$ will now be generalized by introducing the overlap matrix $`q_{ab}=m^am^b`$. This form of the action is different from the one in bm3 ; bmy , where the delta function enforcing the TAP equations was used to eliminate the factor $`J_{ij}m_i^am_j^a`$ in (10). This is a crucial point: this method explicitly breaks the BRST invariance of the action, which is the crucial tool to establish the exact connection with the static free energy. We therefore do not use this method, and keep the whole BRST invariant action nota .

In order to linearize the quadratic terms we introduce the usual Lagrange multipliers, $`m^am^b\lambda _{ab}`$, $`m^ax^bw_{ab}`$, $`\overline{\psi }^a\psi ^bt_{ab}`$ brst1 . After this is done, the integrals in $`x`$ and $`\overline{\psi },\psi `$ become Gaussian and can be performed explicitly. Moreover, the action factorizes and for $`N\mathrm{}`$ we can use the steepest descent method, to get,

$$\beta r\mathrm{\Psi }(\beta ,r)=\underset{n0}{lim}\frac{1}{n}\left[\mathrm{\Sigma }_0+\mathrm{log}\underset{a}{}dm^ae^{(m^a)}\right].$$
(12)
Following bm1 ; bm3 ; brst1  we define,

$`B_{ab}`$ $``$ $`\beta ^2(1q_{aa})\delta _{ab}+t^{ab}`$ (13)
$`\mathrm{\Delta }_{ab}`$ $``$ $`\beta ^2(1q_{aa})\delta _{ab}w^{ab},`$ (14)
and therefore obtain (details will be given elsewhere),

$`\mathrm{\Sigma }_0={\displaystyle \frac{1}{2\beta ^2}}{\displaystyle \underset{ab}{}}(B_{ab}^2\mathrm{\Delta }_{ab}^2){\displaystyle \underset{a}{}}\left(B_{aa}+\mathrm{\Delta }_{aa}\right)(1q_{aa})`$
$`{\displaystyle \underset{ab}{}}\left[{\displaystyle \frac{\beta ^2}{4}}r^2q_{ab}^2+\lambda ^{ab}q_{ab}r\mathrm{\Delta }^{ab}q_{ab}\right]{\displaystyle \frac{1}{2}}\mathrm{log}[(2\pi \beta ^2)^ndetq_{ab}]`$ (15)
$`(m^a)=r{\displaystyle \underset{a}{}}\varphi _0(q_{aa},m^a)+{\displaystyle \underset{ab}{}}\lambda ^{ab}m^am^b+\mathrm{log}det\left({\displaystyle \frac{\delta _{ab}}{1m_a^2}}+B_{ab}\right)`$
$`{\displaystyle \frac{1}{2\beta ^2}}{\displaystyle \underset{ab}{}}\left[\mathrm{tanh}^1m^a{\displaystyle \underset{c}{}}\mathrm{\Delta }^{ac}m^c\right]q_{ab}^1\left[\mathrm{tanh}^1m^b{\displaystyle \underset{c}{}}\mathrm{\Delta }^{bc}m^c\right].`$ (16)
The parameters $`\mathrm{\Delta }_{ab},B_{ab},\lambda _{ab}`$ and $`q_{ab}`$ must be fixed by the saddle point equations, and it is easy to show that $`B_{ab}=0`$ is solution. It is important at this point to consider the role of the supersymmetry. In juanpe  it was noted that action (10) is invariant under a generalization of the BRST supersymmetry brs ; tito ,

$$\delta m_i=ϵ\psi _i\delta x_i=ϵr\psi _i\delta \overline{\psi }_i=ϵx_i\delta \psi _i=0.$$

where $`ϵ`$ is an infinitesimal Grassmann parameter. If we calculate the variation of $`m_i\overline{\psi }_i`$ and $`x_i\overline{\psi }_i`$ juanpe ; brst1 , we obtain the two BRST equations,

$`\overline{\psi }_i\psi _i`$ $`+`$ $`m_ix_i=0`$ (17)
$`r\overline{\psi }_i\psi _i`$ $`+`$ $`x_ix_i=0.`$ (18)
After some algebra, these equations become,

$`\mathrm{\Delta }_{ab}`$ $`=`$ $`\beta ^2q_{ab}r`$ (19)
$`\lambda _{ab}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\beta ^2r^2q_{ab},`$ (20)
and it is possible to show that the remaining saddle point equations are indeed satisfied by (19) and (20). The only saddle point equation we are left with is obtained by doing the variations of (15) and (16) with respect to $`\lambda _{ab}`$. This gives,

$$q_{ab}=m^am^b,$$
(21)
where the average is performed with the distribution $`\mathrm{exp}((m^a))`$. If we use equations (19), (20) and (21) into (15) and (16), and make the change of variable $`m^ah^a=\mathrm{tanh}^1(m^a)`$, we finally obtain,

$$\beta \mathrm{\Psi }(\beta ,r)=\mathrm{log}2+\frac{\beta ^2}{4n}[r\underset{ab}{\overset{n}{}}q_{ab}^2\underset{a}{\overset{n}{}}(1q_{aa})^2]\frac{1}{nr}\mathrm{log}\underset{a}{\overset{n}{}}\frac{dh^a}{\sqrt{2\pi \beta ^2detq_{ab}}}\mathrm{cosh}(h^a)^re^{\frac{1}{2\beta ^2}_{ab}^nh^aq_{ab}^1h^b}.$$
(22)
Expression (22) is different from the one of bm3 ; bmy : the BRST supersymmetry automatically selects one saddle point in the space of parameters, and in so doing it drastically reduces the number of parameters, compared to bm3 ; bmy . The computation of $`\mathrm{\Psi }(\beta ,r)`$ has at this point the same degree of difficulty as the one of the standard free energy $`F(\beta )`$, with just one overlap matrix $`q_{ab}`$ to be fixed variationally. We shall now show that the connections between $`\mathrm{\Psi }(\beta ,r)`$ and $`F(\beta )`$ are in fact much deeper than that.

The general form of the quenched free energy in the SK model is sk ,

$$\beta F(\beta )=\frac{\beta ^2}{4}+\frac{\beta ^2}{2n_s}\underset{\alpha >\beta }{\overset{n_s}{}}Q_{\alpha \beta }^2\frac{1}{n_s}\mathrm{log}\underset{[\sigma ^\alpha ]}{}\mathrm{exp}\left[\frac{\beta ^2}{2}\underset{\alpha \beta }{\overset{n_s}{}}Q_{\alpha \beta }\sigma ^\alpha \sigma ^\beta \right],$$
(23)
where $`Q_{\alpha \beta }`$ is the $`n_s\times n_s`$ overlap matrix, with $`n_s0`$. If in (22) we use the relation,

$$\mathrm{cosh}(h_a)^r=\frac{1}{2^r}\underset{[\tau _a^\mu =\pm 1]}{}e^{h_a_\mu ^r\tau _a^\mu },$$
(24)
we obtain,

$$\beta \mathrm{\Psi }(\beta ,r)=\frac{\beta ^2}{4}+\frac{\beta ^2}{2n}\left[r\underset{a>b}{\overset{n}{}}q_{ab}^2+\frac{r1}{2}\underset{a}{\overset{n}{}}q_{aa}^2\right]\frac{1}{nr}\mathrm{log}\underset{[\tau _a^\mu ]}{}\mathrm{exp}\left[\frac{\beta ^2}{2}\left(\underset{ab}{\overset{n}{}}\underset{\mu \nu }{\overset{r}{}}\tau _a^\mu q_{ab}\tau _b^\nu \underset{a}{\overset{n}{}}\underset{\mu }{\overset{r}{}}q_{aa}\right)\right].$$
(25)
The trace terms in equations (23) and (25) suggest the relation $`n_s=rn`$. Once this identification is done, we can connect the $`\sigma ^\alpha `$ spin variables ($`\alpha =1,\mathrm{},n_s`$), to the $`\tau _a^\mu `$ spin variables ($`a=1,\mathrm{},n;\mu =1,\mathrm{},r`$) in the following way,

$$(\sigma _1,\mathrm{},\sigma _{n_s})=(\tau _1^1,\mathrm{},\tau _1^r,\mathrm{}\mathrm{},\tau _n^1,\mathrm{},\tau _n^r).$$

Let us now assume that the potential $`\mathrm{\Psi }(\beta ,r)`$ (and thus the TAP complexity) is calculated at $`k`$ levels of replica symmetry breaking (RSB) mpv . The TAP overlap matrix $`q_{ab}`$ is then given by,

$$q_{ab}^{(k)}=q_0+\underset{i=1}{\overset{k+1}{}}(q_iq_{i1})\epsilon _{ab}^{(n,y_i)},$$
(26)
with $`y_{k+1}=1`$ and $`\epsilon _{ab}^{(n,1)}=\delta _{ab}`$. The matrices $`\epsilon ^{(n,y_i)}`$ are $`n\times n`$ ultrametric block matrices, equal to one on the diagonal blocks of size $`y_i`$, and zero elsewhere. The variables $`y_i`$ are thus the replica symmetry breaking points. In the TAP approach the diagonal of the overlap matrix, $`q_{aa}=q_{k+1}`$, contains the self-overlap of the states, and for this reason $`y_{k+1}=1`$. There are $`k+1`$ values of the overlap, but only $`k`$ nontrivial breaking points, and thus $`q_{ab}`$ is a $`k`$-RSB matrix. Given this form of $`q_{ab}`$, it is possible to prove (details will be given elsewhere) that,

$$\underset{ab}{\overset{n}{}}\underset{\mu \nu }{\overset{r}{}}\tau _a^\mu q_{ab}^{(k)}\tau _b^\nu \underset{a}{\overset{n}{}}\underset{\mu }{\overset{r}{}}q_{aa}^{(k)}=\underset{\alpha \beta }{\overset{rn}{}}Q_{\alpha \beta }^{(k+1)}\sigma _\alpha \sigma _\beta ,$$
(27)
where $`Q_{\alpha \beta }^{(k+1)}`$ is a standard $`rn\times rn`$ RSB matrix, with $`k+1`$ levels of replica symmetry breaking. More precisely,

$$Q_{\alpha \beta }^{(k+1)}=q_0+\underset{i=1}{\overset{k+1}{}}(q_iq_{i1})\epsilon _{\alpha \beta }^{(rn,ry_i)}.$$
(28)
From this formula we see that the entries of $`Q_{\alpha \beta }^{(k+1)}`$ are the same as $`q_{ab}^{(k)}`$, whereas the $`k+1`$ replica symmetry breaking points $`x_i`$ of $`Q_{\alpha \beta }^{(k+1)}`$ are rescaled by a factor $`r`$, that is $`x_i=ry_i`$. In particular, the largest breaking point $`xx_{k+1}`$ of the static matrix $`Q_{\alpha \beta }^{(k+1)}`$ is given by,

$$x=r.$$
(29)
By inserting relation (27) into (25), we finally obtain,

$$\mathrm{\Psi }(\beta ,r|q_{ab}^{(k)})=F(\beta |Q_{ab}^{(k+1)}),$$
(30)
We have thus proved that the thermodynamic potential $`\mathrm{\Psi }(\beta ,r)`$ calculated at the $`k`$ RSB level is equal to static free energy $`F(\beta )`$ calculated at the $`k+1`$ RSB level. The replica symmetry breaking points of the static matrix $`Q_{\alpha \beta }`$ are simply the ones of the TAP matrix $`q_{ab}`$ rescaled by the parameter $`r`$, and the extra $`k+1`$-th breaking point of $`Q_{\alpha \beta }`$ is equal to $`r`$. This rescaling was first noted in bmy , and later in potters , although the lack of BRST symmetry of those calculations prevented to prove equation (30).

From equation (3), and given the relation between $`\mathrm{\Psi }(\beta ,r)`$ and $`F(\beta )`$, we finally have the general Legendre equation connecting the quenched complexity of the TAP states to the standard static free energy in the SK model,

$$\mathrm{\Sigma }(\beta ,f)=\beta xf\beta xF(\beta ;x),$$
(31)
with the largest breaking point $`x`$ fixed by the equation,

$$f=F(\beta ;x)+x\frac{F(\beta ;x)}{x}.$$
(32)
This result can be summarized as follows: the supersymmetric quenched complexity of the TAP states is the Legendre transform of the static free energy with respect to the largest breaking point $`x`$ of its overlap matrix. This result allows to compute the exact quenched complexity in the SK model, that is the complexity at the full\- RSB level, as the Legendre transform of the full-RSB static free energy leuzzi .

It has been conjectured in riko  that in systems with more than one step of replica symmetry breaking the complexity of clusters at level $`i`$ is given by the Legendre transform of the free energy with respect to the breaking point $`x_i`$. For $`x_i=x_{max}`$ clusters are just states, and our result is recovered. It would be interesting to study whether the conjecture of riko  can be exactly proved within the supersymmetric formalism used here.

We thank A. Crisanti, L. Leuzzi, R. Monasson, A. Montanari, F. Ricci-Tersenghi, and T. Rizzo for some interesting discussions.