text:
# The Multi-Component Nature of the Vela Pulsar Nonthermal X-ray Spectrum

## 1 Introduction

The Vela pulsar (PSR B0833-45) is the strongest $`\gamma `$-ray source in the sky, but it is one of the most difficult pulsars to detect at X-ray energies. This is in part because it is embedded in a very bright X-ray synchrotron nebula providing a large unpulsed background, but also because its pulsed X-ray emission is comparatively weak. Its spectral power peaks in the $`\gamma `$-ray band, a characteristic of most other pulsars having detected emission at high energies, but it falls off in the hard X-ray band.

The Vela pulsar is well-known at high $`\gamma `$-ray energies. Its light curve as measured by CGRO/EGRET (Kanbach et al. 1994) consists of two relatively narrow peaks separated by 0.43 in phase. The first peak lags the radio peak phase by 0.12. The region between the two $`\gamma `$-ray peaks contains a low level of emission. Observations using COMPTEL and OSSE on CGRO (Strickman et al. 1996 and references therein) show that, while the light curve is similar in the low-energy $`\gamma `$-ray band to that at higher energies, the spectrum is considerably harder. The details of this rollover and comparison of light curve behavior are still sketchy due to the relatively poor statistics of the low energy $`\gamma `$-ray results. Since various models of pulsar high-energy emission make predictions concerning the shape and location of this rollover, and the behavior of the spectrum below the rollover, better measurement of pulsed emission in the X-ray band is important.

The first detection of pulsed emission at X-ray energies was made by ROSAT in the 0.1 - 2 keV band (Ogelman et al. 1993). The observed spectrum was consistent with a blackbody. Given the nonthermal nature of pulsed emission at other energies, this component may come from a distinct emission mechanism such as thermal emission from the hot polar cap or neutron star surface. The first detection of a nonthermal X-ray pulsed component, made with our first RXTE observation, is reviewed in the next section.

In this paper, we report the results of our analysis of a 300 ks RXTE Cycle 3 observation. With significantly improved statistics, we are able to clearly confirm the existence of at least two non-thermal spectral components and present phase-resolved spectra.

## 2 RXTE Cycle 1 Results

We detected the Vela pulsar for the first time in the 2 - 30 keV band during a 93 ks RXTE Cycle 1 observation (Strickman, Harding & De Jager 1999 \[SHD99\]) during January 1997. We found pulsed emission with a high degree of significance in each of three broad energy bands: 2 - 8 keV, 8-16 keV and 16-30 keV. The light curve in each of these energy bands showed two peaks separated by about 0.4 in phase, very similar to the light curves seen by CGRO instruments. The two peaks, however, exhibited quite different characteristics. The first peak was very narrow and was phase-aligned with the EGRET $`\gamma `$-ray peak. The second peak was much broader and appeared to shift in phase (away from the first peak) from lower to higher energies. In the 16-30 keV band, this peak was near the phase of the EGRET second peak, but in the 2-8 keV band the peak was closer to the first RXTE peak and near the phase of the second optical peak. Due to limited statistics, we were able to determine the spectrum only of the two main peaks (Peaks 1 and 2) and what was called the “Peak 1 precursor”, a broad region of enhanced emission preceding Peak 1. We found that the spectrum of Peak 1 (photon index $`0.68\pm 0.14`$) was significantly harder than Peak 2 (photon index $`1.17\pm 0.12`$), and both peaks were harder than the Peak 1 precursor (photon index $`1.57\pm 0.29`$).

From these results, we speculated that Peak 2 was a blend of two separate components: a hard component whose spectrum and light curve connect smoothly with the OSSE, COMPTEL and EGRET gamma-ray spectrum and light curves, and a soft component, whose spectrum and light curve appear to be associated more with the optical emission. We also noted that the Peak 1 precursor was at the phase of the radio peak and had the softest spectrum of any RXTE phase components.

## 3 RXTE Cycle 3 Observations and Analysis

The Cycle 3 observations were carried out during April/May and July/August 1998 with a good exposure of 274 ks. Data were collected in GoodXenon event-by-event mode for the entire observation. We epoch-folded the data at the pulsar period using the standard fasebin suite of ftools supplied by the RXTE Science Support Center. The pulsar ephemeris was derived from the same dataset used by the Princeton Pulsar Database (Arzoumanian et al. 1992), but recalculated by D. Nice (private communication, 1999) so as to produce a single ephemeris valid over the entire time span of our observations. We summed the resulting energy dependent light curves into three broad energy bands (2-8 keV, 8-16 keV and 16-30 keV) to study broad energy dependence of light curve features.

In the analysis of our Cycle 1 observation, we used a simple “on-pulse” minus “off-pulse” technique to extract phase resolved spectra. This analysis is inappropriate here since we want to separate spectra from only partially resolved features in the light curve. Instead we have created a model of the light curve shape based on the broad band light curves. The model consists of a constant “background” level with a set of five sinusoidal “peaks” superimposed. The peaks are each of form

$`C(i)`$ $`=`$ $`A(i)|\mathrm{cos}[{\displaystyle \frac{\pi }{2}}(\varphi \varphi __0(i))]|^{\xi (i)},`$ (1)
$`\xi (i)`$ $`=`$ $`{\displaystyle \frac{0.693}{\mathrm{ln}[\mathrm{cos}(\frac{\pi }{2}W(i))]}}`$
where $`A(i)`$, $`\varphi _0(i)`$ and $`W(i)`$ are the amplitude, center phase and full width at half maximum of peak $`i`$ respectively.

We fit this model to the light curve data in 94 native PCA energy channels roughly spanning the 2-30 keV band. Preliminary fits to the broad band light curves indicated that, in most cases, $`\varphi _0(i)`$ and $`W(i)`$ did not vary much with energy. Hence the values of these parameters were either fixed or computed from linear functions of energy derived from the preliminary fits. Hence, the only free parameter in the individual channel light curve fits was $`A(i)`$. For each energy channel, the counts for each peak were determined by integrating the appropriate best fit sinusoid over phase. The resulting spectra were placed in PHA files, where they could be accessed and fit with XSPEC as described in the following section.

## 4 Results

Figure 1 shows broad-band light curves including the RXTE light curves in three energy bands together with those of EGRET, OSSE, Chandra and Optical observations. The RXTE light curves show a narrow peak (Peak 1) at the phase of the EGRET first peak and with a phase width of $`0.035\pm 0.004`$, and a second peak that is now clearly seen to be a blend of two components. The first of these components, which we call Peak 2-soft, is closer to Peak 1 and is the less intense of the two components. It is more distinctly visible in the light curve at low energies. The second component, Peak 2-hard, is further from Peak 1 and is significantly broader than Peak 2-soft, with mean widths in phase of $`0.11\pm 0.02`$ and $`0.04\pm 0.01`$ respectively. There are two other statistically significant peaks that appear in the RXTE light curve: Peak 4 at the radio phase (0.0) with width $`0.03\pm 0.01`$ and a weaker peak, Peak 3, leading the radio peak with width $`0.10\pm 0.04`$. These two peaks made up the broad enhancement which we denoted the Peak 1 precursor in our Cycle 1 data. The various peak center phases do not appear to vary with energy except for Peak 2-hard, which is not consistent with a constant. It varies by 0.08 in phase from the low end to the high end of the RXTE band. This may in part be an artifact of the sinusoid fitting process. In all cases, the peak widths are consistent with no variation with energy over the 2-30 keV band.

To examine phase-resolved spectral behavior, we have performed fits to the spectra derived from each of the sinusoid model peaks using the full RXTE channel resolution. We used the XSPEC WABS\*PEGPWRLW spectrum model, a power law with photoelectric absorption, and have fixed the hydrogen column density to a constant $`1\times 10^{20}\mathrm{atoms}\mathrm{cm}^2`$ for all fits, based on ROSAT results (Ogelman et al. 1993). The resulting best fit photon power-law indices and uncertainties are listed in Table 1, along with $`\chi ^2`$ values. It is clear that a model in which all the power law indices are consistent with a single value can be rejected handily ($`\chi ^2`$ of 32.0 for 4 dof for a model with all indices equal to their mean). Peak 1 contributes the bulk of $`\chi ^2`$ in this case and is obviously harder than the remainder of the spectra. The other four spectra are only marginally inconsistent with fluctuations about their mean ($`\chi ^2`$ of 8.7 for 3 dof). Figure 2 shows the resulting photon spectra of the different fitted peaks in the RXTE light curve, along with data from CGRO, ROSAT and optical observations. The Peak 1 spectrum seems to connect smoothly with the OSSE spectrum, confirming and extending the turnover from the COMPTEL and EGRET spectra. The Peak 2-hard spectrum is much softer than Peak 1 and does not seem to have as smooth a connection to the OSSE spectrum. Peak 2-soft appears softer than Peak 2-hard, but is not well represented by the simple power-law model. Uncertainties in the two Peak 2 components are increased by possible contamination of each other due to the rather arbitrary nature of our peak shape model and the proximity of the two components. However, we have no a priori reason to use any specific model or assume a more complex model shape (e.g. asymmetric). Within these constraints, we believe that these spectra are representative of the component spectra. The Peak 4 (at the radio phase) spectrum is well resolved but soft and weak. Within uncertainties, it is indistinguishable from Peak 2-soft. The Peak 3 (radio precursor) spectrum has not been plotted in Figure 2, because the errors are large, but it also is consistent with the spectra of Peak 2-soft and Peak 4.

Our finding of clearly separate emission components in the phase resolved spectrum is confirmed by examination of the broad-band light curves in Figure 1. The hardest RXTE spectral components, Peak 1 and Peak 2-hard, are phase-aligned with EGRET Peaks 1 and 2. Interestingly, we find no evidence for the interpeak emission component seen by EGRET. This is consistent with EGRET’s finding that this is the hardest emission component in the Vela phase-resolved spectrum. The softest RXTE spectral components, Peak 2-soft, Peak 3 and Peak 4 are all phase-aligned with features in the optical and radio light curves. Peak 2-soft is aligned with the second optical peak. Peak 4, which is strongest in the 2-5 keV band, is at the phase of the radio pulse (phase 0.0/1.0) and with a peak in the optical light curve at the radio phase. Peaks 1 and 2 are roughly in phase with the first two peaks in the Chandra HRC profile. Peak 3 is in the vicinity of the very strong third peak in the Chandra profile.

In summary, then, we have identified separate hard and soft emission components in the Vela pulsar spectrum. An extrapolation of the soft component spectrum to lower energies is roughly consistent with the optical spectrum (Nasuti et al. 1997, Mignani and Caraveo 2001), while the hard component continues trends from higher energy. The hardest component, Peak 1, is phase aligned with features at higher energy, while the soft components, Peaks 2-soft and the radio phase peak, are aligned with features at lower energies. Peak 2-hard is an exception, being not much harder than the soft components but phase aligning with a $`\gamma `$-ray feature. Further observations with improved statistics would allow us to separate the two components of Peak 2 with more confidence.

In Figure 3, we plot our phase-averaged RXTE spectrum for better comparison with the other phase-averaged results, including the recent Chandra observations (Pavlov et al. 2001). Although it seems to smoothly bridge the gap between the ROSAT/Chandra and OSSE spectra, the RXTE phase-averaged spectrum is clearly not a simple power law but a transition between hard and soft spectral components. Pavlov et al. (2001) determined that a two-component model, blackbody plus power law, best fit the Chandra phase-averaged point-source spectrum in the 0.25 - 8 keV band. Our phase-averaged spectrum roughly agrees with the power-law component of the Chandra spectrum in hardness, even though it is not well-represented by a single power-law. Also, the current Chandra result represents total flux from the point source, while our result is “on-pulse” minus “off-pulse,” which explains the difference in normalization and could lead to a variation in shape as well. We will be able to make a better comparison when Chandra phase-resolved spectra are available.

## 5 Discussion

The RXTE Cycle 3 observations have independently confirmed the multicomponent nature of pulsed emission from the Vela pulsar in the energy range 2 - 30 keV as suggested by the Cycle 1 observations. With the improved statistics of the Cycle 3 data, we have been able to separate the broad second peak in the RXTE X-ray light curve into soft (Peak 2-soft) and hard (Peak 2-hard) spectral components which maintain their phase integrity throughout the RXTE energy range. In addition, we have discovered a new feature in the RXTE light curve: a peak (Peak 4) at the phase of the radio pulse with an extremely soft spectrum. There is, in addition, significant emission leading the radio phase (Peak 3). Peaks 1 and 2-hard make up the hard spectral component whose light curve peaks are in phase with those of the gamma-ray light curve, and whose spectrum smoothly connects to the 100 keV - 5 GeV spectrum. The Vela pulsar now joins the Crab pulsar in having phase-aligned pulse profiles and a continuous spectrum from hard X-ray to hard-gamma rays. Unlike Vela, the phase alignment of the Crab pulse profile continues down to the radio band. RXTE Peaks 2-soft, 3 and 4 make up the soft component and each are phase aligned with features in optical and radio light curves. Their spectra, consistent with each other in both flux and spectral index, extrapolate to the optical flux points.

