text:
# Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign

## 1 INTRODUCTION

The uninterrupted 7-day ASCA observation of the TeV blazar Mrk 421 in 1998 has revealed that day-scale X–ray flares seen in previous observations were probably unresolved superpositions of many smaller flares (Takahashi et al., 2000). The nearly–continuous observation allowed not only the possibility to track the individual flares entirely from the rise to decay, but it also enabled quantitative statistical tests of the time series by employing the power spectrum or the structure function (Kataoka et al., 2001; Tanihata et al., 2001).

The main characteristic of blazars is their high flux observed from radio to $`\gamma `$–rays coupled with strong variability and strong polarization. These properties are now successfully explained by the scenario where blazars are active galactic nuclei (AGN), possessing jets aligned close to the line of sight, and accordingly the Doppler-boosted non-thermal emission from the jet dominates other emission components (e.g. Blandford & Königl 1979; Urry & Padovani 1995). This is what makes blazars critical in understanding jets in AGN.

The broadband spectra of blazars consist of two peaks, one in the radio to optical–UV range (and in some cases, reaching to the X–ray band), and the other in the hard X–ray to $`\gamma `$–ray region. The high polarization of the radio to optical emission suggests that the lower energy peak is produced via the synchrotron process by relativistic electrons in the jet. The higher energy peak is believed to be due to Compton up-scattering of seed photons by the same population of relativistic electrons. Several possibilities exist for the source of the seed photons; these can be the synchrotron photons internal to the jet (Jones, O’Dell, & Stein, 1974; Ghisellini & Maraschi, 1989), but also external, such as from the broad emission line clouds (Sikora, Begelman, & Rees, 1994) or from the accretion disk (Dermer, Schlickeiser, & Mastichiadis, 1992; Dermer & Schlickeiser, 1993).

The blazars with peak synchrotron output in the X–ray range also emit strongly in the $`\gamma `$–ray energies, and the brightest of those have been detected in the TeV range with ground-based Cherenkov arrays. These are the so-called “TeV blazars.” In TeV blazars, variability of the synchrotron flux is measured to be the strongest and most rapid in the X–ray band, and thus it provides the best opportunity to study the electrons that are accelerated to the highest energies. In particular, the synchrotron peak is a very important observable in two aspects: first because the flux at the peak represents the total emitted power from the blazar, and secondly because the peak frequency reflects the maximum energy of radiating particles gained in the acceleration process.

Mrk 421 is among the closest known blazars, at redshift of 0.031, and was the first blazar (and also the first extragalactic source) discovered to be a TeV emitter. It was first detected as a weak source by EGRET (Lin et al., 1992), and 9 months later, Whipple detected a clear signal from this object between 0.5 and 1.5 TeV (Punch et al., 1992; Petry et al., 1996). Flux variability on various time scales has been observed, including a very short flare with a duration of $``$1 hour (Gaidos et al., 1996). Ever since, Mrk 421 has been repeatedly confirmed to be a TeV source by various ground-based telescopes. It has also been one of the most studied blazars, and have been a target of several multi-wavelength campaigns (Macomb et al., 1995; Maraschi et al., 1999; Takahashi et al., 2000).

