text:
# Overview of PHENIX Results from the First RHIC Run W.A. Zajc,8 for the PHENIX Collaboration:

## 1 INTRODUCTION

The PHENIX experiment has been designed to measure a broad variety of signals from both heavy ion and polarized proton-proton collisions at RHIC. The pursuit of penetrating probes generated in the early stages of the collision, combined with a program of hadron measurements, provides a detector with unparalleled capabilities to address observables sensitive to all stages of the collision process. This same detector is also very well-suited to the study of gluon and anti-quark contributions to the proton spin.

The PHENIX detector consists of three spectrometers: two muon spectrometers covering the full azimuth for $`1.1<|\eta |<2.4`$ and a central spectrometer consisting of two arms each subtending $`90^o`$ in azimuth and with $`|\eta |<0.35`$. A central magnet provides an axial field, while each muon spectrometer contains a magnet that produces a roughly radial field. The central arms contain three tracking sub-systems: pad chambers (PC), drift chambers (DC) and time-expansion chambers (TEC); two forms of electromagnetic calorimetry (PbSc and PbGl); a time-of-flight hodoscope (TOF) and ring imaging Cerenkov counter (RICH). These sub-systems, together with a set of beam-beam counters (BBC) located in the region $`3<|\eta |<3.9`$, provide superb hadron and electron identification over a broad range of transverse momentum. The muon spectrometers use cathode strip chambers in three stations for tracking (muTr), and five layers of Iarocci tubes interleaved with iron absorber for muon identification (muID). Global event characterization is achieved via a multiplicity and vertex detector (MVD) consisting of silicon strips and pads covering $`|\eta |<2.5`$, and the RHIC-standard Zero-Degree Calorimeters (ZDCs), which detect neutral particles emitted along the beam directions. The front end electronics for all sub-systems are clocked synchronously with the beam crossing frequency of 9.4 MHz. A set of Level-1 triggers derived from various sub-systems is used to initiate readout of the entire detector through a pipelined high bandwidth data acquisition system capable of archiving 20 MB/s.

For the first physics run of RHIC in the summer of 2000, the portions of the PHENIX detector shown in Figure 1 were instrumented. Elements of all sub-systems, with the exception of the muTr, were in place and read out. Small subsets of the MVD and muID front end electronics were connected and tested as part of an engineering run. All other sub-systems were instrumented in fractions ranging from 25% to 100% of their ultimate aperture and were used in the physics results which are presented here. Independent minimum bias triggers were formed using coincidences between the BBC counters and between the ZDC counters. A total of approximately 5M events was recorded at $`\sqrt{s_{NN}}=130`$ GeV. The primary trigger used for most of the results presented below is based on the BBC coincidence with an additional offline requirement that restricts the collision vertex to $`|z|<20`$ cm. Whenever possible, physics quantities are presented as a function of centrality and/or the number of participants $`N_{part}`$. Details on the determination of these quantities are presented in the following section and in the associated references.

## 2 GLOBAL OBSERVABLES

The systematic variation of particle yields and the produced transverse energy with the number of participants reflects the underlying reaction mechanisms. For example, Gyulassy and Wang have emphasized that such a study can discriminate between cascade models and models which incorporate gluon saturation effects. PHENIX has studied the production near mid-rapidity of both charged particles and of transverse energy as a function of centrality. The deposition of energy in the ZDC’s is correlated with that of charge in the BBC’s to provide an unambiguous mapping between these observables and the centrality of the collision,

as shown in Figure 2. A Glauber model is then used to determine the number of participants $`N_P`$ and the number of binary collisions $`N_C`$ for each centrality class. (This methodology is used consistently for all such studies presented in this contribution.)

The number of produced charged particles is determined from the correlation between hits in two layers of Pad Chambers with the vertex location, which provides (on a statistical basis) the charged particle multiplicity distribution in the interval $`|\eta |<0.35`$. After corrections for acceptance, efficiencies, decays and double hits, the charged particle pseudo-rapidity distribution $`dN_{ch}/d\eta `$ is calculated for each centrality bin, scaled by the corresponding number of participant pairs, then plotted versus $`N_P`$, as shown in Figure 3a. Also shown there are comparisons to the model predictions of EKRT, which does not reproduce the trend of the data. (Results very similar to the PHENIX data were reported at this conference by the PHOBOS Collaboration.)

A similar analysis has been performed for transverse energy measured in the PHENIX PbSc calorimeter. A careful treatment of the contributions from produced energy in the aperture, the in-flux from scattering sources, and both the in-flux and out-flux from decays is performed to convert the deposited energy seen in the calorimeter to the equivalent transverse energy. The total transverse energy density per participant pair is then calculated in the same $`N_P`$ bins as used in the multiplicity analysis. The trend, shown in Figure 3b., is essentially identical to that found for charged multiplicity, and indicates that particle and transverse energy production is not simply proportional to the number of participants. The yield may be described as a superposition of terms proportional to participants and binary collisions,

$$\mathrm{Yield}=AN_P+BN_C.$$
(1)
The results of this procedure are given in Table 1. The extent to which the ratio $`B/A`$ is significantly different from zero may be taken as evidence for the role of binary collisions in contributing to production of particles and transverse energy at RHIC.

Also shown in Figure 3 are comparisons of both the charged multiplicity and the transverse energy per participant pair to the distributions measured by the WA98 collaboration at the CERN SPS, clearly demonstrating the increase of both quantities at RHIC. In the case of energy density calculated following the Bjorken prescription, the value of 4.6 $`\mathrm{GeV}/\mathrm{fm}^3`$ is roughly 60% larger than that found at the SPS. It is interesting to note that this increase does not result from an increase in $`dE_T/d\eta /dN_{ch}/d\eta `$ ; the value is very nearly equal to that found at the SPS.

Another global feature of hadronic production is the azimuthal pattern of emission, as parameterized by elliptic flow. Typically such analyses are performed with respect to the reaction plane determined for each event. It is also possible to determine the elliptic flow pattern by measuring the auto-correlation function of particles in the event, that is, by calculating the correlation function $`C(\mathrm{\Delta }\varphi )`$, where $`\mathrm{\Delta }\varphi `$ is the difference in azimuthal angle between two particles from the same event. PHENIX has performed such an analysis using charged particles with $`p_T>200`$ MeV/c in the interval $`|\eta |<0.35`$. The correlation function is calculated using an event-mixing prescription to determine the background distribution, then fit to extract the (assumed positive) $`v_2`$ coefficient:

$$C(\mathrm{\Delta }\varphi )1+2|v_1|^2\mathrm{cos}(\mathrm{\Delta }\varphi )+2|v_2|^2\mathrm{cos}(2\mathrm{\Delta }\varphi )$$
(2)
Standard PHENIX event characterization methods are used to study the dependence of $`v_2`$ versus centrality in each $`p_T`$ bin. A sample of the results are shown in Figure 4. These trends, which are in good agreement with data from STAR and preliminary results from PHOBOS, suggest that the high density matter formed at RHIC efficiently translates the initial spatial asymmetry into a corresponding one in momentum space.

## 3 IDENTIFIED HADRONS

The PHENIX central arm detectors have been designed to provide particle identification over the broadest possible momentum range. The primary tool for charged hadron identification is the time-of-flight difference between the BBC and the highly-segmented TOF hodoscope, which spans $`\mathrm{\Delta }\varphi =45^o`$ in the East spectrometer arm. The overall time resolution $`\sigma 115`$ ps permits unambiguous $`\pi /K`$ separation to at least $`p_T=1.5`$ GeV/c. Normalized minimum bias $`p_T`$ spectra for $`\pi ^\pm `$, $`\mathrm{K}^\pm `$, $`p`$’s and $`\overline{p}`$’s are shown in Figure 5. The shape of the spectra clearly depends on the particle species, with pions having the lowest $`p_T`$ and protons and anti-protons having the largest. The dependence of the local slopes increases with both centrality and with particle mass, consistent with expectations from radial flow.

The ratio of anti-protons to protons at $`y=0`$ is of particular interest, since it is a direct measure of the net baryon content in the central region. The $`\overline{p}/p`$ ratio has been studied as a function of both transverse momentum and centrality, and found to be only weakly dependent on $`p_T`$ and independent of centrality within systematic errors (Figure 6). For minimum bias collisions, the $`\overline{p}/p`$ ratio in the interval 0.8 GeV/c $`<p_T<`$ 3.0 GeV/c is 0.64 $`\pm `$ 0.01(stat.) $`\pm `$ 0.07(sys.). This value (which is consistent with that observed by the other RHIC experiments) together with the measured yields shown in Figure 5 signifies that the central region in heavy ion collisions at RHIC is meson-dominated, as distinct from the baryon-dominated case at the SPS and lower energies, where the ratio is never greater than $`0.1`$.

Charged pion pairs in the central arms have also been used to perform an HBT analysis. Two separate analyses were performed, one using pions identified with the TOF as described above, and the second in which the 700 ps time-of-flight resolution of the PbSc EmCal (EMC) is used to identify pions with transverse momentum $`p_T<0.7`$ GeV/c. In addition to the obvious utility of comparing two independent data sets for this analysis, the larger acceptance of the EMC substantially increases the number of available pairs (by a factor of 5). Correlation functions are calculated in the standard $`q_{Tside},q_{TOut},q_{Long}`$ projections of the relative momentum for both $`\pi ^+\pi ^+`$ and $`\pi ^{}\pi ^{}`$ pairs and are then fit over the full 3D phase space to the form

$$C_2(q_{Tside},q_{TOut},q_{Long})=1+\lambda \mathrm{exp}[(q_{Tside}R_{Side})^2(q_{TOut}R_{Out})^2(q_{Long}R_{Long})^2].$$
(3)
The results, presented in Figure 7 and Table 2, show little if any variation from values obtained at lower energies. Given the significantly higher multiplicities and densities observed at RHIC, this is somewhat puzzling, and perhaps indicates significant dynamic effects on the radii from strong transverse expansion driven by these higher densities. In this context, it should be noted that the values of the radii and $`\lambda `$ (for the TOF analysis), as well as the dependence on pair transverse momentum are all consistent with measurements from the STAR collaboration.

## 4 PARTICLE PRODUCTION AT HIGH TRANSVERSE MOMENTUM

The high center-of-mass energy provided by RHIC, corresponding to values of $`\sqrt{s_{NN}}`$ where hard scattering at the partonic level is observed in $`p`$-$`p`$ and $`p`$-$`\overline{p}`$ collisions, offers the exciting possibility of using perturbative probes amenable to quantitative calculation to explore hot nuclear matter. A first step in this program is the measurement of the transverse momentum spectrum for charged particles and its variation with the number of participants in the collision. PHENIX has performed such an analysis using charged tracks reconstructed with the drift and pad chambers of the central arm. The raw distribution of tracks is formed for six exclusive centrality classes, again using the standard PHENIX event classification described in Section 2. Corrections are made in each centrality bin for the spectral distortions due to momentum resolution and to backgrounds from scattering and decays.

Results are shown in Figure 9 for the interval $`0.4<p_T<5`$ GeV/c, which spans six decades in yield.