The hard X-ray band is a complex transition region in the Vela pulsar spectrum, where at least two non-thermal emission components possibly produced at several different locations in the magnetosphere are simultaneously present. The 2-30 keV Vela spectrum may thus be a “Rosetta Stone” that will provide some valuable clues to deciphering pulsar emission mechanisms. Unfortunately, from the RXTE observations alone the theoretical picture is not yet clear, but there are a few interesting connections to be made with the existing models.

Pulsar high-energy emission models generally have fallen into two main categories. Polar cap models (e.g. Daugherty & Harding 1996, Sturner et al. 1995) advocate particle acceleration and emission in the magnetic polar regions within several stellar radii of the neutron star surface. High energy photon spectra are produced by primary particle curvature or inverse-Compton emission and additional synchrotron and inverse-Compton components from electron-positron pair cascades. In outer gap models (e.g. Cheng, Ho & Ruderman 1986, Romani 1996), particle acceleration takes place in the outer magnetosphere, in vacuum gaps along null charge surfaces. Although the same radiation processes that contribute to polar cap emission also produce the outer-gap high-energy emission, the pair production now requires soft X-ray photons from a hot neutron star rather than a strong magnetic field.

Although the RXTE hard component spectrum connects to the CGRO $`\gamma `$-ray spectrum, the X-ray spectrum is harder, requiring a break around 50-100 keV, a feature already suggested by OSSE (Strickman et al. 1996). Polar cap models predict such a break in the synchrotron spectrum of the pair cascade at the local cyclotron energy, blueshifted by the relativistic motion of the pairs along the magnetic field (Harding & Daugherty 1999). The predicted turnover energy is an increasing function of local field strength and thus provides a measure of the height of the emission above the neutron star surface. With pair Lorentz factors around $`\mathrm{\Gamma }_\pm 20`$, the turnover at 50-100 keV in the Vela hard X-ray spectrum would indicate an emission height of 1-2 stellar radii above the surface. The outer gap model of Romani (1996) predicts a turnover in the gamma-ray curvature spectrum and a transition to synchrotron radiation in the Vela spectrum, but at an energy around 1 MeV.

The RXTE soft component suggests continuous emission from optical to hard X-rays. In polar cap models, this emission cannot be synchrotron radiation since it is well below the cyclotron energy, but inverse Compton scattering radiation would be possible. Zhang & Harding (2000) have argued that inverse-Compton scattering of thermal X-rays by cascade pairs will produce observable soft X-ray emission and the computed spectrum of this radiation (Dyks et al. 2001) is possibly consistent with the Vela RXTE phase-averaged spectrum and is likely to contribute at least some of the soft emission component. However, scattering of much softer photons, possibly radio, is required to produce a non-thermal emission component extending down to the optical band. The fact that one of the RXTE soft component peaks is at the radio phase makes this a potentially attractive prospect for further study. In outer gap models, synchrotron radiation of backflowing particles from the gap accelerator produces non-thermal X-rays (Cheng & Zhang 1999). The predicted spectral index of -1.9 is consistent with the spectral index of the RXTE soft component.

We would like to thank George Pavlov and Divas Sanwal for providing their data on the Chandra pulse profile in advance of publication and for comments on the manuscript. We also thank David Nice for providing the ephemeris from the Princeton Pulsar Database.


format: text

text:
# Parity-Violating Electron Scattering as a Probe of Supersymmetry

\[

## Abstract

We compute the one-loop supersymmetric (SUSY) contributions to the weak charges of the electron ($`Q_W^e`$) and proton ($`Q_W^p`$) using the Minimal Supersymmetric Standard Model (MSSM). These $`q^2=0`$ vector couplings of the $`Z^0`$-boson to fermions will be determined in two fixed-target, parity-violating electron scattering experiments. The SUSY loop contributions to $`Q_W^p`$ and $`Q_W^e`$ can be substantial, leading to several percent corrections to the Standard Model values for these quantities. We show that the relative signs of the SUSY loop effects on $`Q_W^e`$ and $`Q_W^p`$ are correlated and positive over nearly all of the MSSM parameter space, whereas inclusion of R-parity nonconserving interactions can lead to opposite sign relative shifts in the weak charges. Thus, a comparison of $`Q_W^p`$ and $`Q_W^e`$ measurements could help distinguish between different SUSY scenarios.

\]

The search for physics beyond the Standard Model (SM) of electroweak and strong interactions is a primary objective for particle and nuclear physics. Historically, parity-violating (PV) interactions have played an important role in elucidating the structure of the electroweak interaction. In the 1970’s, PV deep inelastic scattering (DIS) measurements performed at the Stanford Linear Accelerator Center (SLAC) confirmed the SM prediction for the structure of weak neutral current interactions. These results were consistent with a value for the weak mixing angle given by $`\mathrm{sin}^2\theta _W1/4`$, implying a tiny $`V`$(electron)$`\times A`$(quark) neutral current interaction. Subsequent PV measurements – performed at both very low scales using atoms as well as at the $`Z^0`$-pole in $`e^+e^{}`$ annihilation – have been remarkably consistent with the results of the SLAC DIS measurement.

More recently, the results of cesium atomic parity-violation (APV)  and deep inelastic $`\nu `$\- ($`\overline{\nu }`$-) nucleus scattering have been interpreted as determinations of the scale-dependence of $`\mathrm{sin}^2\theta _W`$. The SM predicts how this parameter will evolve from its precisely measured value at the $`Z^0`$-pole. The cesium APV and neutrino DIS measurements imply $`2\sigma `$ and $`+3\sigma `$ deviations, respectively, from the predicted evolution of $`\mathrm{sin}^2\theta _W`$ (defined in the $`\overline{\mathrm{MS}}`$ scheme). If conventional atomic or hadron structure effects are ultimately unable to account for these discrepancies, the results of these precision measurements would point to new physics.

In light of this situation, two new measurements involving polarized electron scattering have taken on added interest: PV Möller ($`ee`$) scattering at SLAC and elastic, PV $`ep`$ scattering at the Jefferson Lab (JLab). In the absence of new physics, both measurements could be used to determine $`\mathrm{sin}^2\theta _W`$ at the same scale ($`|Q^2|0.03`$ $`(\mathrm{GeV}/c)^2`$) – falling between the scales relevant to the APV and neutrino DIS measurements – with comparable precision in each case. Any significant deviation from the SM prediction for $`\mathrm{sin}^2\theta _W`$ at this scale would provide striking evidence for new physics, particularly if both measurements report a deviation. On the other hand, agreement would imply that the most likely explanations for the cesium APV and neutrino DIS results are atomic and hadron structure effects within the SM.

In this Letter, we analyze the prospective implications of the parity-violating electron scattering (PVES) measurements for supersymmetry (SUSY). Although no supersymmetric particle has yet been discovered, there exists strong motivation for believing that SUSY is a component of the “new” Standard Model. For example, the existence of low-energy SUSY is a prediction of many string theories; it offers a solution to the hierarchy problem, providing a mechanism for maintaining the stability of the electroweak scale against large radiative corrections; and it results in coupling unification close to the Planck scale. In light of such arguments, it is clearly of interest to determine what insight about SUSY the new PVES measurements might provide.

In the simplest version of SUSY – the Minimal Supersymmetric Standard Model (MSSM) – low-energy precision observables experience SUSY only via tiny loop effects involving virtual supersymmetric particles. In the MSSM, the requirement of baryon minus lepton number ($`BL`$) conservation leads to conservation of the R-parity quantum number, $`P_R=(1)^{2S+3(BL)}`$, where $`S`$ denotes spin. Every SM particle has $`P_R=+1`$ while the corresponding superpartner, whose spin differs by $`1/2`$ unit, has $`P_R=1`$. Conservation of $`P_R`$ implies that every vertex has an even number of superpartners. Consequently, for processes like $`eeee`$, all superpartners must live in loops. Such loops generate corrections – relative to the SM amplitude – of order $`(\alpha /\pi )(M/\stackrel{~}{M})^2`$, where $`M`$ denotes a SM particle mass and $`\stackrel{~}{M}`$ is a superpartner mass. Thus, low-energy experiments must probe an observable with a precision of few tenths of a percent or better in order to discern SUSY loop effects. Low-energy charged current experiments have already reached such levels of precision, and the corresponding implications of these experiments for the MSSM have been discussed elsewhere.

The leading-order SM contribution to the PV $`ee`$ and $`ep`$ asymmetries is governed by the tree-level vector coupling of the $`Z^0`$-boson to these fermions – the so-called “weak charge”. At tree-level in the SM this coupling is suppressed: $`Q_W^p=Q_W^e=14\mathrm{sin}^2\theta _W0.1`$. One-loop SM electroweak radiative corrections further reduce this tiny number, leading to the predictions $`Q_W^e=0.0449`$ and $`Q_W^p=0.0721`$. The fortuitous suppression of these couplings in the SM renders them more transparent to the possible effects of new physics. Consequently, experimental precision of order a few percent, rather than a few tenths of a percent, is needed to probe SUSY loop corrections. As we show below, the possible magnitude of these effects may be as large as the proposed experimental error bars for the $`Q_W^e`$ and $`Q_W^p`$ measurements (8% and 4%, respectively ). Moreover, the relative sign of the effect in both cases is correlated – and positive – over nearly all available SUSY parameter space. To our knowledge, this correlation is specific to the MSSM, making it a potential low-energy signature of this new physics scenario.

The content of the MSSM has been described in detail elsewhere, so we review only a few features here. The particle spectrum consists of the SM particles and the corresponding superpartners: spin-0 sfermions ($`\stackrel{~}{f}`$), spin-$`1/2`$ gluinos ($`\stackrel{~}{g}`$), and spin-$`1/2`$ mixtures of Higgsinos and electroweak gauginos, the neutralinos ($`\stackrel{~}{\chi }_{14}^0`$) and charginos ($`\stackrel{~}{\chi }_{1,2}^\pm `$). In addition, the Higgs sector of the MSSM contains two doublets (“up”- and “down”-types, respectively), whose vacuum expectations $`v_u`$ and $`v_d`$ are parameterized in terms of $`v=\sqrt{v_u^2+v_d^2}`$ and $`\mathrm{tan}\beta =v_u/v_d`$. Together with the SU(2)L and U(1)Y couplings $`g`$ and $`g^{}`$ respectively, $`v`$ is determined from $`\alpha `$, $`M_Z`$, and the Fermi constant extracted from the muon lifetime, $`G_\mu `$, while $`\mathrm{tan}\beta `$ remains a free parameter. The MSSM also introduces a coupling between the two Higgs doublets characterized by the dimensionful parameter $`\mu `$.

Degeneracy between SM particles and their superpartners is lifted by the SUSY-breaking Lagrangian, which depends in general on 105 additional parameters. These include the SUSY-breaking Higgs mass parameters; the gaugino masses $`M_{1,2,3}`$; the left- (right-)handed sfermion mass parameters $`M_{\stackrel{~}{f}_L}^2`$ ($`M_{\stackrel{~}{f}_R}^2`$); and terms which mix $`\stackrel{~}{f}_L`$ and $`\stackrel{~}{f}_R`$ into mass eigenstates $`\stackrel{~}{f}_{1,2}`$. One expects the magnitude of the SUSY-breaking parameters to lie somewhere between the weak scale and $`1`$ TeV. Significantly larger values can reintroduce the hierarchy problem.

Theoretical models for SUSY-breaking mediation provide relations among this large set of parameters, generally resulting in only a few independent parameters at the SUSY-breaking or GUT scale. According to the model-independent analysis of Ref. , however, the superpartner spectrum implied by these models conflicts with the combined constraints of low-energy charged current data, $`M_W`$, and the muon anomalous magnetic moment unless one allows for nonconservation of $`P_R`$. Here, we adopt a similar model-independent approach and do not impose any specific relations among SUSY-breaking parameters. To our knowledge, no other model-independent analysis of MSSM corrections to neutral current observables has appeared in the literature, nor have the complete set of corrections to low-energy PV observables been computed previously (see, e.g. ).