The multi-wavelength campaign of Mrk 421 in 1998 was one of the first opportunities to observe a blazar in the TeV range using several telescopes located in different locations in the world, so as to have as continuous coverage as possible (Aharonian et al., 1999; Piron et al., 2001). Observations in other frequencies included X–ray observations by ASCA, RXTE (Takahashi et al., 2000), and BeppoSAX (Maraschi et al. 1999; Fossati et al. 2000a,b), EUV observations by EUVE, optical observation with BVRI filters organized by the WEBT collaboration (http://www.to.astro.it/blazars/webt/), and radio observations at the Metsähovi Radio Observatory.

In this paper, we present the results of the spectral analysis of X–ray and EUV data during the 1998 April campaign. In particular, we estimate the location of the synchrotron peak in the $`\nu F_\nu `$ spectrum. So far, quantitative studies of the variation of the synchrotron peak have been conducted only for the two brightest blazars, Mrk 421 and Mrk 501. The largest variation of the synchrotron peak energy was observed in Mrk 501 by BeppoSAX, where the peak energy shifted $``$2 orders of magnitude (Pian et al., 1998). Collecting data from different epochs (separated by as long as years), Tavecchio et al. (2001) have shown the relation of the form $`E_{\mathrm{peak}}F_{0.1100\mathrm{k}\mathrm{e}\mathrm{V}}^n`$ with $`n`$2. A similar analysis was done by Fossati et al. (2000b) for the Mrk 421 data obtained by BeppoSAX in 1997 and 1998, which showed a relation of $`E_{\mathrm{peak}}F_{0.110\mathrm{k}\mathrm{e}\mathrm{V}}^n`$ with n=0.55. In this paper we describe a continuous 7-day variation of the synchrotron peak, which allows us to investigate the dynamical change of the synchrotron spectrum during the flares.

We first describe the observations and data analysis in §2. The results of spectral analysis are described in §3; a summary of the results and a discussion of our findings with the emphasis on the temporal evolution of the spectrum are presented in §4.

## 2 OBSERVATION AND DATA ANALYSIS

The week-long multi-frequency campaign of Mrk 421 was carried out in April 1998 (Takahashi et al., 2000). The goal of this paper is to study the evolution of the high energy synchrotron spectrum as a function of time during the campaign, and also the correlation with the TeV flux. This requires the best possible spectral coverage in the X–ray range, for which we assembled all available data collected during the campaign, including: pointings by X–ray observatories ASCA, RXTE, and BeppoSAX, as well as the EUV observations by EUVE.

### 2.1 ASCA

Mrk 421 was observed with ASCA for 7 days continuously (except for the Earth occultation) during 1998 April 23.97 – 30.8 UT (PI: T. Takahashi). The ASCA detectors – two SISs (Solid-state Imaging Spectrometers; Burke et al. 1991; Yamashita et al. 1997) and two GISs (Gas Imaging Spectrometers; Ohashi et al. 1996) were in operation. We refer to Tanihata et al. (2001) to the details of the observation modes and data reduction. The obtained spectra were rebinned such that all bins have the same statistics.

Since late 1994, the efficiency of both SIS detectors below $``$1 keV has been decreasing over time. The problem is believed to be caused by the increased residual dark current level, and also by the decrease in the charge transfer efficiency, although effort is still underway to model the effects. While it has been reported that it may be possible to parameterize the degradation by an additional absorption column as a function of time (Yaqoob, 2000), since our observation concerns the precise shape of the continuum spectrum, we will use only the data obtained by the GIS for spectral analysis. (Note that the X–ray spectra in Figure 2 of Takahashi et al. (2000) show the data obtained by SIS. For this, there is a too strong decrease of inferred flux at the lower energy, which is due to instrumental effects described above.)

### 2.2 RXTE

The RXTE observations (PI: G. Madejski) were mainly coordinated to coincide with the ASCA observation period, but were conducted mostly only during the low-background orbits. The RXTE extends the ASCA bandpass to a higher energy with the instruments PCA (Proportional Counter Array; 2–60 keV) and HEXTE (High Energy X–ray Timing Experiment; 15–200 keV). However, because both the flux from the source and the sensitivities of the detectors decrease towards higher energies, in the analysis below we use only the data from the PCA, in the energy range of 3–25 keV.

All data reduction was performed using the HEASOFT software packages, based on the REX script, provided by the RXTE Guest Observer Facility (GOF). The screening criteria excluded the data when the elevation angle from the Earth’s limb was lower than 10°, and the data acquired during the SAA passages. Since the PCUs are activated during the SAA passage, we used the data only after when the background level dropped to the quiescent level and became stable, which is typically 30 minutes after the peak of a SAA passage.

For the best signal-to-noise, we selected events only from the top xenon layer, and excluded events from other gas layers. For the maximum consistency among the observations, we used only the three detectors which were turned on for all the observations (known as PCU0, PCU1, and PCU2). Since the PCA is not an imaging detector, the blank-sky observations must be used to estimate the background spectrum. For this, we used the background models provided by the RXTE GOF to evaluate the background during each observation. We note that the residual background contamination of the data within the energy bandpass selected by us was always significantly lower than the source counts. The mean net RXTE PCA count rate from Mrk 421 in the 10 - 20 keV range for 3 PCUs as above was at least $`3`$ counts s-1 while the total residuals from blank-sky in this energy range $`\mathrm{𝑝𝑙𝑢𝑠}`$ the fluctuations of the Cosmic X-ray Background were lower than 0.2 counts s-1 (see HEASARC Web page regarding the RXTE PCA background), providing the confidence that an imprecise background subtraction did not affect our results.

### 2.3 BeppoSAX

The BeppoSAX observations (PI: L. Chiappetti) started 3 days before the ASCA observations (Maraschi et al. 1999; Fossati et al. 2000a,b). We use the data obtained by the LECS (Low Energy Concentrator Spectrometer; 0.1–10 keV), which extends the bandpass to lower energy, and also the MECS (Medium Energy Concentrator Spectrometer; 1.3–10 keV). The data analysis was based on the linearized event files, together with the appropriate background event files and response matrices, provided by the BeppoSAX Science Data Center (SDC). Data reduction was performed using the standard HEASOFT software, XSELECT.

We extracted the source photons from a circular region centered at the source with a radius of 8 arcmin for both LECS and MECS. For the background photons, we constructed the spectra by extracting the photons from the same detector regions using blank-sky observations. Each spectrum was rebinned using the grouping templates provided by the BeppoSAX SDC, and we limited the energy range to 0.12–4.0 keV for LECS, and 1.65–10.5 keV for MECS, as suggested in the Cookbook for BeppoSAX NFI Spectral Analysis (1999).

### 2.4 EUVE

The EUVE observations covered the entire campaign of Mrk 421 from April 19.4 to May 1.2 (PI: T. Takahashi). The source photons were extracted in a 12 arcmin aperture. In order to estimate the absolute flux in the EUV range, we normalized it to be consistent with the simultaneous data obtained by the other observatories. Given that the EUV spectrum is difficult to determine within the rather limited bandpass, we used only the observed integrated count rate (cf. Marshall, Fruscione & Carone 1995). For this, we used the period where the observations by all observatories overlapped (MJD 50927.03 - 50927.34; §3.1). During this period, the average EUVE count rate was 0.46$`\pm `$0.01 cts/s. At the same time, the 0.13–0.18 keV flux derived from the combined LECS, MECS, GIS, and PCA spectrum gives $`\nu F_\nu =(1.78\pm 0.04)\times 10^{10}`$ erg cm-2 s-1, where the error represents the statistical error.

An additional uncertainty arises from the fact that the difference in the spectral shape is not taken in to account. This can be estimated from the results in Marshall, Fruscione & Carone (1995). Taking $`N_\mathrm{H}=1.5\times 10^{20}`$ cm-2 for Mrk 421, assuming that the spectrum index varies from 0.5 to 1.0, the flux density would vary from 1363–1394 $`\mu `$Jy for 1 cts/s in the EUVE range. This gives an additional 2.5% error as a systematic error.

In the following analysis, we use the above value for scaling the observed EUVE count rate to flux. Thus,

$$\nu F_\nu (\mathrm{EUVE})=1.78\times 10^{10}\left(\frac{F_{\mathrm{cts}/\mathrm{s}}}{0.46}\right)\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2,$$
(1)
where $`F_{\mathrm{cts}/\mathrm{s}}`$ is the EUVE count rate.

## 3 RESULTS

### 3.1 The Shape of the 0.1–25 keV Spectrum

During April 24, observations by all 3 X–ray satellites, ASCA, BeppoSAX, RXTE, and also EUVE happened to exactly overlap continuously for 26.5 ksec (MJD 50927.03 - 50927.34). This provides us with high quality data over a very broad energy coverage, spanning more than 2 decades, from 0.12 to 25 keV. The net exposure for each instrument during the overlapping 26.5 ks period were 9.1, 13.3, 3.6, and 5.5 ksec for GIS, PCA, LECS, MECS, respectively.

We first fitted the spectrum for each instrument with a power-law model, with a low energy absorption fixed to the Galactic value, $`N_\mathrm{H}=1.5\times 10^{20}`$ cm-2 (Elvis, Wilkes, & Lockman, 1989). The results are shown in the top half of Table 1, suggesting that the fit is not adequate for any of the individual instruments. This certainly suggests that a more complex model, rather than just a simple power-law function, is required to describe the spectrum.

The shape of the synchrotron spectrum is determined from the energy distribution of the radiating electrons, which cuts off at the high energy representing the acceleration limit. As described in the introduction, the high energy end of the synchrotron spectrum in TeV blazars is located in the X–ray range; this is most likely the reason that the high-quality X–ray data cannot be described by a single power-law.

Given the curved spectrum for each instrument, we fitted the combined X–ray spectrum with various functions to model the shape of the observed curvature. We note that our aim is not to develop the shape of the curvature from first principles, but rather, to reproduce the observed data. For this reason, we considered the following models for comparison: (1) A broken power-law, (2) a curved power-law, (3) a power-law with an exponential cutoff, and (4) a quadratic function in $`\mathrm{log}\nu `$-$`\mathrm{log}`$$`\nu F_\nu `$ plane. The curved power-law that we define below is similar to the broken power-law model, but the change in the spectral index is not discrete. The form attempted by us is $`F(E)=KE^{\mathrm{\Gamma }_1}(1+\frac{E}{E_{\mathrm{br}}})^{\mathrm{\Gamma }_2}`$. The form for the exponentially cutoff power law is $`F(E)=KE^{\mathrm{\Gamma }_1}\mathrm{exp}(E/E_{\mathrm{br}})`$, and the quadratic function is described as $`\mathrm{log}\nu F_\nu (E)=log(\nu F_{\nu ,\mathrm{peak}})a(\mathrm{log}E\mathrm{log}E_{\mathrm{peak}})^2`$.

Since there are slight differences in the cross-calibration between the instruments, we allow the normalization to vary for each instrument. For all results below, the flux will be normalized to the GIS value, unless otherwise noted. The $`N_\mathrm{H}`$ is fixed to the Galactic value.

By fitting the spectrum with the models above, we first found that the power-law with an exponential cutoff is not a good representation of the data: the model cutoff is too sharp at high energies. The broken power-law model is also not acceptable, and the curved power-law is significantly better. We also find that the quadratic function fits the spectrum surprisingly well. The fitted spectrum and the residuals to the quadratic model are shown in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign, where all LECS, MECS, GIS, and PCA spectra are shown to be well described with this single function. The results of the spectral fits by the 4 models are summarized in in the bottom half of Table 1.

### 3.2 The Time Resolved Spectra

The EUVE observation covered the whole campaign (April 19.4 to May 1.2), and thus completely overlapped our 7-day ASCA observation. This adds a data point at 0.13-0.18 keV in addition to the ASCA spectrum in the 0.7–7 keV range, providing us the best opportunity to study the continuous spectral variability. In fact this is the longest continuous monitoring of this object so far spanning over the 0.13–7 keV range.

We divided the whole observation into shorter segments with each having a duration of 10 ksec. In order to combine the GIS and EUVE data, the GIS count rate spectra must first be converted to the source flux spectra as the EUVE count rate data converts directly into flux (§2.4). For this, we first fit the GIS spectrum with a curved power-law model with Galactic absorption. This gives the source spectrum in the GIS energy range. Note that this procedure depends on the selected model, but long as it describes the spectrum well, this should introduce the smallest bias. We then combine this with the EUVE data, which gives the combined source spectrum.

In fitting the combined EUVE/GIS spectra, we first started with a single power-law model. In similarity to the results in §3.1, none of the 56 spectra could be fit with a single power-law function, with their average $`\chi _\nu ^2`$ of 7.7, implying the curvature of the spectra for all data sets. Given that a quadratic model fits the spectrum equally well as the curved power-law model, but with fewer parameters (see §3.1), we will assume the quadratic model for all data sets hereafter. In fact, because of the energy gap between the EUVE and GIS, the fit with a curved power-law model resulted in having multiple minima in the $`\chi ^2`$ plane (since the parameters are correlated), and the exact value of the best-fit parameters depended highly on the initial value.

We checked whether the fit was adequate for all data segments. The $`\chi _\nu ^2`$ for each fit (shown in the bottom panel in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign) ranged from 0.71 to 1.45. The number of degrees of freedom ranges from 60 to 187, which gave $`P(\chi ^2)>`$5% for 53 out of the total of 56 data points. For the other 3 data points (corresponding to times 230–240, 270–280, and 450–460 ksec), we found several line or edge-like structures in the spectra, which resulted in the increase of $`\chi _\nu ^2`$. Although it is not clear whether this is due to some systematic effect or if it is a true property of the source, since our aim is to model the continuum spectrum correctly, we excluded the energy ranges where those features are present. This improved the fit significantly, which is shown as the crosses. With this, all the 56 spectra were well described with the quadratic function.

The time evolution of the derived parameters in the fit is shown in the top 3 panels in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign; the peak flux, peak energy, and the curvature parameter $`a`$ (see definition in Table 1). The peak energy is observed to shift between 0.5 keV and 1.2 keV, in which a general trend is found such that the peak energy is higher when the peak flux is higher. This is demonstrated as the filled circles in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign. The fit results are listed in Table 2.

We do note that the actual validity of using a quadratic function in the entire 0.1–25 keV range was checked in detail (including an examination of trends in residuals) only for the epoch when the overlapping BeppoSAX, ASCA, and RXTE observations were available. Since we measured only the EUVE flux below 0.7 keV in the combined EUVE/GIS fit, the derived peak location can depend on the model function we assume. For instance if $`any`$ function was allowed, the peak can be located at any energy within the gap. For this, we emphasize that the derived peak location is the estimated peak, assuming that the spectrum gap can be extrapolated with the same quadratic function, and also that the error bar describes only the statistical error.

We also analyzed the BeppoSAX data during the campaign. Furthermore, in order to make a comparison with data obtained from another campaign, we also analyzed the BeppoSAX data obtained in 1997 (PI: G. Vacanti, L. Chiappetti). We divided the observation into segments of the same 10 ks duration, and fit the combined LECS/MECS spectra with a quadratic function. For the BeppoSAX data, we also fit the spectra with a curved power-law model for comparison. The derived peak flux and peak energy from both models are listed in Table 3. In fitting the combined LECS/MECS spectra, the normalizations are allowed to vary independently, and the values in Table 3 are normalized to the MECS2 value. The fitting results show that although the spectral fit parameter is one less in the quadratic function, both models equally represent the observed spectra. In contrast to the combined EUVE/GIS data, there is no energy gap in the BeppoSAX data, which provides support that our assumed quadratic form is appropriate in estimating the peak location in the EUVE/GIS case. The derived peak flux and peak energy using both functions are similar, but we also note that the quadratic function gives a slightly lower peak energy than the curved power-law function.

The derived peak flux and peak energy assuming the quadratic model for the BeppoSAX observations are shown in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign as the open circles (for 1998) and squares (for 1997). Note that the BeppoSAX data points in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign are multiplied by a normalization factor of 1.07, which we adopted from the combined fit in §3.1 It can be seen that the results during the 1998 campaign lie on a similar slope as the results derived from the EUVE/GIS data during the same campaign.

The relation of peak flux $`\nu F_{\nu ,\mathrm{peak}}`$ and peak energy $`E_{\mathrm{peak}}`$ for the EUVE/GIS data set can be best described with a power-law function of $`E_{\mathrm{peak}}\nu F_{\nu ,\mathrm{peak}}^\alpha `$, where $`\alpha `$=0.72$`\pm `$0.02 (or $`ϵ`$=0.76$`\pm `$0.02 in form $`E_{\mathrm{peak}}F_{0.110\mathrm{k}\mathrm{e}\mathrm{V}}^ϵ`$; see description below). By including also the BeppoSAX data points during 1998, the index becomes $`\alpha `$=0.77$`\pm `$0.02 ($`ϵ`$=0.72$`\pm `$0.02). It can be seen that the data from 1997 lie on a somewhat different slope. The relation derived from the 1997 data alone gives $`\alpha `$=0.96$`\pm `$0.09 ($`ϵ`$=0.79$`\pm `$0.06), slightly steeper than the index derived from the 1998 campaign. Fitting all data together we obtain $`\alpha `$=0.59$`\pm `$0.01 ($`ϵ`$=0.50$`\pm `$0.01). All errors correspond to a 1$`\sigma `$ error.

A similar analysis was performed by Fossati et al. (2000b) for the 1997 and 1998 BeppoSAX data. Fossati et al. (2000b) note a tight relationship between the measured integrated 0.1-10 keV flux $`F_{0.110\mathrm{k}\mathrm{e}\mathrm{V}}`$ and the peak energy $`E_{\mathrm{peak}}`$ such that $`E_{\mathrm{peak}}F_{0.110\mathrm{k}\mathrm{e}\mathrm{V}}^ϵ`$ with $`ϵ`$=0.55$`\pm `$0.05, similar to the index derived from our results including all data sets. For a direct comparison, we fit simultaneously both 1997 and 1998 for BeppoSAX data points only. This gave $`\alpha `$=0.52$`\pm `$0.02 ($`ϵ`$=0.45$`\pm `$0.01). The estimated $`ϵ`$ is slightly smaller than reported in Fossati et al. (2000b), which might be affected by different subdividing of the observations. Nonetheless, it is worth noting that the slope derived from the combined data appears be softer (smaller $`\alpha `$ or $`ϵ`$) than each epoch considered alone.

Furthermore, the relation appeared to also vary between individual flares during the 1998 campaign. In some flares the peak energy clearly increases together with flux, while the correlation is weaker in others. This can be seen clearly in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign, which shows the evolution of the synchrotron peak in each flare. The peak energy has a trend to shift to the higher energy and then return to the initial value together with the flux, but there is no common slope.

### 3.3 Spectral Variations at Flares

In the previous section, we described the spectral variations as the evolution of the location of the peak in the $`\nu F_\nu `$ spectrum. Here, we take advantage of the wide energy range coverage afforded by the use of by EUVE, ASCA, and RXTE data to investigate spectral variations in more detail. Given that the RXTE observations happen mostly during the flare minima rather than at their peaks (see Figure 1 in Takahashi et al. 2000), here we focus on the beginning (rising portion) of the flares.

We divided the observations into 10 ksec segments, similar to the previous section. For each segment, the GIS and PCA spectra are fitted simultaneously with a curved power-law function with an additional constant factor for the normalization. The derived source spectra in the GIS and PCA range are combined with the EUVE data point, which gives the total source spectrum in the energy range of $``$0.1–25 keV. In the top panels of Figures Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 CampaignEvolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign, we show the evolution of the $`\nu F_\nu `$ spectrum derived in this manner for the 2nd, 4th, and 5th flare (starting from 80 ksec, 260 ksec, and 340 ksec, respectively, as defined in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign). Note that the spectra are further rebinned for this Figure. For each flare, we also plot a subsequent spectrum, but those only include the EUVE and GIS data.

To see the evolution more clearly, we calculated the ratio of the rebinned spectra by dividing these “subsequent” spectra by the spectrum measured at the onset of a given flare (integrated over 80–90, 260–270, 340–350 ksec, respectively). This is shown in the bottom panels of Figures Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 CampaignEvolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign. Interestingly, these three flares all showed different features. It can be seen that Flare 2 starts with a hardening at the very highest energy range with a concave shape in the ratio. As time progresses, the whole energy range becomes dominated by a new spectral slope. In contrast to this, in Flare 4 the lower energy photons appear to precede the total flux rise. Flare 5 was similar to Flare 2, but all energy bands increase together, with a larger amplitude at higher energies. We also studied the spectral variation for other flares, but a detailed investigation was rather difficult without a denser temporal coverage by RXTE.

### 3.4 The X–ray and TeV Correlation

Since X–ray and TeV energies both reflect the electrons accelerated to the highest energies in TeV blazars, simultaneous X–ray and TeV observations provide an important tool to probe the emission region of the jet (e.g., Krawczynski et al. 2002, and references therein). On the other hand, since blazars are highly variable on short time scales, exactly simultaneous observations are essential. Thanks to the relatively high brightness of the source during the 1998 campaign, it was possible to generate a significant measurement for the separate TeV pointings every night (Aharonian et al. 1999; Takahashi et al. 2000). Here we provide the X–ray spectrum from the exactly simultaneous time regions with the TeV observations.

For each TeV observation by HEGRA (Aharonian et al., 1999) and Whipple (Takahashi et al., 2000), we collected all existing X–ray data during the TeV exposure, from the ASCA, RXTE, and EUVE observations (the CAT observations mostly overlap with the HEGRA observations). In particular, all Whipple pointings during the ASCA observations happened to coincide with the RXTE observations, providing an X–ray spectrum ranging from 0.1–25 keV, simultaneous with the TeV observation.

We generated the combined X–ray spectrum using the same procedures as described in the previous sections. We then fitted the X–ray spectrum for each period assuming a quadratic function. The X–ray spectra during the HEGRA observations (consisting of the EUVE and ASCA data) are shown in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign(a). The spectrum for each observation is shown as different symbols, and the dashed line shows the modeled quadratic function derived via the fit. The X–ray spectra during the Whipple observations, which consist of the EUVE, ASCA and RXTE data are shown in Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign(b). Note that the spectra are further rebinned in the figures. While each integration time is as short as several hours, it can be seen that the curved shape of the high energy peak of the synchrotron spectrum is clearly resolved, and it is evident that the peak is located in the X–ray range.

The results of fits to each X–ray spectrum concurrent with the TeV pointings are summarized in Table 4. We also list the simultaneous TeV fluxes, in which the HEGRA data were taken by Aharonian et al. (1999). In Figure Evolution of the Synchrotron Spectrum in Mrk 421 during the 1998 Campaign, we show the correlation of the derived location of the synchrotron peak flux with the TeV flux. Since the sensitive energy ranges of the two TeV telescopes are slightly different, we used the integrated flux in Crab units. A good correlation can be seen, showing that the TeV flux is higher when the synchrotron peak flux is higher. A best fit to the relation with a power-law function gives $`F_{\mathrm{TeV}}F_{\mathrm{sy},\mathrm{peak}}^\alpha `$ with $`\alpha `$=1.7$`\pm `$0.3.

## 4 SUMMARY AND DISCUSSION

We have performed detailed analysis of the combined EUVE, ASCA, BeppoSAX, and RXTE data collected during the long look campaign of Mrk 421 in 1998, resulting in the measurement of the time-resolved spectrum in a broad energy range of 0.1–25 keV in segments as short as several ksec. These are among the highest quality spectra with regard to both photon statistics and energy coverage so far for any blazar, providing the precise spectral shape at the highest energy end of the synchrotron spectrum.

We have shown that both curved power-law function and quadratic function in $`\mathrm{log}\nu `$$`\mathrm{log}`$$`\nu F_\nu `$ space reproduce the combined EUVE, BeppoSAX, ASCA, and RXTE spectrum. For the epochs where BeppoSAX data do not exist, we assumed that the spectral gap between $`\mathrm{𝐸𝑈𝑉𝐸}`$ and $`\mathrm{𝐴𝑆𝐶𝐴}`$ can be extrapolated with a quadratic function, and we showed that the energy and the flux at the synchrotron peak in the $`\nu F_\nu `$ spectrum are correlated, showing an overall trend of a higher synchrotron peak energy for higher peak flux, but the details of this relation differed from flare to flare. In particular, within the three flare beginnings, one flare started from a rise in the higher energy while another started from a rise in the lower energy. The relative amplitude of the rise in different energies also differed from flare to flare. An interesting result was shown from the 2nd flare, where it was indicated that this flare started with a hardening at the higher energy. This feature is difficult to describe by a simple cooling or acceleration of a single electron distribution (see e.g. Kataoka 2000), and thus this concave curvature indicates that two different spectral components were likely to co-exist at the beginning of the flare. We remark that such a behavior is observed in only one flare, but the importance is that we observed a flare that could not be explained by a single electron distribution. This suggests an appearance of a new component, which generated the observed flare.

Recently, many theoretical studies has been performed to model the energy dependence, or time lags observed in the day-scale flares in blazars. One approach considers models with a homogeneous emission region (Mastichiadis & Kirk, 1997; Chiaberge & Ghisellini, 1999; Li & Kusunose, 2000). For instance, Mastichiadis & Kirk (1997) have suggested that a sudden increase in the maximum energy of the accelerated electrons could result in the X–ray and TeV flare observed in Mrk 421. Similar results were shown by Kataoka et al. (2000), in reproducing the soft-lag observed in PKS 2155–304. However, most of these models assume a change in some parameter in the emission region (such as the magnetic field, maximum energy of electrons, number density of electrons, etc.). The observational results for Mrk 421 described above lead us to suggest a different scenario.

An alternative scenario has each flare forming as a result of a separate electron distribution. This was suggested in modeling the high state SED of Mrk 501 (Kataoka et al., 1999). They applied the one-zone homogeneous synchrotron self-Compton model to the SED from April 1997, and concluded that a single electron distribution is insufficient to reproduce the observed synchrotron spectrum. This scenario could work for instance when the jet is emitted intermittently from the central engine, and several separate regions generate different emission components. Another viable scenario can be provided by an internal shock model, where the light curve results from a superposition of many flares due to the collisions of shells, which may occur when a faster shell catches up to a slower shell (Ghisellini 2001, Spada et al. 2001; Sikora et al. 2001)

In fact, Tanihata et al. (2003) has shown from simulations that the internal shock model can naturally explain various variability properties observed in TeV blazars. One observational fact is that the X–ray flares always appear to lie on top of an underlying offset-like component. In the internal shock model, flares can be considered as arising from collisions of shells that had the largest difference in the initial velocity. All other collisions generate smaller amplitude, longer flares, which pile up to generate the offset component. In this case, it was shown that the flares will tend to have a higher synchrotron peak energy as compared with that of the offset component.

Assuming that flares are due to an emergence of a new component with a higher synchrotron peak energy than the pre-existing component, the synchrotron peak energy in the observed spectrum would appear as shifting to a higher energy. It is interesting to note that our observations are consistent with this, with the overall trend of the synchrotron peak energy being higher during flares. Furthermore in this case, the spectral evolution would depend on the relation of the new and pre-existing spectra. This would naturally explain the different spectral evolution among different flares.

Finally, we remark again that due to the data gap between the EUVE and ASCA data in our observation, the results concerning the spectrum peak are somewhat tentative. Because of the rapid variability of the source in short time scales, the effective area of the instrument is most important. Long-look observations by the new observatories such as $`XMM`$-$`Newton`$ should provide data allowing a more detailed analysis of the flaring mechanism, leading to the dynamics of the accelerated electrons in blazar jets.

We thank the anonymous referee for numbers of constructive comments to improve the paper. Support for this work was provided by the Fellowship of Japan Society for Promotion of Science for Young Scientists, and by the Department of Energy contract to SLAC no. DE-AC3-76SF00515. We also acknowledge the NASA Chandra grants via SAO grant no. GO1-2113X.


format: text

text:
# 1 Illustration of the concept of overlap in a molecular network. For a pair of paralogs the overlap Ω is defined as the number of common neighbors they have in the network. In the case of transcription network the regulatory overlap Ω_{𝑟𝑒𝑔} counts transcription factors regulating both paralogs, while for the physical interaction network the interaction overlap Ω_{𝑖𝑛𝑡} counts their common binding partners. The pair of paralogs used in this illustration has the overlap Ω=2 out of the total of 5 distinct neighbors of the pair. That corresponds to a normalized overlap of 2/5=0.40.

Upstream Plasticity and Downstream Robustness
in Evolution of Molecular Networks.

Sergei Maslov∗)
Department of Physics, Brookhaven National Laboratory,
Upton, New York 11973, USA
Email: maslov@bnl.gov