In nucleon-nucleon collisions, production of high transverse momentum particles results from hard scattering at the partonic level. Since these cross sections are small, their contribution in nucleus-nucleus collisions, absent collective effects, should scale as the total number of binary collisions. For the most peripheral (80-92% in centrality) bin shown in Figure 9, the shape and yield above $`p_T>2`$ GeV/c are in good agreement with an interpolation of hadron-hadron data scaled by the number of binary collisions ($`N_{coll}=3.7\pm 2`$) calculated for this centrality bin. This is not the case for the most central collisions, which show a deficit at large $`p_T`$ with respect to $`N_{coll}`$ scaling. This is presented in Figure 9, where the scaled ratio

$$\frac{Yield(Central)/N_{coll}(Central)}{Yield(Peripheral)/N_{coll}(Peripheral)}$$
(4)
is plotted for two different values of the peripheral reference set. The suppression at $`p_T>3`$ GeV/c is inconsistent both with enhancements expected in that region for the Cronin effect and also seen in PHENIX data for mid-peripheral data when scaled to the p-p reference distribution. Similar results for such a suppression pattern in unidentified charged particles have been reported by the STAR collaboration.

An independent analysis by PHENIX provides further insight into the nature of the suppression. The superb segmentation $`\mathrm{\Delta }\eta \times \mathrm{\Delta }\varphi =0.01\times 0.01`$ and excellent resolution $`8.2\%/\sqrt{E(\mathrm{GeV})}1.9\%`$ of the PbSc electromagnetic calorimeter are used to extract the $`\pi ^0`$ $`p_T`$ spectrum by reconstruction (on a statistical basis) of their principal $`\pi ^0\gamma \gamma `$ decay mode. Extensive studies based on the mixing of single photon showers into both real and simulated events are used to determine the variation of the reconstruction efficiency with multiplicity and of the background with the event multiplicity.

A detailed simulation is used to estimate the contributions from background and decay particles. Results are shown in Figure 11 for both a peripheral (60-80%) and central (0-10%) sample. As in the unidentified charged particle analysis, the peripheral data are well described with an interpolation of existing hadron-hadron data scaled by the number of binary collisions. (An additional factor derived from p-p data of $`\pi ^0/(h^++h^{})=1/3.2`$ is used to scale the measured yield of unidentified charged particles to neutral pions.) Once again, the central data fall well below the corresponding yield expected if high $`p_T`$ particle production scaled with the number of binary collisions. This is strikingly illustrated by the scaled ratio of Equation 4 for the central to peripheral yields, as shown in Figure 11.

A comparison of the ratio in Figure 9 for unidentified charged particles to that of Figure 11 for $`\pi ^0`$’s suggests that the suppression may be more pronounced in the case of identified (neutral) pions. Should this be established (the present understanding of the systematic errors does not permit a definite conclusion) a consistent description would imply that charged pions are also preferentially suppressed as compared to the unidentified charged particles. Hints of precisely this behavior are seen in Figure 5, where the trends in identified particle spectra suggest that the contributions from protons and anti-protons become comparable to those from pions for $`p_T>2`$ GeV/c. Further work is required to improve both statistical and especially systematic errors to determine if the agreement of these data with “jet quenching” predictions is in fact evidence for an enhanced energy loss mechanism in hot nuclear matter. Clearly essential to that program is a set of detailed measurements of both p-p reference data and proton-nucleus collisions to determine the quantitative value of the Cronin effect at RHIC energies.

## 5 FUTURE PLANS

A major component of the PHENIX physics program is dedicated to measurement of leptonic signals. For example, the $`\pi ^0`$ spectra described in Section 4, while intriguing and valuable in themselves, are also a prelude to the measurement of direct photons. In the central arms, virtual photons and vector mesons can be detected via their decays to $`e^+e^{}`$ pairs, while open charm and bottom production are expected to dominate the production of single electrons for $`p_T>`$2 GeV/c. Measurement of these signals requires $`\pi /e`$ rejection in excess of $`10^3`$ and very careful control of background contributions. Figure 12 shows the first efforts in this program. Charged particles are tracked using the PHENIX drift and pad chambers. Electrons are selected by requiring that at least three PMT’s fire in the RICH, and are further identified via a tight matching cut between tracking momentum and energy deposition in the EmCal.

The power of this combined particle identification from the various PHENIX sub-systems is apparent from Figure 12, where a clean electron spectrum is extracted at a level 2-3 orders of magnitude below that of all charged particles. This approach will form the basis for future PHENIX measurements of vector mesons and open charm at RHIC, and is also applicable to photon measurements via external conversion. The measurement of electrons along with all other PHENIX analyses will benefit greatly from the significantly enhancements made to the detector for Run-2 at RHIC. As shown in Figure 13,

the complete aperture of the central arms will be available (compare to Figure 1), much of the MVD will be instrumented, and an entirely new spectrometer to measure muons will be deployed. These additions in aperture and in capability, coupled with significant upgrades to the data acquisition and triggering system, should result in a hundredfold or more increase in the event sample obtained from Run-2. This will allow PHENIX to explore new signals and observables, to greatly increase statistical precision and understanding of systematic errors on existing analyses, to obtain the vital p-p comparison data, and to begin a program of measurements with polarized protons dedicated to understanding the proton spin.

## 6 ACKNOWLEDGMENTS

We thank the staff of the RHIC project, Collider-Accelerator, and Physics Departments at BNL and the staff of PHENIX participating institutions for their vital contributions. We gratefully acknowledge support from the Department of Energy and NSF (U.S.A.), Monbu-sho and STA (Japan), RAS, RMAE, and RMS (Russia), BMBF and DAAD (Germany), FRN, NFR, and the Wallenberg Foundation (Sweden), MIST and NSERC (Canada), CNPq and FAPESP (Brazil), IN2P3/CNRS (France), DAE (India), KRF and KOSEF (Korea), and the US-Israel Binational Science Foundation.


format: text

text:
# An Evolutionary Scenario for Blazar Unification

## 1 Introduction

During the past decade, significant progress has been made in understanding the physics of active galactic nuclei (AGN) and the connection between different subclasses of AGN. In particular, various lines of evidence indicate that FR II radio galaxies are the parent population of radio-loud quasars, with the relativistic jets of FSRQs being closely aligned to our line-of-sight (for a review, see Urry & Padovani (1995)). Orientation effects are also thought to unify FR I radio galaxies with BL Lac objects. The BL-Lac/FR-I subclasses have lower average luminosities than the FSRQ/FR-II subclasses. Among blazars, which include highly polarized and optically violently variable quasars, FSRQs and BL Lac objects, there appears to be an almost continuous sequence of properties from FSRQs through LBLs to HBLs. This trend is characterized by decreasing bolometric luminosities, a shift of the peak frequencies of their broadband spectral components towards higher values, and a decreasing fraction of power in $`\gamma `$ rays compared with lower-frequency radiation (Sambruna, Maraschi, & Urry, 1996; Fossati et al., 1998).

In the framework of relativistic jet models, the low-frequency (radio – optical/UV) emission from blazars is interpreted as synchrotron emission from nonthermal electrons in a relativistic jet. The high-frequency (X-ray – $`\gamma `$-ray) emission is thought to be produced via Compton upscattering of low frequency radiation by the same electrons responsible for the synchrotron emission. Possible sources of soft seed photons for Compton scattering are the synchrotron photons themselves (Maraschi, Ghisellini, & Celotti, 1992; Bloom & Marscher, 1996) or external photons, presumably dominated by the accretion disk emission which can enter the jet either directly (Dermer, Schlickeiser, & Mastichiadis, 1992; Dermer & Schlickeiser, 1993) or after reprocessing by circumnuclear gas and dust (Sikora, Begelman, & Rees, 1994; Dermer, Sturner, & Schlickeiser, 1997). In addition, a significant contribution to the soft radiation field may also be provided by infrared emission from dust in the vicinity of the AGN (Bla$`\dot{\mathrm{z}}`$ejowski et al., 2000; Arbeiter, Pohl, & Schlickeiser, 2001).

The blazar sequence has been studied by Ghisellini et al. (1998) using a large sample of blazar broadband spectra. They suggest that along the sequence HBL $``$ LBL $``$ FSRQ, an increasing energy density of the external radiation field leads to an increasing amount of Compton cooling. As a result, the maximum energy in the electron distribution decreases, causing the synchrotron and Compton peaks to shift to lower frequencies. This explanation for the blazar sequence was developed further by Georganopoulos, Kirk, & Mastichiadis (2001), who argue that the radiating jet plasma is outside the broad-line scattering region in weak sources and within it in powerful sources. Detailed model fits to the broadband spectra of several blazars also indicate that a decreasing contribution of the external radiation to the seed photon field reproduces the FSRQ $``$ LBL $``$ HBL sequence. While FSRQs generally require a dominant external Compton contribution to produce the observed $`\gamma `$-ray emission (e.g., Sambruna et al. (1997); Mukherjee et al. (1999); Hartman et al. (2001)), HBLs can be successfully fitted with pure synchrotron self-Compton (SSC) models (e.g., Mastichiadis & Kirk (1997); Pian et al. (1998); Petry et al. (2000)). The LBL BL Lacertae appears to be intermediate between those two classes, requiring a non-negligible contribution from external Comptonization to reproduce its $`\gamma `$-ray spectrum (Madejski et al., 1999; Böttcher & Bloom, 2000). Orientation effects certainly play a role in spectral variations among sources within a given subclass, or even to account for the blue blazar subclass (Georganopoulos, 2000), but cannot explain the full range of blazar properties (Sambruna, Maraschi, & Urry, 1996).

It is uncertain whether the different blazar types are connected through an evolutionary sequence, or whether they constitute limited periods of AGN activity on parallel evolutionary paths. An evolutionary sequence has recently been proposed by D’Elia & Cavaliere (2000) and Cavaliere & D’Elia (2001) who consider the efficiency of energy extraction from the accretion flow onto the central supermassive black hole through magnetic fields anchored in an optically thick accretion disk vs. energy extraction from rotational energy of the black hole via the Blandford-Znajek (BZ) mechanism (Blandford & Znajek, 1977). They argue that BL Lac objects can plausibly be powered by the BZ mechanism without a significant contribution from an accretion disk, while FSRQs require a dominant contribution from the accretion disk to power the relativistic jets responsible for the production of the nonthermal broadband spectra. They suggest that a gradual depletion of the central environment in FSRQs by accretion onto the central black hole may be responsible for a transition from the disk-powered jet production (in FSRQs) to the BZ-powered jet production mode (in BL Lac objects). Furthermore, they find that a cosmological BL Lac birth rate proportional to the rate of “deaths” of FSRQs (through depletion of their central environments) is consistent with the cosmological evolution observed for both quasars and BL Lac objects. The masses of the central black holes in radio quasars generally exceed $`3\times 10^8M_{}`$ (Wandel, 1999; Laor, 2000). Similarly detailed mass determinations have so far not been possible for the lineless BL Lacertae objects. The cluster environments of FR II radio galaxies and radio quasars on the one hand, and FR I radio galaxies and BL Lacs on the other hand appear consistent with an evolution of the former into the latter class of objects (Hill & Lilly, 1991; Yee & Ellingson, 1993; Urry & Padovani, 1995). Also, the distribution of central black-hole masses in nearby non-active galaxies appears to extend to larger values than the distribution of radio-quasar black-hole masses (Laor, 2000), consistent with some fraction of present-day giant elliptical galaxies hosting the remnants of earlier quasar activity.