MSSM loop corrections to elementary $`e`$-$`f`$ amplitudes consist of several topologies. Contributions to gauge-boson self energies can be expressed entirely in terms of the oblique parameters $`S`$, $`T`$, and $`U`$ in the limit that $`\stackrel{~}{M}M_Z`$. Since present collider limits allow for fairly light superpartners, however, we do not work in this limit. Consequently, the corrections arising from the photon self-energy ($`\mathrm{\Pi }_{\gamma \gamma }`$) and $`\gamma `$-$`Z^0`$ mixing tensor ($`\mathrm{\Pi }_{Z\gamma }`$) contain a residual $`q^2`$-dependence not embodied by the oblique parameters. Vertex corrections, external leg corrections, and box graphs are process-dependent and cannot be parameterized in any general way. Note that one must consider corrections to both neutral current and charged current amplitudes, since the former are normalized to $`G_\mu `$ and since $`\mathrm{sin}^2\theta _W`$ is calculated in the MSSM using $`\alpha `$, $`M_Z`$, $`G_\mu `$ and radiative corrections to these quantities. We evaluate these corrections using the modified dimensional reduction renormalization scheme ($`\overline{\mathrm{DR}}`$) .

Including these corrections, the general structure for an elementary $`A(e)\times V(f)`$ amplitude is

$$=\frac{G_\mu }{2\sqrt{2}}Q_W^f\overline{e}\gamma _\mu \gamma _5e\overline{f}\gamma ^\mu f$$
(1)
with

$$Q_W^f=\rho _{PV}\left[2T_3^f4\kappa _{PV}Q_f\mathrm{sin}^2\theta _W\right]+\lambda _f.$$
(2)
Here, $`\rho _{PV}=1+\delta \rho _{PV}`$ and $`\kappa _{PV}=1+\delta \kappa _{PV}`$, with $`\delta \rho _{PV}`$ and $`\delta \kappa _{PV}`$ denoting a universal set of corrections to PV amplitudes and $`\lambda _f`$ containing the effects of process-specific vertex, external leg, and box graph effects. It is useful to express the SUSY contributions to $`\delta \rho _{PV}`$ and $`\delta \kappa _{PV}`$ in terms of oblique parameters as well as a residual $`q^2`$-dependence of $`\mathrm{\Pi }_{Z\gamma }`$ and $`\mathrm{\Pi }_{\gamma \gamma }`$:

$`\delta \rho _{PV}^{\mathrm{susy}}`$ $`=`$ $`\widehat{\alpha }T\widehat{\delta }_{VB}^\mu `$ (3)
$`\delta \kappa _{PV}^{\mathrm{susy}}`$ $`=`$ $`\widehat{\alpha }\left({\displaystyle \frac{\widehat{c}^2}{\widehat{c}^2\widehat{s}^2}}\right)T+\widehat{\alpha }\left[{\displaystyle \frac{1}{4\widehat{s}^2(\widehat{c}^2\widehat{s}^2)}}\right]S`$ (6)
$`+{\displaystyle \frac{\widehat{c}}{\widehat{s}}}\left[{\displaystyle \frac{\widehat{\mathrm{\Pi }}_{Z\gamma }(q^2)}{q^2}}{\displaystyle \frac{\widehat{\mathrm{\Pi }}_{Z\gamma }(M_Z^2)}{M_Z^2}}\right]`$
$`+\left({\displaystyle \frac{\widehat{c}^2}{\widehat{c}^2\widehat{s}^2}}\right)\left[{\displaystyle \frac{\mathrm{\Delta }\widehat{\alpha }}{\widehat{\alpha }}}{\displaystyle \frac{\widehat{\mathrm{\Pi }}_{\gamma \gamma }(M_Z^2)}{M_Z^2}}+\widehat{\delta }_{VB}^\mu \right].`$
Here, the hat denotes quantities renormalized in the $`\overline{\mathrm{DR}}`$ scheme; $`\widehat{\delta }_{VB}^\mu `$ are the SUSY vertex, external leg, and box graph corrections to the $`\mu `$-decay amplitude; $`\widehat{s}`$ ($`\widehat{c}`$) is the sine (cosine) of the weak mixing angle in the $`\overline{\mathrm{DR}}`$ scheme defined at the scale $`\mu =M_Z`$; and $`\mathrm{\Delta }\widehat{\alpha }`$ is the SUSY contribution to the difference between the fine structure constant and the electromagnetic coupling renormalized $`\mu =M_Z`$. Note that superpartner loops give $`\widehat{\mathrm{\Pi }}_{Z\gamma }(q^2)q^2`$.

In order to evaluate the potential size of SUSY corrections, we generated set of $`3000`$ different combinations of SUSY-breaking parameters, chosen randomly from flat distribution in mass parameters and $`\mathrm{ln}\mathrm{tan}\beta `$. The former were bounded below by present collider limits and bounded above by 1000 GeV, corresponding to the $`𝒪`$(TeV) naturalness limit. We also restricted $`\mathrm{tan}\beta `$ to lie in the range $`1.4<\mathrm{tan}\beta <60`$ as required by perturbativity of third generation quark Yukawa couplings, and we allowed for left-right mixing among sfermions. In order to avoid unacceptably large flavor-changing neutral currents, we have also assumed no intergenerational sfermion mixing.

For each combination of parameters, we evaluate superpartner masses and mixing angles, which we then use as inputs for computing the radiative corrections. We also separately evaluate the corresponding contributions to the oblique parameters. The latter are tightly constrained from precision electroweak data. We rule out any parameter combination leading to values of $`S`$ and $`T`$ lying outside the present 95% confidence limit contour for these quantities. We note that this procedure is not entirely self-consistent, since we have not evaluated non-universal MSSM corrections to other precision electroweak observables before extracting oblique parameter constraints. As noted in Ref. , where MSSM corrections to $`Z`$-pole observables were evaluated using different models for SUSY-breaking mediation, non-universal effects can be as large as oblique corrections. Nevertheless, we expect our procedure to yield a reasonable estimate of the oblique parameter constraints. Since $`S`$ and $`T`$ do not dominate the low-energy SUSY corrections (see below), our results depend only gently on the precise allowed ranges for these parameters.

For the case of charged current observables, gluino loops can generate substantial corrections when the masses and mixing angles for $`u`$\- and $`d`$-type squarks are not identical. In contrast, gluinos decouple entirely from the one-loop MSSM corrections to semi-leptonic neutral current PV observables. Moreover, MSSM Higgs contributions to vertex, external leg, and box graph corrections are negligible due to the small, first- and second-generation Yukawa couplings. The light Higgs contribution to gauge boson propagators has already been included via the oblique parameters, while the effects of other MSSM Higgs bosons are sufficiently small to be neglected.

In Fig. 1, we plot the shift in the weak charge of the proton, $`\delta Q_W^p=2\delta Q_W^u+\delta Q_W^d`$, versus the corresponding shift in the electron’s weak charge, $`\delta Q_W^e`$, normalized to the respective SM values. The corrections can be as large as $`4\%`$ ($`Q_W^p`$) and $`8\%`$ ($`Q_W^e`$) – roughly the size of the proposed experimental errors. Generally speaking, the magnitudes of $`\delta Q_W^{e,p}`$ decrease as SUSY mass parameters are increased. The largest effects occur when at least one superpartner is relatively light. An exception occurs in the presence of significant mass splitting between sfermions, which may lead to sizable contributions. However, such weak isospin-breaking effects also increase the magnitude of $`T`$, so their impact is bounded by oblique parameter constraints. This consideration has been implemented in arriving at Fig. 1. We also observe that the presence or absence of sfermion left-right mixing affects the distribution of points, but does not significantly change the range of possible corrections. For the situation of no left-right mixing, the points are more strongly clustered near the origin. Thus, while corrections of the order of several percent are possible in either case, large effects are more likely in the presence of left-right mixing.

The shifts $`\delta Q_W^{e,p}`$ are dominated by $`\delta \kappa _{PV}^{\mathrm{susy}}`$. Non-universal corrections involving vertex corrections and wavefunction renormalization experience significant cancellations, while box graphs are numerically suppressed. We find that $`\delta \kappa _{PV}^{\mathrm{susy}}`$ is nearly always negative, corresponding to a reduction in the effective $`\mathrm{sin}^2\theta _W`$ for the PVES experiments. Within $`\delta \kappa _{PV}^{\mathrm{susy}}`$ itself, contributions from the various terms in Eq. (6) have comparable importance, with some degree of cancellation occurring between the effects of $`S`$ and $`T`$. Thus, the oblique parameter approximation gives a rather poor description of the MSSM effects on the weak charges.

As evident from Fig. 1, the relative sign of the corrections to both $`Q_W^p`$ and $`Q_W^e`$ is nearly always the same and nearly always positive. Since $`Q_W^p>0`$ ($`Q_W^e<0`$) in the SM, SUSY loop corrections give $`\delta Q_W^p>0`$ ($`\delta Q_W^e<0`$). This correlation is significant, since the effects of other new physics scenarios can display different signatures. For example, for the general class of $`E_6`$ neutral gauge bosons (with mass $``$ 1000 GeV), the effects on $`Q_W^p`$ and $`Q_W^e`$ also correlate, but $`\delta Q_W^{e,p}/Q_W^{e,p}`$ can have either sign in this case. Leptoquark interactions, in contrast, would not lead to discernible effects in $`Q_W^e`$ but could induce sizable shifts in $`Q_W^p`$ .

As a corollary, we also note that the relative importance of SUSY loop corrections to the weak charge of heavy nuclei probed with APV is suppressed. The shift in the nuclear weak charge is given by $`\delta Q_W(Z,N)=(2Z+N)\delta Q_W^u+(2N+Z)\delta Q_W^d`$. Since the sign of $`\delta Q_W^f/Q_W^f`$ due to superpartner loops is nearly always the same, and since $`Q_W^u>0`$ and $`Q_W^d<0`$ in the SM, a strong cancellation between $`\delta Q_W^u`$ and $`\delta Q_W^d`$ occurs in heavy nuclei. This cancellation implies that the magnitude of $`\delta Q_W(Z,N)/Q_W(Z,N)`$ is generally less than about 0.2% for cesium and is equally likely to have either sign. Since the presently quoted uncertainty for the cesium nuclear weak charge is about 0.6%, the present $`2\sigma `$ deviation from the SM prediction does not substantially constrain the SUSY parameter space.

While agreement between experimental values for $`Q_W^{e,p}`$ and the SM predictions would not produce significant new constraints on the MSSM, a deviation of $`2\sigma `$ or more could help distinguish between the MSSM and other new physics scenarios. For example, a $`+2.5\sigma `$ deviation for $`Q_W^{e,p}`$, along with the cesium APV result, would be consistent with an $`E_6`$ $`Z^{}`$ boson but fairly difficult to accomodate in the MSSM. To illustrate, we plot in Fig. 2 the impact of such a result on the MSSM parameter space. Here, we have assumed no L-R mixing, set $`\mathrm{tan}\beta =10`$, and taken $`\mu =M_1=M_2/2=100`$ GeV. The plot indicates the allowed left-handed slepton and squark masses ($`\stackrel{~}{M}_L`$ and $`\stackrel{~}{M}_Q`$), assuming for simplicity a common mass for all generations. Present collider searches rule out the dark shaded region, while charged current data exclude the light shaded area. The hatched region would be ruled out (at 95% C.L.) by a large, positive $`\delta Q_W^{p,e}`$, yielding an upper bound $`\stackrel{~}{M}_Q300`$ GeV. Should future collider limits exceed this value, then a $`Z^{}`$ would be the favored explanation of of large cesium APV and $`Q_W^{p,e}`$ deviations.

Alternately, one may relax the assumption of $`P_R`$ conservation. Doing so leads to new interactions of the type $`\stackrel{~}{e}\overline{e}\nu `$ and $`\stackrel{~}{d}^{}\overline{e}^cu`$ with coupling strengths $`\lambda _{ijk}`$ and $`\lambda _{ijk}^{}`$ (the subscripts refer to generation number). These interactions give rise to new tree-level exchange of sfermions between SM fermions. As discussed in Ref. , charged current data, $`M_W`$, and the results of cesium APV can be accommodated under this scenario if $`\lambda _{12k}`$ and $`\lambda _{11k}^{}`$ are nonzero. The interior of the truncated ellipse in Fig. 1 shows the possible corrections to $`Q_W^p`$ and $`Q_W^e`$, given the constraints from other electroweak data. We observe that the prospective effects of $`P_R`$ non-conservation are quite distinct from SUSY loops. The shifts $`\delta Q_W^p/Q_Q^p`$ and $`\delta Q_W^e/Q_W^e`$ have opposite signs over most of the allowed region, in contrast to the situation for SUSY loop effects or $`E_6`$ $`Z^{}`$ bosons. Thus, a comparison of results for the two PVES experiments could help determine whether this extension of the MSSM is to be favored over other new physics scenarios.