Kim Sneppen
Nordita, Blegdamsvej 17, 2100 Copenhagen Ø, Denmark
Email: Sneppen@nbi.dk

Kasper Astrup Eriksen
Department of Theoretical Physics, Lund University,
Sölvegatan 14A, SE-223 62 Lund, Sweden.
and

Department of Physics, Brookhaven National Laboratory,
Upton, New York 11973, USA
Email: kasper@thep.lu.se.

∗)Corresponding Author

Abstract

Background:
Gene duplication followed by functional divergence of associated proteins is a major force shaping molecular networks in living organisms. Recent availability of system-wide data for yeast S. Cerevisiae allow us to access the effects of gene duplication on robustness and plasticity of molecular networks.

Results:
We demonstrate that the upstream transcriptional regulation of duplicated genes diverges fast, losing on average 4% of their common transcription factors for every 1% divergence of their amino acid sequences. In contrast, the set of physical interaction partners of their protein products changes much slower. The relative stability of downstream functions of duplicated genes, is further corroborated by their ability to substitute for each other in gene knockout experiments.

Conclusion:
Apparently the upstream regulation of genes evolves much more rapidly than the downstream functions of the associated proteins. This is in accordance with a view where it is regulatory changes that mainly drives evolution. Any evolutionary model has eventually to account for this disparity and we have here quantified its size on a genome wide scale. In this context a very important open question is to what extent our results for duplicated genes within yeast (paralogs) carries over to homologeous proteins in different species (orthologs).