Here we expand on the suggestion of a cosmological evolutionary sequence (D’Elia & Cavaliere, 2000; Cavaliere & D’Elia, 2001) from FSRQs to BL Lac objects by demonstrating that such a sequence plausibly and self-consistently explains the trend in their spectral properties. We employ a combined SSC/external Compton model for relativistic jet emission, as described in §2, and derive parameter values from the observed multifrequency signatures of FSRQs and HBLs. The evolutionary sequence connecting the respective parameter regimes is presented in §3. We discuss our results and suggest a possible evolutionary connection between radio-loud and radio-quiet AGN classes in §4.

## 2 Blazar Model and Parameter Estimation

The radiating plasma is modeled by a spherical emission region of co-moving radius $`R_b`$ that moves outward along the axis of the jet with bulk Lorentz factor $`\mathrm{\Gamma }=(1\beta _\mathrm{\Gamma }^2)^{1/2}`$. Electrons are injected into the plasma blob with a comoving power-law distribution $`Q(\gamma )=Q_0\gamma ^s`$, for electron Lorentz factors in the range $`\gamma _1\gamma \gamma _2`$. The jet power $`L_j`$ corresponding to the energy input of injected particles into the jet is assumed to be proportional to the accretion disk power $`L_D`$, as suggested by the observed proportionality between kinetic jet power and narrow line emission luminosity in radio-loud AGN reported by Rawling & Saunders (1991) and the proportionality between core radio luminosity and broad emission line luminosity in radio-loud AGN found by Celotti, Padovani & Ghisellini (1997). In addition, Xu, Livio & Baum (1999) have found a similar proportionality between the 5 GHz radio luminosity and the \[O III\] 5007 emission line luminosities both for radio-loud and radio-quiet AGN. The jet power is related to the injected particle distribution through $`L_j=f_jL_D=m_ec^2_{\gamma _1}^{\gamma _2}𝑑\gamma \gamma Q(\gamma )`$.

Electrons lose energy through adiabatic expansion, which is treated in terms of an escape time scale $`t_{\mathrm{esc}}\eta R_b/c`$ with $`\eta 1`$, and through radiative losses. These losses are generally dominated by synchrotron and Thomson processes, so they can be parameterized by the functional form $`\dot{\gamma }_{\mathrm{rad}}=\nu _0\gamma ^2`$, where $`\nu _0=4c\sigma _T\left(u_B+u_{\mathrm{sy}}+u_{\mathrm{ext}}\right)/3m_ec^2`$. The terms $`u_B`$, $`u_{\mathrm{sy}}`$, and $`u_{\mathrm{ext}}`$ represent the comoving magnetic field, synchrotron, and external radiation field energy densities, respectively. If $`u_{\mathrm{ext}}`$ is dominated by reprocessed accretion-disk radiation due to circumnuclear gas and dust with a scattering optical depth $`\tau _{\mathrm{repr}}`$ at a characteristic distance $`R_{\mathrm{sc}}`$ from the central source, then $`u_{\mathrm{ext}}\mathrm{\Gamma }^2L_D\tau _{\mathrm{repr}}/(4\pi R_{\mathrm{sc}}^2c)`$. Our approach is similar to the one used by Ghisellini et al. (1998), but differs from their model by our explicitly accounting for adiabatic losses and escape. These processes lead to the formation of a low-energy cutoff in the electron spectrum and force us to consider two distinct radiative regimes of electron cooling and of the equilibrium electron spectra (see below), while the Ghisellini et al. (1998) approach always leads to a $`\gamma ^2`$ electron spectrum for $`1\gamma \gamma _1`$. We neglect pair processes, which generally turn out to be unimportant in the course of detailed spectral modeling of blazars.

Electrons injected with $`\gamma \gamma _{\mathrm{cr}}c/(R_b\nu _0)`$ cool rapidly prior to escape, whereas electrons injected with $`\gamma \gamma _{cr}`$ lose only a small fraction of their energy before they escape. The competition between electron injection, cooling and escape establishes an equilibrium distribution given by