We thank R. Carlini, B. Filippone, S. Heinemeyer, R.D. McKeown, and M. Wise for useful comments. This work was supported in part by the U.S. Department of Energy and National Science Foundation.


text:
# Overlap Distribution of the Three-Dimensional Ising Model

## I Introduction

In this paper we investigate the two replica overlap probability density $`P_L(q)`$ for the three-dimensional ($`3d`$) Ising model. On a $`L^3`$ lattice with periodic boundary conditions $`q`$ is defined by

$$q=\frac{1}{N}\underset{i=1}{\overset{N}{}}s_i^{(1)}s_i^{(2)}\mathrm{with}N=L^3,$$
(1)
where $`s_i^{(1)}`$ and $`s_i^{(2)}`$ are the spins of two copies (replica) of the system at temperature $`T=1/\beta `$. The distribution of the overlap $`q`$ is of major importance in spin-glass investigations , where it plays the role of an order parameter, often called Parisi order parameter.

To our knowledge this quantity has never been investigated for simple spin systems like the $`3d`$ Ising model. One reason is certainly that one has in that situation the magnetization $`m`$ as an explicit order parameter at hand and a description of the critical properties based on the magnetic probability density $`P_L^m(m)`$ is believed to be identical to one based on $`P_L(q)`$, in particular $`q=m^2`$. However, the overlap probability density is an interesting object for study on its own merits and we find remarkable differences between the shapes of $`P_L(q)`$ and $`P_L^m(m)`$. Therefore, we find it worthwhile to have the properties of $`P_L(q)`$ documented for the Ising model, which is by orders of magnitude easier to simulate than spin glasses, since the dynamics is much faster and only one (instead of many) realization needs to be simulated.

In the vicinity of the critical point, by finite-size scaling (FSS) arguments  $`P_L(q)`$ can, in leading order for $`L`$ large, be written as

$$P_L(q)=\frac{1}{\sigma _L}P^{}(q^{})\mathrm{with}q^{}=\frac{q}{\sigma _L}.$$
(2)
Here $`P^{}`$ is a universal, $`L`$-independent function and $`\sigma _L`$ is the standard deviation of $`q`$ with respect to the probability density $`P_L(q[1,+1])`$ (or $`P_L(q[0,1])`$ when appropriate).

A major focus of our investigation is on the tails of the $`P_L(q)`$ distribution, which we control for $`L=36`$ at $`T_c`$ over 500 orders of magnitude by using the multi-overlap MC algorithm . This is also of interest in view of a conjecture by Bramwell et al.  that a variant of extreme order statistics describes the asymptotics of certain probability densities for a large class of correlated systems. Besides the Ising model with some $`T(L)T_c`$ as $`L\mathrm{}`$, their class includes the $`2dXY`$ model in the low temperature phase, turbulent flow problems, percolation models and some self-organized critical phenomena. For large $`L`$ the asymptotic behavior is claimed to be described by an $`L`$-independent curve, which for the overlap variable would read ($`q^{}\mathrm{}`$)

$$P^{}(q^{})=C\mathrm{exp}\left[a\left(q^{}q_{\mathrm{max}}^{}e^{b(q^{}q_{\mathrm{max}}^{})}\right)\right].$$
(3)
Here $`C`$, $`a`$, $`b`$ are constants and $`q_{\mathrm{max}}^{}=q_{\mathrm{max}}/\sigma _L`$, where $`q_{\mathrm{max}}`$ is the position of the maximum of the probability density $`P_L(q)`$ at positive $`q`$. Equation (3) is a variant of Gumbel’s first asymptote , see  for reviews of extreme order statistics.

However, Eq. (3) is in contradiction with the widely accepted large-deviation behavior, based on the proportionality of the entropy with the volume ,

$$P_L(q)\mathrm{exp}\left[Nf(q)\right],$$
(4)
where, for large $`N`$, $`f(q)`$ does not depend on $`N`$. Our data support Eq. (4). Using the multimagnetical approach  a similar study of the tails could be performed for the magnetic probability density $`P_L^m(m)`$, but this is outside the scope of our present paper.

We like to point out that for the overlap distribution of spin glasses the status of Eq. (4) is unclear due to the quenched average. Our previously reported result  demonstrates that for the $`3d`$ Edwards-Anderson Ising spin glass a probability distribution of the form (3) gives an excellent description of the tails of the Parisi overlap distribution. Because of the special nature of its phase transition, Eq. (4) may also be questioned for the $`2dXY`$ model, where the extreme order asymptotics (3) (more precisely a variant of it) with $`a=\pi /2`$ is found in the spin wave approximation . However, the range of validity of this perturbative argument is unclear, at least to us. It may be worthwhile to employ the methods of Ref. , or those of the present paper, to perform a careful numerical investigation of the the $`2dXY`$ model with respect to these questions.

We have performed simulations at and below the critical (Curie) temperature $`T_c`$ of the Ising model phase transition. We approximate $`T_c`$ by the value of Ref.

$$\beta _c=\frac{1}{T_c}=\mathrm{0.221\hspace{0.17em}654}$$
(5)
and present our results for $`T_c`$ in Sect. II. Besides addressing the question of the asymptotic behavior of the overlap distribution, we estimate the critical exponent ratio $`2\beta /\nu `$ from the FSS behavior of the standard deviation $`\sigma _L`$. For the temperatures below the Curie temperature we choose

$$\beta _1=0.232\mathrm{and}\beta _2=0.3.$$
(6)
For $`\beta =0.232`$ multimagnetical results are available , which determine the probability density $`P_L^m(m)`$ of the magnetization over many orders of magnitude.

Our numerical results were obtained with the spin-glass code of the investigations of  by simply choosing all the exchange coupling constants to be equal to $`+1`$. A code which is specialized to the Ising model would be far more efficient. Therefore, we have limited our present simulations to small and medium sized lattices. As the results appear already quite clear, there seems to be no particularly strong reason to push on towards (much) larger systems.

## II Results at the critical temperature

Figure 1 shows our overlap probability density results $`P_L(q)`$ at the critical temperature (5). They rely on a statistics of 32 independent runs (with different pseudo random number sequences) for lattices up to size $`L=30`$ and on 16 independent runs for our largest lattice, $`L=36`$. After calculating the multi-overlap parameters  the following numbers of sweeps were performed per repetition (i.e. independent run) : $`2^{19}`$ for $`L=4`$, $`2^{21}`$ for $`L=6`$, $`2^{22}`$ for $`L=8`$, $`2^{23}`$ for $`L=12,16`$, $`2^{24}`$ for $`L=24`$, $`2^{25}`$ for $`L=30`$ and, again, $`2^{24}`$ for $`L=36`$ (with the present computer program this lattice size became too time consuming to scale its CPU time properly). Not to overload Fig. 1 error bars are only shown for selected values of $`q`$, whereas the lines are drawn from all data. The probability densities are normalized to

$$_1^{+1}𝑑qP_L(q)=1,$$
(7)
and we show only the $`q0`$ part because of the symmetry $`P_L(q)=P_L(q)`$. We cut the range at $`q=0.5`$, because for $`L8`$ the probability densities are almost zero for $`q0.5`$.

Somewhat surprisingly we find the maximum of our $`P_L(q)`$ probabilities at $`q_{\mathrm{max}}=0`$, in contrast to the magnetization where one finds a double peak at $`T_c`$, see for instance numerical data in Ref.  and analytical results of Ref. , both with periodic boundary conditions. In our simulation we kept a time series for the magnetization, which reproduces the expected double peaked histograms at $`T_c`$ (as the accuracy of our magnetization histogram is lower than that of results in the literature, we abstain from giving a figure). Such differences are expected. For instance, while the two low-temperature magnetization values $`m=\pm 1`$ correspond to four overlap configurations, two with $`q=+1`$ and two with $`q=1`$, the inverse is not true. There are altogether $`2^N`$ overlap configurations with $`q=+1`$ and another $`2^N`$ with $`q=1`$.

Figure 2 shows $`\mathrm{ln}\left[P_L(q)\right]`$ versus $`q`$. The ordinate is cut off at $`1000`$, to cover a range with results from at least two lattice sizes. The $`L=36`$ lattice continues to exhibit accurate results down to $`1200`$, thus the data from this system cover $`1200/\mathrm{ln}(10)=521`$ orders of magnitude.

The collapse of the $`P_L(q)`$ functions (2) on one universal curve $`P^{}(q^{})`$ is depicted in Fig. 3. The figure shows some scaling violations, which become rather small from $`L24`$ onwards. The standard deviation $`\sigma _L`$ behaves with $`L`$ according to

$$\sigma _LL^{2\beta /\nu }\left(1+c_2L^\omega +\mathrm{}\right)$$
(8)
Note that the ratio $`\beta /\nu `$ is defined for the magnetization, and by FSS theory  $`\sigma _L^mL^{\beta /\nu }`$ holds for the standard deviation of the magnetization. The factor of two difference in the exponent of Eq. (8) comes from dimensionality. Scaling relations and estimates of the Ising model critical exponents are reviewed in Ref. . In particular, $`2\beta /\nu =d2+\eta `$ holds. Our estimate of $`2\beta /\nu `$ from a four-parameter fit of Eq. (8) to our data for the standard deviation $`\sigma _L`$ is $`2\beta /\nu =1.0293(28)`$ with $`Q=0.31`$ the goodness of fit (see  for the definition of $`Q`$). Restricted to our $`L24`$ lattices the more stable two-parameter fit to the leading behavior of Eq. (8) gives

$$\frac{2\beta }{\nu }=1.030\pm 0.005\mathrm{with}Q=0.36.$$
(9)
The most accurate estimates of the literature  cluster around $`\eta =0.036`$ with an error of a few units in the last digit. Within the conventional statistical uncertainties this is consistent with our $`2\beta /\nu `$ values. The two-parameter fit becomes quickly inconsistent when the smaller lattices with $`L<24`$ are included, with a trend towards smaller values of $`2\beta /\nu `$. Therefore, we conjecture that there will be a slightly increasing trend when larger lattices should become available. Because its larger error bar reflects to some extent systematic uncertainties, we prefer Eq. (9) over the four-parameter fit as our final estimate.

In Fig. 4 we show the logarithm $`\mathrm{ln}[P^{}(q^{})]`$ of the rescaled overlap probability densities and we see a breakdown of scaling for sufficiently large $`q^{}`$. The ordering of the lattices is that the rightmost curve corresponds to the $`L=36`$ lattice. The smaller lattices deviate from it. From the left: First, the $`L=4`$ lattice (not visible), next the $`L=6`$, then $`L=8`$, $`L=12`$, $`L=16`$, $`L=24`$ and last $`L=30`$. The agreement is over a larger and larger range in $`q^{}`$. However, scaled back to $`q`$, it concerns the vicinity of $`q=0`$. To quantify this, we have collected in Table I the $`q^{}`$ and corresponding $`q`$ values at which the deviation

$$\mathrm{}_{36}\mathrm{ln}[P_L^{}(q^{})]\mathrm{ln}[P_{36}^{}(q^{})]\mathrm{ln}[P_L^{}(q^{})]$$

becomes 1/2, a deviation too small to be visible on the scale of Fig. 4. The $`q^{}`$ values are seen to increase, whereas the corresponding $`q`$ values decrease.

It is well known that the requirement of consistency of a universal probability density (2) with the functional form (4) determines the function $`f(q)`$. Namely, to leading order the scaling of the function (2) $`P_L^{}(q^{})`$ implies that

$$L^df\left(q^{}L^{2\beta /\nu }\right)=g(q^{})$$
(10)
is an $`L`$-independent function. Therefore,

$$f(q)q^{d\nu /2\beta }$$
(11)
holds. Note that the non-critical Gaussian behavior is a special case, obtained for $`d\nu /2\beta =2`$.

In contrast to Eqs. (2) and (11), the functional form (4) is expected to hold for all $`q`$, when $`L`$ becomes large. This is easily tested by plotting