Key words: Network, Evolution, Gene duplication, Paralogs, Alternate pathways.

Background

Biological processes are rarely performed by single isolated molecules. Instead, they typically involve a coordinated activity of many molecules forming a neighborhood in biomolecular networks. Changes in molecular networks are thus coupled to evolution of new functions and functional relationships in the organism. Gene duplication is an important source of raw material for the evolutionary development of molecular networks in a given species . Immediately after a duplication event the pair of duplicated genes is thought to be identical in both sequences and functional roles in the cell. However, with time their properties including their position within a network diverge. Here we quantify this divergence in the yeast Saccharomyces Cerevisiae using several recent system-wide data sets. To this end we measure: 1) The similarity of positions of duplicated genes in the transcription regulatory network  given by the number of transcription regulators they have in common; 2) The similarity of the set of binding partners  of their protein products, and their ability to substitute for each other in knock-out experiments . These measures reflect, correspondingly, the upstream and downstream properties of molecular networks around the duplicated genes.

Results

The first measure of divergence of duplicated genes compares sets of their transcriptional regulators. Such a set contains information about different conditions under which the gene is expressed, and thus reflects its functional roles in the cell. To quantify the similarity of transcriptional regulation of a pair of genes we introduce the concept of “regulatory overlap” $`\mathrm{\Omega }_{reg}`$ given by the number of transcription factors that bind to the upstream regions of both these genes (see Fig. 1 for a general illustration).