$$n_{e,eq}(\gamma )=n_{e,eq}^0\{\begin{array}{cc}(\gamma /\gamma _m)^p\hfill & \text{for }\gamma _l\gamma \gamma _m\hfill \\ (\gamma /\gamma _m)^{(s+1)}\hfill & \text{for }\gamma _\mathrm{m}\gamma \gamma _2\hfill \end{array}$$
(1)
Two distinct radiative regimes can be defined. If $`\gamma _1<\gamma _{\mathrm{cr}}`$, most of the injected electrons cool inefficiently and the electron population is in the slow cooling regime. In this case, $`p=s`$, $`\gamma _l=\gamma _1`$, $`\gamma _m=\gamma _{cr}`$, and $`n_{e,eq}^0=Q_0t_{\mathrm{esc}}\gamma _m^s/V_b`$, where $`V_b`$ is the blob volume. When $`\gamma _1>\gamma _{\mathrm{cr}}`$, the electron population is in the fast cooling regime, $`p=2`$, $`\gamma _l=\gamma _{cr}`$, $`\gamma _m=\gamma _1`$, and $`n_{e,eq}^0=Q_0\gamma _m^{(1+s)}/[V_b\nu _0(s1)]`$. The synchrotron, and external-Compton radiation emitted by these relativistic electron populations are evaluated using the formulae given by Dermer, Sturner, & Schlickeiser (1997). The SSC radiation is calculated using the analytic solution of Tavecchio, Maraschi, & Ghisellini (1998).

In the absence of an external Compton component, we recover the SSC model used to fit data from HBLs such as Mrk 501, Mrk 421, or PKS 2155-304 (e.g., Mastichiadis & Kirk (1997); Pian et al. (1998); Tavecchio, Maraschi, & Ghisellini (1998); Bednarek & Protheroe (1999); Petry et al. (2000); Sambruna et al. (2000); Fossati et al. (2000); Kataoka et al. (2000)). The parameters in an SSC model are the Doppler factor $`D=\left[\mathrm{\Gamma }(1\beta _\mathrm{\Gamma }\mathrm{cos}\theta _{\mathrm{obs}})\right]^1`$, the magnetic field $`B`$, $`R_b`$, $`\gamma _1`$, $`\gamma _2`$, the injection luminosity $`L_j`$, and the injection spectral index $`s`$. In principal, the escape time scale parameter $`\eta `$ could also be considered a free parameter, though we assume (see Kataoka et al. (2000)) that $`\eta 1`$, i.e., $`t_{\mathrm{esc}}R_b/c`$. These studies find typical HBL parameters of $`R_b10^{15}`$$`10^{16}`$ cm, $`\gamma _210^6`$, $`B0.1`$ G, $`D`$ 20 – 30, and a rather broad range of injection luminosities $`L_j10^{37}`$ ergs s-1$`10^{41}`$ erg s-1, depending on the various model details and the individual flares that are modeled.

The addition of the external Compton component in FSRQs complicates parameter estimation. However, since detailed modeling of blazar spectra using self-consistent, combined SSC + ERC spectral models is generally a rather time-consuming task, simple formulae to achieve order-of-magnitude estimates for the relevant model parameters are desirable in order to find a convenient starting point for more detailed modeling efforts. In the following, we derive such estimates, although we must caution that some of the observables, in particular those pertaining to the SSC component, are rather uncertain, which may introduce significant scatter in the derived analytical parameter estimates. We assume that the X-ray spectra of FSRQs are dominated by the SSC process, as has been argued by various authors on the basis of spectral fitting as well as variability considerations (e.g., Mukherjee et al. (1999); Hartman et al. (2001); Sikora et al. (2001)), while the $`\gamma `$-ray spectrum is dominated by the external Compton component. The peak photon energies and peak $`\nu F_\nu `$ fluxes, $`\nu F_\nu =10^{10}F_{10}`$ ergs cm-2 s-1, of the Compton and synchrotron components are determined from simultaneous multiwavelength observations. In the following, photon energies are given in dimensionless units $`ϵ=h\nu /m_ec^2=10^nϵ_n`$, and we define $`ϵ_BB/B_{\mathrm{cr}}=2.3\times 10^{14}B(\mathrm{G})`$. The jets in blazars are generally believed to be oriented close to the line of sight where the Doppler beaming factor $`D\mathrm{\Gamma }`$. Thus, for simplicity, we assume $`D=\mathrm{\Gamma }`$ in the following. The model parameters can then be estimated from the observables through the expressions

$$ϵ_{\mathrm{sy}}=\frac{Dϵ_B\gamma _{\mathrm{cr}}^2}{1+z},$$
(2)
$$ϵ_{\mathrm{ERC}}=\frac{D^2ϵ_{}\gamma _{\mathrm{cr}}^2}{1+z},$$
(3)
$$ϵ_{\mathrm{SSC}}=\frac{Dϵ_B\gamma _{\mathrm{cr}}^4}{1+z},$$
(4)
$$\left(\nu F_\nu \right)_{\mathrm{sy}}=\frac{V_bD^4}{4\pi d_L^2}\frac{4}{3}c\sigma _Tu_B\gamma _{\mathrm{cr}}^2n_ef_{\mathrm{sp}}f_{\mathrm{sy}},$$
(5)
$$\left(\nu F_\nu \right)_{\mathrm{ERC}}=\frac{V_bD^4}{4\pi d_L^2}\frac{4}{3}c\sigma _Tu_{\mathrm{ext}}\gamma _{\mathrm{cr}}^2n_ef_{\mathrm{sp}}f_{\mathrm{ERC}},\mathrm{and}$$
(6)
$$\left(\nu F_\nu \right)_{\mathrm{SSC}}=R_b\frac{V_bD^4}{4\pi d_L^2}\frac{4}{3}c\sigma _T^2u_B\gamma _{\mathrm{cr}}^4\left(n_ef_{\mathrm{sp}}\right)^2f_{\mathrm{SSC}}.$$
(7)
The subscripts “sy,” “SSC,” and “ERC” refer to the synchrotron, SSC, and external Compton components, respectively, and $`ϵ_{}`$ is the mean photon energy of the external soft photon field in the stationary frame of the AGN. The factor $`f_{\mathrm{sp}}`$ in equations (5) – (7) is a normalization factor defined through $`_1^{\mathrm{}}𝑑\gamma n_e(\gamma )\gamma ^2=\gamma _{\mathrm{cr}}^2n_ef_{\mathrm{sp}}`$, and $`f_{\mathrm{SSC}}`$ in equation (7) is a correction factor between the $`\nu F_\nu `$ peak value and the total energy output in the SSC component, which differ significantly due to the substantial spectral broadness of the SSC emission. Typically, $`f_{\mathrm{SSC}}0.1`$. For the more strongly peaked synchrotron and ERC components, the corresponding correction factor is $`f_{\mathrm{sy}}f_{\mathrm{ERC}}1`$.

Combining the above estimates, and defining $`d_L=10^{28}d_{28}`$ cm and $`R_{\mathrm{sc}}=10^{18}R_{\mathrm{sc},18}`$ cm, we find

$$D=3.2\sqrt{\frac{ϵ_{\mathrm{ERC},2}ϵ_{\mathrm{sy},7}(1+z)}{ϵ_{,5}ϵ_{\mathrm{SSC},1}}},$$
(8)
$$B=1.4\sqrt{\frac{ϵ_{\mathrm{sy},7}^3ϵ_{,5}(1+z)}{ϵ_{\mathrm{SSC},1}ϵ_{\mathrm{ERC},2}}},\mathrm{G}$$
(9)
$$L_D\tau _{\mathrm{repr}}=3\times 10^{45}\frac{F_{\mathrm{ERC},10}}{F_{\mathrm{sy},10}}\frac{(1+z)R_{\mathrm{sc},18}^2ϵ_{\mathrm{sy},7}^2ϵ_{,5}^2}{ϵ_{\mathrm{ERC},2}^2}\mathrm{erg}\mathrm{s}^1,$$
(10)
$$R_{\mathrm{esc}}=1.2\times 10^{16}\frac{ϵ_{\mathrm{ERC},2}}{(1+z)ϵ_{,5}}\sqrt{\frac{ϵ_{\mathrm{SSC},1}}{ϵ_{\mathrm{sy},7}^5}}\left(1+\frac{F_{\mathrm{ERC},10}}{F_{\mathrm{sy},10}}\right)^1\mathrm{cm},$$
(11)
$$R_b=10^{17}\frac{d_{28}}{\sqrt{1+z}}\sqrt{\frac{F_{\mathrm{sy},10}^2}{F_{\mathrm{SSC},10}}}\sqrt{\frac{ϵ_{,5}ϵ_{\mathrm{SSC},1}^3}{ϵ_{\mathrm{ERC},2}ϵ_{\mathrm{sy},7}^5}}\sqrt{\frac{(f_{\mathrm{SSC}}/0.1)}{f_{\mathrm{sy}}^2}}\mathrm{cm},\mathrm{and}$$
(12)
$$f_j=4(3s)\frac{d_L^2\left(\nu F_\nu \right)_{\mathrm{sy}}}{L_D\gamma _1^{s2}f_{\mathrm{sy}}}\left(\frac{ϵ_{\mathrm{sy}}}{ϵ_{\mathrm{SSC}}}\right)^{(2+s)/2}\left(\frac{ϵ_{}}{ϵ_{\mathrm{ERC}}}\right)^2\left(1+\frac{\left(\nu F_\nu \right)_{\mathrm{ERC}}}{\left(\nu F_\nu \right)_{\mathrm{sy}}}\right).$$
(13)
Note that equation (13) requires an independent estimate of $`\gamma _1`$. Assuming that the optical spectrum is produced through synchrotron emission from strongly cooled electrons, the injection spectral index $`s`$ can be estimated through the optical spectral index $`\alpha _{\mathrm{opt}}`$ from the relation $`s=2\alpha _{\mathrm{opt}}`$.

The requirement that $`\eta =R_{\mathrm{esc}}/R_b1`$ yields a constraint on the external soft photon energy $`ϵ_{}`$ from equations (11) and (12), given by

$$ϵ_{}\frac{2.4\times 10^6}{\left(\sqrt{1+z}d_{28}\right)^{2/3}}\frac{ϵ_{\mathrm{ERC},2}}{ϵ_{\mathrm{SSC},1}^{2/3}}\left(\frac{\sqrt{F_{\mathrm{SSC},10}}}{F_{\mathrm{sy},10}+F_{\mathrm{ERC},10}}\right)^{2/3}\left(\frac{f_{\mathrm{sy}}}{\sqrt{(fssc/10)}}\right)^{2/3}.$$
(14)
However, we point out that the observables $`ϵ_{\mathrm{SSC}}`$ and $`\left(\nu F_\nu \right)_{\mathrm{SSC}}`$ are generally uncertain by at least an order of magnitude since the SSC component does generally not show up as a pronounced bump in the broadband spectra of FSRQs. Also the ERC peak energy $`ϵ_{\mathrm{ERC}}`$ can rarely be determined to an accuracy of better than an order of magnitude. Thus, great caution must be used when drawing conclusions from Eq. (14).

For the example of 3C 279, one of the best-observed EGRET-detected FSRQs, we can estimate typical values of the observables as $`ϵ_{\mathrm{sy}}3\times 10^7`$, $`ϵ_{\mathrm{ERC}}300`$, $`ϵ_{\mathrm{SSC}}0.1`$, $`F_{\mathrm{sy},10}1`$, $`F_{\mathrm{ERC},10}5`$, and $`F_{\mathrm{SSC},10}0.5`$ (see, e.g., Hartman et al. (2001)). For 3C 279, $`z=0.538`$ and $`d_{28}=0.77`$, and we take $`\tau _{\mathrm{repr}}=0.1`$, $`\gamma _1=100`$, and $`s=2.2`$. Equation (14) then implies $`ϵ_{}\stackrel{<}{}1.8\times 10^6`$ assuming $`\eta =1`$ (corresponding to $`E1`$ eV or $`\nu 2\times 10^{14}`$ Hz). Possible implications of this soft seed photon energy estimate will be discussed in §3. Other model parameter estimates based on these observables are $`D=28`$, $`B=2.2`$ G, $`R_{\mathrm{sc}}=1.2\times 10^{18}\sqrt{L_{46}(\tau _{\mathrm{repr}}/0.1)}`$ cm, $`R_b=1.7\times 10^{15}`$ cm, $`f_j=6.4\times 10^4L_{46}^1(\gamma _1/100)^{0.2}`$, implying $`L_j6.4\times 10^{42}(\gamma _1/100)^{0.2}`$ ergs s-1. The high-energy cut-off of both the synchrotron and the Compton components are poorly known in 3C 279 and other FSRQs, so no reliable constraint on the maximum electron energy $`\gamma _2`$ can be derived in the FSRQ case.

The parameter estimates for the FSRQ case quoted above are in good agreement with those found by Ghisellini et al. (1998) for FSRQs in general, and through more detailed spectral modeling of 3C 279 in particular (e.g., Hartman et al. (2001)), which indicates that our analytic estimates may indeed be useful to assess a convenient starting point for more detailed spectral modeling of FSRQs.

## 3 Evolutionary Sequence from FSRQs to HBLs

In the previous section, we obtained typical parameter values for HBLs and FSRQs. These estimates indicate that the most significant differences between model parameters of FSRQs and HBLs are found in the values of the jet power $`L_j`$ and the magnetic field $`B`$, which are both found to decrease by more than an order of magnitude from FSRQs to HBLs. The lack of evidence for a significant contribution of an external Compton component to the $`\gamma `$-ray spectra of HBLs indicates that the energy density of the external soft photon field — parameterized through $`L_D`$, $`\tau _{\mathrm{repr}}`$, and $`R_{\mathrm{sc}}`$ — decreases by $`\stackrel{>}{}3`$ orders of magnitude along the sequence of blazar types. This is in agreement with the results of Ghisellini et al. (1998).

The decreasing energy density of the external soft photon field along the sequence FSRQ $``$ LBL $``$ HBL, in tandem with the decreasing jet power suggests that this is a consequence of a decreasing accretion rate along this sequence. This, in turn, is in very good agreement with the suggestion of D’Elia & Cavaliere (2000) and Cavaliere & D’Elia (2001) of an evolutionary sequence from FSRQs to BL Lacs characterized by such a decline in the accretion rate. We are thus motivated to try to reconstruct the evolutionary sequence FSRQ $``$ LBL $``$ HBL by varying the accretion rate (and other model parameters dependent on the accretion rate) in our specific spectral blazar model described in the previous section. In order to do so, we are starting out with parameters appropriate to reproduce the broadband spectrum of a typical FSRQ. The FSRQ phase constitutes an early phase of blazar evolution in which the central regions of the galaxy are rich in gas and dust, leading to a high accretion rate onto the central, supermassive black hole. At the same time, the circumnuclear material efficiently reprocesses and scatters the accretion-disk radiation, leading to the observed strong optical emission lines in the broad line region and to a high energy density of the external soft photon field in the jet.

Due to the limited supply of matter in the nuclear region, the average density of the circumnuclear material will gradually decrease, leading to a decreasing accretion rate and a decreasing reprocessing efficiency. We parametrize this evolutionary transition by a gradual decline of the accretion rate and assume that the reprocessing efficiency $`\tau _{\mathrm{repr}}`$ is proportional to the accretion rate. This should be a reasonable approximation for the long-term evolution on cosmological time scales, although the actual BLR density and the instantaneous accretion rate may be less strictly correlated on shorter time scales. Furthermore, we assume that the fraction of accreted energy channeled into the relativistic outflow, $`f_j`$, remains constant, as suggested by the observed $`L_j`$ vs. $`L_{\mathrm{BLR}}`$ and $`L_j`$ vs. $`L_{\mathrm{NLR}}`$ correlations (Rawling & Saunders, 1991; Celotti, Padovani & Ghisellini, 1997; Xu, Livio & Baum, 1999), and that the transverse extent of the jet, $`R_b`$, does not change dramatically during this evolution. The magnetic field is chosen to be a constant fraction of the equipartition magnetic field so that $`u_B=0.1u_e`$, which automatically leads to the required reduction of the magnetic field along the evolutionary sequence.

The resulting sequence of broadband spectra, shown in Fig. 1, provides a quantitative explanation for the observed trend of luminosities and peak photon energies in the FSRQ $``$ LBL $``$ HBL sequence (Fossati et al., 1998). It does not account for the spectral properties of specific blazars, which depend on variations of parameter values in different sources and, importantly, the angle between the direction of the jet axis and the direction to the observer. We note that the apparent minimum at $`\nu 10^{20}`$ Hz in the FSRQ spectrum (corresponding to the $`\tau _{\mathrm{repr}}=0.1`$ curve in Fig. 1) might be filled in by an additional contribution from Compton-scattered accretion disk radiation that enters the jet directly, which has been neglected in our semi-analytic model.

Our estimate (Eq. ) for the average energy $`ϵ_{}2\times 10^6`$ of the dominant external soft photon field has interesting implications. Although one has to keep in mind that the analytical estimate for $`ϵ_{}`$ may be rather uncertain due to the uncertainty in the observables pertaining to the SSC component, we found that significantly higher values of $`ϵ_{}`$ did not allow an acceptable representation of a typical FSRQ spectrum with our model. This probably indicates that the external soft photon field is dominated by optical line emission from the partially ionized broad-line regions. The external photon field peaking at $`ϵ_{}`$ is unlikely to be due to Thomson scattering off free electrons in a highly ionized circumnuclear environment, since in that case $`ϵ_{}`$ should reflect the peak energy of the optically thick accretion disk, which would be inconsistent with the evidence for thermal excess emission, peaking at $`2\times 10^{15}`$ Hz, observed in the low $`\gamma `$-ray state of 3C 279 by Pian et al. (1999) and with the peak frequency of $`3\times 10^{15}`$ Hz of the big blue bump seen in 3C 273 (Lichti et al., 1995) which is generally attributed to the emission from the optically thick accretion disk.

## 4 Discussion and Summary

The proposed evolutionary scenario provides a simple physical connection between different blazar subclasses in terms of the depletion of the accretion flow onto supermassive black holes. After an early phase of rapid black-hole growth, the time scale for this depletion is given by the matter free-fall time $`t_{\mathrm{ff}}7\times 10^8m_8T_4^{3/2}\mathrm{yr}`$ at the Bondi-Hoyle accretion radius, where $`m_8`$ is the mass of the central black hole in units of $`10^8M_{}`$ and $`T=10^4T_4`$ K is the average temperature of the circumnuclear material. This is generally shorter than the Hubble time, and gives a plausible time scale to explain a finite duration of blazar-like activity in the nuclei of elliptical galaxies. The scenario proposed here is also in accord with the X-ray faintness of giant elliptical galaxies known to contain supermassive black holes, which may be remnants of past blazar activity. The central sources in these galaxies have been successfully modeled with advection-dominated accretion flow solutions at very low accretion rates, $`l1`$ (Fabian & Rees, 1995; Di Matteo & Fabian, 1997; Di Matteo et al., 2000), suggesting a central region depleted after an earlier period of more efficient accretion.

If this picture is correct, then there should also be objects in an earlier phase of rapid black-hole growth that are accreting at rates $`l1`$. Following the scenario outlined by Sanders et al. (1988), this earlier stage of blazar evolution would comprise merging galaxies, infrared luminous galaxies, and radio-quiet QSOs. Accretion flows near the Eddington limit might produce optically-thick accretion disks extending all the way to the innermost stable orbit. If the phenomenology of Galactic black-hole candidates can be scaled to supermassive black holes accreting at comparable values of $`l`$, then radio jet formation should be quenched in such a mode of accretion as is observed, e.g., in GRS 1915+105 (Mirabel et al., 1998; Feroci et al., 1999; Belloni, Migliari, & Fender, 2000), GX 339-4 (Fender et al., 1999; Corbel et al., 2000), and XTE J1550-564 (Corbel et al., 2001). This idea is supported by recent theoretical work that links the existence and energetics of a central advection-dominated inflow-outflow system to the formation of relativistic jets (Blandford & Begelman, 1999; Becker, Subramanian, & Kazanas, 2001). The ASCA detections of relativistically broadened Fe K$`\alpha `$ emission in MCG-6-30-15 (Tanaka et al., 1995) and IRAS 18325-5926 (Iwasawa et al., 1996) and the rapid X-ray variability seen in IRAS 13224-3809 and PHL 1092 (Boller, 2000), suggest that at least in a number of cases an optically thick accretion disk does extend towards the innermost stable orbit. In most cases the evidence is yet inconclusive, but a similar mode of accretion can generally not be ruled out.

The formation of radio jets in blazars and radio galaxies, if related to the formation of an advection-dominated accretion mode in the inner portions of the accretion flow, is triggered by a decreasing Eddington ratio. The decline of $`l`$ might be due to a combination of a decreasing accretion rate and an increasing black-hole mass. An implication of this scenario is that the masses of the central black holes in galaxies which host BL Lacs should, on average, be greater than the masses of black holes in galaxies which host FSRQs or optical QSOs. Moreover, subclasses at earlier stages in the blazar sequence should exhibit increasingly stronger cosmological evolution. Bade et al. (1998) find evidence for negative cosmological evolution of X-ray selected BL Lac objects, and Stickel et al. (1991) find evidence for positive cosmological evolution in a radio-selected BL Lac sample, in accord with this picture. Our work provides strong support for the suggestion by D’Elia & Cavaliere (2000) and Cavaliere & D’Elia (2001) along parallel lines, but on the basis of general energy considerations and a generic cosmological evolution scenario for quasars evolving into BL Lac objects. We have combined the basic idea of such a cosmological evolution with a self-consistent spectral blazar model and demonstrated that the spectral properties of the different blazar subclasses can be consistently explained by such an evolutionary sequence. The blazar model adopted here is similar to the one previously used by Ghisellini et al. (1998) — but with the important addition of particle escape and adiabatic losses — who had argued that the different blazar subclasses could be unified by an increasing level of external radiation to the soft radiation field in the emitting region.

To summarize, we have proposed an evolutionary scenario linking FSRQs, LBLs, and HBLs through gradual depletion of the circumnuclear environment of a supermassive black hole. As the accretion power declines with time, less gas and dust is available to reprocess accretion-disk radiation and produce an external Compton-scattered component in the blazar jet. Analytic formulae are given to estimate the relevant model parameters for FSRQs. We have calculated a spectral sequence for different blazar subclasses that are related through a one-parameter model defined by the optical depth of the circumnuclear matter. This scenario links radio-loud AGNs with progenitor merger galaxies and QSOs that constitute the high Eddington-ratio limit of the evolutionary sequence. In analogy with the observed anti-correlation between radio and soft X-ray activity in some Galactic black-hole candidates, these progenitor sources would accrete with large Eddington ratios, resulting in an optically-thick accretion disk that quenches relativistic jet production.

We thank the anonymous referee for a very careful review and helpful and constructive comments. The work of MB is supported by NASA through Chandra Postdoctoral Fellowship grant PF 9-10007 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS 8-39073. The work of CD is supported by the Office of Naval Research.


text:
# 1 Introduction

## 1 Introduction

The original treatments of the influence of magnetic and non–magnetic disorder on superconductors applied to classic BCS superconductors , has been reexamined for superconductors with $`d`$–wave symmetry of order parameter . On the other hand, in high temperature superconductors the distance between the chemical potential and Van Hove singularity was found to be relatively small. This has lead to formulation of Van Hove scenario for high temperature superconductors which says that optimal critical temperature is reached when chemical potential passes through the Van Hove singularity in the density of states . However, doping with charge carriers does not only change the density in the system but also smear the density of states eliminating its singularities and, specially for anisotropic superconductors, introduce the electron pair–breaking phenomenon .

## 2 Superconductivity in a Disordered System

We start with the single band Hubbard model with an attractive extended interaction which is described by the following Hamiltonian :

$$H=\underset{ij\sigma }{}t_{ij}c_{i\sigma }^{}c_{j\sigma }+\frac{1}{2}\underset{ij}{}U_{ij}n_in_j\underset{i}{}(\mu \epsilon _i)n_i.$$
(1)
In the above $`n_i=n_i+n_i`$ is the charge on site labeled $`i`$, $`\mu `$ is the chemical potential. Disorder is introduced into the problem by allowing the local site energy $`\epsilon _i`$ to vary randomly from site to site, $`c_{i\sigma }^{}`$ and $`c_{i\sigma }`$ are the Fermion creation and annihilation operators for an electron on site $`i`$ with spin $`\sigma `$, $`t_{ij}`$ is the amplitude for hopping from site $`j`$ to site $`i`$ and finally $`U_{ij}`$ is the attractive interaction ($`U_{ij}<0`$), between electrons on neighbour sites ($`ij`$).

Here we will assume, for simplicity, that the random site energy $`\epsilon _i`$ has a uniform distribution $`\epsilon _i[\delta /2,\delta /2]`$. Following the usual way we shall apply the Coherent Potential Approximation (CPA) to mean that the coherent potential $`𝚺(E)=𝚺(i,i;E)`$ , in a site approximation, is defined by the zero value of an averaged t-matrix $`𝑻(i,i;E)`$.

The linearized gap equation in Hartree-Fock-Gorkov approximation :

$$1=\frac{|U|}{\pi }_{\mathrm{}}^{\mathrm{}}dE\mathrm{tanh}\frac{E}{2k_BT_c}\mathrm{Im}\frac{\overline{G}_{11}^d(E)}{2E\mathrm{Tr}𝚺(E)},$$
(2)
where $`T_c`$ is the critical temperature, $`k_B`$ denotes the Boltzmann constant, $`\overline{G}_{11}^d(E)`$ is an averaged electron Green function which defines the weighted density of states (DOS) of $`d`$-wave electron states $`\overline{N}_d(E)`$:

$$\overline{N}_d(E)=\frac{1}{\pi }\mathrm{Im}\overline{G}_{11}^d(E)=\frac{1}{\pi N}\underset{𝒌}{}\mathrm{Im}\frac{\eta _𝒌^2}{4}\frac{1}{E\mathrm{\Sigma }_{11}(E)\epsilon _𝒌+\mu }.$$
(3)
$`\mathrm{\Sigma }_{11}(E)`$ describes the electron self energy in the normal disordered system and $`\eta _𝒌=2(\mathrm{cos}k_xcosk_y)`$. Examples of projected densities of states $`N_d(E)`$ for disordered 2D system stem are presented in Fig. 1a while the corresponding self energies are plotted in Fig. 1b.

Equation 1 can be rewritten in terms of Matsubara frequencies $`\omega _n`$ as

$$1=\frac{|U|}{2k_BT_c}\underset{n}{}\frac{\overline{G}_{11}^d(ı\omega _n)}{2ı\omega _n\mathrm{Tr}𝚺(ı\omega _n)}$$
(4)
In Fig. 2a we show the density $`N^d(ı\omega )=\mathrm{Im}\overline{G}^d(ı\omega _n)/\pi `$ versus imaginary frequency $`\omega `$, where $`E=ϵ_f+ı\omega `$ for Fermi energy $`ϵ_f`$ chosen at Van Hove singularity ($`ϵ_f=0`$, Fig. 1a). In the standard treatments  it was assumed to be constant. However in our case, due to the Van Hove singularity, it depends strongly on $`\omega `$. The corresponding projected density $`N^d(ı\omega )`$ can be roughly approximate by a simple formula $`N_d(ı\omega )=c/(\omega +b)`$ where $`c`$ and $`b`$ are constants (Fig. 2a). Approximating also $`\mathrm{\Sigma }(i\omega )=ı|\mathrm{\Sigma }_0|\mathrm{sgn}(\omega )`$ we get from Eq. 2 the analytic pair-breaking formula:

$$\psi \left(\frac{1}{2}\right)\psi \left(\frac{1}{2}+\frac{b}{2\pi T_{c0}}\right)=\psi \left(\frac{1}{2}+\rho _c\right)\psi \left(\frac{1}{2}+\rho _c+\frac{b}{2\pi T_c}\right),$$
(5)
where $`\rho _c=|\mathrm{\Sigma }_0|/(2\pi T_c)`$ is a pair–breaking parameter and $`T_{c0}`$ is the critical temperature for a clean superconductor. Note that for large $`b`$ ($`b\mathrm{}`$) Eq. 4 transforms into the standard Abrikosov-Gorkov equation . Fig. 2b presents the comparison between these analytic formulae as well as the numerical results obtained from the gap equation (Eq. 2). Both analytic and numerical results (Fig. 2b) show that Van Hove singularity make the superconducting system more robust in presence the disorder.

## 3 Conclusions

Analyzing the effect of disorder on disordered $`d`$–wave superconductor we have found the additional influence of the Van Hove singularity. In the presence of a weak disorder $`\delta <t`$ we observe only small change of the projected density of states $`\overline{N}_d`$ (Fig. 1a) present in the gap equation (Eqs. 2,4) In the same time we observe the rapid degradation of $`T_c`$, which is connected with the pair-breaking effect (Eq. 2). Interestingly, Van Hove singularity modifies the standard Abrikosov-Gorkov formula, originally obtained for a constant density of states, increasing the critical value of $`|\mathrm{\Sigma }_0|`$ which destroys the superconducting phase (Fig. 2b). That result is in qualitative agreement with the experimental results on Zn substitutions .

## Acknowledgements

The author would like to thank Prof. K.I. Wysokinski, Prof. B. L. Györffy and Dr. J.F. Annett for helpful discussions.


text:
# Novel Theory for Topological Structure of Vortices in a Bose-Einstein condensate

## I $`\varphi `$Mapping Topological Current Theory for $`\times \stackrel{}{V}`$

The basic field of condensate wave-function $`\psi (\stackrel{}{x})`$ is a section of complex line bundle, i.e., a section of $`2`$-dimensional real vector bundle on $`R^3`$:

$$\psi (\stackrel{}{x})=\varphi ^1(\stackrel{}{x})+i\varphi ^2(\stackrel{}{x}).$$
(8)
Following $`\varphi `$-mapping theory a $`2`$-dimensional unit vector is defined as $`n^a=\frac{\varphi ^a}{\varphi }(a=1,2),`$ where $`\varphi ^2=\varphi ^a\varphi ^a=\psi ^{}\psi .`$Substituting these formulas into Eq.(2), it is easy to find out $`V^i=\frac{\mathrm{}}{m}ϵ_{ab}n^a_in^b,`$ and the curl of $`\stackrel{}{V}`$ can be expressed in terms of $`n^a`$:

$$(\times \stackrel{}{V})^i=\frac{\mathrm{}}{m}ϵ^{ijk}ϵ_{ab}_jn^a_kn^b.$$
(9)
Using $`_in^a=\frac{_i\varphi ^a}{\varphi }+\varphi ^a_i\frac{1}{\varphi }`$ and the Green function relation in $`\varphi `$-space, $`\frac{}{\varphi ^a}\frac{}{\varphi ^a}\mathrm{ln}\varphi =2\pi \delta ^2(\stackrel{}{\varphi }),`$ one can directly prove a novel expression for $`\times \stackrel{}{V}`$:

$$(\times \stackrel{}{V})^i=\frac{h}{m}\delta ^2(\stackrel{}{\varphi })D^i(\frac{\varphi }{x})=\frac{h}{m}j^i,$$
(10)
where $`D^i(\frac{\varphi }{x})=\frac{1}{2}ϵ^{ijk}ϵ_{ab}_j\varphi ^a_k\varphi ^b`$ is the Jacobian vector, and

$$j^i=\frac{1}{2\pi }ϵ^{ijk}ϵ_{ab}_jn^a_kn^b=\delta ^2(\stackrel{}{\varphi })D^i(\frac{\varphi }{x})$$
(11)
is just a simple $`2`$-dimensional case of the $`N`$-dimensional $`\varphi `$-mapping topological current . This formula including $`\delta ^2(\stackrel{}{\varphi })`$ to describe the singularities of $`\stackrel{}{\varphi }`$ is just the precise topological expression for $`\times \stackrel{}{V}`$ that Landau and Feynman expected to find out long time ago. Therefore an important conclusion is reached: $`\times \stackrel{}{V}=0,iff\stackrel{}{\varphi }0;\times \stackrel{}{V}0,iff\stackrel{}{\varphi }=0`$.

The implicit function theory shows that , under the regular condition $`\stackrel{}{D}\left(\frac{\varphi }{x}\right)0`$, the general solutions of

$$\varphi ^1(x,y,z)=0,\varphi ^2(x,y,z)=0$$
(12)
can be expressed as

$$x=x_j(s),y=y_j(s),z=z_j(s),(j=1,2,\mathrm{},N)$$
(13)
which represent $`N`$ isolated singular strings $`L_j`$ with parameter $`s`$. These strings are just known as the vortex lines.

In $`\delta `$-function theory , one can prove

$$\delta ^2(\stackrel{}{\varphi })=\underset{j=1}{\overset{N}{}}\beta _j_{L_j}\frac{\delta ^3(\stackrel{}{x}\stackrel{}{x}_j(s))}{\left|D(\frac{\varphi }{u})\right|_{\mathrm{\Sigma }_j}}𝑑s,$$
(14)
where $`D(\frac{\varphi }{u})_{\mathrm{\Sigma }j}=(\frac{1}{2}ϵ^{jk}ϵ_{mn}\frac{\varphi ^m}{u^j}\frac{\varphi ^n}{u^k}),`$ and $`\mathrm{\Sigma }_j`$ is the $`j`$th planer element transversal to $`L_j`$ with local coordinates $`(u^1,u^2)`$. The positive integer $`\beta _j`$ is the Hopf index of $`\varphi `$-mapping. Meanwhile it can be proved that the direction vector of $`L_j`$ is

$$\left(\frac{d\stackrel{}{x}}{ds}\right)_{x_j}=[\stackrel{}{D}(\frac{\varphi }{x})/D(\frac{\varphi }{u})_{\mathrm{\Sigma }_j}]_{x_j}.$$
(15)
Then from Eqs.(14) and (15) we find the important inner topological structure of $`\times \stackrel{}{V}`$:

$$\times \stackrel{}{V}=\frac{h}{m}\underset{j=1}{\overset{N}{}}\beta _j\eta _j_{L_j}\frac{d\stackrel{}{x}}{ds}\delta ^3(\stackrel{}{x}\stackrel{}{x}_j(s))𝑑s,$$
(16)
where the positive integer $`\beta _j`$ is the Hopf index of $`\varphi `$-mapping, and $`\eta _j`$ is the Brouwer degree, $`\eta _j=\pm 1.`$ And the winding number of $`\stackrel{}{\varphi }`$ around $`L_j`$ is $`W_j=\beta _j\eta _j.`$ Therefore the vorticity of vortex line $`L_j`$ is $`\mathrm{\Gamma }_j=_{\mathrm{\Sigma }_j}\times \stackrel{}{V}𝑑\stackrel{}{s}=\frac{h}{m}W_j,`$ where $`\mathrm{\Sigma }_j`$ is the $`j`$th planer element transversal to $`L_j`$; and the total vorticity on a surface $`\mathrm{\Sigma }`$ should be

$$\mathrm{\Gamma }=_\mathrm{\Sigma }\times \stackrel{}{V}𝑑\stackrel{}{s}=\frac{h}{m}\underset{j=1}{\overset{N}{}}W_j.$$
(17)
We stress that there are no hypothesis in the deduction above. Eqs.(10) and (16) are called the differential forms of the quantization condition, which cannot be derived from the single-valued principle of wave function, and are more essential than the integral form (Eq.(17)).

For the GL theory, Eqs.(4) and (5) lead to

$$\stackrel{}{A}+\lambda ^2\times \stackrel{}{B}=\frac{mc}{e}\stackrel{}{V},$$
(18)
where $`\lambda `$ is the penetration depth, $`\lambda ^2=\frac{1}{\rho }\frac{mc^2}{4\pi e^2}`$. In London approximation, $`\rho `$ and therefore $`\lambda `$, are treated as constants; hence when noticing $`\stackrel{}{B}=\times \stackrel{}{A}`$, $`\stackrel{}{B}=0`$ and Eq.(10), we find a topological equation for $`\stackrel{}{B}`$:

$$\stackrel{}{B}\lambda ^2^2\stackrel{}{B}=\frac{mc}{e}\times \stackrel{}{V}=\frac{hc}{e}\delta ^2(\varphi )\stackrel{}{D}(\frac{\varphi }{x}).$$
(19)
This formula directly leads to

$$\stackrel{}{B}\lambda ^2^2\stackrel{}{B}=\mathrm{\Phi }_0\underset{j=1}{\overset{N}{}}W_j_{L_j}\frac{d\stackrel{}{x}}{ds}\delta ^3(\stackrel{}{x}\stackrel{}{x}_j(s))𝑑s,$$
(20)
where $`\mathrm{\Phi }_0`$=$`\frac{hc}{e}`$ is the unit flux quantum. We see that in simple case $`W_j=1`$, the above equation is just the so-called modified London equation . This expression says that, when the condensate wave function $`\psi `$ has no zero values, $`\stackrel{}{\varphi }0`$, i.e., $`\delta ^2(\stackrel{}{\varphi })=0`$, and $`\stackrel{}{B}\lambda ^2^2\stackrel{}{B}=0,`$ which just corresponds to the Meissner state; while in the case of mixed state, $`\stackrel{}{\varphi }`$ possesses $`N`$ isolated zeros, $`\delta ^2(\stackrel{}{\varphi })0`$, thus a type-II superconductor is penetrated by an array of $`N`$ vortices, with each one carrying a quantum flux proportional to the winding number $`W_j`$.

## II Decomposition of $`U(1)`$ Gauge Potential and The Vortex with Quantized Flux

The decomposition theory of gauge potential in $`SO(N)`$ and $`SU(N)`$ gauge theories is now playing a more and more important role in theoretical studies, because it virtually inputs topological  and other important informations  to the gauge potential. In the theory of superconductivity, $`\psi `$ is a condensate wave function describing the charged continuum, so the covariant derivative in $`U(1)`$ gauge theory is introduced to describe the interaction between $`\psi `$ and the electromagnetic field:

$$D_i\psi =_i\psi i\frac{e}{\mathrm{}c}A_i\psi ,(i=1,2,3)$$
(21)
where $`A_i`$ is the magnetic gauge potential vector. The complex conjugate of $`D_i\psi `$ is $`D_i^{}\psi ^{}=_i\psi ^{}+i\frac{e}{\mathrm{}c}A_i\psi ^{}`$. And the magnetic field tensor is given by

$$f_{ij}=_iA_j_jA_i.$$
(22)
Multiplying $`D_i\psi `$ with $`\psi ^{}`$ and $`D_i^{}\psi ^{}`$ with $`\psi `$ respectively, we can deduce the decomposition formula for $`U(1)`$ gauge potential:

$$A_i(\psi )=\frac{\mathrm{}c}{2ie}\frac{1}{\psi ^{}\psi }[(\psi ^{}_i\psi _i\psi ^{}\psi )(\psi ^{}D_i\psi D_i^{}\psi ^{}\psi )].$$
(23)
The above expression $`A_i=A_i(\psi )`$ means that the magnetic gauge potential possesses an inner structure in terms of charged condensate wave function $`\psi `$ and $`\psi ^{}`$. The inner structure of $`A_i(\psi )`$ with Eq.(22) gives a theory that in superconductivity how the stationary motion of condensate wave function creates an intrinsic magnetic field. This is the important physical meaning of decomposition of $`U(1)`$ gauge potential in quantum mechanics.

Furthermore it has been proved that the covariant derivative part $`[(\psi ^{}D_i\psi D_i^{}\psi ^{}\psi )]`$ corresponds to the gradient of a phase factor: $`(_i\lambda )`$ . Thus this covariant derivative part contributes nothing to the field tensor $`f_{ij}`$, so it can be ignored, and

$$A_i(\psi )=\frac{\mathrm{}c}{2ie}\frac{1}{\psi ^{}\psi }(\psi ^{}_i\psi _i\psi ^{}\psi ).$$
(24)
It should be emphasized that the above $`U(1)`$ gauge potential decomposition theory together with $`\varphi `$-mapping theory has been successfully used to study many other topological problems in physics .

In $`\varphi `$-mapping theory $`A_i(\psi )`$ can be rewritten in terms of $`n^a`$ as $`A_i(\psi )=\frac{\mathrm{}c}{e}ϵ_{ab}n^a_in^b,`$ and $`f_{ij}`$ becomes $`f_{ij}=2\frac{\mathrm{}c}{e}ϵ_{ab}_in^a_jn^b.`$ Therefore the intrinsic magnetic field vector from $`A_i(\psi )`$ is expressed as

$$B_i(\psi )=\frac{1}{2}ϵ_{ijk}f_{jk}=\mathrm{\Phi }_0\frac{1}{2\pi }ϵ_{ijk}ϵ_{ab}_jn^a_kn^b.$$
(25)
Using Eq.(10) we have $`B_i(\psi )=\mathrm{\Phi }_0\delta ^2(\stackrel{}{\varphi })D^i(\frac{\varphi }{x}),`$ which gives the topological structure of intrinsic magnetic field $`B_i(\psi )`$; $`B_i(\psi )`$ does not matter with the external magnetic field. As before, the zero points of $`\stackrel{}{\varphi }(x)`$, i.e., the singular vortex lines in superconductivity contribute to intrinsic magnetic field as

$$B_i(\psi )=\mathrm{\Phi }_0\underset{j=1}{\overset{N}{}}W_j_{L_j}\frac{dx^i}{ds}\delta ^3(\stackrel{}{x}\stackrel{}{x}_j(s))𝑑s.$$
(26)
This leads to an important phenomenon that, the magnetic flux coming from the stationary motion of $`\psi `$ itself is quantized

$$\mathrm{\Phi }=_\mathrm{\Sigma }\stackrel{}{B}(\psi )𝑑\stackrel{}{s}=\mathrm{\Phi }_0\underset{j=1}{\overset{N}{}}W_j,$$
(27)
and each singular vortex line $`L_j`$ carries a magnetic flux $`\mathrm{\Phi }_j=W_j\mathrm{\Phi }_0.`$

The above decomposition theory of $`U(1)`$ gauge potential naturally arrives at the conclusion that, in superconductivity continuum, the $`N`$ isolated singular vortices are just $`N`$ isolated topological elementary excitations carrying with magnetic fluxoid, while their quantum numbers are characterized by topological numbers $`W_j=\beta _j\eta _j`$. We see that the $`\varphi `$-mapping topological theory in this paper is independent of concrete physical models, that gives a profound understanding to the nature of the creation of the vortex lines and the flux quantization in BEC.

Comparing Eq.(24) with Eq.(2) it directly follows a simple relation between $`V_i`$ and $`A_i`$

$$V_i=\frac{e}{mc}A_i,$$
(28)
which is just the London’s assumption . We stress that, the essence and the significance of this relation are not truly realized until now the inner structure of gauge potential is revealed and therefore the stationary motion of condensate wave function is naturally related to the intrinsic magnetic field.

## III Spatial Bifurcation of Vortex Lines

The Solution (13) of Eq.(12) is based on the condition $`\stackrel{}{D}\left(\varphi /x\right)0`$. When it fails, i.e.,

$$\stackrel{}{D}\left(\frac{\varphi }{x}\right)=0$$
(29)
at some points (marked as $`\stackrel{}{r}_j^{}`$) along $`L_j`$, the functional relationship between coordinate $`x`$ and $`z`$, or $`y`$ and $`z`$ is not unique in the neighborhood of $`\stackrel{}{r}_j^{}`$, because the direction of the zero line expressed by

$$\frac{dx}{dz}=D^1\left(\frac{\varphi }{x}\right)/D^3\left(\frac{\varphi }{x}\right)|_{\stackrel{}{r}_j^{}},\frac{dy}{dz}=D^2\left(\frac{\varphi }{x}\right)/D^3\left(\frac{\varphi }{x}\right)|_{\stackrel{}{r}_j^{}}$$
(30)
is indefinite at $`\stackrel{}{r}_j^{}`$. Hence this very point $`\stackrel{}{r}_j^{}`$ is called a bifurcation point of the two-component vector in $`3`$-dimensional space.

According to the $`\varphi `$-mapping theory, the Taylor expansion of the solution of Eq.(12) in the neighborhood of $`\stackrel{}{r}_j^{}`$ can be generally expressed as : $`A(xx_j^{})^2+2B(xx_j^{})(zz_j^{})+C(zz_j^{})^2+\mathrm{}=0,`$ where $`A,B`$ and $`C`$ are constants. This leads to

$$A(\frac{dx}{dz})^2+2B\frac{dx}{dz}+C=0orC(\frac{dz}{dx})^2+2B\frac{dz}{dx}+A=0.$$
(31)
The solutions of Eq.(31) give different branches of the zero lines, i.e., the vortex lines at bifurcation points. In following four main cases in the branch process are simply discussed (the detailed deduction and figures may be found in Ref ):

Case 1 ($`A0`$): For $`\mathrm{\Delta }=4(B^2AC)>0`$, from Eq.(31) we get two different spatial directions at the bifurcation point

$$\frac{dx}{dz}_{1,2}=\frac{B\pm \sqrt{B^2AC}}{A}.$$
(32)
This is the intersection of two vortex lines of different directions.

Case 2 ($`A0`$): For $`\mathrm{\Delta }=4(B^2AC)=0`$, we get only one direction at the point

$$\frac{dx}{dz}_{1,2}=\frac{B}{A},$$
(33)
which includes three sub-cases: (a) Two vortex lines tangentially intersect; (b) Two vortex lines merge into one line; (c) One vortex line splits into two lines.

Case 3 ($`A=0,C0`$): For $`\mathrm{\Delta }=4(B^2AC)>0`$, from Eq.(31) we have

$$\frac{dz}{dx}_{1,2}=\frac{B\pm \sqrt{B^2AC}}{C}=0,\frac{2B}{C}.$$
(34)
There are two sub-cases: (a) Three vortex lines merge into one line; (b) One vortex line splits into three lines.

Case 4 ($`A=C=0`$): Eq.(31) gives respectively

$$\frac{dx}{dz}=0,\frac{dz}{dx}=0.$$
(35)
This case shows that two curves normally intersect at the bifurcation point, which is similar to case 3.

It should be noted that, noticing the continuity of topological current $`\stackrel{}{j}`$ from Eq.(11) ($`_ij^i=0`$), at the bifurcation point the sum of the topological charge of final vortex line(s) is required to be equal to that of the initial line(s) for a fixed index $`j`$: $`_f\beta _{j_f}\eta _{j_f}=_i\beta _{j_i}\eta _{j_i}.`$

This work was supported by the National Natural Science Foundation and the Doctor Education Fund of Educational Department of the People’s Republic of China.


text:
# Studies of spin-orbit scattering in noble-metal nanoparticles using energy level tunneling spectroscopy

## Abstract

The effects of spin-orbit scattering on discrete electronic energy levels are studied in copper, silver, and gold nanoparticles. Level-to-level fluctuations of the effective $`g`$-factor for Zeeman splitting are characterized, and the statistics are found to be well-described by random matrix theory predictions. The strength of spin-orbit scattering increases with atomic number and also varies between nanoparticles made of the same metal. The spin-orbit scattering rates in the nanoparticles are in order of magnitude agreement with bulk measurements on disordered samples.

Spin-orbit (SO) interactions are central to determining the energy-level structure and dynamics of atomic nuclei, large-atomic-number atoms, and a wide variety of solid-state systems. They have important effects on electronic band structures , they influence the symmetry class of random-matrix theories applied to quantum systems , and they are put to use in proposals for controllable manipulation of electron spins . A particularly fundamental method for probing SO interactions in solid-state systems is to examine individual quantum energy levels in nanoparticles, both because SO matrix elements can be measured directly and because SO interactions produce strong modifications of the effective $`g`$-factors for spin-Zeeman splitting. Previously, measurements of three individual energy levels in particular gold nanoparticles showed that SO effects reduce the $`g`$-factor to $`g0.280.45`$ , in sharp contrast with Al particles, which ordinarily give $`g2`$ . Adding 4% gold to an aluminum nanoparticle also increased the SO scattering, thereby reducing the $`g`$-factors in the particle to $`g0.50.8`$ . In this paper we report systematic studies of SO effects in noble-metal nanoparticles covering a range of atomic number ($`Z`$) – copper, silver, and gold. We characterize the mesoscopic fluctuations, present within single particles, of the $`g`$-factors and SO matrix elements $`H_{\mathrm{so}}`$, and find that the statistics for the g-factors are well-described by recent predictions of random matrix theory (RMT) . From the values of the g-factors and the mean energy-level spacing $`\delta `$ , we extract SO scattering rates for each material, and we find order of magnitude agreement with expectations based on weak localization and electron spin resonance measurements.

SO interactions cause the electronic energy levels in a metal grain to be not purely spin-up or spin-down, and as a consequence they reduce the $`g`$-factor for spin-Zeeman splitting below the free-electron value of 2. In low-$`Z`$ nanoparticles for which SO scattering is sufficiently weak (with $`H_{\mathrm{so}}\delta `$) perturbation theory can be used to calculate the $`g`$-factor for a given electronic state $`n`$ :

$$g_n=2\left(12\underset{mn}{}\frac{|\mathrm{\Psi }_m|H_{\mathrm{so}}|\mathrm{\Psi }_n|^2}{(E_nE_m)^2}\right).$$
(1)
Here the $`\mathrm{\Psi }_m`$ denote the unperturbed (pure-spin) single-electron eigenstates with energy $`E_m`$. It is expected that the SO matrix elements will vary between different pairs of unperturbed energy levels because they depend on the precise nature of the wavefunctions, which are chaotic and fluctuating . The local energy level spacing will also vary. As a result, $`g_n`$ should exhibit mesoscopic fluctuations for different levels, $`n`$. For heavier elements, perturbation theory ceases to be valid. However, for large SO-interaction strengths, Brouwer et al. and Matveev et al. have recently employed RMT to predict the statistical distribution $`P(g)`$ of $`g`$-factors in a nanoparticle . In , the strength of SO interactions is tuned continuously, allowing numerical studies of the full crossover from the Guassian orthogonal ensemble (GOE) to the Gaussian symplectic ensemble (GSE) for the energy-level statistics. In this crossover, it is found that the mean value of $`P(g)`$ shifts to smaller g-factors and the variance of the distribution exhibits a maximum for $`g`$ 1.2.

Our tunneling samples are fabricated in a vertical geometry through a 3–10 nm bowl-shaped hole etched in a Si3N4 membrane  (see inset, Fig. 1). One electrode is made by evaporating 1750 Å of Al on top of the membrane to fill the bowl. Following oxidation in 50 mTorr of O2 for 3 minutes to form a tunnel barrier, 5–20 Å of Cu, Ag, or Au is evaporated on the lower side of the device to form particles with radii of 3–5 nm. A second tunnel barrier is then formed by the evaporation of 11 Å of Al2O3. This is followed by 1750 Å of Al to make the lower tunneling electrode. The samples are screened at 4.2 K to select those whose current-voltage ($`I`$-$`V`$) curves exhibit Coulomb-staircase steps, indicating tunneling through a single nanoparticle. Then more detailed $`I`$-$`V`$ curves are measured in a dilution refrigerator, with all electrical leads passing through copper-powder filters to minimize the effects of RF radiation.

Figure 1 shows $`I`$-$`V`$ and $`dI/dV`$-$`V`$ curves for a Cu

device (Cu$`\mathrm{\#}`$1). As $`V`$ is increased, $`I`$ increases in discrete steps due to tunneling through well-resolved energy levels within the nanoparticle . The resonances in the differential conductance directly correspond to the energy level spectrum of the nanoparticle. In an applied magnetic field ($`H`$) most resonances in $`dI/dV`$ divide into 2 peaks, allowing measurements of the energy-difference for spin-Zeeman splitting between Kramers-pair eigenstates, $`\mathrm{\Delta }E=g\mu _BH`$. In order to determine the effective $`g`$-factors, we must convert the magnitude of the splitting from the measured value in voltage to energy, correcting for capacitive division of the voltage across the two tunnel junctions according to $`\mathrm{\Delta }E=\mathrm{\Delta }V/(1+C_1/C_2)`$ where $`C_1/C_2`$ is the capacitance ratio for the two tunnel junctions. This ratio can be determined by measuring the shift in $`V`$ for the resonances as the aluminum leads are changed from superconducting to normal in small magnetic fields, or by comparing the voltage for tunneling through the same quantum state for positive and negative $`V`$ . The low-field $`g`$-factors for sample Cu$`\mathrm{\#}`$1 fall in the range from 1.30$`\pm `$0.04 to 1.82$`\pm `$0.05.

In the weak SO-scattering limit where perturbation theory is valid, the effects of SO interactions on the magnetic field evolution of neighboring energy levels (shown in detail in Fig. 1(c) and in broader view in Fig. 2(a)) can be modeled by diagonalizing a simple $`2\times 2`$ Hamiltonian matrix with off-diagonal matrix elements $`H_{\mathrm{so}}`$ coupling the eigenstates. The minimum energy difference at the avoided level crossing directly gives 2$`H_{\mathrm{so}}`$ between those two states. The solid lines in Fig. 1(c) correspond to solutions of $`2\times 2`$ Hamiltonians with $`H_{\mathrm{so}}`$=76 $`\mu `$eV for the avoided crossing near 4.9 meV, and $`H_{\mathrm{so}}`$=134 $`\mu `$eV

for the avoided crossing near 5.4 meV. These matrix elements are in the size range expected based on the measured $`g`$-factors and the mean level spacing $`\delta =`$ 0.70 meV. If the level spacing were uniform and equal to this value, and all the matrix elements $`H_{\mathrm{so}}`$ were $``$ 130 $`\mu `$eV, then Eq. (1) would predict $`g=26.6(H_{\mathrm{so}}/\delta )^2`$ 1.77. If energy-level fluctuations are included, the estimate for the typical $`g`$-factor is somewhat smaller . In accord with Eq. (1), we have seen that in a given sample the states which happen to have smaller-than-average energy spacings from neighboring levels generally have smaller-than-average $`g`$-factors.

For the heavier elements Ag and Au, the g-factors become smaller with increasing atomic number. Figure 2 shows color scale plots of $`dI/dV`$ for increasing magnetic field in Cu, Ag, and Au samples. We find $`g`$-factors in the range 1.3 to 1.9 for Cu, 0.25 to 1.1 for Ag, and 0.05 to 0.19 in Au for the samples shown. Table 1 lists the measured mean values of the g-factors and their standard deviations for all the samples we have measured with more than three resolvable resonances.

The most direct way to compare the statistics of the measured $`g`$-factors with theoretical predictions is by plotting the integrated probability distribution of the $`g`$-factors for each sample. In Fig. 3 the points correspond to the g-factors from several Au, Ag, and Cu samples, while the solid lines are the theoretical predictions from Brouwer et al. with the SO-scattering strength parameter $`\lambda `$ adjusted to minimize the least squares error between the theoretical and experimental distributions. Our results are in good qualitative agreement with the RMT predictions. In particular, for a single value of the adjustable parameter $`\lambda `$, the theoretical distributions account reasonably for both the mean value of $`g`$ and the

standard deviation for each sample.

We do observe significant differences in the average value of the $`g`$-factor for different nanoparticles composed of the same material, e.g., for the two Cu samples depicted in Fig. 3 and for the Ag samples #1 and #2. This is true even though the mean level spacings for the two Cu samples are nearly identical. Because of these differences, we cannot meaningfully consolidate the data from many samples to improve the quality of our statistics in making comparisons to the RMT predictions. The physical origin of SO interactions at low temperatures within noble metal nanoparticles is expected to be scattering from defects or the nanoparticle surface . Sample-to-sample variations in the strength of SO interactions are therefore reasonable in nanoparticles, if the defect concentration or the surface-to-volume ratio may vary, for instance due to size and shape differences between particles.

A few of the energy levels shown in Fig. 2 do not exhibit Zeeman splitting. We observe this only in samples where we can identify the tunneling transitions as corresponding to changes from an odd to an even number of conduction electrons within the particle. We make this identification by observing the lowest energy tunneling resonance. If the nanoparticle initially has an odd number of electrons, then the first energy level accessible for tunneling will have only one spin state available, and therefore the lowest energy transition does not undergo Zeeman splitting (Fig. 2(a,c,d)) . In other words, the tunneling transition occurs via a spin-singlet state. Similarly, we can explain the lack of observed splitting in some higher energy resonances as due to tunneling via excited spin-singlet states. (See the states near 7 meV in Fig. 2(c) and 0.3 meV in Fig. 2(d).) It is also possible that the Zeeman splitting may sometimes be sufficiently small for a level in Ag or Au that we cannot resolve it. Our resolution limit is approximately $`g`$ 0.03. From the theoretical probability distributions of Fig. 3, about 0.02% of the states in Ag and 1.4% of the states in Au can be expected to have g$`<`$0.03.

In some samples, background-charge noise and nonequilibrium effects can interfere with accurate determinations of the statistics of the $`g`$-factors. For instance, as a function of the voltage across the device or as a function of time, the background charge affecting the nanoparticle can sometimes take on an ensemble of different values, shifting the energy levels that one is trying to measure. When this happens, current steps due to background-charge changes may be misinterpreted as due to quantum states, or the same quantum state may be measured several times at different values of voltage. In one Cu sample not shown in this paper, we measured six apparently different resonances all with an identical $`g`$-factor of 1.50, which we ascribe to this effect. This behavior was not seen in any other samples measured in this study. If several closely-spaced levels are observed in a sample, ways to diagnose the cause as background noise include the presence of negative differential conductance even with normal-state electrodes, the presence of apparent level crossings rather than avoided crossings, and identical $`g`$-factors for several resonance peaks.

An important question for understanding the properties of nm-scale devices is whether the strength of SO scattering is consistent with expectations based on measurements made by other techniques on larger samples, or whether there is new physics at small scales. Because effective $`g`$-factors measured by different experimental techniques can probe physically distinct quantities , we will make the comparison based on estimates for SO scattering times ($`\tau _{\mathrm{so}}`$) that we can extract from our data. We will also make comparison only to samples with mean free paths for electron scattering comparable to the size of our nanoparticles, since this scattering is the probable source of SO-induced spin relaxation . In the limit of strong SO scattering in a nanoparticle, the relationship between the $`g`$-factor and $`\tau _{\mathrm{so}}`$ has been solved analytically :

$$g^2=\frac{3}{\pi \mathrm{}}\tau _{\mathrm{so}}\delta .$$
(2)
For the Au samples, using the values of $`\delta `$ listed in Table 1, we find $`1/\tau _{\mathrm{so}}=(28)\times 10^{12}s^1`$. The Ag and Cu samples fall outside the strong-scattering regime. However by fitting to the measured $`g`$-factors the numerical distributions calculated by Brouwer et al., we can still determine $`\tau _{\mathrm{so}}`$ for each sample from the single fitting parameter $`\lambda `$ defined in : $`1/\tau _{\mathrm{so}}=\lambda ^2\delta /(\pi \mathrm{})`$. By this method, the SO scattering rates determined for the two Ag samples are quite different, $`1/\tau _{\mathrm{so}}=3\times 10^{11}s^1`$ and $`2\times 10^{12}s^1`$, while for Cu $`1/\tau _{\mathrm{so}}=(28)\times 10^{11}s^1`$.

For Cu, we can also use perturbation theory to check the rate determined from the RMT approach. The value of $`\tau _{\mathrm{so}}`$ within perturbation theory is

$$1/\tau _{\mathrm{so}}=\frac{2\pi |H_{\mathrm{so}}|^2}{\mathrm{}\delta }.$$
(3)
For $`H_{\mathrm{so}}130\mu `$eV and $`\delta `$=0.7 meV, we find $`1/\tau _{\mathrm{so}}2\times 10^{11}s^1`$, in agreement with the rate for Cu#1 listed in Table 1.

These values for $`\tau _{\mathrm{so}}`$ can be compared to the SO scattering rates measured previously by weak localization, and from the linewidths of conduction electron spin resonance (CESR) in metal foils. From weak localization measurements on quench-condensed noble-metal films, Bergmann  found $`1/\tau _{\mathrm{so}}=5.4\times 10^{12}s^1`$ for Au, $`2.6`$-$`3.1\times 10^{11}s^1`$ for Ag (depending on the disorder) and $`3.2\times 10^{11}s^1`$ for Cu. The values for Au and Cu agree to within a factor of three with all the rates we extract, while the weak-localization rate for Ag coincides with sample Ag#2 but not Ag#1. The CESR comparisons are less direct, because these measurements are performed on crystalline foils with thicknesses on the order of microns or more. In these experiments it is generally found that the spin-flip rate scales inversely with sample thickness due to the dominant effect of surface scattering . The measured values include spin-flip rates of $`7\times 10^9s^1`$ for Au foils one hundred microns thick , $`1.5\times 10^8s^1`$ for a 5-$`\mu `$m-thick Ag foil , and $`2.6\times 10^8s^1`$ for a 1.2 $`\mu `$m Cu film . Assuming that the rates would scale as 1/thickness down to a thickness of 5 nm, we arrive at estimates of $`1\times 10^{14}s^1`$ for Au, $`1\times 10^{11}s^1`$ for Ag, and $`6\times 10^{10}s^1`$ for Cu. These SO scattering rates are also in rough order-of-magnitude agreement with the values obtained from the $`g`$-factor analysis.

In conclusion, we have measured the $`g`$-factor distributions in Cu, Ag, and Au nanoparticles using tunneling spectroscopy. The $`g`$-factors within a given sample exhibit mesoscopic fluctuations, with distributions in good accord with random-matrix-theory predictions. We have observed significant differences in both the $`g`$-factor distributions and $`\tau _{\mathrm{so}}`$ between different nanoparticles composed of the same metal, illustrating the importance of the specific sample structure in the origin of SO scattering. The SO scattering rates we extract for the nanometer-scale particles are in rough order-of-magnitude agreement with extrapolations from previous CESR and weak localization measurements.

We thank Piet Brouwer and Xavier Waintal for discussions and the use of their simulation curves. This work was supported by the Packard Foundation and the NSF (DMR-0071631). JRP acknowledges the support of a NSF Graduate Research Fellowship. Sample fabrication was performed at the Cornell node of the National Nanofabrication Users Network.