$$f_L(q)=\frac{1}{N}\mathrm{ln}\left[P_L(q)\right],$$
(12)
as is done in Fig. 5, and seeing if $`f_L(q)`$ is $`L`$-independent up to $`O(1/N)`$ terms as it should. Not to obscure the behavior by too large symbols, the lines are plotted without error bars. A fit of the scaling form (11) to our $`f_{36}(q)`$ data for $`q<0.2`$, $`f(q)=0.0000490.073q^{d\nu /2\beta }`$ with $`2\beta /\nu =1.030`$ from Eq. (9), is also included in the figure.

We see excellent convergence towards an $`L`$-independent function, where the higher lying curves correspond to the smaller lattices ($`L=4`$ being the one on top). However, the scaling behavior (11) only holds in the vicinity of $`q=0`$. To make this quantitatively more precise, we subtract the function $`f_{36}(q)`$ from the others and plot the difference in Fig. 6 (at selected values of $`q`$ we now include barely visible error bars by plotting $`f_L(q)\pm \mathrm{}f_L(q)f_{36}(q)`$). In the large volume limit (for both of two lattices) the difference $`|f_{L_1}(q)f_{L_2}(q)|`$ should be bounded by a constant

$$\left(\frac{1}{N_1}\frac{1}{N_2}\right).$$

We see in Fig. 6 that the $`L=24`$ and 30 curves fall almost on the $`L=36`$ one (the zero line). Note that the figure is cut off at $`q=0.62`$. The number of sweeps needed to propagate the system over the full admissible $`q`$-range scales in the multi-overlap ensemble at least proportional to the system size $`N`$. Aiming at a comparable statistics for all system sizes, the required computer time thus grows at least proportional to $`N^2`$. Therefore and because of numerical problems with the floating point representation caused by the extreme smallness of $`P_L(q)`$ for $`q1`$ when $`L`$ is large, we restricted the overlap simulations to $`q[0.7,+0.7]`$ for the lattice sizes $`L=16,24`$, and 30, and to $`q[0.62,+0.62]`$ for the $`L=36`$ lattice. Nevertheless the smallest values of $`P_L(q)`$ we sampled were those of the $`L=36`$ lattice.

On the basis of Eq. (4) the plots of Figs. 5 and 6 had to be expected. The conjecture (3) of Bramwell et al.  appears to be ruled out for the Ising model. Namely, when Eq. (3) (with $`q_{\mathrm{max}}^{}=0`$) and Eq. (12) are both valid in some region of $`bq^{}=bL^{2\beta /\nu }`$, one finds to leading order

$$bqL^{2\beta /\nu }=d\mathrm{ln}(L)+\mathrm{ln}[f(q)],$$
(13)
and $`f(q)`$ can only be $`L`$-independent if $`b`$ is not a constant, but depends on $`q`$ and $`L`$.

## III Results below the critical temperature

Below the critical temperature we made 16 independent runs per lattice size with the following numbers of sweeps per repetition: $`2^{16}`$ for $`L=4`$, $`2^{17}`$ for $`L=6`$, $`2^{18}`$ for $`L=8`$, $`2^{19}`$ for $`L=12`$, $`2^{20}`$ for $`L=16(\beta =0.232)`$ and $`2^{21}`$ for $`L=16(\beta =0.3)`$. The overlap probability densities $`P_L(q)`$ for $`\beta =0.232`$ and $`\beta =0.3`$ are shown together in Fig. 7. Clearly, the peaks moved away from zero and are now at $`q_{\mathrm{max}}=0.3408`$ ($`L=16,\beta =0.232`$) and $`q_{\mathrm{max}}=0.8237`$ ($`L=16,\beta =0.3`$). In Figs. 8 and 9 we show the logarithms of these probability densities at $`\beta =0.232`$ and $`\beta =0.3`$, respectively. The scales in these two figures are chosen to accommodate all the $`P_L(0)`$ data, but not their tails, which continue down to much lower values. For the largest ($`L=16`$) difference between $`P_L(0)`$ and the maximum of $`P_L(q)`$, we see that it increases from about four orders of magnitude at $`\beta =0.232`$ to about 65 orders of magnitude at $`\beta =0.3`$.

Note that the most likely $`P_L(0)`$ configurations are those where one replica stays around magnetization $`m=0`$ and the other around the maximum of the magnetic probability density at positive or negative magnetization $`m`$. It follows that the ratio

$$R_L=\frac{P_L(0)}{P_L(q_{\mathrm{max}})},$$
(14)
where $`P_L(q_{\mathrm{max}})`$ is the maximum of $`P_L(q)`$, is related to the interface tension $`F_s`$ according to the formula introduced in  for the ratio $`P_L^m(0)/P_L^m(q_{\mathrm{max}})`$,

$$R_L=CL^p\mathrm{exp}\left[2L^2F_s\right]+\mathrm{}.$$
(15)
Here $`C`$, $`p`$ are constants and $`p=1/2`$ in the one-loop capillary wave approximation or one-loop $`\mathrm{\Phi }^4`$-theory , compare the discussion in Ref. . Two-loop $`\mathrm{\Phi }^4`$-theory is considered in Ref. . The correction is large and it appears that a reliable estimate of $`p`$ does not exist.

To determine the interface tension one may first calculate the lattice size dependent effective interface tensions

$$F_{s,L}=\frac{1}{2L^2}\mathrm{ln}R_L,$$
(16)
and then make an extrapolation of $`F_{s,L}`$ for $`L\mathrm{}`$. Table II collects our $`F_{s,L}`$ results where the error bars with respect to the last digits are given in parentheses. For the sake of easy comparison, we list also some $`F_{s,L}`$ results of Ref.  at $`\beta =0.232`$, obtained by applying the definitions (14) and (16) to the probability density of the magnetization. It is notable that the $`F_{s,L}`$ estimates from the overlap densities increase monotonically in the listed range of lattice sizes, whereas the $`F_{s,L}`$ estimates from the magnetization show a more complex behavior: Up to $`L=12`$ they decrease, then they turn around to increase and the increase has been followed  up to lattices of size $`L=32`$.

We pursue a similar fitting strategy as in Ref. . As there, it turns out that our data do not really support fits to more than two parameters and that including the capillary wave term with the one-loop theoretical coefficient $`p`$ does not lead to any improvements of the goodness $`Q`$ of the fits. In essence we are left with fits to the leading, likely effective, correction

$$F_{s,L}=F_s+\frac{a_1}{L},$$
(17)
and, due to our small lattice sizes, finite-size corrections are so big that the best we can do is a fit of the interface tensions from the largest two lattices, $`L=12`$ and 16. This yields the estimates

$$F_s=0.03002(24)\mathrm{at}\beta =0.232,$$
(18)
and

$$F_s=0.30583(68)\mathrm{at}\beta =0.3,$$
(19)
which are (under the circumstances of our limited system sizes) in reasonably good agreement with results of Hasenbusch and Pinn (HP) , for a review see Ref. . Again, our error bars are purely statistical and do not reflect systematic errors due to our small lattice sizes.

Our result (18) at $`\beta =0.232`$ is lower than the multimagnetical estimate of Ref. . This moves into the right direction and indicates that the resolution of the inconsistency between the multicanonical and the HP estimate at $`\beta =0.232`$, discussed in the paper by Zinn and Fisher , has its origin in the complex finite-size scaling behavior of $`F_{s,L}`$ estimates from the magnetization, which could be resolved by simulating larger systems. It is notable that this difficulty of the extrapolation appears to be limited to a small neighborhood of $`\beta =0.232`$, as the multimagnetical $`F_s`$ estimates  at $`\beta =0.227`$ and $`\beta =0.2439`$ are perfectly consistent with HP, see Fig. 1 of Ref. .

## IV Summary and Conclusions

In Sect. II we have analyzed the critical behavior of the overlap variable $`q`$. In essence agreement with the standard scaling picture is found, but with some new insights. In particular, we exhibit in Table I that scaling appears to be confined to a small $`q`$ (but still large $`q^{}=qL^{2\beta /\nu }`$) neighborhood. It may be worthwhile to check whether the magnetic probability distribution, for which comparable simulations are easier to perform, exhibits a similar behavior. Further, we find support in favor of standard large deviations (4), instead of the form (13) derived from Gumbel’s first asymptote (3).

Below the critical point, in Sect. III, we estimate interface free energies from our overlap probability densities. The results are smoother than those from the probability density of the magnetization  and tend to reconcile discrepancies noted in Ref. . But, like at the critical point, considerably larger lattices would be needed to reach high precision results.

###### Acknowledgements.

We would like to thank Gernot Münster und Uwe-Jens Wiese for helpful correspondence. This work was in part supported by the US Department of Energy under contract DE-FG02-97ER41022. The numerical simulations were performed on the Compaq SC256 computer of CEA in Grenoble under grant p526 and on the T3E computer of NIC in Jülich under grant hmz091.


text:
# Evidence for a Low-Spin to Intermediate-Spin State Transition in LaCoO3

## Abstract

We present measurements of the magnetic susceptibility and of the thermal expansion of a LaCoO3 single crystal. Both quantities show a strongly anomalous temperature dependence. Our data are consistently described in terms of a spin-state transition of the Co3+ ions with increasing temperature from a low-spin ground state ($`t_{2g}^{\mathrm{\hspace{0.17em}6}}e_g^0`$) to an intermediate-spin state ($`t_{2g}^{\mathrm{\hspace{0.17em}5}}e_g^1`$) without (100 K - 500 K) and with ($`>500`$K) orbital degeneracy. We attribute the lack of orbital degeneracy up to 500 K to (probably local) Jahn-Teller distortions of the CoO6 octahedra. A strong reduction or disappearance of the Jahn-Teller distortions seems to arise from the insulator-to-metal transition around 500 K.

Transition-metal oxides have fascinating physical properties as e.g. high-temperature superconductivity in the cuprates or colossal magnetoresistance in the manganites. Their properties are often governed by a complex interplay of charge, magnetic, structural, and orbital degrees of freedom. Moreover, for a given oxidation state some transition metals display different spin states as it is the case in various cobalt oxides. Quite recently a class of layered cobalt compounds with the chemical composition REBaCo2O5+δ (RE = rare earth) has attracted considerable interest. These compounds show a broad variety of ordering phenomena and other transitions, e.g. (antiferro- and/or ferro-) magnetic order, charge and/or orbital order, metal-insulator transitions or spin-state transitions vogt00a ; suard00a ; moritomo00a ; martin97a ; respaud01a ; wu01a ; wang01a ; roy02a ; akahoshi99a . For TlSr2CoO5 it has been proposed that a metal-insulator transition is driven by a spin disproportionation, which consists of an alternating ordering of Co3+ ions in an intermediate-spin state (IS: $`t_{2g}^{\mathrm{\hspace{0.17em}5}}e_g^1`$; S = 1) and in a high-spin state (HS: $`t_{2g}^{\mathrm{\hspace{0.17em}4}}e_g^2`$; S = 2) doumerc99a ; foerster01a .

The occurrence of Co3+ in different spin states is known since the 1950s from LaCoO3 jonker53a ; bhide72a , which transforms with increasing temperature from a non-magnetic insulator to a paramagnetic insulator around 100 K and shows an insulator-to-metal transition around 500 K. But even for this rather simple pseudo-cubic perovskite the nature of these transitions is still unclear. The ground state is usually attributed to the low-spin configuration (LS: $`t_{2g}^{\mathrm{\hspace{0.17em}6}}e_g^0`$; S = 0) and the paramagnetic behavior above 100 K to the thermal population of an excited state. However, the question whether the excited state has to be identified with the HS or the IS state is subject of controversial discussions. Early publications often assume a LS/HS scenario asai94a ; itoh94a ; yamaguchi96a . In order to explain the insulating nature up to 500 K an ordering of LS and HS Co3+ ions has been proposed which vanishes at the insulator-to-metal transition senaris95a ; raccah67a . Yet the presence of a HS configuration below 400 K has been questioned on the basis of X-ray absorption and photoemission experiments abbate94a . Alternative descriptions of LaCoO3 favoring a LS/IS scenario potze95a ; saitoh97b ; asai98a ; yamaguchi97a ; kobayashi00b  are mainly based on the results of LDA+U calculations korotin96a , which propose that due to a strong hybridization between Co-$`e_g`$ levels and O-$`2p`$ levels the IS state is lower in energy than the HS state. Within this scenario the occurrence of orbital order and its melting have been proposed in order to explain the insulating nature below 500 K and the insulator-to-metal transition, respectively korotin96a . Up to now there is no experimental evidence, neither for an orbitally ordered state, nor for a HS/LS superstructure.