The information about duplicated genes used in this study was extracted from the list of all pairs of paralogous (evolutionary related) proteins found in the yeast genome , while the system-wide data for its transcription regulatory network was taken from the Ref.  (see Methods for more details). Fig. 2A shows the distribution of the regulatory overlap for different values of the percent identity (PID) of amino acid sequences of paralogous proteins.

From this figure one can see that the regulatory overlap decreases as a function of the PID. While multiple overlaps dominate the distribution for PID $``$ 80%, at lower values of PID they disappear in favor of smaller overlaps.

To quantify the average rate of loss of the regulatory overlap $`\mathrm{\Omega }_{reg}`$, Fig 2B shows its average value as a function of PID. The regulatory overlap in this plot is normalized by the ancestral connectivity of a gene, estimated as the total number of distinct transcription factors that are involved in regulation of at least one of the pair of proteins (see Fig 1). One interesting feature of Fig. 2B is that even pairs of proteins whose amino acid sequences are 100% identical to each other have only 30% overlap in their upstream regulation. The overlap becomes 60% when measured in units of the smaller among the two numbers of regulatory inputs of a pair of paralogs. The second feature of Fig. 2B is a gradual decline of the average regulatory overlap over the whole range of sequence similarities. The data in Fig. 2B may be fitted with an exponential decay with a rate corresponding to a 4% chance of losing a given common regulator of a pair for every 1% decrease in their amino acid sequence identity. Thus already at PID$`=`$80% half of the regulatory overlap present at PID$`=`$100% is lost. The decline in the regulatory overlap presented in Fig. 2A,B is in accord with a recently published analysis of similarity of microarray expression patterns  of paralogs. In fact, due to a more direct information about the gene regulation contained in the dataset of Ref. , our analysis extends the gradual decline to much lower PID than could be detected  from the microarray data. After we submitted this manuscript for the first time another study has observed a rapid decline in the overlap of upstream motivs . As this study was carried out as a function of Ks the long time behaviour can not be extracted.

We now consider the second measure of divergence, thus concentrating on downstream functional properties of duplicated genes. Such properties are in part reflected in their set of physical interaction partners of their protein products. The similarity of physical interaction neighborhoods of a pair of proteins can be quantified as the “interaction overlap” $`\mathrm{\Omega }_{int}`$ given by the number of proteins that bind to both of them (See Fig. 1). In our study we use the system-wide information about protein-protein physical interactions obtained by combining two high throughput two-hybrid datasets . Fig 3A shows the average value of the interaction overlap $`\mathrm{\Omega }_{int}`$ between duplicated genes. Again it decreases with decreasing PID, reflecting gradual loss/change of physical interaction partners of proteins in the course of evolution.

A similar analysis, but as a function of $`K_s`$ – the number of silent nucleotide substitutions per site – was previously reported by Wagner . In agreement with that study, we find that gene duplicates are more likely to share interaction partners than one expects by pure chance alone (the horizontal line in the Fig. 3A). Our set of yeast paralogs contains 366 pairs with both proteins present in the combined set of . Out of these pairs 138 (38%) share at least one interaction partner. We find that the decrease in average overlap becomes systematic only for PID$`<60`$%, while above 60% it remains roughly constant in both Uetz , Ito , and combined datasets (Fig. 3A). Unlike in Fig. 2B, the overlap shown in Fig 3A is not normalized by the total number of distinct interaction partners of a paralogous pair. The corresponding normalized plot confirms the above main conclusions.

