text:
# On the statistical significance of the GZK feature in the spectrum of ultra high energy cosmic rays

## 1 Introduction

The presence or lack of a feature in the spectrum of ultra-high energy cosmic rays (UHECRs) is key in determining the nature of these particles and their sources. Astrophysical proton sources distributed homogeneously in the universe produce a feature in the spectrum due to the production of pions off the cosmic microwave background. This feature, consisting of a rather sharp suppression of the flux, occurs at energies above $`7\times 10^{19}`$ eV, as a result of the threshold in the production of pions in the final state of a proton-photon inelastic interaction. This important prediction was made independently by Greisen and by Zatsepin and Kuzmin . The resulting spectral feature is now known as the GZK cutoff (or feature, as we prefer to call it). Alternative models for UHECR sources that involve new physical processes may evade the presence of this feature (see, e.g., ). Recent reviews on the origin and propagation of the ultra-high energy cosmic rays can be found in , while a recent review of the observations can be found in .

The detection of cosmic ray events with energy above $`E_{GZK}7\times 10^{19}`$ eV does not necessarily imply that the GZK feature is not present: what characterizes the presence of the GZK feature is the relative number of events above and below $`E_{GZK}`$ when both sides of the spectrum can be accurately determined. The steep injection spectra required to fit the observations below $`E_{GZK}`$ imply that only a handfull of events above $`10^{20}`$ eV can be detected during the operation time of experiments such as AGASA and HiRes. This makes the identification of the GZK feature by these experiments extremely difficult. The problem is exacerbated by the fluctuations due to the discreteness of the process of photo-pion production, as will be discussed below. These uncertainties need to be considered when attempting a determination of the best fit injection spectrum of the particles, and in order to quantify the statistical significance of the presence or absence of the GZK feature in the observed spectrum.

The discrepancy between the results of the two largest experiments has generated much debate. AGASA  reports a higher number of events above $`E_{GZK}`$ than expected while HiRes reports a flux consistent with the GZK feature. Here we investigate in detail the statistical significance of this discrepancy as well as the significance of the presence or absence of the GZK feature in the data. We find that neither experiment has the necesseray statistics to establish if the spectrum of UHECRs has a GZK feature. In addition, a systematic error in the energy determination of the two experiments seems to be required in order to make the two sets of observations compatible in the low energy range, $`10^{18.5}10^{19.6}`$ eV, where enough events have been detected to make the measurements reliable. The combined systematic errors in the energy determination is likely to be $`30\%`$. If we decrease the AGASA energies by $`15\%`$ while increasing HiRes energies also by $`15\%`$, the two experiments predict compatible fluxes at energies below $`E_{GZK}`$ and at energies above $`E_{GZK}`$ the fluxes are within $`2\sigma `$ of each other. In this case, the best fit injection spectrum has a spectral index of $`2.6`$ but a determination of the GZK feature has very low significance. The detection or non-detection of the GZK feature in the cosmic ray spectrum remains open to investigation by future generation experiments, such as the Pierre Auger project  and the EUSO  and OWL experiments.

This paper is planned as follows: in §2 we describe our simulations. In §3 we illustrate the present observational situation, limiting ourselves to AGASA and HiResI, and compare the data to the predictions of our simulations. We conclude in §4.

## 2 UHECR spectrum simulations

We assume that ultra-high energy cosmic rays are protons injected with a power-law spectrum in extragalactic sources. The injection spectrum is taken to be of the form

$$F(E)dE=\alpha E^\gamma \mathrm{exp}(E/E_{\mathrm{max}})dE$$
(1)
where $`\gamma `$ is the spectral index, $`\alpha `$ is a normalization constant, and $`E_{\mathrm{max}}`$ is the maximum energy at the source. Here we fix $`E_{\mathrm{max}}=10^{21.5}`$ eV, large enough not to affect the statistics at much lower energies. As shown in , the induced spectrum of a uniform distribution of sources in space is almost indistinguishable from a distribution with the observed large scale structure in the galaxy distribution. Based on this result, we assume a spatially uniform distribution of sources and do not take into account luminosity evolution in order to avoid the introduction of additional parameters.

We simulate the propagation of protons from source to observer by including the photo-pion production, pair production, and adiabatic energy losses due to the expansion of the universe. In each step of the simulation, we calculate the pair production losses using the continuous energy loss approximation given the small inelasticity in pair production ($`2m_\mathrm{e}/m_\mathrm{p}10^3`$). For the rate of energy loss due to pair production at redshift $`z=0`$, $`\beta _{\mathrm{pp}}(E,z=0)`$, we use the results from . At a given reshift $`z>0`$,

$$\beta _{\mathrm{pp}}(E,z)=(1+z)^3\beta _{\mathrm{pp}}((1+z)E,z=0).$$
(2)
Similarly, the rate of adiabatic energy losses due to redshift is calculated in each step using