In this paper we present a study of the thermal expansion $`\alpha `$ and of the magnetic susceptibility $`\chi `$ of a LaCoO3 single crystal. The combined analysis of $`\alpha `$ and $`\chi `$ gives clear evidence for a thermal population of the IS state without orbital degeneracy. The lack of orbital degeneracy could arise from orbital order as proposed in Ref. korotin96a  or it can be interpreted as a consequence of Jahn-Teller (JT) distortions of the CoO6 octahedra with Co3+ in the JT-active IS state. Above 600 K, our $`\chi (T)`$ analysis suggests the presence of the IS state with orbital degeneracy, which may arise from a suppression or strong reduction of the JT distortion due to the insulator-to-metal transition.

The crystal used in this study was cut from a large single crystal ($`l`$$``$$`8`$ cm; $`ø`$$``$$`6`$ mm) grown by the floating-zone technique in an image furnace. The crystal is strongly twinned as usual for distorted perovskites. The magnetization has been measured by a SQUID magnetometer in the temperature range from 2 K up to 300 K in an applied field of 50 mT and by a Faraday balance in the temperature range from 200 K up to 1000 K in a field of 1 T. A high-resolution measurement of the linear thermal expansion $`\alpha `$$`=`$$`\frac{1}{L}\frac{L}{T}`$ has been performed using a capacitance dilatometer from 4 K to 180 K.

The magnetic susceptibility of our LaCoO3 crystal (Fig. 1) agrees well with that found in previous studies bhide72a ; itoh94a ; yamaguchi96a ; saitoh97b . The maximum around 100 K signals the spin-state transition of the Co3+ ions. For higher temperatures $`\chi (T)`$ shows (i) a Curie-like decrease up to about 500 K, (ii) a temperature-independent plateau between 500 K and 600 K and (iii) again a Curie-like decrease above 600 K. The increase of $`\chi (T)`$ below 30 K arises most probably from a Curie contribution due to magnetic impurities and/or oxygen non-stoichiometry. In order to obtain the Curie susceptibility $`\chi ^C`$ of Co3+ we subtract a term $`(P/T+\chi _0)`$ from the raw data. Here $`P/T`$ represents the impurity contribution and $`\chi _0=\chi _{\mathrm{dia}}+\chi _{\mathrm{vV}}`$ the sum of the diamagnetic contribution of the core electrons and of the paramagnetic van Vleck susceptibility of Co3+. A fit of the low-temperature data gives $`P`$$`=`$$`0.02`$ emuK/mole and $`\chi _0`$$`=`$$`6.510^4`$ emu/mole. The magnitude of $`P`$ allows to estimate an impurity content of less than 1% polarons . Our value of $`\chi _0`$ is close to those observed e.g. in Ref. itoh94a  in LaCoO3 or in ZnCo2O4 with Co3+ in the LS state miyatani66 .

The susceptibility of a system with a non-magnetic ground state and a magnetic excited state reads

$$\chi ^C(T)=\frac{N_\mathrm{A}g^2\mu _\mathrm{B}^2}{3k_\mathrm{B}T}\frac{\nu S(S+1)(2S+1)e^{\mathrm{\Delta }/T}}{1+\nu (2S+1)e^{\mathrm{\Delta }/T}}.$$
(1)
Here $`N_\mathrm{A}`$ is the Avogadro number, $`\mu _\mathrm{B}`$ the Bohr magneton and $`k_\mathrm{B}`$ the Boltzmann constant, $`\mathrm{\Delta }`$ denotes the energy splitting of the two states, $`g`$ is the Landé factor, $`S`$ the spin and $`\nu `$ the orbital degeneracy of the excited state. For simplicity we consider a purely ionic model for LaCoO3 and the spin-only values for the magnetic moments, where the ground state of Co3+ is the LS state with S = 0 and the electronic configuration $`t_{2g}^{\mathrm{\hspace{0.17em}6}}e_g^0`$. The excited state is either the IS (S = 1 and $`t_{2g}^{\mathrm{\hspace{0.17em}5}}e_g^1`$) or the HS (S = 2 and $`t_{2g}^{\mathrm{\hspace{0.17em}4}}e_g^2`$) state. In the HS (IS) state the $`t_{2g}^{}`$ ($`t_{2g}^{}`$ and $`e_g^{}`$) levels are only partially filled. Therefore the HS state consists of 3 orbital states, which are degenerate in a cubic crystal field leading to $`\nu `$$`=`$$`3`$. The IS state contains 6 orbital states but even in a cubic field these states are split into two orbital triplets separated by an energy of about 1 eV due to the coulomb interaction within the 3d shell sugano . Thus, an orbital degeneracy $`\nu `$$`=`$$`3`$ is expected for the IS state, too. The orbital degeneracy is lifted by lower symmetries and the crystal symmetry of LaCoO3 is only rhombohedral. It is, however, unclear whether this can be observed in $`\chi ^C(T)`$. First, the rhombohedral distortion of LaCoO3 is very small. Second, the crystal field of the Co3+ levels is mainly determined by the surrounding O2- ions and it is not clarified whether the CoO6 octahedra are distorted, i.e. whether the local symmetry of the Co3+ ions is less than cubic. From neutron scattering only one Co-O distance is reported thornton86a , whereas optical data give evidence for different Co-O distances yamaguchi97a . This question is further discussed below. Our combined analysis of $`\chi `$ and $`\alpha `$ will allow an unambiguous decision which of the four scenarios, LS/IS or LS/HS with or without orbital degeneracy yields the appropriate description of the spin-state transition of LaCoO3.

The inset of Fig. 1 shows an attempt to describe $`\chi ^C(T)`$ within a LS/HS scenario. Note that Eq. 1 has only $`\mathrm{\Delta }`$ as a free parameter (The $`g`$-factor may be varied only to some extent). In order to reproduce the strong increase of $`\chi ^C`$ below 100 K one has to use an energy splitting $`\mathrm{\Delta }`$$``$$`290`$ K but in this case the calculated $`\chi ^C(T)`$ for $`T`$$`>`$$`100`$ K is much larger than the experimental data (solid line in the inset of Fig. 1). In this calculation we have set $`g`$$`=`$$`2`$ and $`\nu `$$`=`$$`1`$. Assuming an orbital degeneracy $`\nu `$$`=`$$`3`$ of the HS state even increases the discrepancy between experimental and calculated data. There are some possibilities to improve the description within a LS/HS scenario. One is to use $`g`$$``$$`1.1`$ (and $`\mathrm{\Delta }`$$``$$`190`$ K) but such a small $`g`$-factor is very unlikely chargeorder . Another one is to introduce an antiferromagnetic nearest-neighbor coupling $`J_{\mathrm{AF}}`$ as has been done in Ref. yamaguchi96a . On a mean field level this leads to

$$\chi ^{MF}(T)=\frac{\chi ^C(T)}{1+\frac{k_\mathrm{B}}{N_\mathrm{A}g^2\mu _\mathrm{B}^2}zJ_{\mathrm{AF}}\chi ^C(T)}$$
(2)
where $`\chi ^C(T)`$ is given by Eq. 1 and $`z`$ is the number of nearest neighbors. A fit according to Eq. 2 for $`T`$$``$$`400`$ K is shown by the dashed line in the inset of Fig. 1. By setting $`g`$$`=`$$`2`$, $`\nu `$$`=`$$`1`$ and $`z`$$`=`$$`6`$ we obtain $`\mathrm{\Delta }`$$``$$`220`$ K and $`J_{\mathrm{AF}}`$$``$$`56`$ K (a fit for $`\nu `$$`=`$$`3`$ yields $`\mathrm{\Delta }270`$ K and $`J_{\mathrm{AF}}`$$``$$`62`$ K). Such a strong antiferromagnetic coupling is, however, in clear contradiction to neutron scattering experiments which give evidence for a weak ferromagnetic coupling asai94a . From these arguments we conclude that a LS/HS scenario does not give a consistent description of the magnetic properties of LaCoO3 up to about 500 K.

In contrast to the LS/HS scenario, a good description is found for $`\chi ^C(T)`$ within a LS/IS scenario. This is shown by the solid line in the main panel of Fig. 1. By setting $`\nu =1`$ the fit yields $`\mathrm{\Delta }180`$ K and a reasonable $`g2.1`$. When an orbital degeneracy $`\nu =3`$ of the IS state is assumed, the quality of the fit becomes worse (not shown in Fig. 1). Thus, the fit of $`\chi ^C(T)`$ favors an orbitally non-degenerate IS state below 500 K. Orbital degeneracy might become important above 600 K as indicated by the dashed line in Fig. 1, which is obtained by ’switching on’ the orbital degeneracy (setting $`\nu =3`$) and leaving the other parameters unchanged ($`\mathrm{\Delta }180`$ K and $`g2.1`$). Obviously, the dashed line is quite close to the experimental $`\chi ^C(T)`$ for $`T>600`$ K.

In Fig. 2 we show the linear thermal expansion of LaCoO3 and of La0.82Sr0.18CoO3. Our high-resolution data confirm previous results obtained by neutron diffraction asai94a  but allow a more detailed analysis. Whereas $`\alpha `$ of La0.82Sr0.18CoO3 shows a weak monotonous increase with temperature as expected for ordinary solids the thermal expansion of LaCoO3 is highly unusual: It is rather large and has a pronounced maximum around 50 K. In view of the spin-state transition that occurs in LaCoO3 (but not in La0.82Sr0.18CoO3) a straightforward interpretation of the anomalous behavior of $`\alpha `$ can be given. In the LS state of Co3+ all electrons occupy $`t_{2g}`$ levels whereas in the IS and HS state $`e_g`$ levels are also occupied. Since the $`e_g`$ states are oriented towards the surrounding negative O2- ions a population of Co3+ in the IS (or HS) state causes an additional widening of the lattice. Simply speaking, the ionic radius of Co3+ depends on its spin state and increases with increasing number of electrons in the $`e_g`$ states ($`r_{\mathrm{Co}^{3+}}^{\mathrm{HS}}>r_{\mathrm{Co}^{3+}}^{\mathrm{IS}}>r_{\mathrm{Co}^{3+}}^{\mathrm{LS}}`$).

Note that it is not possible to detect a sharp anomaly of $`\alpha `$ at a characteristic temperature in agreement with specific heat data which also do not show such an anomaly stolen97a . That means, the spin-state transition in LaCoO3 is not a phase transition in the thermodynamic sense. This justifies the description of the spin-state transition by a thermal population of an excited magnetic state from the LS state which remains the state of the lowest energy (as assumed in Eq. 1). The (additional) relative length change due to the spin-state transition is proportional to the thermal population of the excited state, and the anomalous thermal expansion is given by its temperature derivative, i.e. by

$$\mathrm{\Delta }\alpha (T)=d\frac{\mathrm{\Delta }\nu (2S+1)e^{\mathrm{\Delta }/T}}{T^2(1+\nu (2S+1)e^{\mathrm{\Delta }/T})^2}.$$
(3)
The product $`\nu (2S+1)`$ is the total degeneracy of the excited IS (HS) state and $`d`$ is determined by the different Co–O bond lengths for Co3+ in the IS (HS) and in the LS state, respectively.

The anomalous thermal expansion $`\mathrm{\Delta }\alpha `$ of LaCoO3 shown in the right panel of Fig. 2 is obtained by subtracting $`\alpha `$ of La0.82Sr0.18CoO3 from the raw data. With respect to $`\mathrm{\Delta }\alpha `$ the different scenarios only differ by the total degeneracy $`\nu (2S+1)`$, which amounts to 3 and 5 for the IS and HS state without and to 9 and 15 for the IS and HS state with orbital degeneracy, respectively. In Fig. 2 we only show the fits for a LS/IS scenario without ($`\nu (2S+1)`$ = 3; solid line) and for a LS/HS scenario with orbital degeneracy ($`\nu (2S+1)`$ = 15; dashed line). The two other fits are lying between the solid and the dashed curve. The fit for $`\nu (2S+1)`$ = 3 gives the best description of the experimental $`\mathrm{\Delta }\alpha (T)`$ and yields $`\mathrm{\Delta }`$$``$$`185`$ K in good agreement with $`\mathrm{\Delta }`$$``$$`180`$ K from the fit of $`\chi ^C`$. The deviations above 50 K may arise from the uncertainty in background determination and/or a temperature dependence of $`\mathrm{\Delta }`$. Depending on the model the values of $`d`$ vary between $`0.66`$ % and $`0.38`$ % (see table 1). They are much smaller than the difference of the Co–O bond lengths of 3% obtained from the sums of the tabulated ionic radii shannon76  of O2- and Co3+ in the LS and the HS state ($`1.89`$ Å and $`1.95`$ Å), respectively, giving further evidence against a LS/HS and for a LS/IS scenario in LaCoO3, because the Co–O bond length of the IS state is expected to lie between those of the LS and the HS state of Co3+.