An alternative way to quantify the extent of the divergence/redundancy of downstream functions of a pair of duplicate genes is to examine the viability of a null-mutant lacking one of them. Gu, et al.  recently analyzed the fraction of yeast genes essential for the survival of the cell separately for singleton genes that lack a duplicated partner (paralog) in the yeast genome and the genes that retained at least one such partner. It was found that the fraction of essential genes is approximately 4 times higher among singleton genes than among ones protected by a highly similar paralog. It was also demonstrated that the protective role of a paralog persists down to rather low levels of the amino acid sequence similarity (PID). In Fig. 3B we confirm these findings using a different set of lethal and viable null-mutants of Ref.  as well as demonstrate that the magnitude of the effect is the strongest among nuclear proteins, where the largest fraction of essential proteins resides. Notice that the lethality (especially that of nuclear proteins) shows a dramatic change at PID around 60%, indicating that paralogs with a higher level of similarity can typically substitute for each other.

Discussion

Having presented different measures of the upstream and the downstream divergence of duplicated genes we are now in a position to discuss them in a wider context. Comparing Fig. 2B to Fig. 3A,B one concludes that changes in the upstream regulation of duplicated genes happen more readily than changes in their downstream function. The overlap in the set of binding partners (Fig. 3A) and the ability of duplicates to substitute for each other (Fig. 3B) remain virtually constant down to PID of 60%, while the average regulatory overlap at this PID has dropped to about 20% of its maximum (Fig. 2B). Thus our results indicate that duplicated genes would still have the ability to partially substitute for downstream functions of each other at the time when the repertoire of their regulatory connections has substantially changed. Such genes would be less constrained in evolving new functions , and thus would contribute to a greater evolutionary plasticity of the network.

Conclusions

The evolution of a biological species organism modifies it on multiple levels ranging from sequences of individual molecules, to their coordinated activity in the cell (molecular networks), all the way up to the phenotype of the organism itself. While its manifestations both on the level of sequences and phenotypes are well documented, the data needed to quantify the evolutionary changes on the level of molecular networks have appeared only very recently. System-wide studies such as high-throughput two hybrid assays of protein-protein interactions , the large scale study of transcriptional regulation , and the whole genome assay of the viability of null-mutants in yeast  have allowed us to go beyond describing particular cases of evolution of molecular networks and look at its large scale dynamics.

Our results show that genetic regulations of duplicated proteins in yeast change faster than both their amino acid sequences and their protein interactions partners. It is tempting to extend this observation to pairs of homologous proteins in different species (orthologs) that diverged from each other as a result of a speciation (as opposed to gene duplication) event. This would help to explain how species with very similar gene contents can evolve novel properties on a relatively short timescale. However, such an inter-species comparison of molecular networks has to wait for a completion of large-scale studies of closely related model organisms.

Methods:

As a source of information about duplicated genes we use the set  consisting of 4443 pairs of paralogous yeast proteins. This set was obtained by blasting all yeast proteins against each other with a conservative E-value cutoff of $`10^{30}`$. We curated this dataset by removing 72 known transposable elements, http://genome-www.stanford.edu/Saccharomyces, and their homologs (108 proteins all together). That left us with 2739 pairs of paralogous yeast proteins formed by 1891 proteins (about 30% of the genome) with at least one homolog in the yeast genome. Pairs of paralogous proteins in this set are characterized by a fairly broad distribution of the percent identity (PID) of their amino acid sequences in the interval from 20% to 100%, as shown in Fig. 4.

It is worthwhile to note that while the number of silent substitutions $`K_s`$ per site in a pair of duplicated genes is commonly used as a proxy of the time elapsed since the duplication event , the PID (or the number of non-silent substitutions per site $`K_a`$=1-PID/100) is rather a very crude estimate of the extent of their functional similarity.

The system-wide data describing the transcription regulatory network of yeast was taken from the Ref. , which reports the in-vivo study of binding between 106 transcription factors to the upstream regulatory regions of genes encoding all 6270 of yeast proteins. Since the number of transcriptional regulators in this dataset is quite large, the probability that by pure chance the same transcription factor would be incorrectly detected among upstream regulators both duplicated genes is relatively small. Thus the contribution of false positives to the regulatory overlap $`\mathrm{\Omega }_{reg}`$ is insignificant. On the other hand, false positives significantly affect the average number of regulatory inputs of each individual proteins. Thus the presence of a considerable fraction of false positives would manifest itself in sensitivity of the normalized regulatory overlap with respect to the P-value cutoff of the data. We failed to find any such dependence when investigating regulatory overlap for different P-values (data now shown): While the average number of regulations per gene decreases six-fold (from 2 to 0.33) when the P-value cutoff is lowered from $`10^2`$ to $`10^4`$, both the initial drop and the rate of the exponential decay fit in Fig. 2B remains virtually unchanged. This suggests that false positives are not a significant part of the experimental dataset of Ref.  at least up to $`10^2`$, and validates the robust nature of the parameters extracted from the Fig. 2B.

As a source of information about binding partners of yeast proteins we combine two high throughput two-hybrid datasets: the core Ito et al. set  (806 interactions among 797 proteins) and the extended Uetz et al. dataset  downloaded from the website of this group (1446 interactions among 1340 proteins) The resulting dataset consists of a total of 1734 proteins joined by 2111 non-redundant interactions. Using this combined dataset we found that even 100% identical proteins share on average only 30% of their binding partners. However, unlike in the case of the upstream regulation, the set of interaction partners of a protein in a two-hybrid experiment is essentially determined by its amino acid sequence. We attribute the 30% overlap in the set of binding partners of identical proteins to false positives/negatives inevitably present in high-throughput two-hybrid experiments. The presence of false positives/negatives also manifests itself in the fact that two independent system-wide two-hybrid experiments ,  have only 141 interacting pairs in common.

The system-wide data on viability of null-mutants used in our study was obtained from Ref.  in which 1103 essential (non-viable null-mutants) and 4678 non-essential (viable null-mutants) yeast proteins were reported. Actual lists of viable and non-viable null-mutants as discovered in Ref.  were downloaded from the Saccharomyces Genome Database http://genome-www.stanford.edu/Saccharomyces.

Authors Contributions

All authors contributed to both the ideas and writing the manuscript in close collaboration. All authors read and approved the manuscript.

Acknowledgment

Work at Brookhaven National Laboratory was carried out under Contract No. DE-AC02-98CH10886, Division of Material Science, U.S. Department of Energy. Two of us (K.E and K.S.) thank the Institute for Strongly Correlated and Complex Systems at Brookhaven National Laboratory for financial support during visits when part of this work was completed. S.M. and K.S. acknowledge the support of the NSF grant PHY99-07949 (work at the KITP, University of California at Santa Barbara). We thank John Little for critically reviewing the manuscript.


text:
# Silent Phase Qubit Based on 𝑑-Wave Josephson Junctions

## Abstract

We report on design and fabrication of a new type of flux qubit that capitalizes on intrinsic properties of submicron YBCO grain boundary junctions. The operating point is protected from the fluctuations of the external fields, already on the classical level; the effects of external perturbations are absent, to the second or third order, depending on the character of the coupling. We propose an experiment to observe quantum tunneling and Rabi oscillations in the qubit. Estimate of the decoherence due to fluctuations of the external flux is presented.

Over the last few years a series of experiments expt ; mooij ; vion  provided a conclusive evidence of quantum superposition in meso- and macroscopic superconductors. The task at hand is scaling up of the system, with two goals in mind: to probe how far quantum superposition can be pushed into macroscopic world, and to develop an element base for quantum computing, which only becomes viable on the scale 10–100 qubits.