$$\beta _{\mathrm{rsh}}(E,z)=H_0\left[\mathrm{\Omega }_M(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}$$
(3)
with $`H_0=75\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

The photo-pion production is simulated in a way similar to that described in ref. . In each step, we first calculate the average number of photons able to interact via photo-pion production through the expression:

$$N_{\mathrm{ph}}(E,\mathrm{\Delta }s)=\frac{\mathrm{\Delta }s}{l(E,z)},$$
(4)
where $`l(E,z)`$ is the interaction length for photo-pion production of a proton with energy $`E`$ at redshift $`z`$ and $`\mathrm{\Delta }s`$ is a step size, chosen to be much smaller than the interaction length (typically we choose $`\mathrm{\Delta }s=100\mathrm{kpc}/(1+z)^3`$).

In Fig. 1 we plot the interaction length for photopion production used in  (solid thin line), and in  (triangles). The dashed line is the result of our calculations (see below), which is in perfect agreement with the results of . The apparent discrepancy at energies below $`10^{19.5}`$ eV with the prediction of Ref.  is only due to the fact that we consider only microwave photons as background, while in  the infrared background was also considered. For our purposes, this difference is irrelevant as can be seen from the loss lengths plotted in Fig. 1. The rightmost thick solid line is the loss length for photopion production , while the other thick solid line is the loss length for pair production. In the present calculations, we do not use the loss length of photopion production which is related to the interaction length through an angle averaged inelasticity. Unlike what was done in , in the current simulations we evaluate the inelasticity for each single proton-photon scattering using the kinematics, rather than adopting an angle averaged value.

We calculate the interaction length, $`l(E)`$, as:

$$l(E)^1=d\epsilon n(\epsilon )_1^{+1}d\mu \frac{1\mu \beta }{2}\sigma (s)$$
(5)
where $`n(\epsilon )`$ is the number density of the CMB photons per unit energy at energy $`\epsilon `$, $`\beta `$ is the velocity of the proton, $`\mu `$ is the cosine of the interaction angle, and $`\sigma (s)`$ is the total cross section for photo-pion production for the squared center of mass energy

$$s=m^2+2\epsilon E\left(1\mu \beta \right).$$
(6)
For the calculation of the interaction length we adopt the $`\mathrm{p}+\gamma \mathrm{hadrons}`$ cross section given in , considering only the photons of the cosmic microwave background as targets. The calculated interaction length (see dashed line in Fig. 1) is in good agreement with the interaction length calculated in  and in .

Once the interaction length is known, we then sample a Poisson distribution with mean $`N_{\mathrm{ph}}(E,\mathrm{\Delta }s)`$, to determine the actual number of photons encountered during the step $`\mathrm{\Delta }s`$. When a photo-pion interaction occurs, the energy $`ϵ`$ of the photon is extracted from the Planck distribution, $`n_{ph}(ϵ,T(z))`$, with temperature $`T(z)=T_0(1+z)`$, where $`T_0=2.728`$ K is the temperature of the cosmic microwave background at present. Since the microwave photons are isotropically distributed, the interaction angle, $`\theta `$, between the proton and the photon is sampled randomly from a distribution which is flat in $`\mu =\mathrm{cos}\theta `$. Clearly only the values of $`ϵ`$ and $`\theta `$ that generate a center of mass energy above the threshold for pion production are considered. The energy of the proton in the final state is calculated at each interaction from kinematics. The simulation is carried out until the statistics of events detected above some energy reproduces the experimental numbers. By normalizing the simulated flux by the number of events above an energy where experiments have high statistics, we can then ask what are the fluctuations in numbers of events above a higher energy where experimental results are sparse. The fits are therefore most sensitive to the energy regions below $`E_{GZK}`$ and give a good estimate of the uncertainties in the present experiments for energies above $`E_{GZK}`$. In this way we have a direct handle on the fluctuations that can be expected in the observed flux due to the stochastic nature of photo-pion production and to cosmic variance.

The simulation proceeds in the following way: a source distance is generated at random from a uniform distribution in a universe with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`\mathrm{\Omega }_m=0.3`$. In a Euclidean universe, the flux from a source would scale as $`r^2`$ where $`r`$ is the distance between the source and the observer. On the other hand, the number of sources between $`r`$ and $`r+dr`$ would scale as $`r^2`$, so that the probability that a given event has been generated by a source at distance $`r`$ is independent of $`r`$: sources at different distances have the same probability of generating any given event. In a flat universe with a cosmological constant, this is still true provided the distance $`r`$ is taken to be

$$r=c_{t_g}^{t_0}\frac{dt}{R(t)},$$
(7)
where $`t_g`$ is the age of the universe when the event was generated, $`t_0`$ is the present age of the universe, and $`R(t)`$ is the scale factor of the universe.

Once a source distance has been selected at random, a particle energy is also assigned from a distribution that reflects the injection spectrum, chosen as in Eq. (1). This particle is then propagated to the observer and its energy recorded. This procedure is repeated until the number of events above a threshold energy, $`E_{th}`$ is reproduced. With this procedure we can assess the significance of results from present experiments with limited statistics of events. There is an additional complication in that the aperture of the experiment usually depends on energy. This is taken into account by allowing the event to be detected or not detected depending upon the function $`H(E)`$ that describes the energy dependence of the aperture.

We only study the spectrum above $`10^{18.5}`$ eV where the flux is supposed to be dominated by extragalactic sources. For this energy range, we focus on the experiments that have the best statistics: AGASA and HiResI. For the AGASA experiment the exposure is basically flat above $`10^{19}`$ eV, while for HiRes the exposure is as plotted in Fig. (2) . For AGASA data, the simulation is stopped when the number of events above $`E_{th}=10^{19}`$ eV equals 866. For HiRes this number is 300. Note that while for AGASA the number of detected events actually corresponds to the generated events, for HiRes the number of detected events requires a correction due to the energy dependence of the aperture $`H(E)`$. This correction allows one to reconstruct the observed spectrum. The statistical error in the energy determination is accounted for in our simulation by generating a detection energy chosen at random from a Gaussian distribution centered at the arrival energy $`E`$ and with width $`\mathrm{\Delta }E/E=30\%`$ for both experiments.

Our simulations reproduce well the predictions of analytical calculations, in particular at the energies where energy losses may be approximated as continuous. In Fig. 3, we compare the results of our simulation with analytical calculations carried out as in ref. . The curves plotted in the figure are the so called modification factors, defined in  for three different values of the source redshift ($`z=0.002`$, $`z=0.02`$ and $`z=0.2`$ from top to bottom). The differential injection spectrum is taken to be a power law $`E^{2.1}`$. The points with error bars are the results of our simulations with $`2\times 10^6`$ particles produced by sources at the redshifts given above. The agreement between the simulated and analytical calculations at low energies is clear. At the energies where photo-pion production becomes important, simulations predict a slightly larger flux than analytical calculations. This is a well known result, and is due to the discreteness of photo-pion energy losses, that allow some particles to reach the detector without appreciable losses.

In Fig. 4 we compare the results of our simulations (points with error bars) with analytical calculations of the diffuse flux of UHECRs (continuous line) expected if sources with no luminosity evolution inject UHECRs with a spectrum $`E^{2.7}`$ . The agreement is evident.

## 3 AGASA versus HiResI

The two largest experiments that measured the flux of UHECRs report apparently conflicting results. The data of AGASA and HiResI on the flux of UHECRs multiplied as usual by the third power of the energy are plotted in Fig. 5 (the squares are the HiResI data while the circles are the AGASA data).

The figure shows that HiResI data are systematically below AGASA data and that HiResI sees a suppression at $`10^{20}`$ eV that resembles the GZK feature while AGASA does not.

We apply our simulations here to the statistics of events of both AGASA and HiResI in order to understand whether the discrepancy is statistically significant and whether the GZK feature has indeed been detected in the cosmic ray spectrum. The number of events above $`10^{19}`$ eV, $`10^{19.6}`$ eV and $`10^{20}`$ eV for AGASA and HiResI are shown in Table 1. For reasons that will be clear below, we also show in the same table the number of events above the same energy thresholds, in the case that AGASA has a systematic error that overestimates the energy by $`15\%`$ while HiResI systematically underestimates the energy by $`15\%`$.

In order to understand the difference, if any, between AGASA and HiRes data we first determine the injection spectrum required to best fit the observations. In order to do this, we run 400 realizations of the AGASA and HiRes event statistics, as reported in the column labelled data in Table 1.

The injection spectrum is taken to be a power law with index $`\gamma `$ between 2.3 and 2.9 with steps of 0.1. For each injection spectrum we calculated the $`\chi ^2`$ indicator (averaged over 400 realizations for each injection spectrum). The errors used for the evaluation of the $`\chi ^2`$ are due to the square roots of the number of observed events. As reported in Table 2, the best fit injection spectrum can change appreciably depending on the minimum energy above which the fit is calculated. In these tables we give $`\chi _e^2(N)`$, where $`N`$ is the number of degrees of freedom, while the subscript, $`e`$, is the logarithm of $`E_{th}`$ (in eV), the energy above which the fit has been calculated. The numbers in bold-face correspond to the best fit injection spectrum. These fits are dominated by the low energy data rather than by the poorer statistics at the higher energies.

If the data at energies above $`10^{18.5}`$ eV are taken into account for both experiments, the best fit spectra are $`E^{2.8}`$ for AGASA and $`E^{2.6}`$ for HiRes. If the data at energies above $`10^{19}`$ eV are used for the fit, the best fit injection spectrum is $`E^{2.6}`$ for AGASA and between $`E^{2.7}`$ and $`E^{2.8}`$ for HiRes. If the fit is carried out on the highest energy data ($`E>10^{19.6}`$ eV), AGASA prefers an injection spectrum between $`E^{2.5}`$ and $`E^{2.6}`$, while $`E^{2.8}`$ or $`E^{2.9}`$ fit better the HiRes data in the same energy region. Note that the two sets of data uncorrected for any possible systematic errors require different injection spectra that change with $`E_{th}`$.

In order to quantify the significance of the detection or lack of the GZK flux suppression, we report in Tables 3 and 4 the mean number of events above the indicated energy threshold, $`N(E>E_{th})`$, for different injection specta. In parenthesis , we show the discrepancy between the observations as in Table 1 and the simulations in number of standard deviations, $`\sigma `$, where $`\sigma ^2=N(E>E_{th})^2N(E>E_{th})^2`$. From Tables 3 and 4 we can see that while HiResI is consistent with the existence of the GZK feature in the spectrum of UHECRs, AGASA detects an increase in flux, but only at the $`2.5\sigma `$ level for the best fit injection spectra.

A more graphical representation of the uncertainties involved are displayed in Figs. 6 for AGASA and 7 for HiRes, for two choices of the injection spectrum. These plots show clearly the low level of significance that the detections above $`E_{GZK}`$ have with low statistics. The large error bars that are generated by our simulations at the high energy end of the spectrum are mainly due to the stochastic nature of the process of photo-pion production: in some realizations some energy bins are populated by a few events, while in others the few particles in the same energy bin do not produce a pion and get to the observer unaffected. The large fluctuations are unavoidable with the extremely small statistics available with present experiments. On the other hand, the error bars at lower energies are minuscule, so that the two data sets (AGASA and HiResI) cannot be considered to be two different realizations of the same phenomenon. Instead, systematic errors in at least one if not both experiments are needed to explain the discrepancies at lower energies.

Taking into account the (theoretical) error bars in the analysis makes the significance of the presence or absence of the GZK feature much weaker than the consistency checks shown in Tables 3 and 4. In order to estimate this effect, we proceed in the following way: we calculate the expected number of events above some threshold with its corresponding standard deviation ($`\sigma _{sim}`$), as determined by the fluctuations in the flux simulation. The observed number of events above the same threshold is also known with the error bar $`\sigma _{obs}`$. The discrepancy between the two is now calculated using the error $`\sigma _{tot}=(\sigma _{sim}^2+\sigma _{obs}^2)^{1/2}`$. Our results are summarized in Tables 3 for AGASA and 4 for HiRes. The numbers with error bars are the simulated expectations, while the discrepancy between simulations and observations, calculated as described above is in square brackets, in units of $`\sigma _{tot}`$. It becomes clear that the effective discrepancy between predictions and the AGASA data is at the level of $`2.12.5\sigma `$. Therefore a definite answer to the question of whether the GZK feature is there or not awaits a significant improvement in statistics at high energies.

As seen in Fig. 5, the difference between the AGASA and HiResI spectra is not only in the presence or absence of the GZK feature: the two spectra, when multiplied by $`E^3`$, are systematically shifted by about a factor of two. This shift suggests that there may be a systematic error either in the energy or the flux determination of at least one of the two experiments. Possible systematic effects have been discussed in  for the AGASA collaboration and in  for HiResI. At the lower end of the energy range the offset may be due to the notoriously difficult determination of the AGASA aperture at threshold while the discrepancies at the highest energies is unclear at the moment. In any case, a systematic error of $`15\%`$ in the energy determination is well within the limits that are allowed by the analysis of systematic errors carried out by both collaborations.

In order to illustrate the difficulty in determining the existence of the GZK feature in the observed data in the presence of systematic errors, we split the energy gap by assuming that the energies as determined by the AGASA collaboration are overestimated by $`15\%`$, while the HiRes energies are underestimated by the same factor. The number of events above an energy threshold is again reported in Table 1, in the column labelled $`15\%`$. In this case the total number of events above $`10^{19}`$ eV is reduced for AGASA from 866 to 651, while for HiResI it is enhanced from 300 to 367. We ran our simulations with these new numbers of events and repeat the statistical analysis described above. The values of $`\chi ^2`$ obtained in this case are reported in Table 5.

For AGASA, the best fit injection spectrum is now between $`E^{2.5}`$ and $`E^{2.6}`$ above $`10^{19}`$ eV and above $`10^{19.6}`$ eV (the $`\chi ^2`$ values are identical). For the HiRes data, the best fit injection spectrum is $`E^{2.6}`$ for the whole set of data, independent of the threshold. It is interesting to note that the best fit injection spectrum appears much more stable after the correction of the $`15\%`$ systematics has been carried out. Moreover, the best fit injection spectra as derived for each experiment independently coincides for the corrected data unlike the uncorrected case. This suggests that combined systematic errors in the energy determination at the $``$ 30% level may in fact be present.

The expected numbers of events with energy above $`10^{19.6}`$ eV and above $`10^{20}`$ eV with the deviation from the data are reported in Tables 6 and 7: while HiResI data remain in agreement with the prediction of a GZK feature, the AGASA data seem to depart from such prediction but only at the level of $`1.8\sigma `$. The systematics in the energy determination further decreased the significance of the GZK feature, such that the AGASA and HiResI data are in fact only less than $`2\sigma `$ away from each other.

We can use the same procedure illustrated above to evaluate the effect of the error bars in the simulated data compared to the data corrected by 15%. The results are reported in square brackets in Tables 6 (for AGASA) and 7 (for HiRes), showing that the effective discrepancy between expectations (with uncertainties due to discrete energy losses and cosmic variance) and AGASA data is even smaller, only at the $`1.5\sigma `$ level. In Fig. 8, we plot the simulated spectra for injection spectrum $`E^{2.6}`$ and compare them to observations of AGASA (upper plot) and HiRes (lower plot).

## 4 Conclusions

We considered the statistical significance of the UHECR spectra measured by the two largest experiments in operation, AGASA and HiRes. The spectrum released by the HiResI collaboration seems to suggest the presence of a GZK feature. This has generated claims that the GZK cutoff has been detected, reinforced by data from older experiments . However, no evidence for such a feature has been found in the AGASA experiment. We compared the data with theoretical predictions for the propagation of UHECRs on cosmological distances with the help of numerical simulations. We find that the very low statistics of the presently available data hinders any statistically significant claim for either detection or the lack of the GZK feature.

A comparison of the spectra obtained from AGASA and HiResI shows a systematic shift of the two data sets, which may be interpreted as a systematic error in the relative energy determination of about $`30\%`$. If no correction for this systematic shift is carried out, the AGASA data are best fit by an injection spectrum $`E^\gamma `$ with $`\gamma =2.6`$ for energies above $`10^{19}`$ eV. The fit steepens to $`\gamma =2.8`$ when considering events down to $`10^{18.5}`$ eV. For HiResI, the best fit is between $`\gamma =2.7`$ and $`\gamma =2.8`$ for events with energy above $`10^{19}`$ eV and $`\gamma =2.6`$ for events above $`10^{18.5}`$ eV. With these best fits to the injection spectrum the AGASA data depart from the prediction of a GZK feature by $`2.6\sigma `$ for $`\gamma =2.6`$. The HiRes data are fully compatible with the prediction of a GZK feature in the cosmic ray spectrum. The fit to the data with energy above $`10^{19}`$ eV is probably less susceptible to contamination by a possible galactic flux. In this case the AGASA departure from the expected GZK feature is, as stressed above, at the level of about $`2.6\sigma `$. Taking into account the uncertainties derived from the simulations, attributed to the discreteness of the photo-pion production and to cosmic variance, this discrepancy becomes even less significant ($`2.3\sigma `$). It is clear that, if confirmed by future experiments with much larger statistics, the increase in flux relative to the GZK prediction hinted by AGASA would be of great interest. This may signal the presence of a new component at the highest energies that compensates for the expected suppression due to photo-pion production, or the effect of new physics in particle interactions (for instance the violation of Lorentz invariance or new neutrino interactions).

Identifying the cause of the systematic energy and/or flux shift between the AGASA and the HiRes spectra is crucial for understanding the nature of UHECRs. This discrepancy has stimulated a number of efforts to search for the source of these systematic errors including the construction of hybrid detectors, such as Auger, that utilize both ground arrays and fluorescence detectors. A possible overestimate of the AGASA energies by $`15\%`$ and a corresponding underestimate of the HiRes energies by the same amount would in fact bring the two data sets in agreement in the region of energies below $`10^{20}`$ eV. In this case both experiments are consistent with a GZK feature with large error bars. The AGASA excess is at the level of $`1.7\sigma `$ ($`1.4\sigma `$ if the observational uncertainties are also taken into account). Interestingly enough, the correction by $`15\%`$ in the error determination implies that the best fit injection spectrum becomes basically the same for both experiments ($`E^{2.6}`$).

With the low statistical significance of either the excess flux seen by AGASA or the discrepancies between AGASA and HiResI, it is inaccurate to claim either the detection of the GZK feature or the extension of the UHECR spectrum beyond $`E_{GZK}`$ at this point in time. A new generation of experiments is needed to finally give a clear answer to this question. In Fig. 9 we report the simulated spectra that should be achieved in 3 years of operation of Auger (upper panel) and EUSO (lower panel). The error bars reflect the fluctuations expected in these high statistics experiments for the case of injection spectrum $`E^{2.6}`$. (Note that the energy threshold for detection by EUSO is not yet clear.) It is clear that the energy region where statistical fluctuations dominate the spectrum is moved to $`10^{20.6}`$ eV for Auger, allowing a clear identification of the GZK feature. The fluctuations dominated region stands beyond $`10^{21}`$ eV for EUSO.

## Acknowledgments

We thank Douglas Bergman for providing the number of events in different bins for the HiRes results. We also thank Venya Berezinsky, James Cronin, Masahiro Teshima, Mario Vietri and Alan Watson for a number of helpful comments. This work was supported in part by the NSF through grant AST-0071235 and DOE grant DE-FG0291-ER40606 at the University of Chicago.


format: text

text:
# Josephson nanocircuit in the presence of linear quantum noise

## Abstract

We derive the effective hamiltonian for a charge-Josephson qubit in a circuit with no use of phenomenological arguments, showing how energy renormalizations induced by the environment appear with no need of phenomenological counterterms. This analysis may be important for multiqubit systems and geometric quantum computation.

Keywords: Josephson effect; quantum computation; decoherence

Josephson junction based nanocircuits have been proposed for the implementation of quantum gates and quantum coherent behavior has been recently observed in charge-based devices (charging energy $`E_C`$ larger than the Josephson energy $`E_J`$). Decoherence in these devices is due to several sources, as fluctuations of the circuit, backaction of the measuring apparatus or noise due to background charges in the substrate. Fluctuations of the circuit are modeled by coupling the system to an environment of harmonic oscillators which mimics the external impedences (see Fig.1). An effective hamiltonian $`H_{eff}`$ is considered, which represents a spin-boson model (or a multistate version), the central variable being the charge $`Q`$ in the island and the environment being fixed in a phenomenological way. The environment produces decoherence and energy shifts, which may in principle be large. In dissipative quantum mechanics shifts are usually treated either by introducing counterterms or by writing $`H_{eff}`$ in terms of renormalized quantities. We present here a model for the electromagnetic environment and we derive a multistate $`H_{eff}`$ using no phenomenological argument. This has two motivations. First in principle bare circuit parameters are well defined and tunable, so we want to know precisely how this reflects on $`H`$. Second the role of induced shifts, which is minor in the devices of Refs., may be crucial in various situations (e.g. geometric quantum computation, dynamics of registers and error correction devices).

We consider the Cooper pair box of Fig.1. The external impedance is modeled by a suitable $`LC`$ transmission line and the Lagrangian of the system is

$$L=\underset{i=1,2}{}\frac{C_i\dot{\varphi }_i^2}{2}V_J(\varphi _1)+\underset{\alpha }{}\left[\frac{C_\alpha \dot{\varphi }_\alpha ^2}{2}\frac{\varphi _\alpha ^2}{2L_\alpha }\right]$$

where $`\dot{\varphi }`$ are voltage drops and the Josephson energy is $`V_J(\varphi _1)=E_J\mathrm{cos}(2e\varphi _1/\mathrm{})`$. The environment is fully specified by the elements $`C_\alpha `$ and $`L_\alpha `$. The circuit is introduced by the constraint $`\dot{\varphi }_1+\dot{\varphi }_2+_\alpha \dot{\varphi }_\alpha =V_x`$, which allows to eliminate one variable and to write

$`L`$ $`=`$ $`{\displaystyle \frac{C_\mathrm{\Sigma }\dot{\eta }^2}{2}}V_J\left(2e(\kappa _2\mathrm{\Phi }\eta )/\mathrm{}\right)+L_b`$ (1)
$`L_b`$ $`=`$ $`{\displaystyle \frac{1}{2}}C_e\dot{\mathrm{\Phi }}^2+{\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{C_\alpha }{2}}\dot{\varphi }_\alpha ^2{\displaystyle \frac{1}{2L_\alpha }}\varphi _\alpha ^2\right]`$ (2)
where $`\eta =\kappa _2\varphi _2\kappa _1\varphi _1`$, $`C_\mathrm{\Sigma }=_iC_i`$, $`C_e=C_1C_2/C_\mathrm{\Sigma }`$, $`\kappa _{1,2}=C_{1,2}/C_\mathrm{\Sigma }`$, and $`\dot{\mathrm{\Phi }}=V_x_\alpha \dot{\varphi }_\alpha `$.

We next diagonalize $`L_b`$ and obtain the form

$$L_b=\underset{\alpha }{}\left[\frac{m_\alpha }{2}\dot{x}_\alpha ^2\frac{m_\alpha \omega _\alpha ^2}{2}x_\alpha ^2\right]C_eV_x\underset{\alpha }{}d_\alpha \dot{x}_\alpha .$$

In order to determine the parameters, notice that $`L_b`$ describes a series $`C_eZ`$ circuit, with $`_\beta d_\beta \dot{x}_\beta =_\alpha \dot{\varphi }_\alpha =V_Z`$. By comparing the linear response with the known classical dynamics of $`V_Z`$ (as explained in Refs.) we determine the spectral density ($`\omega >0`$)

$$J^{}(\omega )=\underset{\alpha }{}\frac{\pi d_\alpha ^2\delta (\omega \omega _\alpha )}{2m_\alpha \omega _\alpha }=\mathrm{Re}\left[\frac{Z(\omega )/\omega }{1+i\omega Z(\omega )C_e}\right]$$

We stress that $`L_b`$ is quadratic therefore the procedure above is an exact way to perform the diagonalization, which uses classical circuit theory as a tool. It’s validity for quantum harmonic oscillators is guaranteed by the Ehrenfest theorem.

To get rid of the $`\mathrm{\Phi }`$ in the potential term in Eq.(1), we perform a (canonical) transformation $`\chi =\eta \kappa _2\mathrm{\Phi }`$ and obtain the total Lagrangian $`L=L_a+L_b`$ where

$$L_a=\frac{C_\mathrm{\Sigma }}{2}\left(\dot{\chi }+\kappa _2V_x\kappa _2\underset{\alpha }{}d_\alpha \dot{x}_\alpha \right)^2+V_J(2e\chi /\mathrm{})$$

One can verify that the variable canonically conjugated to $`\chi `$ is the charge $`Q`$ in the island. The Hamiltonian corresponding to $`L`$ reads

$`H`$ $`=`$ $`{\displaystyle \frac{Q^2}{2C_1}}+Q\kappa _2{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{d_\alpha }{m_\alpha }}p_\alpha E_J\mathrm{cos}\left({\displaystyle \frac{2e}{\mathrm{}}}\chi \right)`$ (4)
$`+{\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{p_\alpha ^2}{2m_\alpha }}+{\displaystyle \frac{1}{2}}\omega _\alpha ^2x_\alpha ^2\right]+C_eV_x{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{d_\alpha }{m_\alpha }}p_\alpha `$
where $`p_\alpha `$ are conjugated to $`x_\alpha `$ and we used the relation $`_\alpha d_\alpha ^2/m_\alpha =1/C_e`$. The system variable $`Q`$ is coupled with the momenta $`p_\alpha `$ of the environment.

A further canonical transformation of the environment ($`\stackrel{~}{x}_\alpha =p_\alpha /(m_\alpha \omega _\alpha )`$, $`\stackrel{~}{p}_\alpha =m_\alpha \omega _\alpha x_\alpha `$) yields a hamiltonian where $`Q`$ is coupled with the coordinates

$`H_{eff}`$ $`=`$ $`{\displaystyle \frac{Q^2}{2C_1}}+V_J\left({\displaystyle \frac{2e}{\mathrm{}}}\chi \right)+\kappa _2Q{\displaystyle \underset{\alpha }{}}c_\alpha \stackrel{~}{x}_\alpha `$ (5)
$`+`$ $`{\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{\stackrel{~}{p}_\alpha ^2}{2m_\alpha }}+{\displaystyle \frac{m_\alpha \omega _\alpha ^2}{2}}\stackrel{~}{x}_\alpha ^2\right]+C_eV_x{\displaystyle \underset{\alpha }{}}c_\alpha \stackrel{~}{x}_\alpha `$ (6)
where $`c_\alpha =d_\alpha \omega _\alpha `$ and we introduce the spectral density $`J(\omega )=\pi _\alpha \delta (\omega \omega _\alpha )c_\alpha ^2/(2m_\alpha \omega _\alpha )=\omega ^2J^{}(\omega )`$.

If we isolate the dc bias $`V_x(t)=V_x+\delta V_x(t)`$ and redefine $`\stackrel{~}{x}_\alpha +c_\alpha C_eV_x/(m_\alpha \omega _\alpha ^2)x_\alpha `$ we finally obtain

$`H_{eff}`$ $`=`$ $`{\displaystyle \frac{Q^2}{2C_1}}\kappa _2V_xQ+V_J\left({\displaystyle \frac{2e}{\mathrm{}}}\chi \right)+Q\kappa _2{\displaystyle \underset{\alpha }{}}c_\alpha x_\alpha `$ (7)
$`+`$ $`{\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{p_\alpha ^2}{2m_\alpha }}+{\displaystyle \frac{m_\alpha \omega _\alpha ^2}{2}}x_\alpha ^2\right]+C_e\delta V_x(t){\displaystyle \underset{\alpha }{}}c_\alpha x_\alpha `$ (8)
Notice that the capacitance $`C_1`$ (and not $`C_\mathrm{\Sigma }`$) enters the $`Q^2`$ term and if we put $`c_\alpha =0`$ we do not obtain the Cooper pair box hamiltonian. This is correct because the environment represents global fluctuations of the circuit, not only of $`Z`$. Notice that a static $`Q`$ shifts the equilibrium positions of the oscillators and also produces a $`Q`$-dependent shift of the zeroes of their energies. We can single out the corresponding term and focus on it. In the situations described in Refs. this term should be canceled by introducing a counterterm, because we have informations only about renormalized effective parameters of the model. In our case we have information about the bare parameters therefore there is no reason to cancel the $`Q`$-dependent shift of the zero of the oscillator energies. It can be reabsorbed in the charging energy if we write the oscillator hamiltonian using the shifted values $`x_\alpha +c_\alpha \kappa _2Q/(m_\alpha \omega ^2)`$, which produces the extra term $`\kappa _2^2Q^2_\alpha c_\alpha ^2/(2m_\alpha \omega _\alpha ^2)=Q^2/(2C_e)`$ and

$`H_{eff}`$ $`=`$ $`{\displaystyle \frac{Q^2}{2C_\mathrm{\Sigma }}}\kappa _2V_xQ+V_J\left({\displaystyle \frac{2e\chi }{\mathrm{}}}\right)+C_e\delta V_x{\displaystyle \underset{\alpha }{}}c_\alpha x_\alpha `$ (9)
$`+`$ $`{\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{p_\alpha ^2}{2m_\alpha }}+{\displaystyle \frac{m_\alpha \omega _\alpha ^2}{2}}\left(x_\alpha +{\displaystyle \frac{c_\alpha \kappa _2Q}{m_\alpha \omega ^2}}\right)^2\right]`$ (10)
which reduces to the non dissipative form by letting $`c_\alpha =0`$. This is a convenient starting point for a weak coupling analysis also because a static shift in the equilibrium points of the oscillators has no effect even if it is large.


text:
# First-Order Phase Transition in a Reaction-Diffusion Model With Open Boundary: The Yang-Lee Theory Approach

## 1 Introduction

It is known that one-dimensional driven diffusive systems can exhibit many interesting collective phenomena such as jamming and spontaneous symmetry breaking . Applications of such models are divers and include the kinetics of biopolymerization  and traffic flow , but to give just a few examples.
A powerful framework for studying equilibrium phase transitions in classical statistical mechanics is the Yang-Lee theory of phase transition . Here we shall briefly explain the key points of this theory. Consider a simple model consisting of $`L`$ spins in thermal contact with a heat reservoir. The hallmark of a phase transition in this system is the appearance of nonanalyticities in its free energy which can be defined in terms of the partition function of the system $`Z_L`$ as $`f=\frac{1}{L}\mathrm{ln}Z_L`$. In the thermal equilibrium the partition function (which is a function of temperature $`T`$ in our example) determines the statistical properties of the system. However, since the partition function has no real and positive zeros in the complex-$`T`$ plane, there is no scope for a phase transition in a finite system. As we increase the number of spins $`L`$ to infinity, the partition function zeros might accumulate towards a point $`T_0`$ on the real axis so there is the possibility of a phase transition at this point. The density of zeros near the accumulation point determines the order of the phase transition. At a first-order phase transition there is a nonzero density of zeros at the critical point whereas the density of zeros decays to zero at a continuous transition. Although we explained the theory of partition function zeros with reference to a system described by a canonical partition function, we should note that the theory still holds when one is working in grand canonical ensembles. In this case the fugacity (or chemical potential) plays similar role as the temperature $`T`$. More generally one can look for the zeros of the partition function in the complex plane of any of its intensive variables .
The lack of such a theory in the non-equilibrium statistical mechanics context has recently encouraged people to examine whether a generalized version of it can still work for out-of-equilibrium systems. Quite similar to the equilibrium statistical mechanics one can define a generalized partition function as sum over the steady state weights. These weights can be calculated using different methods such as matrix product formalism in which the stationary probability distribution function is written in terms of the products of non-commuting operators . Then we look at the zeros of the partition function in the thermodynamic limit in the complex plane of any intensive variable of the model. The zeros might divide the complex plane into several regions where the free energy of the system is analytic. Each region corresponds to a different phase and as we explained above, the order of transition depends on how the zeros accumulate towards the positive real axis. Several out of equilibrium models have been proposed and studied using the generalized version of the Yang-Lee theory -, and it seems that the results obtained from application of this theory accord themselves well with the known results obtained from other approaches.
In this paper we will introduce a one-dimensional reaction-diffusion model with open boundary and apply the Yang-Lee theory thereto. This model has already been studied with reflecting boundaries  while the open boundary problem has not been solved. We shall investigate the effects of the open boundary condition on the phase diagram and also the phase transition points of this model. Our work will also be another example (with more complicated reaction rules) that shows the applicability of the Yang-Lee theory to study the critical behaviors of the out of equilibrium systems.
In chapter (2) we define the model. The chapter (3) is devoted to the application of the matrix product formalism in order to find the steady states weights and therefore the partition function of the model. Using this formalism allows us to find the stationary weights of the model exactly and calculate many interesting quantities such as the density profile of particles on the chain and also the correlation functions. It turns out that our model has a first-order phase transition between a low-density and a high-density phase provided that there is a constraint on the reaction rates. In chapter (4) we shall apply the Yang-Lee theory to our model and calculate the grand canonical partition function of the system and its line of zeros in the thermodynamic limit. The line of zeros is a circle and the density of zeros is constant over it all. In the Yang-Lee theory language this signifies a first-order phase transition. Finally in section (5) we shall conclude and explain the generalizations.

## 2 The Model

Consider a one-dimensional lattice of length $`L`$. From the left boundary classical particles are injected with rate $`\alpha `$ if the target site is empty. They can also be extracted from there with rate $`\beta `$, provided that it is occupied by a particle. In the bulk of the chain particles diffuse to the left and right. When two of them meet, they can merge into a single particle. Similarly, a single particle can split into two particles. All of these processes take place preferentially to the leftward direction. Specifically, the reaction rules in the bulk of the chain are

$$\begin{array}{cccc}\text{diffusion to the left:}\hfill & & \mathrm{}+AA+\mathrm{}\hfill & \text{with rate}q\hfill \\ \text{diffusion to the right:}\hfill & & A+\mathrm{}\mathrm{}+A\hfill & \text{with rate}q^1\hfill \\ \text{coagulation at the left:}\hfill & & A+AA+\mathrm{}\hfill & \text{with rate}q\hfill \\ \text{coagulation at the right:}\hfill & & A+A\mathrm{}+A\hfill & \text{with rate}q^1\hfill \\ \text{decoagulation to the left:}\hfill & & \mathrm{}+AA+A\hfill & \text{with rate}\mathrm{\Delta }q\hfill \\ \text{decoagulation to the right:}\hfill & & A+\mathrm{}A+A\hfill & \text{with rate}\mathrm{\Delta }q^1\hfill \end{array}$$
(1)
in which $`A`$ and $`\mathrm{}`$ represent an occupied and an empty site respectively. At the left boundary we have

$$\begin{array}{cccc}\text{injection at the first site:}\hfill & & \mathrm{}A\hfill & \text{with rate}\alpha \hfill \\ \text{extraction at the first site:}\hfill & & A\mathrm{}\hfill & \text{with rate}\beta .\hfill \end{array}$$
(2)
It is seen that (1) consists of transitions which involve only the configurations of nearest neighboring sites. For $`q>1`$ the particles have a tendency to move in the leftward direction. In the following section we will use the Matrix Product Formalism to find the steady states weights of this model.

## 3 Application of the Matrix Product Formalism

The time evolution operator of our model can be written as

$$H=\underset{j=1}{\overset{L1}{}}h_{j,j+1}+h_1^{(L)}$$

in which

$$h_{j,j+1}=^{(i1)}h^{(Li1)}$$

and

$$h_1^{(L)}=h^{(L)}^{(L1)}$$

where $``$ is a $`2\times 2`$ identity matrix, $`h`$ is a $`4\times 4`$ matrix for the bulk interactions and $`h^{(L)}`$ is a $`2\times 2`$ matrix for particle input and output at the left boundary. In a basis $`(\mathrm{}\mathrm{},\mathrm{}A,A\mathrm{},AA)`$ the bulk evolution operator $`h`$ and $`h^{(L)}`$ read

$$h=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& (\mathrm{\Delta }+1)q& q^1& q^1\\ 0& q& (\mathrm{\Delta }+1)q^1& q\\ 0& \mathrm{\Delta }q& \mathrm{\Delta }q^1& q+q^1\end{array}\right),h^{(L)}=\left(\begin{array}{cc}\alpha & \beta \\ \alpha & \beta \end{array}\right).$$
(3)
In the following we show that the stationary probability distribution of any configuration of our model can be calculated exactly using the matrix product formalism (MPF). Let us first briefly review the MPF. According to this approach, for models with Hamiltonian $`H=_{j=1}^{L1}h_{j,j+1}+h_1^{(L)}+h_L^{(R)}`$, the stationary probability distribution $`P(𝒞)`$ of any configuration $`𝒞`$ is assumed to be of the form

$$P(𝒞)=\frac{1}{Z_L}W|\underset{i=1}{\overset{L}{}}(\tau _iD+(1\tau _i)E)|V$$
(4)
($`\tau _i=0`$ if the site $`i`$ is empty and $`\tau _i=1`$ if it is occupied by a particle) with the following property

$$HP(𝒞)=0.$$

The factor $`Z_L`$ in (4) is a normalization factor and can easily be obtained using the normalization condition $`_𝒞P(𝒞)=1`$. The matrices $`D`$ and $`E`$ are square matrices and stand for the presence of a particle and an empty site. These matrices beside the vectors $`W|`$ and $`|V`$ satisfy the following algebra

$$\begin{array}{c}h\left[\left(\begin{array}{c}E\\ D\end{array}\right)\left(\begin{array}{c}E\\ D\end{array}\right)\right]=\left(\begin{array}{c}\overline{E}\\ \overline{D}\end{array}\right)\left(\begin{array}{c}E\\ D\end{array}\right)\left(\begin{array}{c}E\\ D\end{array}\right)\left(\begin{array}{c}\overline{E}\\ \overline{D}\end{array}\right),\hfill \\ \\ <W|h^{(L)}\left(\begin{array}{c}E\\ D\end{array}\right)=<W|\left(\begin{array}{c}\overline{E}\\ \overline{D}\end{array}\right),\hfill \\ \\ h^{(R)}\left(\begin{array}{c}E\\ D\end{array}\right)|V>=\left(\begin{array}{c}\overline{E}\\ \overline{D}\end{array}\right)|V>.\hfill \end{array}$$

The auxiliary matrices $`\overline{E}`$ and $`\overline{D}`$ are also square matrices. For the present model with $`h^{(R)}=0`$, one can easily see from (3) and (3) that the corresponding algebra is

$$\begin{array}{c}[C,\overline{C}]=[E,\overline{E}]=0\hfill \\ \\ \overline{E}CE\overline{C}=(q+q\mathrm{\Delta }+q^1)ECq(1+\mathrm{\Delta })E^2q^1C^2\hfill \\ \\ \overline{C}EC\overline{E}=(q^1+q^1\mathrm{\Delta }+q)CEq^1(1+\mathrm{\Delta })E^2qC^2.\hfill \end{array}$$
(5)
with

$$\begin{array}{c}W|((\alpha +\beta )E+\overline{E}\beta C)=W|\overline{C}=0\hfill \\ \\ \overline{E}|V=\overline{C}|V=0\hfill \end{array}$$
(6)
in which we have defined $`C:=D+E`$ and $`\overline{C}:=\overline{D}+\overline{E}`$. Traditionally, one should find either a representation for the algebra (5) and (6)/ or calculate the stationary weights (4) directly from the algebra without using any representations. In the following we will show that a finite-dimensional representation of the algebra (5) and (6) exists under special constraint on the parameters $`\alpha `$, $`\beta `$, $`q`$ and $`\mathrm{\Delta }`$. One can easily check that for $`q^21+\mathrm{\Delta }`$ the following matrices and vectors

$$\begin{array}{c}C=\left(\begin{array}{cc}1+\mathrm{\Delta }& 0\\ 0& q^2\end{array}\right),\overline{C}=0,E=\left(\begin{array}{cc}1& \lambda \\ 0& q^2\end{array}\right),\overline{E}=\left(\begin{array}{cc}\frac{q^21}{q}\mathrm{\Delta }& \frac{\mathrm{\Delta }}{q}\lambda \\ 0& 0\end{array}\right)\hfill \\ \\ |V=\left(\begin{array}{c}\frac{\lambda }{q^21}\\ 1\end{array}\right),W|=(\frac{q\mathrm{\Delta }(q^2q\beta 1)}{\lambda (\beta +\beta \mathrm{\Delta }q\mathrm{\Delta })},1)\hfill \end{array}$$
(7)
provide a two-dimensional representation for the algebra (5) and (6) provided that $`\alpha =(q^1q+\beta )\mathrm{\Delta }`$. For all rates to be positive we assume $`q>1`$ and $`\beta qq^1`$. For $`q^2=1+\mathrm{\Delta }`$, $`C`$ is not diagonalizable and we find another representation

$$\begin{array}{c}C=\left(\begin{array}{cc}1& \lambda \\ 0& 1\end{array}\right),\overline{C}=0,E=\left(\begin{array}{cc}\frac{1}{q^2}& \lambda \\ 0& 1\end{array}\right),\overline{E}=\left(\begin{array}{cc}\frac{(q^21)^2}{q^3}& \frac{q^21}{q}\lambda \\ 0& 0\end{array}\right)\hfill \\ \\ |V=\left(\begin{array}{c}\frac{q^2\lambda }{q^21}\\ 1\end{array}\right),W|=(\frac{q^3\beta q^22q+q^1+\beta }{\lambda (q^3\beta q^2q+\beta )},1).\hfill \end{array}$$
(8)
In (7) and (8), $`\lambda `$ is a free parameters. Using the representations (7) and (8) one can easily calculate the density profile of the particles in the stationary state and also their density-density correlations. The mean density of particle on the chain at site $`i`$, $`\rho _i`$, can be written in terms of the matrices $`C`$ and $`E`$ and the vectors $`|V`$ and $`W|`$

$$\rho _i=\frac{W|C^{i1}(CE)C^{Li}|V}{W|C^L|V}.$$
(9)
By using (7) for $`q^21+\mathrm{\Delta }`$ and after some algebra we find

$`\rho _i=`$ (10)
$`{\displaystyle \frac{\mathrm{\Delta }(1+\mathrm{\Delta })^1(1q^2+q\beta )(\mathrm{\Delta }(1+\mathrm{\Delta })^Lq^{2i}(q^21)(1+\mathrm{\Delta })^iq^{2L})q^{2i+1}}{q\mathrm{\Delta }(q^21)(q^{2L}(1+\mathrm{\Delta })^L)+\beta (\mathrm{\Delta }q^2(1+\mathrm{\Delta })^Lq^{2L}(q^21)(1+\mathrm{\Delta }))}}.`$
Two different phases are apparent: $`q^2>1+\mathrm{\Delta }`$ and $`q^2<1+\mathrm{\Delta }`$. In the thermodynamic limit $`L\mathrm{}`$ we find two different expressions for the mean density of particles at site $`i`$. In the first phase, which is called the *low-density phase*, we obtain

$$\rho _i=\frac{q^2\mathrm{\Delta }(q^1q+\beta )}{(1+\mathrm{\Delta })(\beta (1+\mathrm{\Delta })q\mathrm{\Delta })}e^{i/\xi }\text{for}q^2>1+\mathrm{\Delta }.$$
(11)
As can be seen, the density of particles at the left boundary has the largest value and will quickly drop to zero in the bulk of the chain with characteristic length $`\xi `$ where $`\xi ^1=|\mathrm{ln}\frac{q^2}{1+\mathrm{\Delta }}|`$.

For the second phase, the *high-density phase*, we obtain

$$\rho _i=\frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}\frac{q^21}{1+\mathrm{\Delta }}e^{(Li)/\xi }\text{for}q^2<1+\mathrm{\Delta }$$
(12)
where $`\xi ^1=|\mathrm{ln}\frac{q^2}{1+\mathrm{\Delta }}|`$. In high-density phase the density of particles is nearly constant in the bulk of the chain equal to $`\frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}`$ and drops to zero near the right boundary. For $`q^2=1+\mathrm{\Delta }`$ we use (8) to find the following expression for the mean density of particles at site $`i`$ in the thermodynamic limit

$$\rho (x)=\frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}(1x)\text{for}q^2=1+\mathrm{\Delta }$$
(13)
in which $`x:=\frac{i}{L}`$ and $`0x1`$. In Figure (1) we have plotted the density profile of the particles on a chain of length $`L=100`$ for different values of $`\mathrm{\Delta }`$ obtained from the Monte carlo simulation. One can readily check that the Monte Carlo results accord with the exact analytical results (11)-(13). In order to see the phase transition more clearly let us define the density of particles averaged over the entire system as an order parameter

$$\rho :=\frac{1}{L}\underset{i=1}{\overset{L}{}}\rho _i.$$
(14)
Plot of (14) as a function of $`\mathrm{\Delta }`$ for $`q=3`$, $`\beta =3`$ and two different values of system length $`L=200,2000`$ is shown in the Figure (2). It can be seen that as the length of the chain approaches the thermodynamic limit $`L\mathrm{}`$, a discontinuity appears in $`\rho `$ at the transition point $`\mathrm{\Delta }_c=q^21`$ which illustrates a first-order phase transition in the system.

Using the matrix representation of the algebra (5) and (6) we can also calculate two-point density correlation functions of this model in each phase exactly. The connected two-point density-density correlation function is defined as

$$\rho _i\rho _j_c=\rho _i\rho _j\rho _i\rho _j\text{for}i<j.$$
(15)
As we saw in (9), the density of particles on the chain at an arbitrary site $`i`$ can be written in terms of the operators $`C`$ and $`E`$ and the boundary vectors $`W|`$ and $`|V`$. The expression $`\rho _i\rho _j`$ can also be written in terms of the same operators and the vectors as

$$\rho _i\rho _j=\frac{W|C^{i1}(CE)C^{ji1}(CE)C^{Lj}|V}{W|C^L|V}.$$
(16)
By using the matrix representation of the algebra we find

$$\rho _i\rho _j=\frac{q\mathrm{\Delta }^2(q^2q\beta 1)}{Z_L(\beta +\beta \mathrm{\Delta }\mathrm{\Delta }q)}[\frac{\mathrm{\Delta }(1+\mathrm{\Delta })^{L2}}{q^21}q^{2(Lj)}(1+\mathrm{\Delta }^{j2})]$$
(17)
in which $`Z_L=W|C^L|V`$ is given by (19). Now we look for the thermodynamic limit of (15) in each phase. Using (10), (15) and (17) we obtain the following expressions for the thermodynamic limit of the connected two-point correlation function of particles

$$\rho _i\rho _j_c=(\frac{\mathrm{\Delta }}{1+\mathrm{\Delta }}\rho _i)\rho _j.$$
(18)
One should use (11) or (12) for $`\rho _i`$ in each phase. This expression is valid for both $`q^2>1+\mathrm{\Delta }`$ and $`q^2<1+\mathrm{\Delta }`$ and has also a simple explanation; namely, in the low-density phase the connected two-point correlation function (15) is non-zero only in the close vicinity of the left boundary where $`i<j\xi `$. For $`j\xi `$, $`\rho _i\rho _j_c0`$. Thus, in the whole region $`\xi j<L`$, particles do not feel any correlations. In the high-density phase both $`i`$ and $`j`$ in (15) should be near the last site of the chain far from the left boundary, otherwise $`\rho _i\rho _j_c0`$.

## 4 Application of the Yang-Lee Theory

Now we will examine the Yang-Lee theory. We define the partition function of the system $`Z_L`$ as the sum of the stationary states weights given by (4)

$$Z_L=\underset{𝒞}{}P(𝒞)=\underset{\{\tau _i=0,1\}}{}W|\underset{i=1}{\overset{L}{}}(\tau _iD+(1\tau _i)E)|V=W|C^L|V.$$

Using the matrix representation of the operator $`C`$ and the vectors $`|V`$ and $`W|`$ in (7) one obtains

$$Z_L=q^{2L}+\frac{\mathrm{\Delta }q(q^2q\beta 1)(1+\mathrm{\Delta })^L}{(q^21)(\beta +\beta \mathrm{\Delta }q\mathrm{\Delta })}.$$
(19)
We look for the zeros of our partition function (19) in the complex-$`q`$ plane at fixed $`\mathrm{\Delta }`$ and $`\beta `$ in the thermodynamic limit $`L\mathrm{}`$. In Figure (3) we have plotted the numerical estimates of these roots for $`\mathrm{\Delta }=8`$, $`\beta =3`$ and $`L=100`$. It can be seen that the roots lie on a circle of radius $`q_c=\sqrt{1+\mathrm{\Delta }}=3`$, as we had predicted above. The roots also intersect the real-$`q`$ axis at an angle $`\frac{\pi }{2}`$ which is the sign of first-order phase transition.

Beside the numerical estimates, one can obtain the line of zeros and also their density near the positive real $`q`$ axis analytically. Following  we define the extensive part of the free energy as

$$g(q,\mathrm{\Delta },\beta )=\underset{L\mathrm{}}{lim}\frac{1}{L}\mathrm{ln}Z_L$$
(20)
in which $`Z_L`$ is given by (19). Now the line of zeros can be obtained from

$$Reg_1=Reg_2.$$
(21)
Here $`g_1`$ and $`g_2`$ are the thermodynamic limits of $`g`$ on the left and the right hand side of the phase transition point. By using (19), (20) and (21) and after some calculations one finds the line of zeros

$$x^2+y^2=(1+\mathrm{\Delta }).$$
(22)
In (22) the parameters $`x`$ and $`y`$ are defined as the real and the imaginary parts of $`q`$ respectively. This result is in quite agreement with the numerical estimate of zeros given in Figure (3). In order to find the density of zeros $`\mu `$ near the positive real $`q`$ axis, we use

$$2\pi \mu (s,\mathrm{\Delta },\beta )=\frac{}{s}Im(g_1g_2)$$
(23)
where $`s`$ is the arc length of the line of zeros which is measured from the critical point and increases along with the line of zeros. The functions $`g_1`$ and $`g_2`$ have the same definitions as mentioned above. After some algebra we find $`\mu =\frac{1}{\pi \sqrt{1+\mathrm{\Delta }}}`$. It is seen that the density of zeros is constant everywhere on the circle beside in the vicinity of the positive real $`q`$ axis; therefore, it confirms the existence of a first-order phase transition in our model.

## 5 Concluding Remarks

In this paper we introduced a one-dimensional out-of-equilibrium model consisting of one class of particles which diffuse, coagulate and decoagulate on an open chain. The left boundary of the chain is assumed to be open so that the particles can enter or leave the system. Exact calculations using the MPF show that on a specific plane in the parameters space a first-order phase transition takes place between a low-density and a high-density phase. The density profile of the particles in each phase is obtained exactly. In the low-density phase the density profile of particles far from the left boundary drops exponentially to zero with the length scale $`\xi =|\mathrm{ln}\frac{q^2}{1+\mathrm{\Delta }}|^1`$. We should mention that in the stationary state there exists only one characteristic length scale $`\xi `$, while the same model without injection and extraction of particles has three different length scales . In the high-density phase the density of particles is nearly constant in the bulk; however, it drops to zero near the right boundary with the same length scale $`\xi `$. It is shown that the steady state of our model can be written in terms of the superposition of shocks . In order to answer the question whether we can apply the classical Yang-Lee theory to the out-of-equilibrium systems such as the one introduced here, we have seen that the application of this theory gives same results obtained from the MPF approach. We calculated the line of partition function zeros exactly. It intersects the positive real $`q`$ axis at right angles at $`q_c=\sqrt{1+\mathrm{\Delta }}`$. The density of zeros is uniform at the cross point, which means a first-order phase transition at critical point. Our model, which is more complex than those studied before, is another example for the generality of the classical Yang-Lee picture of phase transitions in one-dimensional out-of-equilibrium models. It can be shown that the MPF also works when the particles are allowed to leave the chain from the right boundary; however, there should be still a constraint on the reaction rates. The whole phase diagram of this model which can be studied at the mean field level or using the Monte-Carlo simulation, is still an open problem.

## References


text:
# Modeling Supply Chains and Business Cycles as Unstable Transport Phenomena

## I Modelling Supply Networks

The production units: We will investigate a system with $`u`$ production units (e.g. machines or factories) $`a,b,c\{1,2,\mathrm{},u\}`$ producing $`p`$ products $`i,j,k\{1,2,\mathrm{},p\}`$. The respective production process is characterized by parameters $`n_b^j`$ and $`f_b^i`$: In each production step, production unit $`b`$ requires $`n_b^j`$ products (educts) $`j\{1,\mathrm{},p\}`$ and produces $`f_b^i`$ products $`i\{1,\mathrm{},p\}`$. The number of production steps of production unit $`b`$ per unit time shall be $`P_b(t)`$. It is the product of the processing (departure) rate $`\mu _b`$ and the probability $`p_b`$ of the production unit being occupied:

$$P_b(t)=\mu _b(t)p_b(\rho _b,c_b,s_b).$$
(1)
Herein,

$$\rho _b=\lambda _b/\mu _b$$
(2)
denotes the channel utilization, $`\lambda _b`$ the feeding rate, $`c_b`$ the number of parallel channels of the production unit, and $`s_b`$ the storage capacity. Changing the number $`c_b`$ of channels or the storage capacity $`s_b`$ is costly. For the time being, we will therefore assume $`c_b`$ to be constant. The expression for $`p_b(\rho _b,c_b,s_b)`$ is well-known for the stationary state of many queuing systems . For example, for a $`M/M/1:(\mathrm{}/\text{FIFO})`$ process (i.e. one channel with first-in-first-out serving, unlimited storage capacity, Poisson-distributed arrival times and exponentially distributed service intervals), one finds

$$p_b(\rho _b,1,\mathrm{})=\mathrm{min}(\rho _b,1),$$
(3)
so that, in this particular case, we have the relation

$$P_b(t)=\mu _b(t)\mathrm{min}[\rho _b(t),1]=\mathrm{min}[\lambda _b(t),\mu _b(t)].$$
(4)
Feeding rates: The feeding rate $`\lambda _b`$ is not just the product of the concentrations of educts required for a product, as that would be the case for chemical reactions. Instead, the feeding rate is determined by the minimum arrival rate $`A_{ab}^j`$ of required educts $`j`$, divided by the number $`n_b^j`$ of educts required for one production step. (The difference is comparable to the Probabilistic AND and the Fuzzy AND in Fuzzy Logic.) As the arrival rates $`A_{ab}^j`$ are given by the input buffer flows $`\nu _b^jI_b^j`$ (where $`\nu _b^j`$ is the rate of getting educt $`j`$ from the the input buffer into the production unit $`b`$ and $`I_b^j`$ is the number of educts stored in the input buffer), the resulting feeding rate is

$$\lambda _b(t)=\underset{j}{\mathrm{min}}\left(\frac{\nu _b^j(t)I_b^j(t)}{n_b^j}\right).$$
(5)
Note that, due to stochastic variations (which can cause queuing effects), efficient production is related to $`\lambda _br\mu _b`$ with $`1>r0.7`$ . Therefore, it would be reasonable to use transport rates $`\nu _b^j`$ according to $`\nu _b^j=r\mu _bn_b^j/I_b^j`$. However, usually there are capacity constraints $`V_b^j`$ in the transport of product $`j`$ in production unit $`b`$, which may be reflected by a function $`U(x,V)`$ such as

$$U(x,V)=\frac{x}{1+x/V}\{\begin{array}{cc}x\hfill & \text{ for }xV\hfill \\ V\hfill & \text{ for }xV.\hfill \end{array}$$
(6)
Therefore, we will make the specification

$$\nu _b^j(t)=U(\frac{r\mu _b(t)n_b^j}{I_b^j(t)},V_b^j(t)).$$
(7)
Input and output buffers: Let us assume, each production unit $`b`$ has input buffers for required educts $`i`$ and output buffers (a warehouse) for the products. We will assume the input buffers are filled with $`I_b^i(t)`$ educts $`i\{1,\mathrm{},p\}`$ and the output buffers with $`B_b^i(t)`$ products waiting to be delivered. If $`D_{ab}^i(t)`$ denotes the delivery flow of products $`i`$ from production unit $`a`$ to $`b`$, the change of an input buffer stock with time is given by the conservation equation

$$\frac{dI_b^i}{dt}=\underset{a}{}D_{ab}^i(t)n_b^iP_b(t),$$
(8)
as $`_aD_{ab}^i(t)`$ is the number of products $`i`$ delivered from various sources (production units) $`a`$ and $`n_b^iP_b`$ is the number of educts $`i`$ used up for production per unit time. Analogously, the dynamics of an output buffer stock is determined by the equation

$$\frac{dB_b^j}{dt}=f_b^jP_b(t)\underset{c}{}D_{bc}^j(t),$$
(9)
as $`f_b^jP_b`$ is the number of newly generated products $`j`$, and $`_cD_{bc}^j(t)`$ are deliveries to other production units $`c`$. The delivery flow $`D_{ab}^i`$ is given by the delivery rate $`\nu _{ab}^i`$ times the available products $`B_a^i`$ in the buffer:

$$D_{ab}^i(t)=\nu _{ab}^i(t)B_a^i(t).$$
(10)
The delivery flows are adapted to the order flows $`O_{ab}^i`$. Ideally, one would have the relation $`D_{ab}^i=O_{ab}^i`$, but due to capacity constraints $`V_{ab}^i`$, we will again assume

$$\nu _{ab}^i(t)=U(\frac{O_{ab}^i(t)}{B_a^i(t)},V_{ab}^i(t)),$$
(11)
which implies $`D_{ab}^iO_{ab}^i`$.

Adaptation of capacities: One important aspect of production is the adaptation of production capacities. This requires usually a lot of time compared to the other time scales of production. The adaptation time $`T_b`$ of the processing rate $`\mu _b`$ is much greater than the other adaptation times, so that we will assume the relation

$$\frac{d\mu _b}{dt}=\frac{1}{T_b}\left[W(\underset{j}{}\frac{f_b^j}{B_b^j},\mathrm{})\mu _b\right].$$
(12)
Herein, $`W(y,\mathrm{})`$ is some control function reflecting the strategy by the production management to adapt the processing rate $`\mu _b`$, e.g. as a function of the output buffer stocks $`B_b^j`$ or other variables. Normally, it will increase with decreasing output buffer stocks $`B_b^j`$, but the function $`W(y,\mathrm{})`$ saturates due to financial, spatial, or technological limitations and inefficiencies in the processing of high order flows. In this study, we will assume a function of the form

$$W(1/z,\mathrm{})=A\frac{1+Bz}{1+Cz+Dz^2}$$
(13)
with $`z=1/y`$ and suitable parameters $`A`$, $`B`$, $`C`$, and $`D`$.

Of course, the transportation capacities $`V_b^j`$ and $`V_{ab}^i`$ are adapted as well. One may set up separate equations for this, but for simplicity, we will here assume the proportionality relations

$$V_b^j(t)=\widehat{V}_b^j\mu _b(t)\text{and}V_{ab}^i(t)=\widehat{V}_{ab}^i\mu _b(t).$$
(14)
Order flows, delivery networks, and price dynamics: The flow of orders is basically given by the flow $`P_bn_b^i`$ of educts required for the production by unit $`b`$. If the proportions $`q_{ab}^i`$ with $`_aq_{ab}^i=1`$ reflect orders from different producers, we have

$$O_{ab}^i=q_{ab}^iP_bn_b^i.$$
(15)
The fractions $`q_{ab}^i`$ characterize the delivery or supply network. One can imagine various scenarios. For example, the fractions could be modelled by an evolutionary selection equation with a selection rate $`\nu _q`$ :

$$\frac{dq_{ab}^i}{dt}=\nu _q\left(F_{ab}^i(t)\underset{a^{}}{}F_{a^{}b}^i(t)q_{a^{}b}^i(t)\right)q_{ab}^i(t).$$
(16)
It makes sense to relate the fitness $`F_{ab}^i`$ to the inverse price $`p_{ab}^i`$ of product $`i`$, which producer $`a`$ takes from production unit $`b`$:

$$F_{ab}^i(t)=1/p_{ab}^i(t).$$
(17)
There are many pricing strategies, but the common element seems to be that the price changes according to the law of supply and demand. Therefore, one conceivable specification would be

$$p_{ab}^i(t)=p_i^0\frac{O_{ab}^i(t)}{D_{ab}^i(t)},$$
(18)
where $`p_i^0`$ is the price when supply $`D_{ab}^i`$ and demand $`O_{ab}^i`$ agree.

Other specifications are, of course, possible as well, as the above formulas partly depend on the strategies of the human decision makers involved.

## II Markets and Business Cycles

We will now introduce a simplified and aggregate description of an economy. For this, we will summarize a whole economic sector by one production unit $`b`$ and replace input and output buffers by markets for the products $`i`$ of the $`u`$ different economic sectors. To describe this, we define the macroscopic stock level of products $`i`$ in the market by

$$N_i(t)=\underset{b}{}I_b^i(t)+\underset{b}{}B_b^i(t).$$
(19)
According to the balance equations (8) and (9), we find

$$\frac{dN_i}{dt}=\underset{b}{}(f_b^in_b^i)P_b(t).$$
(20)
Furthermore, in the arguments of the remaining mathematical relations, we have to replace input buffer stocks $`I_b^i`$ and output buffer stocks $`B_b^i`$ by the respective number $`N_i`$ of products in the market, as the market replaces the buffers. In this way, we obtain the feeding rates

$$\lambda _b(t)=\underset{j}{\mathrm{min}}\left(\frac{\nu _b^j(t)N_j(t)}{n_b^j}\right),$$
(21)
which enter $`P_b`$ according to equation (1), the transport rates

$$\nu _b^j(t)=U(\frac{r\mu _b(t)n_b^j}{N_j(t)},V_b^j),$$
(22)
and the processing rates

$$\frac{d\mu _b}{dt}=\frac{1}{T_b}\left[W(\underset{j}{}\frac{f_b^j}{N_j(t)},\mathrm{})\mu _b\right].$$
(23)
For $`P_b\lambda _b`$, one would have to solve this equation together with

$$\frac{dN_i}{dt}=\underset{b}{}(f_b^in_b^i)\mu _b(t)\underset{j}{\mathrm{min}}\left[U(\frac{rn_b^j}{N_j(t)},\widehat{V}_b^j)\frac{N_j(t)}{n_b^j}\right].$$
(24)
For a linear supply chain with $`n_b^i=\delta _{b,i+1}`$ and $`f_b^i=\delta _{b,i}`$ (where $`\delta _{k,l}=1`$, if $`k=l`$, otherwise 0), this leads to

$$\frac{dN_i}{dt}=V_i(t)N_{i1}(t)V_{i+1}(t)N_i(t)$$
(25)
with

$$V_i(t)=\mu _i(t)U(\frac{r}{N_{i1}(t)},\widehat{V}_i^{i1}).$$
(26)
Interestingly, Eqs. (23) and (25) basically agree with the traffic flow model by Hilliges and Weidlich , where (25) is analogous to the equation for the vehicle density and (23) corresponds to the velocity equation. These equations are known to behave linearly unstable with respect to perturbations of $`N_i(t)`$, if $`T_b`$ exceeds a certain threshold which depends on the maximum slope $`dW/dN_i`$ . This is the reason why one can frequently observe an instability of supply chains, called the “bullwhip effect” (e.g. in the “beer distribution game” ). A thorough analysis shows that the Hilliges-Weidlich model would, in fact, be more suitable for the description of supply chains than of traffic: First, their velocity equation does not contain a convection term, as it should be the case for traffic. Second, their equations normally do not show metastability , i.e. a region in which traffic flow breaks down, when a perturbation exceeds a critical amplitude, but where smaller perturbations fade away. For this reason, a supply chain of the above kind is expected to behave either stable or linearly unstable, no matter how large the perturbation amplitude is. The most interesting point, however, is the reaction of the system to a periodic perturbation in the consumption rate $`V_{p+1}(t)`$, as the resulting variations in the stock levels are synchronized and much slower (see Fig. 1). For this reason, they may explain business cycles as self-organized phenomenon.

## III Summary and Outlook

In this contribution, a theory of supply networks has been sketched, which may open a new area of econophysics . This theory is developed to help understand the dynamical phenomena, breakdowns, instabilities and inefficiencies of production processes and supply networks. Future work will have to address questions such as the relevance of the network structure for the resulting dynamics, possible control strategies, the role of the market and pricing mechanism, etc. This specifically concerns the choice of the factor $`r`$, the specification of the functions $`U`$, $`W`$, and the equations suggested in the paragraph on delivery networks and price dynamics. Here, we have focussed on an application to markets, i.e. macroeconomics. For this, we have replaced input and output buffers by markets for different products \[see Eq. (19)\]. One may view this as some kind of micro-macro link, as we have started with equations for single firms (a microeconomic level), e.g. Eqs. (8), (9), and ended up with equations (20) for market sectors. Assuming a sectoral structure of economics, one can relate the resulting equations with the Hilliges-Weidlich model, which has originally been developed for traffic flow. These equations can describe the “bullwhip effect” due to their linear instability in a certain regime of operation. The underlying mechanism is the slow adaptation of the processing rate to changes in the order flows or stock levels in the market. Interestingly, the resulting oscillations in the stock levels of the different products $`i`$ have a characteristic frequency, which can be much lower than the underlying fast variations in the consumption rate. These oscillations synchronize among different economic sectors and may explain business cycles as a self-organized phenomenon with slow dynamics. As the variables in the model are operational and measurable, the model can be tested and calibrated with empirical data. Apart from some equations which were not further applied in this study, most of the proposed model equations were conservation equations, equations given by the product flows, or relations derived with stochastic concepts used in queuing theory. They reflect the transport and interaction of products, so that the physics of driven many-particle systems and of complex systems can make some significant contributions to the new multidisciplinary field of self-organization phenomena in production and supply networks.


text:
# Untitled Document

An alternative to the variation of the fine structure constant

M. E. DE SOUZA

Departamento de Física, Universidade Federal de Sergipe, Campus Universitário, 49100-000 São Cristóvão, Sergipe, Brazil - e-mail: mdesouza@ufs.br

Considering that quasars evolve into normal galaxies the average expansion velocity of matter in quasars (or galaxies) can be estimated approximately. Such velocity makes us to see the light of quasars slightly Doppler shifted which explains the misleading variation of the fine structure constant.

There have been reports1-4 in the literature proposing a time variation of the fine structure constant. They are based on observations of absorption lines of gas clouds seen against background quasars. Thus the clouds are quasar absorbers. This supposed time variation of $`\alpha `$ motivated the recent and very important article of Davies5 which proposes that the velocity of light may vary with time. This would cause profound changes in the foundatins of physics.

This work presents an alternative explanation for the phenomenon observed in the gas clouds which maintains $`\alpha `$ as a fundamental constant of nature, that is, independent of time. It is based on a proposal for the evolution of galaxies from quasars. The remarkable work of Márquez et al.$`^6`$has shown that very high redshift elliptical galaxies appear to harbor quasars. They have also shown that such galaxies are very small and have diameters smaller than 1kpc. All the studied objects (about 15 quasars) have extended structures of ionized gas around them (this fact had already been presented by other researchers). They have found other galaxies in the fields of the studied objects only a few kpc away from them. Some of the quasars present asymmetric radio sources with collimated one-sided jets of extended ionized gas. This means that galaxies are born as quasars that become galaxies by means of the shedding of matter (ionized gas) from their cores. The same kind of phenomenon has been observed in galaxies. Very recent data$`^7`$of NGC 6240, which is considered a typical protogalaxy, show that “approximately 70% of the total radio power at 20cm originates from the nuclear region ($``$1.5kpc), of which half is emitted by two unresolved (R$``$30pc) cores and half by a diffuse component. Nearly all of the other 30% of the total radio power comes from an arm-like region extending westward from the nuclear region”. A very important property of many quasars is their brightness which can vary from night to night. This flickering may have its origin in the outward motion of large quantities of matter from their cores. This brightness variability is also present in Seyfert galaxies which are powerful sources of infrared radiation. Many of them are also strong radio emitters. For example, over a period of a few months, the nucleus of the Seyfert galaxy M77(or NGC1068) switches on and off a power output equivalent to the total luminosity of our galaxy8. It is also worth noting that the nuclei of Seyfert galaxies are very bright and have a general starlike appearance. Researchers have found that some Seyfert galaxies exhibit explosive phenomena8. For example, M77 and NGC4151 expel huge amounts of gas from their nuclei. The spectra of both galaxies show strong emission lines, just as quasars. Shaver et al.$`^9`$ have found that there is almost no quasar for $`z<0.5`$. This clearly shows that quasars evolve into galaxies.

This ejection of matter is quite common in active galaxies also. Let us mention just some of them. NGC 2992 presents a jet-like structure and a circum-nuclear ring10. Falcke and Biermann$`^{11}`$ report that there is a large scale emission-like jet going outward from the core of NGC 4258 with a mass of about $`4\times 10^{35}`$kg and with a kinetic power of approximately $`10^{42}`$ergs/s and expansion velocity of 2000km/s. This is of the same order of supernovae velocities. It is well known that BL Lacertae objects are powerful sources of radio waves and infrared radiation. They share with quasars the fact of exhibiting a starlike appearance and of showing short-term brightness fluctuations. As some quasars do, they also have a nebulosity around the bright nucleus. Researchers$`^{12}`$ have managed to obtain the spectrum of their nebulosity. The spectrum of the nebulosity is strikingly similar to the spectrum of an elliptical galaxy(M32s spectrum, in this case). In terms of the evolution above described they are simply an evolutionary stage of a quasar towards becoming a more quiet galaxy. Radio galaxies share with BL Lacertae objects many of the properties of quasars. As Heckman et al.$`^{12}`$ have shown, in the middle and far infrared (MFIR) quasars are more powerful sources of MFIR radiation than radio galaxies. Also, there have been investigations showing that the emission from the narrow-line region(NLR) in radio-loud quasars is stronger than in radio galaxies of the same radio power13,14,15. Goodrich and Cohen$`^{16}`$ have studied the polarization in the broad-line radio galaxy 3C 109. After the intervening dust is taken into account the absolute V-magnitude of this galaxy becomes $`26.6`$ or brighter, which puts it in the quasar luminosity range. The investigators suggest that “many radio galaxies may be quasars with their jets pointed away from our direct line of sight”. It has also been established that radio galaxies are found at intermediate or high redshifts and that they are clearly related to galactic evolution because as the redshift increases cluster galaxies become bluer on average, and contain more young stars in their nuclei. This is also valid for radio galaxies: the higher the redshift, the higher their activity. All these data show that a radio galaxy is just an evolutionary stage of an active galaxy towards becoming a quiet galaxy, i.e., it is just a stage of the slow transformation by means of an overall expansion of a quasar into a normal galaxy.

In the light of the above considerations the nuclei of old spirals must exhibit a moderate activity. This is actually the case. The activity must be inversely proportional to the galaxys age, i.e., it must be a function of luminosity. The bluer they are, the more active their nuclei must be. As discussed above there must also exist a relation between this activity and the size of the nucleus(as compared to the disk) in spiral galaxies. Our galaxy has a mild activity at its center. Most of the activity is concentrated in a region called Sagittarius A, which includes the galactic center. It emits synchroton and infrared radiations. Despite its large energy output Sagittarius A is quite small, being only about 40 light years in diameter. Besides Sagittarius A our galaxy exposes other evidences showing that in the past it was a much more compact object: a) Close to the center, on opposite sides of it, there are two enormous expanding arms of hydrogen going away from the center at speeds of 53km/s and 153km/s; b) Even closer to the center there is the ring called Sagittarius B2 which is expanding at a speed of 110km/s(8,17). It is worth noting that the speeds are low(as compared to the velocities of relativistic electrons from possible black holes). This phenomenon is not restricted to our galaxy. Recent high-resolution molecular-line observations of external galaxies have revealed that galactic nuclei are often associated with similar expanding rings18.

This paper will not deal with the mechanism of the expulsion of matter from the center of quasars and will be restricted to showing what may be happening in the case of the clouds absorption lines. It is well known that the size of a quasar is smaller than 1kpc and that a typical galaxy has a radius of about 50.000 light-years $`5\times 10^{20}m`$ and, as the age of the Universe is about $`10^{17}s`$, the average expansion velocity of matter in galaxies since the time when they were quasars is of the order of

$$<V_e>5\times 10^3m/s.$$
(1)
This figure is just a rough estimate. For very high redshift quasars the expansion speed was at least an order of magnitude higher as it is inferred from the data on shedding of matter from the cores of quasars and galaxies. Therefore, the light of quasars should be Doppler shifted to higher frequencies for the light that is coming from points close to the line of sight. This means that the quasar light that reaches a gas cloud is slightly Doppler shifted due to the expansion of matter with respect to a rest-frame sitting in the center of the quasar. The shift is given by

$$\frac{\mathrm{\Delta }\nu }{\nu _0}=\frac{<V_e>}{c}1.67\times 10^5$$
(2)
which is of the same order of magnitude of the shifts attributed to variation of the fine structure constant around $`z=1`$ in the work of Webb et al.4. Take notice that if we attribute the above result to a variation of the fine structure function we obtain

$$\frac{\mathrm{\Delta }\alpha }{\alpha }=\frac{1}{2}\frac{\mathrm{\Delta }\nu }{\nu _0}=0.835\times 10^5$$
(3)
which is a figure of the order of magnitude of that found in Webb et al.4. Thus, comparing results from quasars at two different redshifts $`z_1`$ and $`z_2`$ (with $`z_1>`$ $`z_2`$), since more active quasars have higher redshifts, it is expected $`<V_e>(z=z_1)`$ to be larger than $`<V_e>(z=z_2)`$ and, of course, $`|\mathrm{\Delta }\nu (z_1)|>|\mathrm{\Delta }\nu (z_2)|`$. For very high redshift quasars $`\frac{\mathrm{\Delta }\nu }{\nu _0}`$ may even reach $`10^4`$ or $`10^3`$ which is in line with the upper limit of Cowie et al.19 according to whom $`|\frac{\mathrm{\Delta }\alpha }{\alpha }|<3.5\times 10^4`$ for $`z3`$.

Therefore, the observations of the variation of the fine structure constant may reveal, actually, how quasars evolve towards becoming quiet, normal galaxies. And thus, the foundations of physics continue being solid as ever!

References

1. Marciano, W. Phys. Rev. Lett. 52, 489 (1984)

2. Barrow, J. D. Phys. Rev. D 35, 1805 (1987)

3. Damour, T. and Polyakov, A. M. Nucl. Phys. B423, 532 (1994)

4. Webb, J. K. et al. Phys. Rev. Lett. 82, 884 (1999)

5. Davies, P. C. W. et al. Nature 418, 603 (2002)

6. Márquez, I. et al. astro-ph/9810012.

7. Colbert, E. J. M. et al. in The Radio Emission from the Ultra-Luminous Far-Infrared Galaxy NGC 6240,

astro-ph/9405046.

8. Kaufmann,III, W. J. in Galaxies and Quasars (W.H.Freeman and Company, San Francisco, 1979).

9. Shaver, P. A. et al.  astro-ph9801211.

10. Chapman, S. C. et al. astro-ph/9810250.

11. H. Falcke, H., Biermann, P. L. astro-ph/9810226.

12. Heckman, T. M. et al.  Ap.J., 391, 39 (1992).

13. Baum, S, Heckman, T. M. Astrophys. J 336, 702 (1989).

14. Jackson, N, Browne I. Nature, 343, 43 (1990).

15. Lawrence, A. Mon. Not. R. Astr. Soc., 1992.

16. Goodrich, R. W., Cohen, M. H. Astrophys. J. 391, 623 (1992).

17. Sofue, Y. Astro. Lett. Comm. 28, 1 (1990).

18. Nakai, N. et al. Pub. Astr. Soc. Japan 39, 685 (1987).

19. Cowie, L. L., Songalia, A. Astrophys. J. 453, 596 (1995)


text:
# Spin relaxation in (110) and (001) InAs/GaSb superlattices

## Abstract

We report an enhancement of the electron spin relaxation time ($`T_1`$) in a (110) InAs/GaSb superlattice by more than an order of magnitude (25 times) relative to the corresponding (001) structure. The spin dynamics were measured using polarization sensitive pump probe techniques and a mid-infrared, subpicosecond PPLN OPO. Longer $`T_1`$ times in (110) superlattices are attributed to the suppression of the native interface asymmetry and bulk inversion asymmetry contributions to the precessional D’yakonov Perel spin relaxation process. Calculations using a nonperturbative 14-band nanostructure model give good agreement with experiment and indicate that possible structural inversion asymmetry contributions to $`T_1`$ associated with compositional mixing at the superlattice interfaces may limit the observed spin lifetime in (110) superlattices. Our findings have implications for potential spintronics applications using InAs/GaSb heterostructures.

Semiconductor heterostructures based on the InAs/GaSb/AlSb system have gained attention recently for potential applications in the rapidly growing field of semiconductor-based spintronics,Cartoixa:2001 ; Hammar:2002 ; Gilbert:2002 ; Heida:1998  in which the focus is to develop novel electronic devices and circuits that exploit the spin property of the electron.Wolf:2001 ; ALS:book  By utilizing some of the unique characteristics of the InAs/GaSb/AlSb system, including the strong spin orbit interaction in InAs and GaSb, the ability to form a variety of interface types (type-I, type-II staggered, and type-II broken gap), and the high electron mobility of InAs, it may be possible to realize a host of novel high-speed spin-sensitive electronic devices. A crucial factor in the design of any spin-sensitive electronic device is the electron spin relaxation time (T1), which must be sufficiently long to allow for transport and$`/`$or processing of the spin-polarized electrons. In InAs/GaSb/AlSb heterostructures, which contain bonds between semiconductor constituents with no atom in common, the nature of the interfaces can have a dramatic influence on the electron spin dynamics. The electron spin relaxation time was recently measured in (001) InAs/GaSb superlattices,Olesberg:2001  where it was found that the asymmetric potential at the interfaces, referred to as native interface asymmetry (NIA),ALS:book ; Olesberg:2001 ; Lau:2002 ; Vervoort:1999 ; Krebs:1996  strongly dominates the precessional spin relaxation process. In addition to affecting the band gap and intervalence band absorption properties of these superlattices,Lau:2002  NIA leads to spin relaxation times $``$ 1 ps, which are more than an order of magnitude shorter than expected considering bulk inversion asymmetry (BIA)Lau:2001  contributions alone.

The interface contribution to spin relaxation in these no-common-atom superlattices may be removed through (110) growth, since the mixed anion-cation interface planes in this orientation lead to a symmetric interface potential. (See Fig. 1.) Since NIA strongly dominates spin relaxation in (001)-oriented InAs/GaSb superlattices, the electron spin lifetime is expected to be considerably longer in (110) InAs/GaSb superlattices. Additionally, suppression of the BIA contribution in (110) heterostructures has been predictedDK:1986  and recently verified experimentally in GaAs/AlGaAs quantum wells.Ohno:1999  This suppression has not been seen in any other material system. Here we observe it in InAs/GaSb superlattices.

We measured spin relaxation times in (110) and (001) InAs/GaSb superlattices using a mid-infrared femtosecond optical parametric oscillator and polarization sensitive differential transmission experiments. We observed an increase in $`T_1`$ by more than an order of magnitude for (110) growth. Theoretical calculations using a nonperturbative 14-band $`𝐤𝐩`$ nanostructure model agree well with experiment. Due to the suppression of NIA and BIA contributions to spin relaxation in (110) InAs/GaSb superlattices, our calculations indicate that $`T_1`$ is highly sensitive to asymmetry arising from compositional mixing at the superlattice interfaces, an effect that may limit $`T_1`$ in (110) InAs/GaSb superlattices.

The InAs/GaSb superlattices were grown by molecular beam epitaxy on (001) or (110)-oriented InAs substrates (+/- 0.1 degree) in a Fisons VG-80 machine equipped with shuttered EPI group-III evaporators and shuttered EPI valved group V cracker cells. The superlattice structures contain 2.1 nm thick InAs layers and 3.7 nm thick GaSb layers, with a total of 85 periods. The (110) and (001) superlattices were grown in consecutive runs to minimize differences in extrinsic sample characteristics. X-ray diffraction measurements indicate high order resolution of the SL peaks, and features of the RHEED pattern are highly streaked, indicating good quality superlattice growth. The room temperature continuous wave photoluminescence signal from the superlattices was of similar intensity and line width for the (001) and (110) structures, and the carrier recombination time, measured by time-resolved photoluminescence up-conversion experiments, was found to be similar in the two superlattices (1-2 ns), indicating consistency of sample quality.

$`T_1`$ times were measured at low temperature (115 K) using pump probe experiments involving 200 fs, circularly polarized mid-infrared pulses (4.0-4.6 $`\mu `$m). The laser source for these measurements is an optical parametric oscillator (OPO) that relies on parametric down conversion in a 1 mm thick periodically poled LiNbO3 (PPLN) crystal.Ebrahimzadeh:2001  The OPO is synchronously pumped by a Titanium sapphire oscillator, and provides 20 to 100 mW of average idler output power over the spectral tuning range (2.7-4.6 $`\mu `$m). The output pulses have a full-width at half maximum bandwidth of 18 meV.

The confinement-induced splitting of the band edge heavy-hole and light-hole states is calculated to be $``$ 130 meV for the (001) and (110) InAs/GaSb superlattices considered here. In this case, circularly polarized idler pulses tuned close to the band gap will excite only the heavy-hole to conduction band transition, leading to a 100% spin-polarized distribution of electrons. The band gap of the (110) superlattice is slightly lower ($``$ 40 meV) than the corresponding (001) structure. Band structure calculations for these superlattices suggest that this difference is due to lower electron and hole confinement energies in the (110) superlattice in addition to the differing influence of strain for different directions of superlattice growth. In order to compare measured spin lifetimes for identical electron kinetic energies, the OPO was tuned to 35 meV above the band gap for each structure. The limited tuning range of the OPO, together with the low band gaps of the InAs/GaSb superlattices studied \[240 meV (110) and 280 meV (001) at 77 K\], restricted our experiments to temperatures between 77 K and 115 K, over which the spin lifetimes were not found to vary significantly.

The optically-injected spin-polarized electrons generate absorption bleaching of the associated interband transition through state filling. By monitoring the transmission of a weaker, delayed probe pulse that has the same (SCP) or opposite (OCP) circular polarization as the pump pulse, $`T_1`$ may be extracted. The differential transmissivity was measured versus the time delay between the pump and probe pulses using a liquid N2-cooled InSb detector and lock-in detection methods. In our experiments, care was taken to ensure that the pump beam impinges on the superlattice sample at normal incidence. In this case, the probe beam interrogates the relaxation time of electron spins polarized in the growth direction. The optically injected carrier density was estimated using the measured pump pulse fluence and the calculated absorption spectrum for each structure to be 1$`\times `$1016 cm-3 for (110) and 3$`\times `$1016 cm-3 for (001), with an uncertainty of $`\pm `$ 40 %.

Fig. 2 shows the results of polarization-sensitive differential transmission measurements on InAs/GaSb superlattices with (001) (Fig. 2(a)) and (110) (Fig. 2(b)) directions of growth. Note the time scales on the x-axes in Fig. 2(a) and Fig. 2(b), which differ by more than an order of magnitude. The large negative pulse-width-limited feature appearing in both data sets originates from two-photon absorption in the InAs substrates. This feature provides a convenient marker for zero time delay, but may otherwise be ignored for the purpose of this discussion. The loss of spin polarization in the optically-injected carrier distribution is indicated by the convergence of the differential transmission signals for the SCP and OCP polarization geometries. ¿From a comparison of the data in Fig. 2(a) and Fig. 2(b), it is clear that spin relaxation occurs on a much longer time scale in the (110) InAs/GaSb superlattice compared to the (001) structure. The OCP signal in Fig. 2(b) grows from $`0`$, indicating that the optically excited carriers are initially fully spin polarized, and that no detectable Coulomb screening effects are present. In the data for the (001) superlattice in Fig. 2(a), the non-zero bleaching signal in the OCP data immediately following the two-photon absorption feature is evidence of spin decay during pulse overlap. $`T_1`$ times are obtained by performing a single exponential fit to the difference between the SCP and OCP curves. Fits were restricted to delay values $``$ 350 fs to avoid the two-photon absorption feature. From the data in Fig. 2, we find $`T_1`$ values of 700 fs and 18 ps for the (001) and (110) superlattices respectively, indicating that the spin relaxation time in InAs/GaSb superlattices is strongly enhanced with (110) growth.

The much longer spin lifetime in the (110) InAs/GaSb superlattice relative to the corresponding (001) structure is primarily attributed to the elimination of the NIA contribution to spin relaxation. Although a larger $`T_1`$ is expected in (110) superlattices due to suppression of BIA,DK:1986  for the short period no-common-atom superlattices investigated here the more crucial distinction arises from the nature of the interfaces. Fig. 1 contains a schematic diagram of the interfaces in InAs/GaSb superlattices for conditions of growth along the (001) \[Fig. 1(a)\] and (110) \[Fig. 1(b)\] directions. For the (001) superlattice, the potential at the interfaces is asymmetric due to tetrahedral bonding between the bulk constituents (InAs and GaSb) that share no atom in common. This effect is referred to as NIA. As in the case of BIA, NIA generates a pseudomagnetic field that serves to relax the electron spins. It has been shown previouslyOlesberg:2001  that NIA strongly dominates the electron spin decay process in a similar (001)-oriented InAs/GaSb superlattice grown on a GaSb substrate. The subpicosecond spin lifetime extracted from the data of Fig. 2(a) is consistent with these earlier studies.Olesberg:2001  In contrast, for an InAs/GaSb superlattice grown in the (110) direction, NIA is absent because the interfaces are mixed, containing equal numbers of both anions and cations (see Fig. 1(b)). In this case, the interface potential is symmetric. Since spin decay in the (001) superlattice is strongly dominated by NIA, with its elimination in the (110) structure the spin lifetime is expected to sharply increase, consistent with the results in Fig. 2(b).

In order to interpret these findings, we have performed calculations of spin relaxation times due to the precessional D’yakonov Perel mechanism, which has been found to dominate in III-V semiconductors and their heterostructures for temperatures above 77 K.ALS:book ; Olesberg:2001 ; Lau:2001 ; OpticalOrientationBook ; Kikkawa:1998 ; Tacheuchi:1997 ; Terauchi:1999 ; Boggess:2000  We employ a nonperturbative $`𝐤𝐩`$ nanostructure model solved in a fourteen-band restricted basis set, in which BIA contributions are included naturallyLau:2001  and NIA contributions are introduced as previously described.Olesberg:2001  Parameters governing the interface asymmetry were taken from reference 7. Neutral impurities were taken as the dominant source of momentum scattering. The spin lifetime varies inversely with momentum scattering time, $`T_1\tau _p^1`$. We find good agreement with experiment using $`\tau _p`$ = 190 fs, which is within the typical range (100-200 fs) for the III-V semiconductors.

The calculated spin relaxation time for the (001) superlattice is $`T_1`$ = 715 fs in good agreement with the experimental value of 700 fs. NIA was found to strongly dominate the electron spin relaxation process for the (001) structure, consistent with earlier findings.Olesberg:2001  For the (110) InAs/GaSb superlattice, our calculations indicate that, with the elimination of NIA and the strong suppression of the remaining BIA contribution,DK:1986  the spin lifetime is expected to be $``$ 600 ps. The much shorter $`T_1`$ observed experimentally (18 ps) indicates that additional contributions to the spin relaxation process are present and substantially dominate over BIA. The electron-hole exchange interaction (the Bir-Aronov-Pikus mechanismBAP:1976 (BAP)) is not a favorable candidate to account for the measured spin lifetime in the (110) superlattice due to the small electron-hole overlapWagner:1993 ; Vinattieri:1993  characteristic of these broken-gap InAs/GaSb superlattices (roughly 20% compared to GaAs/AlGaAs quantum wells) and because the timescale of spin decay (18 ps) is much shorter than that found previously for the BAP interaction.Ohno:1999 ; OpticalOrientationBook  The Elliott-Yafet mechanism, which increases in importance with decreasing band gap,OpticalOrientationBook ; Chasalviel:1975  may limit $`T_1`$ in (110) InAs/GaSb superlattices, although a detailed theoretical treatment for this mechanism in III-V semiconductor superlattices is not presently available.

Because of the large spin-orbit interaction and small band gap in InAs/GaSb superlattices, a substantial contribution to the spin decay process could arise from a small degree of structural inversion asymmetry (SIA) associated with compositional mixing at the superlattice interfaces. In order to assess the sensitivity of $`T_1`$ to such effects, the spin lifetime was calculated assuming that the interface on one side of the InAs layers is not compositionally abrupt, corresponding to a *single monolayer* of In0.5Ga0.5As0.5Sb0.5. Such a situation may arise due to differing growth kinetics at the normal and inverted interfaces, an effect which is commonly present in InAs/GaSb superlattices.Feenstra:1994 ; Steinshnider:2000 ; Lew:1997 ; Wang:1995  With the inclusion of the mixed composition monolayer, the calculated spin lifetime in the (110) InAs/GaSb superlattice was reduced by a factor of 20, yielding a value of 29 ps, in line with the experimental value of 18 ps. This dramatic reduction in $`T_1`$ occurs because, although the pseudomagnetic field associated with BIA is primarily in the growth direction for (110)-oriented superlattices, SIA introduces a field component in the plane of the superlattice that serves to relax the growth-direction polarized electron spins injected in these optical experiments. In contrast, for the (001)-oriented superlattice including such compositional mixing has a much smaller influence on the calculated spin lifetimes because, for (001) superlattices, both the psuedomagnetic fields associated with BIADK:1986  and NIA have large in-plane components. Detailed STM characterization of the present (110) superlattice would help in the conclusive identification of the physical process limiting $`T_1`$.

In summary, we observe an enhancement of the electron spin relaxation time in InAs/GaSb superlattices by more than an order of magnitude (25 times) with (110)-oriented growth using polarization sensitive pump probe techniques and a mid-infrared, subpicosecond PPLN OPO. This enhancement is attributed to suppression of the NIA and BIA contributions to the spin decay process in (110) superlattices. Our findings have implications for potential spintronics applications using InAs/GaSb heterostructures.

This research is supported by the DARPA MDA972-01-C-0002, DARPA/ARO DAAD19-01-1-0490 and the Natural Sciences and Engineering Research Council of Canada.