The separate fits of both $`\chi ^C(T)`$ and $`\mathrm{\Delta }\alpha (T)`$ already favor a LS/IS scenario without orbital degeneracy but much more convincingly this conclusion is obtained by a scaling behavior between both quantities. From Eqs. 1 and 3 a straightforward calculation yields

$$C\mathrm{\Delta }\alpha (T)=\frac{(T\chi ^C(T))}{T}$$
(4)
with $`C`$ = $`\frac{N_\mathrm{A}g^2\mu _\mathrm{B}^2}{3k_\mathrm{B}}\frac{S(S+1)}{d}`$. As shown in Fig. 2 this scaling behavior is well fulfilled by the experimental data. The scaling unambiguously reveals that the anomalous thermal expansion of LaCoO3 arises from the spin-state transition. Moreover, the scaling factor $`C^{exp}`$ = $`172`$ emuK/mole agrees almost perfectly with $`C`$ = $`167`$ emuK/mole expected for a LS/IS scenario without orbital degeneracy (see table 1). In contrast, the other three scenarios yield strongly different scaling factors ranging from 252 to 870 emuK/mole. Thus, our data clearly exclude descriptions within a LS/IS scenario with orbital degeneracy saitoh97b ; asai98a  and also within a LS/HS scenario with or without orbital degeneracy asai94a ; itoh94a ; yamaguchi96a ; senaris95a .

Within an ionic picture the LS/IS scenario in LaCoO3 is surprising. The energies of the different spin states of Co3+ are determined by the balance of the Hund’s rule coupling of parallel spins and the crystal-field splitting $`\mathrm{\Delta }_{CF}`$ between the $`t_{2g}`$ and the $`e_g`$ levels. The ground state is either the HS (small $`\mathrm{\Delta }_{CF}`$) or the LS (large $`\mathrm{\Delta }_{CF}`$) but never the IS state. Moreover, the IS is expected to lie at least about 1 eV above the ground state sugano . In order to explain the much smaller value observed experimentally ($`\mathrm{\Delta }`$$``$$`185`$ K $``$$`0.016`$ eV) the energy of the IS state has to be lowered relative to the LS and HS states. This can arise from a hybridization between the Co-$`e_g`$ and the O-$`2p`$ levels as has been found in LDA+U band-structure calculations korotin96a . We note that an additional stabilization of the IS state can arise from a JT distortion. The IS state is strongly JT-active because of the partially filled $`e_g`$ level. The HS state is only weakly JT-active since the gain of JT energy is much less in a partially filled $`t_{2g}`$ level, and the LS is not JT-active at all. Thus, a JT distortion favors the IS state and it easily explains the lifting of the orbital degeneracy. Remarkably, this picture can also account for the behavior of $`\chi ^C`$ above 500 K. Due to the insulator-to-metal transition there are delocalized charge carriers. The corresponding charge fluctuations usually suppress (or weaken) the JT distortions. Then an orbital degeneracy of $`3`$ has to appear above 500 K as suggested by our $`\chi ^C(T)`$ analysis. In addition, an increase of the energy gap $`\mathrm{\Delta }`$ may be expected because of the loss of JT energy. However, the analysis of $`\chi ^C(T)`$ above 600 K does not allow to determine a high-temperature value of $`\mathrm{\Delta }`$ with sufficient accuracy. This latter point deserves further clarification.

In summary, our combined study of the magnetic susceptibility and the thermal expansion shows that the spin-state transition in LaCoO3 is consistently described by a thermal population of the intermediate-spin state. The intermediate-spin state has no orbital degeneracy up to about 500 K. This may arise from (local) Jahn-Teller distortions of the CoO6 octahedra. We analyzed our data within a simple ionic model of LaCoO3 but we stress that the experimentally observed scaling between $`\mathrm{\Delta }\alpha `$ and $`(\chi T)/T`$ is model-independent and may serve as a sensitive test of more sophisticated models.

Valuable discussions with M. Braden, A. Freimuth, M. Haverkort, Z. Hu, E. Müller-Hartmann, D. Khomskii, and L.H. Tjeng are acknowledged. This work was supported by the Deutsche Forschungsgemeinschaft through SFB 608.


text:
# Formation and Interaction of Membrane Tubes

\]submitted to Phys. Rev. Lett., February 8, 2002