It was suggested, that use of high-Tc cuprates in qubits would dispose of the need to apply fine-tuned external fields to keep it in the operating point, due to the time-reversal symmetry breaking in systems with DD junctionsBlatter ; zagoskin . On the level of a few qubits this is not a major advantage, compared to the relative difficulty of fabrication and threat of extra decoherence from nodal quasiparticles and zero-energy states (ZES) in cuprates. Therefore the research was concentrated on conventional superconductors, where all of the aforementioned successes were achieved.

Nevertheless, the logical progression of research from single qubit to qubit-qubit coupling and further to qubit registers is bringing us to the point where intrinsic bistability of high-Tc qubits will become a major advantage. In the last few years, we saw the development of reliable fabrication of bistable submicron high-Tc structuresIlichev01 ; Ivanov ; Lindstrom . The danger of nodal quasiparticles and ZES seems to be overestimatedAminSmirnov ; Fominov . Finally, putting aside any technological motivation, it is tempting to probe high-Tc compounds for the same kind of macroscopic quantum coherence as was observed in low-Tc superconductors, just to see how and whether these phenomena depend on, e.g., symmetry of the ground state of the system.

The “silent qubit” (for the etimology, see below) can be viewed as a SQUID formed out of high-Tc film, with two DD-grain boundary junctions forming a mesoscopic island (Fig. 1a,c). The crystal lattice and $`d`$-wave order parameter orientations on the two sides of the grain boundary are chosen in such a way that each junction has doubly degenerate ground state Ilichev01 ; amin .

To simplify the analysis, we keep only two harmonics in the current-phase relations of the DD junctions: $`I_i=I_{ci}\mathrm{sin}\varphi _iI_{ci}^{}\mathrm{sin}2\varphi _i`$, where $`I_{ci}`$ is the Josephson current, and $`\varphi _i`$ is the phase difference across the $`i^{\mathrm{th}}`$ junction ($`i=1,2`$). This approximation successfully reproduces the observed behaviour of similar devices in classical regimeIlichev01 ; Lindstrom . The flux quantization condition binds the phases $`\varphi _1`$ and $`\varphi _2`$ to the total flux through the loop $`\mathrm{\Phi }`$ via $`\varphi \varphi _1+\varphi _2=2\pi \mathrm{\Phi }/\mathrm{\Phi }_0`$, where $`\mathrm{\Phi }_0=h/2e`$ is the flux quantum. Introducing the superconducting phase of the island $`\theta (\varphi _1\varphi _2)/2`$, the free energy of the qubit, threaded by an external magnetic flux $`\mathrm{\Phi }_x`$, in units of the Josephson energy $`E_J`$ ($`E_1`$, where $`E_i=I_{ci}\mathrm{\Phi }_0/2\pi `$), is given by

$$𝒰(\theta ,\varphi )=\frac{(\varphi \varphi _x)^2}{2\beta _L}_\varphi \left[\mathrm{cos}\theta \frac{\stackrel{~}{\alpha }_\varphi }{4}\mathrm{cos}(2\theta )\right]+\stackrel{~}{𝒰}(\theta ,\varphi ).$$
(1)
Here $`\alpha _i=2I_{ci}^{}/I_{ci},`$ $`_\varphi =(1+\eta )\mathrm{cos}(\varphi /2),\stackrel{~}{\alpha }_\varphi =\left[(\alpha _1+\eta \alpha _2)\mathrm{cos}\varphi \right]/\left[(1+\eta )\mathrm{cos}\left(\varphi /2\right)\right],`$ $`\eta =E_2/E_1,`$ and $`\varphi _x=2\pi \mathrm{\Phi }_x/\mathrm{\Phi }_0.`$ The dimensionless self-inductance of the loop $`\beta _L2\pi LI_{c1}/\mathrm{\Phi }_0`$ is considered negligible ($`\beta _L0`$); then $`\varphi \varphi _x`$ and the first term in (1) can be dropped.

The last term in (1):

$$\stackrel{~}{𝒰}(\theta ,\varphi )=\left[\eta 1+(\alpha _1\eta \alpha _2)\mathrm{cos}\frac{\varphi }{2}\mathrm{cos}\theta \right]\mathrm{sin}\frac{\varphi }{2}\mathrm{sin}\theta ,$$

is zero when the two junctions are identical (i.e. $`\eta =1`$ and $`\alpha _1=\alpha _2`$). In this case the potential energy minima are at $`\theta =\pm \mathrm{arccos}(1/\stackrel{~}{\alpha }_\varphi )`$ when $`\stackrel{~}{\alpha }_\varphi >1`$ and zero otherwise. The $`\pm `$ signs correspond to the states on the right and left sides of the potential well (i.e. $`|\pm `$), respectively. For $`\stackrel{~}{\alpha }_\varphi >1`$ the potential has two minima, which are degenerate at any external flux $`\varphi _x.`$ The current induced in the loop by the external flux does not depend on the state of the qubit, which justifies the moniker “silent”. The potential profile (1) is similar to the one of persistent current qubit mooij , but here we have only one independent phase (as $`\beta _L0`$), and the problem becomes one-dimensional. The barrier between the potential minima is flux-dependent: $`W=\mathrm{cos}(\varphi _x/2)(\stackrel{~}{\alpha }_\varphi +\stackrel{~}{\alpha }_\varphi ^12).`$

In the general case ($`\alpha _1\alpha _2,`$ $`\eta 1`$) the two minima are only degenerate when $`\varphi _x=0`$. Now a state-dependent persistent current flows in the loop even in zero external field; in units of $`I_{c1}`$, $`_0^\pm =\pm [\eta \sqrt{\stackrel{~}{\alpha }_0^21}/(1+\eta )\stackrel{~}{\alpha }_0^2](\alpha _2\alpha _1)`$, where $`\stackrel{~}{\alpha }_0\stackrel{~}{\alpha }_{\varphi =0}`$.

An intermediate regime ($`\alpha _1=\alpha _2`$, $`\eta 1`$) is most interesting. It takes place when junctions have similar current-phase dependences, but different critical currents (i.e. widths), and should be expected if the junctions are close enough to each other (Fig. 1). At $`\varphi _x=0`$, the equilibrium value of $`\theta `$ is the same for both junctions, and there is no spontaneous current: $`_0^\pm =0`$. At finite $`\varphi _x`$, the induced currents differ for the two states of the qubit, but the difference is of higher order in $`\varphi _x`$, keeping the qubit silent.

Expanding the free energy to the third order in $`\varphi _x`$, we find for the minima:

$$𝒰_{\mathrm{min}}^\pm =A_0+A_2\varphi _x^2\pm A_3\varphi _x^3$$
(2)
where $`A_i`$ are explicit, but cumbersome, functions of $`\stackrel{~}{\alpha }_0`$ and $`\eta `$. As expected, there is no first order dependence on $`\varphi _x`$ ($`A_1=0`$). The second order term in (2) does not depend on the state of the qubit. The first state-dependent term in the minimum energy of the system is $`O(\varphi _x^3)`$ (Fig. 2). Therefore small fluctuations of $`\varphi _x`$ do not affect the degeneracy of the states. The difference between the energy minima grows as the external flux is increased until the point at which the potential barrier vanishes altogether, and the minimum with higher energy disappears (the jumps in Fig. 2).

The current in the loop is found from $`_\varphi ^\pm =_{\varphi _x}𝒰_{\mathrm{min}}^\pm =2A_2\varphi _x\pm 3A_3\varphi _x^2`$. This current generates a state-dependent magnetic flux $`\delta \varphi \varphi \varphi _x=\beta _L_\varphi ^\pm `$. Note that the state-dependent contribution to the induced flux is $`O(\varphi _x^2).`$ For finite self-inductances $`\beta _L`$, $`\delta \varphi `$ has to be calculated self-consistently. Figure 3 shows the result of such a calculation. With the parameters chosen, an external flux close to $`0.2\mathrm{\Phi }_0`$ generates an additional flux ($`0.005\mathrm{\Phi }_0`$) through the loop, which is of the same order as the estimate based on the above expansion.

Tunneling between the potential minima occurs due to the uncertainty relation between the charge $`Q`$ of the island and its superconducting phase $`\theta `$. The tunneling matrix element is approximately given by ($`\mathrm{}=k_\mathrm{B}=1`$) $`\mathrm{\Delta }\omega _p(\varphi _x)e^{\sqrt{\zeta W(\varphi _x)E_J/E_c}}`$, where $`\zeta `$ is a constant of the order of 1. The coefficient $`\omega _p(\varphi _x)\sqrt{\omega _p^+\omega _p^{}}`$ is determined by the frequencies $`\omega _p^\pm `$ of small oscillations in the right and left potential minima, respectively. This dependence follows from the expression of $`\mathrm{\Delta }`$ as the matrix element between the lowest energy states in the two wells. In the case of $`\stackrel{~}{\alpha }_01`$, it is only valid qualitativelyTZA , but enough for the present analysis. Due to the symmetry of the potential profile when $`\alpha _1=\alpha _2`$, the linear dependence on $`\varphi _x`$ cancel and we are left with $`\omega _p(\varphi _x)=\sqrt{E_JE_C(\alpha _0\alpha _0^1)}(1\kappa \varphi _x^2),`$ where $`E_C=e^2/2C`$ is the charging energy, $`C`$ is the effective capacitance of the junctions, and $`\kappa `$ is a dimensionless coefficient of $`O(1)`$. Fluctuations of $`\varphi _x`$ influence $`W`$ and therefore $`\mathrm{\Delta }`$. Expanding the Josephson potential near the origin, we obtain the tunnelling barrier $`W=𝒰_{\mathrm{max}}𝒰_{\mathrm{min}}^\pm =(B_0A_0)+\left(B_2A_2\right)\varphi _x^2+O(\varphi _x^3)`$, where $`B_0=\left(\eta +1\right)\left(1\stackrel{~}{\alpha }_0/4\right),`$ $`B_2=\left[(\stackrel{~}{\alpha }_01)/4\right]\left[\eta /(\eta +1)\right]`$. Again there is no dependence on $`\varphi _x`$ in the lowest order.

Truncating the Hilbert space of the qubit to the two lowest energy states, one can write the effective Hamiltonian of the system as $`H=\frac{1}{2}\mathrm{\Delta }(\varphi _x)\sigma _x+\frac{1}{2}ϵ(\varphi _x)\sigma _z`$, where $`ϵ(\varphi _x)E_JA_3\varphi _x^3`$. All single qubit operations can be realized by applying controlled flux $`\varphi _x.`$ Note that the qubit only leaves the operating point when a finite external flux is applied. Unlike the earlier qubit designsmooij ; vion  this point is protected from external flux fluctuations, already on the classical level \[cf. Eq. (2)\].

Readout of the quantum state can in principle be achieved using a dc-SQUID to measure the magnetic flux generated by the qubit when it is biased by an external magnetic field. In this work the readout SQUID was fabricated onchip using a high-$`T_c`$ grain boundary SQUID (cf. Fig. 1c); alternatively it can be made of aluminum. Although coherent oscillations in an all-aluminum persistent current qubit have previously been observed using a similar readout designmooij , the measured decoherence time appeared to be rather short. Subsequent results of Il’ichev et al.Rabi  have demonstrated a longer decoherence time in a measurement scheme which couples the qubit to a high quality tank circuit. Here, we calculate the decoherence time of our silent qubit in the same setupGrajcar ; Rabi  (see Fig. 1b). The circuit is assumed to have a resonance frequency $`\omega _T`$ and the damping rate $`\gamma _T`$. The mutual inductance is $`M=k\sqrt{L_TL}`$, $`k1`$ being the coupling coefficient, and $`L`$ and $`L_T`$ are the inductance of the qubit and tank respectively.

Without limiting the generality of our approach, we can ascribe all the dephasing and dissipation in the qubit to its interaction with the fluctuating flux $`\varphi _x(t)`$, created by the tank current $`I_T(t)`$. These fluctuations are characterized by correlator $`K(t,t^{})=I_T(t),I_T(t^{})`$, spectral density $`K(\omega )=\omega (\gamma _T\omega _T^2/L_T)[\mathrm{coth}(\omega /2T)+1]/[(\omega ^2\omega _T^2)^2+\gamma _T^2\omega ^2],`$ and dispersion $`I_T^2=K(t,t)`$ $`=(\omega _T/2L_T)\mathrm{coth}(\omega _T/2T).`$ At $`\varphi _x=0`$, the qubit Hamiltonian becomes $`H=(\mathrm{\Delta }/2)\sigma _x\lambda _3I_T^3\sigma _z\lambda _2I_T^2\sigma _x`$, where $`\lambda _3`$ and $`\lambda _2`$ are coupling coefficients depending on the qubit parameters. In the Bloch-Redfield approximationsmirnov , we calculate the energy relaxation rate $`\mathrm{\Gamma }=30\lambda _3^2I_T^2^2(\gamma _T/\mathrm{\Delta }^3)(\omega _T^2/L_T)\mathrm{coth}(\mathrm{\Delta }/2T)`$, together with a dephasing time of the qubit $`\gamma ^1,`$ where $`\gamma =\mathrm{\Gamma }/2+\gamma _0`$, $`\gamma _0=(16\pi /3)\lambda _2^2(\gamma _T^2T^3/L_T^2\omega _T^4),`$ if $`T\omega _T,`$ and $`\gamma _0=\lambda _2^2Q_T(\omega _T/L_T^2)[\mathrm{coth}^2(\omega _T/2T)1]`$ at temperatures $`T\omega _T.`$

Using the experimental data of Ref. TZA,  ($`I_c=0.5`$ $`\mu `$A, $`I_c^{}=0.6`$ $`\mu `$A, $`C10`$ fF), we find $`E_C2`$ GHz, $`\omega _p/2\pi 40`$ GHz, and $`\mathrm{\Delta }/2\pi 1.6`$ GHz. For $`\beta _L0.01`$, which is the value used in Fig. 3, the inductance of the loop is of the order of $`10pH`$. To estimate the contributions of the cubic and quadratic terms to qubit dephasing and dissipation, we chose the following parameters: $`\eta =2,\alpha _1=\alpha _2=2.4`$, and $`E_J=1.66\times 10^{22}J`$. If the tank frequency $`\omega _T/2\pi =10`$ MHz, its quality factor $`Q_T=2000,`$ and the coupling coefficient $`k1/33`$, then contribution of quadratic flux fluctuations to dephasing rate is small, so that the dephasing time due to qubit coupling to the tank is $`\gamma _0^120`$ ms at temperatures of order $`10`$ mK, while the contribution of the cubic fluctuations to the dephasing and relaxation rates is totally negligible. It means that at the operating point the silent qubit is practically decoupled from the fluctuations caused by the controlling circuits. The dominant source of decoherence is from the nodal quasiparticles at the junction, considered in Ref. AminSmirnov, , which may reduce the decoherence time to about 1–100 ns.

In conclusion, a new type of flux qubit, based on specific properties of submicron YBCO grain boundary junctions, is proposed and fabricated. The symmetry of the device provides an operating point, which is stable and protected against the external field fluctuations.

We are grateful to A. Golubov, E. Il’ichev, A. Maassen van den Brink, Y. Nakamura, and A. Shnirman for stimulating discussions.