(\[)

## Abstract

We show that the formation of membrane tubes (or membrane tethers), which is a crucial step in many biological processes, is highly non-trivial and involves first order shape transitions. The force exerted by an emerging tube is a non-monotonic function of its length. We point out that tubes attract each other, which eventually leads to their coalescence. We also show that detached tubes behave like semiflexible filaments with a rather short persistence length. We suggest that these properties play an important role in the formation and structure of tubular organelles.

Biological membranes (such as the endoplasmic reticulum, the Golgi apparatus, the inner mitochondrial membrane, or the plasma membrane) often form highly dynamic tubular networks tube . The formation and transport of membrane tubes (with tens of nanometers in diameter) are thought to involve motor proteins that are able to grab the membrane and pull on it as they move along the filaments of the cytoskeleton motion . Nanotubes (or tethers) can also be pulled out by various experimental techniques (such as hydrodynamic flow waugh82 , micropipettes evans96 , or optical tweezers raucher99 ), and very recently Roux et al. roux02  have managed to set up a minimal in vitro experimental system, in which tubes have been pulled by kinesin motor proteins. The physics of long tubular membranes is relatively simple and is well understood. However, the initial formation of nanotubes from planar or large spherical pieces of membrane is a subtle process, the understanding of which is crucial to the study of various biological processes (involving tube formation or simply membrane pulling). In addition to tube formation we also study the interaction of membrane tubes pulled from the same membrane, and make a few comments on the mechanical properties of detached membrane tubes.

Living cells maintain the surface tension of most of their membranes at a constant level by keeping lipid reservoirs (at fixed chemical potentials) raucher99 . A constant pressure in closed vesicles is also maintained via osmosis. Therefore, for our study we choose the ensemble in which the surface tension $`\sigma `$ and the inside pressure $`p`$ (relative to the outside) are fixed rather than the surface area of the membrane $`A`$ or the volume of the vesicle $`V`$.

With the bending term included the free energy of the membrane can be written as helfrich

$$=\frac{\kappa }{2}(2H)^2dA+\sigma ApVfL,$$
(1)
where $`\kappa `$ is the bending rigidity and $`H`$ is the mean curvature of the membrane. The membrane is pulled in the $`Z`$ direction with a point force $`f`$, and the end-to-end distance of the membrane in this direction is denoted by $`L`$.

For a tube of length $`L`$ and radius $`R`$ the free energy (at $`p=0`$) can be written as $`_{\mathrm{tube}}=[\kappa /(2R^2)+\sigma ]2\pi RLfL`$. To minimize $`_{\mathrm{tube}}`$ the surface tension acts to reduce the radius, while the bending rigidity works against this. The balance between the two sets the equilibrium radius $`R_0`$ and force $`f_0`$, which can be calculated by taking $`_{\mathrm{tube}}/R=0`$ and $`_{\mathrm{tube}}/L=0`$:

$$R_0=\sqrt{\frac{\kappa }{2\sigma }},f_0=2\pi \sqrt{2\sigma \kappa }.$$
(2)
For typical values of $`\kappa 40`$ pNnm and $`\sigma 0.05`$ pN/nm one finds $`R_020`$ nm and $`f_012.6`$ pN.

For simplicity we consider only axisymmetric surfaces with the $`Z`$ axis being the symmetry axis (Fig. 1a upper inset). Such surfaces can be parametrized by the angle $`\psi (S)`$, where $`S`$ is the arclength along the contour. The coordinates $`R(S)`$ and $`Z(S)`$ depend on $`\psi (S)`$ through

$$\dot{R}=\mathrm{cos}\psi \text{and}\dot{Z}=\mathrm{sin}\psi ,$$
(3)
and the mean curvature can be expressed as

$$2H=\dot{\psi }+(\mathrm{sin}\psi )/R.$$
(4)
For axisymmetric surfaces one can derive the so-called general shape equation from the free energy (1) by variational methods axisymmetric ; firstint :

$`\stackrel{\dot{}\dot{}\dot{}}{\psi }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{\psi }^3{\displaystyle \frac{2\mathrm{cos}\psi }{R}}\ddot{\psi }+{\displaystyle \frac{3\mathrm{sin}\psi }{2R}}\dot{\psi }^2+{\displaystyle \frac{3\mathrm{cos}^2\psi 1}{2R^2}}\dot{\psi }`$ (5)
$`+\overline{\sigma }\dot{\psi }{\displaystyle \frac{\mathrm{cos}^2\psi +1}{2R^3}}\mathrm{sin}\psi +{\displaystyle \frac{\overline{\sigma }}{R}}\mathrm{sin}\psi \overline{p},`$
where $`\overline{\sigma }=\sigma /\kappa =1/(2R_0^2)`$ and $`\overline{p}=p/\kappa `$. Taking the first integral of this equation leads to firstint

$`\ddot{\psi }\mathrm{cos}\psi `$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{\psi }^2\mathrm{sin}\psi {\displaystyle \frac{\mathrm{cos}^2\psi }{R}}\dot{\psi }+{\displaystyle \frac{\mathrm{cos}^2\psi +1}{2R^2}}\mathrm{sin}\psi `$ (6)
$`+\overline{\sigma }\mathrm{sin}\psi {\displaystyle \frac{1}{2}}\overline{p}R{\displaystyle \frac{\stackrel{~}{f}}{R}},`$
with the integration constant $`\stackrel{~}{f}=f/(2\pi \kappa )=f/(f_0R_0)`$.

For nearly flat membranes ($`\psi 1`$), after the parameter change $`[\psi (S),R(S)]\psi (R)`$, one can expand the shape equation (6) in powers of $`\psi `$ and keep only the terms up to linear order:

$$R^2\psi ^{\prime \prime }+R\psi ^{}(R^2\overline{\sigma }+1)\psi =\stackrel{~}{f}R\overline{p}R^3/2.$$
(7)
Because of the parameter change, the primes denote derivations with respect to $`R`$ (which, in linear order of $`\psi `$, are identical to derivations with respect to $`S`$). The general solution of this differential equation is

$$\psi (R)=\frac{\stackrel{~}{f}}{\overline{\sigma }}\frac{1}{R}+\frac{1}{2}\frac{\overline{p}}{\overline{\sigma }}R+c_1I_1(R\sqrt{\overline{\sigma }})+c_2K_1(R\sqrt{\overline{\sigma }}),$$
(8)
where $`I_i(x)`$ and $`K_i(x)`$ are modified Bessel functions, and $`c_1`$ and $`c_2`$ are integration constants. $`I_1(R\sqrt{\overline{\sigma }})`$ diverges exponentially for $`R\mathrm{}`$. Because for a big vesicle we expect the shape to converge to that of a sphere, $`c_1`$ must vanish. At $`R=0`$ the divergence of the $`1/R`$ term must be canceled by the $`K_1`$ term, leading to $`c_2=\stackrel{~}{f}/\sqrt{\overline{\sigma }}`$. Integrating $`\psi (R)`$ with respect to $`R`$ gives the shape of the membrane in this linear approximation:

$$Z_{\mathrm{lin}}(R)=Z_0\frac{2R_0f}{f_0}\left[\mathrm{ln}(\frac{R}{\sqrt{2}R_0})+K_0(\frac{R}{\sqrt{2}R_0})\right]\frac{R^2}{2R_{\mathrm{ves}}},$$
(9)
where the integration constant $`Z_0`$ serves as a reference coordinate, and we have expressed $`\sigma `$, $`\kappa `$, and $`p`$ in terms of $`R_0`$, $`f_0`$, and vesicle radius $`R_{\mathrm{ves}}=2\sigma /p`$. The last term is a trivial contribution, describing a spherical vesicle under tension $`\sigma `$ and pressure $`p`$. The second term, which is proportional to $`f`$, is the linear response, and describes the deformation of the vesicle. The quantity between the brackets converges to $`[\mathrm{ln}(2)\gamma ]`$ for $`R0`$, where $`\gamma =0.577\mathrm{}`$ is the Euler constant. For $`p=0`$ and large $`R`$ the logarithmic term dominates, which corresponds to a catenoid, the well known shape of a soap film in cylindrical geometry under zero pressure.

Because the pressure makes only a trivial contribution, and has a negligible effect on tube formation from big vesicles (plug $`pRpR_0=2\sigma R_0/R_{\mathrm{ves}}\sigma `$ into the shape equations), we neglect it from now on, and consider a piece of (initially flat) membrane that spans a ring of radius $`R_{\mathrm{ring}}`$ located at $`Z=0`$.

For large deformations (or pulling forces) the linear approximation breaks down, and we have to solve Eq. (5) together with Eq. (3) numerically. Note that Eq. (6) could also be used, but it is numerically less stable. We start solving the differential equations from the ring, where we impose a zero curvature (or free hinge) boundary condition. Thus, the four initial parameters at $`Z=0`$ are as follows: (i) $`R=R_{\mathrm{ring}}`$; (ii) $`\psi =\mathrm{arcsin}[f/(2\pi \sigma R_{\mathrm{ring}})]\epsilon `$, where the small deviation $`\epsilon `$ from the catenoid shape is chosen (with a shooting and matching technique) such that the contour line reaches the $`Z`$ axis; (iii) $`\dot{\psi }=(\mathrm{sin}\psi )/R`$ ensures zero mean curvature \[see Eq. (4)\]; (iv) finally, $`\ddot{\psi }`$ is determined from Eq. (6).

The results of the numerical solution can be seen in Fig. 1. The main part of Fig. 1a shows the shape (contour line) of the membrane for $`R_{\mathrm{ring}}=20R_0`$ and different values of $`L`$. For small deformations (left line) Eq. (9) gives a very good approximation of the shape. In this linear regime the size of the deformation is approximately $`L_{\mathrm{lin}}=Z_{\mathrm{lin}}(0)Z_{\mathrm{lin}}(R_{\mathrm{ring}})`$.

For larger deformations a tube emerges in the middle, and the linear approximation fails. However, far from the tubular part, where $`\psi `$ is small, the approximation is still valid (Fig. 1a lower inset). Thus, it is convenient to define the tube as the piece between $`L_{\mathrm{lin}}`$ and $`L`$, and the base as the rest of the membrane between 0 and $`L_{\mathrm{lin}}`$. This way, the dependence of the size of the total deformation $`L`$ on the ring radius $`R_{\mathrm{ring}}`$ is absorbed in the size of the base $`L_{\mathrm{lin}}`$, and the length of the tube $`L_{\mathrm{tube}}=LL_{\mathrm{lin}}`$ (i.e., the deviation from the linear approximation) becomes independent of $`R_{\mathrm{ring}}`$.

With this definition it is enough to determine the $`f`$-$`L`$ curve for one particular ring size (e.g., $`20R_0`$), from which the universal (ring size independent) $`f`$-$`L_{\mathrm{tube}}`$ curve can be calculated (Fig. 1b inset). Because only the $`f`$-$`L`$ curves have real physical meaning, one can easily calculate them from the $`f`$-$`L_{\mathrm{tube}}`$ curve for any ring radius ($`R_0`$) by simply adding $`L_{\mathrm{lin}}(f,R_{\mathrm{ring}})`$ to $`L_{\mathrm{tube}}`$. The most intriguing feature of the $`f`$-$`L`$ curves (shown in Fig. 1b for three different values of $`R_{\mathrm{ring}}`$) is their non-monotonicity. The force first grows linearly, in accordance with the linear approximation (9), and converges to $`f_0`$ for large $`L`$. But in between it overshoots by $`0.13f_0`$, and then oscillates about $`f_0`$ with an exponentially decaying amplitude. This oscillation results in infinitely many intervals with negative slopes, which, in the $`f`$-ensemble (where $`f`$ is the control parameter rather than $`L`$), are mechanically unstable and represent an infinite series of first order shape transitions at $`f_0`$. Because $`L_{\mathrm{lin}}`$ is a monotonically increasing function of both $`f`$ and $`R_{\mathrm{ring}}`$, the main peak becomes an overhang for large rings ($`>20000R_0`$), leading to a first order transition even in the $`L`$-ensemble.

For a nearly cylindrical section, like the tubular part of the membrane, the shape equations can be expanded in powers of $`U(Z)=R(Z)R_0`$. After changing parameters $`[\psi (S),Z(S)]\psi (Z)`$, the expansion of Eq. (5) up to first order in $`U(Z)`$ and for $`p=0`$ reduces to bozic01

$$R_0^4U^{\prime \prime \prime \prime }=U.$$
(10)
Here the primes denote derivations with respect to $`Z`$. The solution of this equation is the sum of two exponentially decaying oscillations from the two ends ($`L_{\mathrm{lin}}`$ and $`L`$) of the tube:

$`{\displaystyle \frac{U(Z)}{R_0}}=`$ $`a_1\mathrm{exp}\left({\displaystyle \frac{ZL_{\mathrm{lin}}}{\sqrt{2}R_0}}\right)\mathrm{cos}({\displaystyle \frac{ZL_{\mathrm{lin}}}{\sqrt{2}R_0}}+\alpha _1)`$
$`+`$ $`a_2\mathrm{exp}\left({\displaystyle \frac{LZ}{\sqrt{2}R_0}}\right)\mathrm{cos}({\displaystyle \frac{LZ}{\sqrt{2}R_0}}+\alpha _2),`$ (11)
where the integration constants converge to $`a_10.746`$, $`a_20.726`$, $`\alpha _10.347`$, and $`\alpha _23.691`$ (determined by numerical fitting) as the length of the tube increases. These oscillations can be seen in Fig. 2 as overshootings from both the base and the tip of the tube. They can be understood intuitively by noticing that the smaller (larger) mean curvature near the base (tip) makes the effect of the bending rigidity on the tube radius less (more) pronounced, leading to a smaller (larger) $`R`$.

The expansion of Eq. (6) up to the first non-vanishing order in $`U(Z)`$ gives us the deviation of the force:

$$\frac{ff_0}{f_0}=R_0^2U^{}U^{\prime \prime \prime }\frac{1}{2}R_0^2U^{\prime \prime 2}+\frac{1}{2}\frac{1}{R_0^2}U^2.$$
(12)
Plugging Eq. (11) into Eq. (12) yields

$$\frac{ff_0}{f_0}=2a_1a_2\mathrm{exp}\left(\frac{L_{\mathrm{tube}}}{\sqrt{2}R_0}\right)\mathrm{cos}(\frac{L_{\mathrm{tube}}}{\sqrt{2}R_0}+\alpha _1+\alpha _2),$$
(13)
which is an exponentially decaying oscillation as a function of $`L_{\mathrm{tube}}`$. This explains the observed non-monotonic behavior of the $`f`$-$`L`$ curves, and gives a very good fit to the $`f`$-$`L_{\mathrm{tube}}`$ curve (Fig. 1b inset) with the same values of $`a_i`$ and $`\alpha _i`$ as determined above for the shape oscillations. In a different ensemble, Heinrich et al. heinrich99  also observed similar oscillations for axially strained vesicles.

The surface energy density $`\rho =(\kappa /2)(2H)^2+\sigma `$ of the membrane is plotted in Fig. 2. It grows from $`\sigma `$ to $`2\sigma `$ as we enter the tube from the base, and diverges at the tip. Because at the tip the angle $`\psi `$ is small again, Eq. (8) can be used to describe the surface there, leading to a logarithmic divergence of the mean curvature podgornik95 ; bozic01 .

In experiments, bundles of tubes can often be observed roux02 . To study the interaction between two tubes, let us start with a planar membrane that spans a ring and is pulled perpendicularly by two point forces ($`f_1`$ and $`f_2`$) at a distance $`d`$. If $`d`$ is large enough ($`R_0dR_{\mathrm{ring}}`$), both protrusions (except for the vicinity of the points of pulling) can be described by the leading logarithmic term of the linear approximation (9). Their superposition can be used to calculate the $`d`$-dependence of the free energy of the membrane. Note that $`fL`$ in Eq. (1) has to be replaced by $`f_1L_1+f_2L_2`$, where $`L_1`$ and $`L_2`$ are the sizes of the deformations. Straightforward calculation results in an attractive potential $`(d)=\mathrm{const}+2R_0\mathrm{ln}(d)f_1f_2/f_0`$ between the two deformations. The same attraction has been found between membrane-bound adhesion molecules bruinsma98 . For the attraction between two tubes replace both $`f_1`$ and $`f_2`$ by $`f_0`$.

To see what happens when two tubes get close, we performed numerical energy minimization with the Surface Evolver program brakke92 . We defined the initial topology as a piece of membrane spanning a ring of radius $`12R_0`$ at $`Z=0`$, with two cylindrical protrusions in the positive $`Z`$ direction. We applied zero curvature boundary condition at the ring and reflecting boundary condition at $`Z=24R_0`$ (Fig. 3). To make sure that the two tubes do not coalesce we placed a very narrow ($`R_0`$) cylindrical obstacle between the tubes, perpendicularly to them, and at a distance $`h`$ from the ring (not shown in Fig. 3). The obstacle simply pinned down the middle of the membrane at $`Z=h`$. We then measured the pinning force $`f_\mathrm{p}`$ (exerted by the obstacle on the membrane) as a function of $`h`$. Fig. 3 shows that $`f_\mathrm{p}`$ is always negative (i.e., pushes the membrane), meaning that there is no energy barrier against coalescence, it occurs smoothly. This is consistent with the experiments of Evans et al. evans96 .

Tubes can be detached via fission and then transported individually in the cell as long prolate-like vesicles (keeping their original radius $`R_0`$). Because of volume and area constraints their energy contains only the bending term: $`=(\kappa /2)(2H)^2dA`$. Although with such a high area-to-volume ratio the energetically most favorable shape is the stomatocyte, these prolates are metastable seifert91 . They are also very flexible. If we consider a prolate as a rod and bend it with a curvature $`C`$ ($`1/R_0`$), its energy, for symmetry reasons, increases quadratically with $`C`$. This energy increase can be written as $`(\kappa /2)(\lambda C^2)2\pi R_0L`$, where $`L`$ is the length of the prolate, and the factor $`\lambda =\mathrm{sin}^2\varphi =1/2`$ indicates that during the integration of the energy around a cross-section (parametrized by $`\varphi `$), only the out-of-plane component of the bending counts. Thus, the bending stiffness of a prolate-like vesicle is $`\kappa \pi R_0`$, and its persistence length $`\kappa \pi R_0/(k_\mathrm{B}T)`$ is in the order of a few hundred nanometers.

Although these prolates are thicker than the microtubules, their persistence length is much closer to that of a DNA. This is because their wall is a two-dimensional fluid. Due to this fluidity they do not even resist twisting, which makes them ideal semiflexible filaments.

We close this Letter by discussing some of the biological consequences of our theoretical results. The major biological relevance of the first force peak of the $`f`$-$`L`$ curves is that in order to form a tube, the motor proteins must be able to provide a force that is 13% larger than what is needed to pull a long tube. At a pulling force of $`f=f_0`$ the peak corresponds to a $`2.1\kappa 21k_\mathrm{B}T`$ high energy barrier, which is practically insurmountable for such a big and slow object as a growing tube. Thus, tube formation works on an all-or-nothing basis: motors can pull out tubes only if they are strong enough to overcome the major force peak.

Although at the tip of the tubes the size of the lipids ($`0.5`$ nm) represents a natural cutoff length-scale for the divergence of the energy density, it goes up to tens of $`\sigma `$. In terms of membrane rupture, the energy density acts as an effective surface tension, meaning that the most likely place for rupture to occur is the tip of the tubes. So if biological systems want to avoid rupture, they either have to protect the tips or distribute the pulling forces at larger areas. That could be achieved, e.g., by utilizing cap proteins or lipid rafts.

We have shown that without external pinning tubes coalesce smoothly. Thus, to explain bundle formation, other physical effects that could prevent the coalescence of tubes (e.g. adhesion between the tubes and the cytoskeleton) must be taken into account.

###### Acknowledgements.

We thank A. Roux, G. Cappello, P. Bassereau, and B. Goud for their helpful discussions. This work was supported by the EU Marie Curie Fellowship HPMF-CT-1999-00124.