text:
# Diversity of discrete breathers observed in a Josephson ladder

## Abstract

We generate and observe discrete rotobreathers in Josephson junction ladders with open boundaries. Rotobreathers are localized excitations that persist under the action of a spatially uniform force. We find a rich variety of stable dynamic states including pure symmetric, pure asymmetric, and mixed states. The parameter range where the discrete breathers are observed in our experiment is limited by retrapping due to dissipation.

Nonlinearity and lattice discreteness lead to a generic class of excitations that are spatially localized on the scale comparable to the lattice constant. These excitations, also known as discrete breathers, have recently attracted a lot of interest in theory of nonlinear lattices . It is believed that discrete breathers might play an important role in the dynamics of various physical systems consisting of coupled nonlinear oscillators. It has been even said that that for discrete nonlinear systems breathers might be as important as solitons are for the continuous nonlinear media.

There have been several recent experiments that reported on generation and detection of discrete breathers in diverse systems. These are low-dimensional crystals , anti-ferromagnetic materials , coupled optical waveguides , and Josephson junction arrays . By using the method of low temperature scanning laser microscopy, we have recently reported the first direct visualization of discrete breathers .

In this paper we present new measurements of localized modes in Josephson ladders. Using the same method as in our first experiment, we study an even more tightly coupled lattice of Josephson junctions and observe a rich diversity of localized excitations that persist under the action of a spatially uniform force.

A biased Josephson junction behaves very similar to its mechanical analog that is a forced and damped pendulum. An electric bias current flowing across the junction is analogous to a torque applied to the pendulum. The maximum torque that the pendulum can sustain and remain static corresponds to the critical current $`I_\mathrm{c}`$ of the junction. For low damping and bias below $`I_\mathrm{c}`$, the junction allows for two states: the superconducting (static) state and the resistive (rotating) state. The phase difference $`\phi `$ of the macroscopic wave functions of the superconducting islands on both sides of the junction plays the role of an angle coordinate of the pendulum. According to the Josephson relation, a junction in a rotating state generates dc voltage $`V=\frac{1}{2\pi }\mathrm{\Phi }_0\frac{d\phi }{dt}`$, where $`\mathrm{}`$ is the time average. By connecting many Josephson junctions by superconducting leads one gets an array of coupled nonlinear oscillators.

We perform experiments with a particular type of Josephson junction array called Josephson ladder. Theoretical studies  of these systems have predicted the existence of spatially localized excitations called rotobreathers. Rotobreathers are $`2\pi `$-periodic solutions in time that are exponentially localized in space.

Measurements are performed on linear ladders (with open boundaries) consisting of Nb/Al-AlOx/Nb underdamped Josephson tunnel junctions . An optical image Fig. 1(a) shows a linear ladder that is schematically sketched in Fig. 1(b). Each cell contains 4 small Josephson junctions. The size of the hole between the superconducting electrodes which form the cell is about $`3\times 3`$ $`\mu `$m2. The distance between the Josephson junctions is about $`24`$ $`\mu `$m. The number of cells $`N`$ in the ladder is 10. Measurements presented in this paper have been performed at $`5.2`$ K. The bias current $`I_\mathrm{B}`$ was uniformly injected at every node via thin-film resistors $`R_\mathrm{B}=32`$ $`\mathrm{\Omega }`$. Here we define vertical junctions ($`\mathrm{𝐽𝐽}_V`$) as those in the direction of the external bias current, and horizontal junctions ($`\mathrm{𝐽𝐽}_H`$) as those transverse to the bias. The ladder voltage was read across the central vertical junction. The damping coefficient $`\alpha =\sqrt{\mathrm{\Phi }_0/(2\pi I_\mathrm{c}CR_{\mathrm{sg}}^2)}`$ is the same for all junctions as their capacitance $`C`$ and sub-gap resistance $`R_{\mathrm{sg}}`$ scale with the area and $`C_\mathrm{H}/C_\mathrm{V}=R_{\mathrm{sgV}}/R_{\mathrm{sgH}}`$. The damping $`\alpha `$ in the experiment can be controlled by temperature and its typical values vary between $`0.1`$ and $`0.02`$.

There are two types of coupling between cells in a ladder. The first is the inductive coupling between the cells that is expressed by the self-inductance parameter $`\beta _\mathrm{L}=2\pi LI_{\mathrm{cV}}/\mathrm{\Phi }_0`$, where $`L`$ is the self-inductance of the elementary cell. The second is the nonlinear Josephson coupling via horizontal junctions. The ratio of the horizontal and vertical junction areas is called the anisotropy factor and can be expressed in terms of the junction critical currents $`\eta =I_{\mathrm{cH}}/I_{\mathrm{cV}}`$. If this factor is equal to zero, the vertical junctions are decoupled and operate independently one from another. On the other hand, if this factor goes to infinity the ladder behaves like a parallel one-dimensional array and no rotobreather can exist since no magnetic flux can enter through the horizontal junctions. Thus, it is an important challenge to increase the anisotropy as far as possible. In this work we present measurements for the highest anisotropy factor ($`\eta =0.56`$) studied up to now. In contradiction to the existing theoretical prediction  for this anisotropy value, in the studied parameter range with moderate dissipation we find a rich diversity of localized excitations.

To briefly introduce the role of parameters we like to quote the equations of motion for our system (see Ref.  for details):

$`\ddot{\phi }_l^\mathrm{V}+\alpha \dot{\phi }_l^\mathrm{V}+\mathrm{sin}\phi _l^\mathrm{V}=\gamma {\displaystyle \frac{1}{\beta _\mathrm{L}}}(\mathrm{\Delta }\phi _l^\mathrm{V}+\stackrel{~}{\phi }_{l1}^\mathrm{H}\phi _{l1}^\mathrm{H}),`$ (1)
$`\ddot{\phi }_l^\mathrm{H}+\alpha \dot{\phi }_l^\mathrm{H}+\mathrm{sin}\phi _l^\mathrm{H}={\displaystyle \frac{1}{\eta \beta _\mathrm{L}}}(\phi _l^\mathrm{H}\stackrel{~}{\phi }_l^\mathrm{H}+\phi _l^\mathrm{V}),`$ (2)
$`\ddot{\stackrel{~}{\phi }}_l^\mathrm{H}+\alpha \dot{\stackrel{~}{\phi }}_l^\mathrm{H}+\mathrm{sin}\stackrel{~}{\phi }_l^\mathrm{H}={\displaystyle \frac{1}{\eta \beta _\mathrm{L}}}(\phi _l^\mathrm{H}\stackrel{~}{\phi }_l^\mathrm{H}+\phi _l^\mathrm{V}),`$ (3)
where $`\phi _l^\mathrm{V},\phi _l^\mathrm{H},\stackrel{~}{\phi }_l^\mathrm{H}`$ are the phase differences across the $`l`$th vertical junction and its right upper and lower horizontal neighbors, $`\phi _l=\phi _{l+1}\phi _l`$, $`\mathrm{\Delta }\phi _l=\phi _{l+1}+\phi _{l1}2\phi _l`$, and $`\gamma =I_\mathrm{B}/I_{\mathrm{cV}}`$ is the normalized bias current.

In order to generate a discrete breather in a ladder we used the technique described in Ref. . We used two extra bias leads for the middle vertical Josephson junction (cf. Fig. 1(b)) to apply a local current $`I_{\mathrm{local}}>I_{\mathrm{cV}}`$. This current switches the vertical junction in the resistive state. At the same time, being forced by magnetic flux conservation, one or two horizontal junctions on both sides of the vertical junctions also switch to the resistive state. After that $`I_{\mathrm{local}}`$ is reduced and, simultaneously, the uniform bias $`I_\mathrm{B}`$ is tuned up. In the final state, we keep the bias $`I_\mathrm{B}`$ smaller than $`I_{\mathrm{cV}}`$ and have reduced $`I_{\mathrm{local}}`$ to zero. By changing the starting value of $`I_{\mathrm{local}}`$ it is possible to get more than one vertical junction rotating. When trying the slightly different current sweep sequence described in Ref.  we generated mainly states with many rotating vertical junctions.

In all measurements presented below, we have measured the dc voltage across the middle vertical junction as a function of the homogeneous bias current $`I_\mathrm{B}`$. We used the method of low temperature scanning laser microscopy  to obtain electrical images from the dynamic states of the ladder. The laser beam locally heats the sample and changes the dissipation in an area of several micrometers in diameter. If the junction at the heated spot is in the resistive state, a voltage change will be measured. By scanning the laser beam over the whole ladder we can visualize the rotating junctions.

Various measured ladder states are shown in Fig. 2 as 2D gray scale maps. The gray scale corresponds to the measured voltage response during the laser scanning. In Fig. 2(a) we present the simplest of observed states, an asymmetric single-site rotobreather, where one vertical Josephson junction and the upper adjacent horizontal Josephson junctions are in the resistive state. On the right side of the plot we show the corresponding schematic view with arrows marking the rotating junctions. In the case of Fig. 2(b) two vertical junctions and two horizontal junctions are rotating that corresponds to an asymmetric two-site breather. An asymmetric three-site breather and an asymmetric four-site breather are shown in Fig. 2(c) and (d) respectively. Here, the voltage of the whirling horizontal and vertical junctions are equal due to magnetic flux conservation. The magnetic flux enters the ladder through one horizontal junction, goes through the vertical ones and leaves the ladder through another horizontal junction. When a single magnetic flux quantum is passing through a Josephson junction its phase $`\phi `$ is changing by $`2\pi `$. In contrast to our previous measurements of an annular ladder ,here in ladders with open boundaries we observe asymmetric breathers as frequently as symmetric breathers. The $`I_\mathrm{B}`$-$`V`$ characteristics of these localized states is presented in Fig. 3(a). Particular states indicated on the plot are stable along the measured curves. The more junctions are whirling, the higher is the measured resistance.

Another type of discrete breathers observed in our experiment is shown in the Fig. 2(e)–(h). Figure 2(e) illustrates a symmetric single-site breather with one vertical and all four adjacent horizontal junctions are in the resistive state. Though we call this states symmetric, the upper and the lower horizontal junctions may in general have different voltages. The simplest voltage-symmetric case would be when each horizontal junction voltage is equal to half of the voltage of the whirling vertical junction . We also observed a variety of multi-site breathers of this type. A two-site, a three-site, and a four-site breather is shown in Fig. 2 (f), (g), and (h), respectively. The corresponding $`I_\mathrm{B}`$-$`V`$ characteristics over the stability range of symmetric breathers are presented in Fig. 3(b).

As mentioned above, Figs. 3(a) and (b) show the $`I_\mathrm{B}`$-$`V`$ characteristics with $`I_{\mathrm{local}}=0`$ of asymmetric and symmetric rotobreathers, respectively. The voltage $`V`$ is always recorded locally on the middle vertical junction which was initially excited by the local current injection. The vertical line on the left side corresponds to the superconducting (static) state. The rightmost (also the bottom) curve accounts for the spatially uniform whirling state (all vertical junctions rotate synchronously and horizontal junctions are not rotating). An electrical image of this state is shown in Fig. 1(c). The upper branches in Fig. 3 represent various localized states. The uppermost branch corresponds to an single-site breather, the next lower branch to a two-site breather, and so on.

It is easy to show that the voltage $`V_\mathrm{r}`$ at which the vertical and horizontal junctions switch back to the superconducting state is the same. By using $`R_{\mathrm{sgV}}/R_{\mathrm{sgH}}=\eta `$ and $`I_{\mathrm{cH}}/I_{\mathrm{cV}}=\eta `$ we get with $`I_\mathrm{r}I_\mathrm{c}`$ that $`V_\mathrm{r}=R_{\mathrm{sgV}}I_{\mathrm{rV}}=R_{\mathrm{sgH}}I_{\mathrm{rH}}`$. For the uniform state we observed $`V_\mathrm{r}0.6`$ mV. At about the same voltage, the whirling horizontal and vertical junctions for the asymmetric breather are retrapped to the static state as can be seen in Fig. 3(a). The symmetric breathers show different behavior. Assuming $`V_\mathrm{H}V_\mathrm{V}/2`$, we can expect that already at the voltage $`V=2V_\mathrm{r}1.2`$ mV the horizontal junctions should trap into the superconducting state. But if the voltages of the top and bottom horizontal junctions are not equal while decreasing $`I_\mathrm{B}`$, the retrapping current in the junction with the lower voltage is reached earlier and the retrapping occurs at a higher measured voltage.

We found that in order to explain the retrapping current of the breather states the bias resistors $`R_\mathrm{B}`$ have to be taken into account. Initially, these resistors have been designed for providing a uniform bias current distribution in the ladder. Assuming that the Josephson junctions in the superconducting state are just short circuits we get the electric circuits which are shown as insets in Fig. 4. The total current injected in the breather region is calculated by using the resistance $`R`$ of vertical junctions measured in the spatially-uniform state. Thus we derive the retrapping current as

$$(N+1)I_\mathrm{B}=\left[\frac{(NM+1)}{(1+\delta )R_\mathrm{B}}+\frac{(N+1)[M+(2\delta )\eta ]}{MR}\right]V,$$
(4)
where $`M`$ is the number of whirling vertical junctions. The parameter $`\delta `$ is equal to $`0`$ for asymmetric breathers and to $`1`$ for symmetric breathers. The calculated retrapping current is compared with the experiment in Fig. 4. For the asymmetric breathers the agreement is fairly good. We believe that the larger discrepancy for the symmetric breathers is due to the assumption $`V_\mathrm{H}=V_\mathrm{V}/2`$ that we imposed in derivation of the retrapping current.

One interesting feature in our experiment is the dynamical behavior of the local states in the linear ladder which differs from the previous observations in the annular ladder  with the anisotropy parameter $`\eta =0.44`$. In Fig. 5 this feature is illustrated for 3-site breathers. After the creation of a symmetric rotobreather (middle curve) we are lowering the uniform bias current $`I_\mathrm{B}`$ until the retrapping current of one of the horizontal junctions at each side is reached. Here we observe switching to the corresponding asymmetric 3-site rotobreather state. By further lowering of the bias current we get to the point where both the vertical and the horizontal junctions reach their retrapping currents and the whole ladder goes into the superconducting state. If, instead of lowering the bias current, we start increasing it another switching point is observed where the previous symmetric breather state is recovered. In contrast to this behavior, in the annular ladder  the observed lower instability of the symmetric breather led to an increase of voltage. Moreover, this switching turned the number of rotating junctions to increase, with few more vertical junctions starting to rotate.

In addition to the simplest hierarchic states described above we also found a large variety of more complex localized dynamic states. Some examples are presented in Fig. 6. Fig. 6(a)-(c) are sections from the middle region (M) \[cf. Fig. 1(b)\] of the measured linear ladder. Fig. 6(a) shows an asymmetric 6-site breather with top and bottom horizontal junctions whirling on its sides. Fig. 6(b) illustrates a 3-site rotobreather which is symmetric on one side and asymmetric on the other. Another very peculiar localized state is presented in Fig. 6(c). The rightmost vertical junction is rotating with a lower frequency then the other four vertical junctions on the left. This can be understood from the appearance of an extra rotating horizontal junction in the interior of this state. We can interpret it as a 4-site breather coupled to a single-site breather. All these localized states would be topologically forbidden in the case of an annular ladder as the magnetic flux inside the superconducting circuit should not accumulate, but they are not forbidden in a ladder with open boundaries. The last two pictures (d) and (e) are taken near the border of the ladder. We find here truncated asymmetric and symmetric 6-site breathers. The marginal vertical junction does not require any horizontal junctions to rotate.

In summary, we presented observation of a large variety of spatially-localized dynamic rotobreather states in Josephson ladders with open boundaries. These states can be excited by the local current injection and supported by the uniform current bias. We observe both symmetric and asymmetric rotobreather states predicted in Ref. , as well as more complex mixed states. We believe that the region of existence and stability of rotobreathers does sensitively depend on dissipation, discreteness, and anisotropy parameter of the ladder.


format: text

text:
# Properties of cosmologies with dynamical pseudo Nambu–Goldstone bosons

## I Introduction

Scalar fields have played a central role in models of the very early universe for the past 20 years. In the past few years attention has turned to models in which a scalar field plays a dynamical role at late times, rather than simply being frozen in as a static relic vacuum energy. Such models, which have been dubbed “quintessence” models , could in principle provide a dynamical solution to the cosmological constant problem – namely the question of why the magnitude of the vacuum energy at the present epoch is so much smaller than one might naïvely expect from particle physics models such as various supergravity theories. A dynamical ‘solution’ of the cosmological constant problem would amount to a demonstration that a particular dynamical evolution of the scalar quintessence field is a natural consequence of the cosmological field equations without fine-tuning of parameters, given some reasonable physical assumptions about the initial conditions.

The most notable recent observational evidence which has driven the theoretical interest is the measurement of the apparent magnitude-redshift relationship using type Ia supernovae (SNe Ia) . These results have been interpreted, in the context of a cosmological model containing pressureless dust and a cosmological constant, $`\mathrm{\Lambda }`$, as evidence that the universe is undergoing accelerated expansion at the present epoch (see  and references therein). The validity of this conclusion is currently open to some doubt, however. In particular, a recent analysis by Riess et al.  indicates that the sample of type Ia supernovae shows a possible evolution in rise times from moderate ($`z0.3`$) to large ($`z1`$) redshifts. Although the statistical significance of this result has been diminished – from the $`5.8\sigma `$ level  to the $`1.5\sigma `$ level  – upon a more rigorous treatment of the uncertainties in the data , it remains true that while a systematic evolution in the rise times of the supernovae is not conclusively ruled in, neither is it conclusively ruled out.

Given that an evolution in the shape of the light curves of the supernovae measured in their rest frame remains a real possibility, it would not be surprising if the peak luminosity – which is the effective standard candle used – were also to evolve. Riess et al.  conclude that the type Ia supernovae data could conceivably be explained entirely within the context of an open Friedmann-Robertson-Walker universe together with a reasonable astrophysical evolution model, e.g., a consequence of a time variation of the abundances of relevant heavy elements in the environment of the white dwarf supernovae progenitors. Detailed astrophysical modelling – see, e.g.,  – should hopefully eventually resolve the issue, although at this stage the difference between our theoretical understanding and the observations remains quite substantial .

In many recent papers it has been commonly assumed that the dynamical scalar field, $`\varphi `$, should obey an effective equation of state $`P_\varphi w\rho _\varphi `$ with $`1<w<0`$, at the present epoch, in order to obtain a cosmological acceleration, i.e., a negative deceleration parameter $`q_0`$. Indeed, the condition that $`1<w<0`$ is often taken as a defining characteristic of “quintessence” . The broad picture in this cosmological scenario is that the universe is currently in the early stage of an epoch of inflationary expansion. The motivation for this is that one could then hope to have a model cosmology in which observations such as the type Ia supernovae apparent magnitude–redshift relation could be explained by a cosmological acceleration in a similar fashion to models with a cosmological constant, but with the possibility of explaining why the magnitude of the vacuum energy density and the energy density in ordinary pressureless matter, $`\rho _m`$, are comparable at the present epoch – the so-called “cosmic coincidence problem” .

One attractive feature of homogeneous isotropic cosmological models with dynamical scalar fields is that many of them possess “cosmological scaling solutions” , namely solutions which at late times have energy density components which depend on the cosmic scale factor according to $`\rho a^{m_1}`$ and $`\rho _\varphi a^{m_2}`$ simultaneously, and which act as attractors in the phase space. If $`m_2<m_1`$, which is the case, for example, for simple power-law potentials with inverse powers  $`V(\varphi )\varphi ^\alpha `$, or for certain power-law potentials with positive powers  then the scalar field dominates at late times, producing a quintessence-dominated cosmology with accelerated expansion at late times. If $`m_1=m_2`$, which is the case for exponential potentials , then the scaling solutions are “self-tuning”  – i.e., the dynamics of the scalar field follows that of the other dominant energy component, with a dependence $`\rho _\varphi a^4`$ in the radiation–dominated era and a dependence $`\rho _\varphi a^3`$ in the matter–dominated era.

Even if the ultimate late-time behaviour properties of the solutions is not precisely “self-tuning” in the above sense, models such as those with inverse power law potentials can still effectively act as “tracking solutions”, since for a wide range of initial conditions, the solutions rapidly converge to a common, cosmic evolutionary track . Thus there are a number of ways in which one might hope to solve the “cosmic coincidence problem”, though in practice a degree of tuning of the parameters has been necessary in all models studied to date.

In this paper, we consider a form of quintessence, an ultra-light pseudo Nambu-Goldstone-boson (PNGB)  which is still relaxing to its vacuum state. From the viewpoint of quantum field theory PNGB models are the simplest way to have naturally ultra-low mass, spin-$`0`$ particles and hence perhaps the most natural candidate for a presently-existing minimally-coupled scalar field. The effective potential of a PNGB field $`\varphi `$ can be taken to be of the form

$$V(\varphi )=M^4[\mathrm{cos}(\varphi /f)+1],$$
(1)
where the constant term is to ensure that the vacuum energy vanishes at the minimum of the potential. This potential is characterized by two mass scales, a purely spontaneous symmetry breaking scale $`f`$ and an explicit symmetry breaking scale $`M`$.

The effective PNGB mass is $`m_\varphi M^2/f`$. To obtain solutions with $`\mathrm{\Omega }_\varphi 1`$, the energy scales are essentially fixed  to values $`M10^3\text{eV}`$, interestingly close to the neutrino mass scale for the MSW solution to the solar neutrino problem, and $`fm_{PL}10^{19}\text{GeV}`$, the Planck scale. Since these two energy scales have values which are reasonable from the viewpoint of particle physics, one might hope to explain the coincidence that the vacuum energy is dynamically important at the present epoch.

The cosmology of PNGB models has already been extensively studied in the literature . In particular, a number of constraints have been placed on the parameters $`M`$ and $`f`$ by various sets of observational data . Most recently, Frieman and Waga  have set bounds based on the SNe Ia data of Riess et al. (hereafter R98) on the one hand, and gravitational lensing surveys on the other. Comparing these bounds is of interest, since the SnIa data have been interpreted as favouring a cosmological constant, whereas gravitational lensing data has been used to place upper bounds on $`\mathrm{\Lambda }`$ . They therefore provide complementary tests of the parameter spaces of models with a non-trivial vacuum energy.

In this paper it is our intention to critically study these bounds. Firstly, we will consider how the bounds are affected by the initial value of the scalar field at the beginning of the matter–dominated epoch. Secondly, we wish to investigate how such bounds might be affected in the case of the PNGB model if the observed apparent faintness of type Ia supernovae is at least partly due to an intrinsic evolution of the sources over cosmological time scales, which in view of the results of  would appear to be a very real possibility. The reason for focusing on the PNGB models in such an investigation is suggested by the fact that whereas many quintessence models have been singled out in the literature, perhaps somewhat artificially, simply because they have the property of yielding an accelerated expansion, many different possibilities arise in the PNGB case. Indeed, at very late times, the apparent magnitude–redshift relation for PNGB models ultimately coincides with that of the Einstein–de Sitter model, even though the density of ordinary matter can be low in the PNGB cosmologies. The requirement that the faintness of type Ia supernovae is entirely due to their cosmological distances places rather strong restrictions on the values of the parameters $`M`$ and $`f`$ , because it requires us to exist at an epoch of the PNGB cosmologies which is still quite far removed from our ultimate destiny. If these restrictions are relaxed because of evolutionary effects, then it is quite plausible that other regions of the parameter space of the PNGB models become viable alternatives. Since PNGB cosmologies could therefore still solve the “missing energy problem”, even if the evidence for a cosmological acceleration proves to be ephemeral, we believe it is important to investigate this possibility quantitatively.

We will begin the paper with a qualitative analysis of the solutions, to provide some general insights which will help to guide our quantitative discussion. Although these properties are no doubt already known, to the best of our knowledge an analysis of the phase space of the solutions has never been presented in the literature. Having completed this analysis in Sec. II we will go on to discuss a number of issues relating to numerical integration in Sec. III, and relate the properties of the solutions found numerically to the exact analysis of Sec. II. In Sec. IV we present the main analysis of the constraints imposed on the $`(M,f)`$ parameter space, allowing for the possibility of evolution of peak luminosities in the type Ia supernova sources. Bounds from gravitational lensing statistics are updated in Sec. V, and the implications of our results are discussed at greater length in Sec. VI.

## II Phase-Space Analysis

We will begin by performing an analysis of the differential equations governing the cosmological evolution in a manner similar to previous studies in inflationary and quintessential models .

The classical action for gravity coupled to a scalar field $`\varphi `$ has the form

$$S=d^4x\sqrt{g}\left[\left(\frac{R}{2\kappa ^2}\frac{1}{2}g^{\mu \nu }_\mu \varphi _\nu \varphi V(\varphi )\right)+\right],$$
(2)
where $`\kappa `$ is the Planck constant, $`R`$ is the Ricci scalar, $`gdetg_{\mu \nu }`$, and $`L`$ is the Lagrangian density of non-relativistic matter and radiation. For simplicity, we assume $`\varphi `$ is minimally coupled to the curvature, and we work in units in which $`\mathrm{}=c=1`$.

Consider a spatially-flat Friedmann-Robertson-Walker (FRW) universe containing a fluid with barotropic equation of state $`P_\gamma =(\gamma 1)\rho _\gamma `$, where $`\gamma `$ is a constant, $`0\gamma 2`$, such as radiation ($`\gamma =4/3`$) or dust ($`\gamma =1`$). There is a self-interacting scalar field with the PNGB potential energy density (1) evolving in this universe. The total energy density of this homogeneous scalar field is $`\rho _\varphi =\dot{\varphi }^2/2+V(\varphi )`$. The governing equations are given by

$`\dot{H}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{2}}\left(\rho _\gamma +P_\gamma +\dot{\varphi }^2\right),`$ (3)
$`\dot{\rho _\gamma }`$ $`=`$ $`3H(\rho _\gamma +P_\gamma ),`$ (4)
$`\ddot{\varphi }`$ $`=`$ $`3H\dot{\varphi }{\displaystyle \frac{dV}{d\varphi }},`$ (5)
subject to the Friedmann constraint

$$H^2=\frac{\kappa ^2}{3}\left(\rho _\gamma +\frac{1}{2}\dot{\varphi }^2+V\right),$$
(6)
where $`\kappa ^28\pi G`$, $`H=\dot{a}/a`$ is the Hubble parameter, and an overdot denotes ordinary differentiation with respect to time $`t`$.

We may rewrite the Friedmann constraint as

$$\mathrm{\Omega }_\gamma +\mathrm{\Omega }_\varphi =1$$
(7)
where

$`\mathrm{\Omega }_\gamma `$ $`=`$ $`{\displaystyle \frac{\kappa ^2\rho _\gamma }{3H^2}},`$ (8)
$`\mathrm{\Omega }_\varphi `$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{3H^2}}\left({\displaystyle \frac{1}{2}}\dot{\varphi }^2+V\right)`$ (9)
are the ratios of the energy densities of the barotropic matter and the quintessence field as fraction of the critical density respectively.

In contrast to the case of the cosmologies with an exponential potential  where the dynamics can be reduced to a 2-dimensional autonomous phase plane, for the system (3) - (6) the simplest phase space appears to be 3-dimensional in the full four-dimensional phase space.

There are two alternative choices of variables which are useful to describe the dynamics, which we will discuss in turn.

### A Field Variables

The first choice is to simply use the Hubble parameter, H, the scalar field and its first derivative as the elementary variables. These are of course simply the variable $`u`$, $`v`$, $`w`$ of Frieman and Waga  up to an overall scaling. By defining

$$I\frac{\varphi }{f}=\frac{\kappa \varphi }{};J\kappa \dot{\varphi };$$
(10)
we therefore obtain the system

$`\dot{H}`$ $`=`$ $`{\displaystyle \frac{3\gamma }{2}}H^2+{\displaystyle \frac{\gamma m^2}{2}}(\mathrm{cos}I+1){\displaystyle \frac{(2\gamma )}{4}}J^2,`$ (11)
$`\dot{I}`$ $`=`$ $`{\displaystyle \frac{J}{}},`$ (12)
$`\dot{J}`$ $`=`$ $`3HJ+{\displaystyle \frac{m^2}{}}\mathrm{sin}I,`$ (13)
where for notational simplicity we define $`m^2=\kappa ^2M^4`$, and $`=\kappa f`$. so that $``$ is dimensionless, while $`m`$ has dimensions inverse time. The constraint equation becomes

$$\kappa ^2\rho _\gamma =3H^2m^2(\mathrm{cos}I+1)\frac{1}{2}J^2.$$
(14)
From Eq. (4), it follows that $`\dot{\rho _\gamma }=0`$ if $`\rho _\gamma =0`$. Therefore trajectories do not cross the 2-dimensional $`\rho _\gamma =0`$ surface, which is a hyperboloid in the variables $`H`$, $`\mathrm{cos}(I/2)`$, and $`J`$. Physical trajectories with $`\rho _\gamma >0`$ are forced to lie within the volume of the $`H`$, $`I`$, $`J`$ phase space bounded by the $`\rho _\gamma =0`$ surface.

The only critical points of the system (11)-(13) at finite values of $`H`$, $`I`$, $`J`$ occur at

1. $`\mathrm{C}_{1\pm }`$ at $`H=\pm H_1`$, $`I=0`$ mod $`2\pi `$, $`J=0`$; and
2. $`\mathrm{C}_2`$ at $`H=0`$, $`I=\pi `$ mod $`2\pi `$, $`J=0`$,

where

$$H_1\sqrt{\frac{2}{3}}m$$
(15)
Both of these points in fact lie on the $`\rho _\gamma =0`$ surface. Furthermore, this surface intersects the $`H=0`$ plane only at the isolated points $`\mathrm{C}_2`$. The $`H>0`$ and $`H<0`$ subspaces are thus physically distinct, and the $`H<0`$ subspace simply correspond to the time-reversal of the $`H>0`$ subspace. Therefore we can take $`H>0`$ without loss of generality.

The pattern of trajectories close to the $`\rho _\gamma =0`$ surface can be ascertained by continuity to the $`\rho _\gamma =0`$ solutions, even though the latter are not physical. The $`\rho _\gamma =0`$ subspace is obtained, for example, by regarding (14) as a quadratic equation for $`H`$, and using the solution to eliminate $`H`$, thereby obtaining a 2-dimensional system for $`I`$ and $`J`$ given by (12) and (13).

We plot the resulting $`H>0`$ pattern of trajectories in Fig. 1 for values of $`I[0,2\pi )`$. Since the potential, $`V(\varphi )`$, is periodic the same pattern of trajectories repeats itself as we extend $`I`$ to $`\pm \mathrm{}`$, with trajectories crossing from one “cell” to another at the cell boundaries.

The trajectories which occupy the lower half of Fig. 1 are obtained from those in the upper half by the symmetry $`I2\pi I`$, $`JJ`$ of the differential equations (11) - (14). Physically this simply corresponds to the scalar field rolling from the maximum to the right of the minimum as opposed to the one to the left.

An analysis of small perturbations about the critical points $`\mathrm{C}_1`$ and $`\mathrm{C}_2`$ yield eigenvalues

1. $`\lambda =\sqrt{6}m\gamma ,\frac{m}{2}\left(\sqrt{6}\pm \sqrt{6+4/^2}\right)`$ at $`\mathrm{C}_{1+}`$,
2. $`\lambda =0,\pm \frac{im}{}`$ at $`\mathrm{C}_2`$.

Thus $`\mathrm{C}_{1+}`$ attracts a two-dimensional bunch of trajectories but is a saddle point with respect to trajectories lying in the $`\rho _\gamma =0`$ surface, as is evident from Fig. 1. The 2-dimensional bunch of trajectories which approach $`\mathrm{C}_{1+}`$ are found to correspond to inflationary solution with $`a\mathrm{exp}(\sqrt{2/3}mt)`$ as $`t\mathrm{}`$ and $`\varphi `$ const $`=2n\pi `$, $`n`$ Z. The possible role of scalar fields with PNGB potentials in driving an inflationary expansion of the early universe has been discussed in .

The point $`\mathrm{C}_2`$ is a degenerate case, in particular with regard to perturbations orthogonal to the $`\rho _\gamma =0`$ surface, (i.e. into the surface $`\rho _\gamma >0`$ region), for which the eigenvalue is zero. It is a centre with respect to the trajectories lying in the $`\rho _\gamma =0`$ surface, and when perturbations of higher order are considered it becomes a stable spiral point in the $`\rho _\gamma =0`$ surface as can be seen in Fig. 1. Since there is a degeneracy, however, an alternative choice of phase space variables is desirable. We will defer a discussion of the late time behaviour of the solution near $`\mathrm{C}_2`$ to Sec. II B.

The points $`\mathrm{C}_{1\pm }`$ correspond to models with a scalar field sitting at the maximum of the potential, whereas $`\mathrm{C}_2`$ corresponds to the scalar field sitting at the bottom of the potential well. The separatrices in Fig. 1 which join $`\mathrm{C}_{1\pm }`$ to $`\mathrm{C}_2`$ correspond to the field rolling from the maximum to the minimum. It would appear from Fig. 1 trajectories which spiral into $`\mathrm{C}_2`$ become arbitrarily close to the separatrix at late times.

The separatrices which join the points $`\mathrm{C}_{1\pm }`$ to points at infinity correspond to solutions for which the scalar field reaches the top of the potential hill as $`t\pm \mathrm{}`$, (e.g., the rightmost trajectory in Fig. 1). Finally, there are also straight line separatrices parallel to the $`H`$-axis at each of the points $`\mathrm{C}_{1\pm }`$, extending from $`H=\pm H_1`$ to infinity, which represent solutions with a static scalar field sitting on top of the potential hill.

To examine the critical points at infinity it is convenient to transform to spherical polar coordinates $`r`$, $`\theta `$, and $`\varphi `$ by defining

$`H`$ $`=`$ $`r\mathrm{cos}\theta ,`$ (16)
$`I`$ $`=`$ $`r\mathrm{sin}\theta \mathrm{sin}\varphi ,`$ (17)
$`J`$ $`=`$ $`r\mathrm{sin}\theta \mathrm{cos}\varphi ,`$ (18)
and to bring the sphere at infinity to a finite distance from the origin by the transformation $`r=\rho /(1\rho )`$, $`0\rho 1`$ .

Although the trajectories on the sphere at infinity do not represent physical cosmologies, it is useful to plot them since the form of the trajectories which lie just within the sphere will be similar. On the sphere $`\rho =1`$ we find

$`{\displaystyle \frac{d\theta }{d\xi }}`$ $`=`$ $`\mathrm{sin}\theta \mathrm{cos}^2\theta \left[3\mathrm{sin}^2\varphi +{\displaystyle \frac{2\gamma }{4}}(\mathrm{tan}^2\theta \mathrm{cos}^2\varphi 6)\right],`$ (19)
$`{\displaystyle \frac{d\varphi }{d\xi }}`$ $`=`$ $`{\displaystyle \frac{3}{2}}\mathrm{cos}\theta \mathrm{sin}2\varphi ,`$ (21)
where $`\xi `$ is a new time coordinate defined by $`d\xi =rdt`$. The resulting integral curves are plotted in Fig. 2. By (14) the projection of the physical region $`\rho _\gamma >0`$ onto the sphere at infinity leads to the condition

$$\mathrm{cot}^2\theta >\frac{1}{6}\mathrm{cos}^2\varphi .$$
(22)
Values of $`\theta `$ and $`\varphi `$ which violate this inequality lie in the shaded region.

The critical points on the sphere at infinity are

1. $`\mathrm{A}_{1\pm }^{\mathrm{}}`$, $`\mathrm{A}_{2\pm }^{\mathrm{}}`$: four points at

$`(\theta ,\varphi )\{(\pm \mathrm{tan}^1\sqrt{6},0),(\pm \mathrm{tan}^1\sqrt{6},\pi )\}`$ or

$`H=\pm \mathrm{}`$, $`I/H=0`$, and $`J/H=\pm \sqrt{6}`$.
2. $`\mathrm{B}_\pm ^{\mathrm{}}`$: two points at

$`H=\pm \mathrm{}`$, $`I/H=0`$, and $`J/H=0`$.

Since the projection onto spherical polar coordinates (16) - (18) is degenerate at the north and south poles $`\theta =0,\pi `$, these points are excluded from the chart $`(\theta ,\varphi )`$ but can be included using an alternative hemispherical projection (see Fig. 3).
3. $`\mathrm{C}_{1,2}^{\mathrm{}}`$: two points at

$`(\theta ,\varphi )\{(\pi /2,\pi /2),(\pi /2,3\pi /2)\}`$ or

$`I=\pm \mathrm{}`$, $`H/I=0`$, and $`J/I=0`$.

An analysis of small perturbations shows that the points $`\mathrm{A}_{1+,2+}^{\mathrm{}}`$ are repellors in all directions of the phase space (cf., Table I), while $`\mathrm{A}_{1,2}^{\mathrm{}}`$ are attractors. This therefore represents the most ‘typical’ early behaviour of solutions. For $`H>0`$, $`\mathrm{A}_{1+,2+}^{\mathrm{}}`$ correspond to the limit $`t0`$. We find that $`H\frac{1}{3t}`$ or $`at^{1/3}`$, while $`\kappa \varphi \pm \sqrt{2/3}\mathrm{ln}t`$, for these solutions. The points $`\mathrm{A}_{1,2}^{\mathrm{}}`$ with $`H<0`$ represent the time-reversed solutions.

The point $`\mathrm{B}_+^{\mathrm{}}`$ ($`\mathrm{B}_{}^{\mathrm{}}`$) repels (attracts) a 2-dimensional bunch of trajectories travelling to (from) finite values of $`H`$, $`I`$, $`J`$, but is a saddle point with respect to directions on the sphere at infinity. The points are found to correspond to $`t0`$ with $`H\frac{2}{3\gamma t}`$ or $`at^{2/3\gamma }`$ while $`\kappa \varphi t^n`$, $`n>0`$.

The points $`\mathrm{C}_{1,2}^{\mathrm{}}`$ are the projection of the points $`\mathrm{C}_{1\pm }`$ and $`\mathrm{C}_2`$ into the sphere at infinity. The degenerate eigenvalues simply reflect the degeneracy of the projection.

$`\mathrm{B}_+^{\mathrm{}}`$ acts as a repellor for trajectories with $`\dot{\varphi }0`$, $`\varphi `$ const as $`t0`$. As shown in Fig. 3, trajectories are driven towards $`\mathrm{B}_+^{\mathrm{}}`$ before they reach $`\mathrm{C}_1^{\mathrm{}}`$. This is consistent with the property that when $`H`$ is large ($`3Hm_\varphi `$), the field evolution is over-damped by the expansion, and the field is effectively frozen to its initial value ($`\dot{\varphi }0`$).

### B Energy Density Variables

In view of the eigenvalue degeneracy encountered above, we can alternatively choose to represent the system by the variables $`H`$, $`x`$, and $`y`$, where

$`x`$ $``$ $`{\displaystyle \frac{\kappa \dot{\varphi }}{\sqrt{6}H}}={\displaystyle \frac{J}{\sqrt{6}H}},`$ (23)
$`y`$ $``$ $`{\displaystyle \frac{\kappa \sqrt{V}}{\sqrt{3}H}}={\displaystyle \frac{\sqrt{2}m\mathrm{cos}(I/2)}{\sqrt{3}H}}.`$ (24)
These are the same variables used by Copeland, Liddle and Wands  for the model with an exponential potential. As above, we will consider $`H>0`$ only. The field equations then take the form

$`\dot{H}`$ $`=`$ $`{\displaystyle \frac{3}{2}}H^2\mu ,`$ (25)
$`\dot{x}`$ $`=`$ $`\pm {\displaystyle \frac{my}{}}\sqrt{1{\displaystyle \frac{3y^2H^2}{2m^2}}}+{\displaystyle \frac{3}{2}}Hx(\mu 2),`$ (26)
$`\dot{y}`$ $`=`$ $`{\displaystyle \frac{mx}{}}\sqrt{1{\displaystyle \frac{3y^2H^2}{2m^2}}}+{\displaystyle \frac{3}{2}}Hy\mu ,`$ (27)
where

$$\mu (x,y)\gamma (1y^2)+(2\gamma )x^2.$$
(28)
We note that in these variables

$$x^2+y^2=\mathrm{\Omega }_\varphi $$
(29)
which is why we have adopted the terminology “energy density variables”. The physical region of the phase space will be constrained to lie within the cylinder $`x^2+y^21`$ since $`\mathrm{\Omega }_\varphi 1`$. In the case of the exponential potential analysed in Ref. , one of the differential equations decoupled, and the dynamics was effectively described by a phase plane with trajectories bounded by the circle $`x^2+y^2=1`$. In the present case, however, no such simplification arises.

The physical region of phase space is further restricted by the requirement that $`y^2\frac{2}{3}m^2/H^2H_{1}^{}{}_{}{}^{2}/H^2`$, which is equivalent to $`\mathrm{cos}^2\left(I/2\right)1`$ in terms of the field variables. For values of $`H>H_1`$ each $`H=\text{const}`$ slice of the cylinder $`x^2+y^21`$ is cut off in the $`y`$-direction above and below the $`y=\pm H_1/H`$ lines. Thus the “fundamental cell” of the phase space can be considered to be a cylinder for $`0HH_1`$, capped by a horn for $`H>H_1`$, which tapers off to a line segment $`1x1`$ on the $`x`$-axis as $`H\mathrm{}`$ (see Fig. 4). In fact, the phase space consists of an infinite number of copies of the fundamental cell of Fig. 4 as a result of the periodic structure of the potential. These cells, $`𝒞_n`$, can be labelled by an integer, $`n`$, with the variable $`I`$ lying in the range $`2n\pi I<2(n+1)\pi `$ for each $`n`$. For cells with even $`n`$ the dynamics is described by (25)–(28) with the upper sign in (26) and (27), while for odd $`n`$ one must take the lower sign.

Within the horn portion of each phase space cell the motion of most trajectories is roughly circular in $`H=\text{const}`$ slices, in a clockwise sense in even cells, $`𝒞_{2n}`$, and anti-clockwise in odd cells $`𝒞_{2n+1}`$. However, trajectories can cross from one cell to another along the $`y=\pm H_1/H`$ boundaries of their horns, which correspond to the surfaces $`I=2n\pi `$ in terms of the field variables. For even cells, $`𝒞_{2n}`$, trajectories join the cell $`𝒞_{2n}`$ along the $`y=H_1/H`$ surface and the cell $`𝒞_{2n+1}`$ along the $`y=H_1/H`$ surface. Below the $`H=H_1`$ plane solutions cannot cross from one cell to another, but remain confined within the cylinder $`x^2+y^21`$.

When $`H=0`$ we see from (25) - (28) that $`\dot{x}=my/`$ and $`\dot{y}=mx/`$, so that trajectories which lie in the $`H=0`$ surface are purely concentric circles. Since $`\dot{a}=0`$ in the $`H=0`$ plane, these do not represent physically interesting cosmologies, but by continuity the behaviour of the trajectories just above the plane will be of a spiral nature.

The origin $`H=x=y=0`$ is in fact the critical point corresponding to $`\mathrm{C}_2`$. The nature of the critical point is altered by the change of variables, however. In particular, whereas the eigenvalues for linear perturbations are unchanged, when higher order corrections are considered the point is no longer always an asymptotically stable spiral as was the case in Fig. 1.

Asymptotically as $`t\mathrm{}`$ we have

$`x=A(t)\mathrm{sin}{\displaystyle \frac{mt}{}},`$ (30)
$`y=A(t)\mathrm{cos}{\displaystyle \frac{mt}{}},`$ (31)
where the amplitude is governed by the equation

$$\frac{d}{dt}(A^2)=3HA^2\left(\gamma 2\mathrm{sin}^2\frac{mt}{}\right)(1A^2).$$
(32)
The nature of $`\mathrm{C}_2`$ is now found to depend on $`\gamma `$:

1. If $`\gamma <1`$ we find that over a cycle the average value of the r.h.s. of (32) is negative and $`A^2`$ decreases so that $`\mathrm{C}_2`$ is an asymptotically stable spiral. Furthermore to leading order $`H\frac{2}{3\gamma t}`$ as $`t\mathrm{}`$ or $`at^{2/3\gamma }`$ and $`A(t)`$ is given by

$$A(t)=Bt^{(\gamma 1)/\gamma }\mathrm{exp}\left[\frac{1}{\gamma }_t^{\mathrm{}}\frac{dt^{}}{t^{}}\mathrm{cos}\left(\frac{mt^{}}{}\right)\right]$$
(33)
with $`B`$ constant. The late-time attractor has $`\mathrm{\Omega }_\varphi =0`$ and $`\mathrm{\Omega }_\gamma =1`$
2. If $`\gamma >1`$ then over a cycle the average value of the r.h.s. of (32) is positive and $`A^2`$ increases until it reaches a limit cycle $`A^2=1`$, i.e. $`x^2+y^2=1`$ or $`\mathrm{\Omega }_\varphi =1`$, $`\mathrm{\Omega }_\gamma =0`$. In this case

$$A(t)=1\overline{B}t^{2(1\gamma )}\mathrm{exp}\left[2_t^{\mathrm{}}\frac{dt^{}}{t^{}}\mathrm{cos}\left(\frac{mt^{}}{}\right)\right].$$
(34)
3. If $`\gamma =1`$, which corresponds physically to an ordinary matter-dominated universe, then an intermediate situation obtains. Essentially any of the concentric circles in the $`H=0`$ plane of Fig. 4 can be approached asymptotically giving a universe for which $`\mathrm{\Omega }_\varphi \alpha _1`$ and $`\mathrm{\Omega }_\varphi 1\alpha _1`$ where $`\alpha _1`$ is a constant in the range $`0<\alpha _1<1`$, which depends on the initial conditions and the parameters $`m`$ and $``$.

We observe that for all values of $`\gamma `$,

$$\rho _\varphi \frac{1}{a^3}$$
(35)
at late times. Since $`\rho _\gamma a^{3\gamma }`$, the three different late time behaviours can thus be understood as a consequence of the scalar field either decreasing more rapidly than the barotropic fluid ($`\gamma <1`$), less rapidly ($`\gamma >1`$), or at the same rate ($`\gamma =1`$). The scalar field thus eventually dominates if $`\gamma >1`$, while the barotropic fluid dominates if $`\gamma <1`$. In the interesting critical case of dust filled models ($`\gamma =1`$) both the scalar field and ordinary matter are of cosmologically significant density at late times.

One simplification that is often made in studying quintessence models is to assume that at late times the quintessence field obeys an equation of state

$$P_\varphi =\left(\gamma _\varphi 1\right)\rho _\varphi $$
(36)
with $`\gamma _\varphi `$ effectively constant. While such an assumption is justified in the case of models with a slowly varying scalar field, it does not apply in the present case. In particular, since the effective barotropic index of the scalar field is

$$\gamma _\varphi =\frac{2x^2}{x^2+y^2},$$
(37)
we see that $`\gamma _\varphi \mathrm{sin}^2\left(mt/\right)`$ at late times, so that it remains truly variable, varying from $`0`$ to $`2`$ over each cycle. It follows from (26)–(29) that the scalar energy density parameter obeys the equation

$$\dot{\mathrm{\Omega }}_\varphi =3H\mathrm{\Omega }_\varphi \left(1\mathrm{\Omega }_\varphi \right)\left(\gamma \gamma _\varphi \right)$$
(38)
so that $`\dot{\mathrm{\Omega }}_\varphi 0`$ as $`t\mathrm{}`$, which accords with the late time properties of the solutions observed above.

## III Numerical Integration

We will now consider the extent to which models based on the PNGB potential are constrained by the latest observational evidence. This work extends the studies previously undertaken by various authors . In the most recent analysis, Frieman and Waga have compared the constraints imposed by the high redshift supernovae luminosity distance on the one hand and gravitational lensing bounds on the other. The two measures provide tests which are potentially in opposition. Here we will perform a similar analysis for the PNGB models, but also taking into account the possibility of luminosity evolution which has not been considered in previous studies . We are therefore considering the model of the previous section with $`\gamma =1`$.

To proceed it is necessary to integrate the equations numerically. To do this we introduce the dimensionless variables

$`u`$ $`=`$ $`{\displaystyle \frac{J}{H_0}}={\displaystyle \frac{\kappa \dot{\varphi }}{H_0}},`$ (39)
$`v`$ $`=`$ $`\mathrm{\Omega }_{m0}^{1/3}(1+z),`$ (40)
$`w`$ $`=`$ $`I={\displaystyle \frac{\kappa \varphi }{}}.`$ (41)
where $`\mathrm{\Omega }_{m0}`$ is the fractional energy density of matter at the present epoch, $`t_0`$. The dynamical system then becomes

$`u^{}`$ $`=`$ $`3{\displaystyle \frac{H}{H_0}}u+{\displaystyle \frac{m^2}{H_0^2}}\mathrm{sin}w,`$ (42)
$`v^{}`$ $`=`$ $`v{\displaystyle \frac{H}{H_0}},`$ (43)
$`w^{}`$ $`=`$ $`{\displaystyle \frac{u}{}},`$ (44)
where the Hubble parameter is defined implicitly according to

$$\frac{H}{H_0}\left[v^3+\frac{1}{6}u^2+\frac{m^2}{3H_0^2}(\mathrm{cos}w+1)\right]^{\frac{1}{2}},$$
(45)
and prime denotes a derivative with respect to the dimensionless time parameter $`\tau H_0t`$. Only the variable $`v`$ differs from those used by Frieman and Waga . Our reason for making the choice (40) is that it allows us to integrate the Friedman equation directly rather than a second order equation (11) which follows from the other equations by virtue of the Bianchi identity. This may possibly lead to better numerical stability since it is not necessary to implement the Friedmann constraint separately.

We begin the integration at initial values of $`u`$, $`v`$, and $`w`$ chosen to correspond to initial conditions expected in the early matter-dominated era. The integration proceeds then until the r.h.s. of (45) is equal to 1, thereby determining the value of the present epoch, $`t_0`$, to be the time at which $`H=H_0`$. We are then also able to determine $`\mathrm{\Omega }_{m0}`$, since according to (40)

$$\mathrm{\Omega }_{m0}=v^3(t_0).$$
(46)
The choice of appropriate initial conditions has been previously discussed . In particular, since the Hubble parameter is large at early times, it effectively acts as a damping term in Eq. (3), driving the scalar field to a state with $`\dot{\varphi }=0`$ initially, i.e. $`u=0`$. We take $`v=1101`$ initially, which in view of relation (40) and the fact that $`\mathrm{\Omega }_{m0}0.11`$ corresponds to the early matter dominated era $`1100z\genfrac{}{}{0pt}{}{<}{}3000`$. Results of the integration do not change significantly if $`v`$ is altered to values within the same order of magnitude.

The initial value of the scalar field variable, $`w_iw(t_i)`$, can lead to some variability in predicted cosmological parameters at the present epoch. A few different values of $`w(t_i)`$ have been considered by different authors . However, the only systematic studies of bounds in the $`M`$,$`f`$ parameter space have been performed  for one particular initial value, $`w_i=1.5`$. One must bear in mind that such bounds are also dependent on $`w_i`$, the value of which is not greatly restricted. Given that we are starting with $`u(t_i)0`$, so that the kinetic energy of the scalar field is initially negligible, the only physical restriction on the value of $`w_i`$ comes from the requirement that the scalar field should be sufficiently far from the minimum of the potential, $`V(\varphi )`$, that $`\mathrm{\Omega }_\varphi (t_i)`$ is small. Thus will ensure that $`w_i`$ is consistent with a scalar field that has emerged from the radiation dominated era with $`\mathrm{\Omega }_\varphi `$ sufficiently small that is consistent with bounds set by primordial nucleosynthesis, and by structure formation models. This still leaves considerable latitude for the choice of $`w_i`$, however.

In Figs. 5 and 6 we display contour plots of $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ in the $`M`$,$`f`$ parameter space for two values $`w_i=1.5`$ and $`w_i=0.5`$. Similar figures have been given by Frieman and Waga  in the $`w_i=1.5`$ case, although our resolution is somewhat better. As $`w_i`$ decreases the contour plots do not change significantly in terms of their overall features, but contours with equivalent values shift to lower values of the $`f`$ parameter. For example, for large values of $`M`$ the $`\mathrm{\Omega }_{\varphi 0}=0.7`$ contour lies at a value $`f2.05\times 10^{18}\text{GeV}`$ if $`w_i=1.5`$, while the same contour lies at $`f0.94\times 10^{18}\text{GeV}`$ if $`w_i=0.5`$.

The other principal feature of the plots 5 and 6, which was not commented on in Ref. , is the wave-like properties of the contours at larger values of $`M`$. These features can be readily understood by considering the corresponding plots of the deceleration parameter, $`q_0`$, which is defined in terms of $`u(t_0)`$, $`v(t_0)`$, and $`w(t_0)`$ by

$`q_0`$ $`=`$ $`{\displaystyle \frac{3}{2}}v^3(t_0)+{\displaystyle \frac{1}{2}}u^2(t_0)1,`$ (47)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(1+3\mathrm{\Omega }_{\varphi 0}6y_{0}^{}{}_{}{}^{2}\right).`$ (48)
We display contour plots of $`q_0`$ in the $`M`$,$`f`$ parameter space in Fig. 7. Essentially, as $`M`$ increases for roughly fixed $`f`$ for sufficiently large $`M`$, the value of $`q_0`$ oscillates over negative and positive values, from a minimum of $`q_0=\frac{1}{2}\left(13\mathrm{\Omega }_{\varphi 0}\right)`$ to a maximum $`q_0=\frac{1}{2}\left(1+3\mathrm{\Omega }_{\varphi 0}\right)`$ about a mean of $`q_0=0.5`$. This corresponds to the scalar field, $`\varphi `$, having undergone more and more oscillations by the time of the present epoch. The minimum value of $`q_0`$ is attained when $`\dot{\varphi }=0`$ instantaneously, while the maximum value of $`q_0`$ is attained when $`\varphi `$ is instantaneously passing through the minimum of its potential.

For smaller values of $`M`$ to the left of the plots the scalar field has only relatively recently become dynamical, whereas for larger values of $`M`$, the scalar can already have undergone several oscillations by the time of the present epoch, particularly if $`f`$ is small. This variation can be understood in terms of the asymptotic period of oscillation of solutions which approach $`\mathrm{C}_1`$, which by Eqs. (30), (31) is

$$t_a=2\pi /m=2\pi f/M^2.$$
(49)
The period $`t_a`$ is shorter for larger $`M`$, or for smaller $`f`$. Since the final $`f`$ values plotted in the $`w_i=0.5`$ case are a factor of 2 smaller than the $`w_i=1.5`$ case, this also explains why points with the same value of $`M`$ have undergone more oscillations up to the present epoch for the smaller value of $`w_i`$.

The value of $`H_0t_0`$ oscillates as $`M`$ increases for roughly fixed $`f`$, according to whether the universe has been accelerating or decelerating in the most recent past, with more rapid variation for parameter values with shorter asymptotic periods, $`t_a`$.

The wiggles in the $`\mathrm{\Omega }_{\varphi 0}`$ contours are a residual effect of the oscillation of the scalar field as $`\mathrm{\Omega }_\varphi `$ settles down to a constant value according to (38). The variation in the value of $`\mathrm{\Omega }_{\varphi 0}`$ with the value of $`w_i`$ can be understood from the fact that for smaller values of $`w_i`$ the scalar field begins its evolution in the matter-dominated era closer to the critical points, $`\mathrm{C}_{1+}`$, corresponding to the maximum of the potential, $`V(\varphi )`$. For fixed $`M`$ and $`f`$ the period of quasi-inflationary expansion is therefore longer, and the present value of $`\mathrm{\Omega }_{\varphi 0}`$ larger.

## IV Constraints from High-redshift type IA Supernovae

Empirical calibration of the light curve - luminosity relationship of type Ia supernovae provides absolute magnitudes that can be used as distance indicators. Since the luminosity distance of a light source, $`d_L`$, is defined by

$$\frac{H_0d_L}{c}=(z+1)_\tau ^{\tau _0}(z+1)𝑑\tau ,$$
(50)
with $`\tau =H_0t`$, it is convenient to define an additional variable

$$r=_{\tau _0}^\tau v𝑑\tau ,$$
(51)
which is proportional to the comoving coordinate distance

$$r_{\text{comoving}}=\frac{r}{H_0a(t_0)v(t_0)}.$$
(52)
We then adjoin a differential equation

$$r^{}=v$$
(53)
to the differential equations (39)-(41) when performing the numerical integration. In terms of $`r`$ and $`v`$, the luminosity distance is then determined according to

$$\frac{H_0d_L}{c}=\frac{v}{\mathrm{\Omega }_{m0}^{2/3}}(rr(t_0)),$$
(54)
and can be used in the appropriate distance modulus to compare with the supernovae data.

Frieman and Waga  have recently considered constraints on PNGB models from supernovae data (without source evolution) using the data set of Riess et al. . We will perform a similar analysis, but we will make use of the largest available data set, namely the 60 supernovae published by Perlmutter et al.  (hereafter P98). Of these, 18 low redshift SNe Ia were discovered and measured in the Calán-Tololo survey , and the Supernova Cosmology Project discovered 42 new SNe Ia at redshifts between 0.17 and 0.83. The peak magnitudes of the supernovae are corrected using the “stretch factor” light curve fitting method . This method is based on fitting a time-stretched version of a single standard template to the observed light curve. The stretch factor is then used to estimate the absolute magnitude.

### A Models without evolution

We use the stretch-luminosity corrected effective $`B`$-band peak magnitude in Table 1 of P98 as the absolute magnitude and denote it as $`m_B^{\text{eff}}`$. The distant modulus of a supernova is defined as $`m_B^{\text{eff}}+_0`$, where $`_0`$ is the fiducial absolute magnitude which has not been given in the literatures. A fitting method to obtain $`_0`$, using only the 18 Calán/Tololo supernovae, can be found in . Perlmutter et al. calculate the confidence regions by simply integrating over $`_0`$. In this paper, we will perform an analytic marginalization over $`_0`$ and obtain the marginal likelihood.

Similarly to the $`\chi ^2`$ statistic used by Riess et al. , we define the quadratic form

$$\chi ^2(M,f)=\underset{i=1}{\overset{60}{}}\frac{(m_{B,i}^{\text{eff}}+_0\mu _i)^2}{\sigma _i^2}$$
(55)
where $`\mu _i(z_i;H_0,M,f)`$ is the predicted distance modulus for model parameters $`(M,f)`$, and

$$\sigma _i^2=\sigma _{m_{B,i}^{\text{eff}}}^2+\left(\frac{5\mathrm{log}\mathrm{e}}{z_i}\sigma _{z,i}\right)^2.$$
(56)
The predicted distance modulus is

$$\mu _i=5\mathrm{log}d_L(z_i)+25,$$
(57)
if the luminosity distance, $`d_L`$, as defined by (50) or (54), is given in units of Megaparsecs.

In order to perform analytic marginalization over $`H_0`$ as well as over $`_0`$ we separate out the $`H_0`$ and $`_0`$ dependence from $`\mu _i`$ into a quantity, $`\nu `$, which we define by

$$\mu _i_0=g_i+\nu $$
(58)
where $`g_i(z_i;M,f)`$ depends implicitly only on $`M`$ and $`f`$. We then follow the statistical procedures adopted by Drell, Loredo and Wasserman  and marginalize over $`\nu `$ using a flat prior that is bounded over some range $`\mathrm{\Delta }\nu `$.

The marginal likelihood is

$$(M,f)=\frac{1}{\mathrm{\Delta }\nu }𝑑\nu \mathrm{e}^{\chi ^2/2}=\frac{s\sqrt{2\pi }}{\mathrm{\Delta }\nu }\mathrm{e}^{Q/2}$$
(59)
where

$$Q(M,f)=\underset{i=1}{\overset{60}{}}\frac{(m_{B,i}^{\text{eff}}g_i\widehat{\nu })^2}{\sigma _i^2}$$
(60)
and $`\widehat{\nu }(M,f)`$ is the best-fit value of $`\nu `$ given $`M`$ and $`f`$, with conditional uncertainty $`s`$. We identify $`1/s^2`$ as the coefficient of $`\nu ^2`$ in (55),

$$\frac{1}{s^2}=\underset{i=1}{\overset{60}{}}\frac{1}{\sigma _i^2},$$
(61)
and $`\widehat{\nu }`$ as $`s^2`$ times the coefficient of $`\nu `$, viz.

$$\widehat{\nu }(M,f)=s^2\underset{i=1}{\overset{60}{}}\frac{m_{B,i}^{\text{eff}}g_i}{\sigma _i^2}.$$
(62)
The quadratic form (60) is what one would obtain by calculating the “maximum likelihood” for $`M`$ and $`f`$. Since $`s`$ is independent of $`M`$ and $`f`$, it follows from Eq. (59) that the marginal likelihood is proportional to the maximum likelihood in this case.

Fig. 8 shows the 68.3% and 95.4% joint credible regions for $`M`$ and $`f`$, and is a direct analogue of Fig. 1 of Ref. , where a similar analysis was performed on 37 supernovae given in R98. The position of the region of parameters which are included at both the 68.3% and 95.4% confidence levels is broadly similar to that obtained from the R98 data . Although Frieman and Waga  did not include the parameter region $`M>0.004h\text{eV}`$, we have redone their analysis on the R98 data and find that the parameter values to the right of Fig. 8 which are admitted at the 95.4% level but excluded at the 68.3% level for the P98 dataset, (labelled region II in Fig. 8), are in fact excluded also at the 95.4% confidence level if the 37 supernovae of the R98 dataset are used. It is possible that this discrepancy has its origin in the different techniques used by Riess et al.  to determine the distance moduli. Possible systematic discrepancies in the “stretch factor” method of P98 versus the “multi-colour light curve” and “template fitting methods” of R98 have been discussed in some detail in Ref. .

The importance of the $`2\sigma `$ included parameter region to right of Fig. 8 diminishes, however, if one compares it with Figs. 5 and 6, since it largely corresponds to parameter values with $`\mathrm{\Omega }_{\varphi 0}\genfrac{}{}{0pt}{}{>}{}0.9`$, which can be discounted by dynamical measurements of $`\mathrm{\Omega }_{m0}`$, where $`\mathrm{\Omega }_{m0}+\mathrm{\Omega }_{\varphi 0}=1`$. Furthermore, the few allowed values below the $`\mathrm{\Omega }_{\varphi 0}=0.9`$ contour in this part of the parameter space have unacceptably small values for the age of the Universe, $`H_0t_0`$.

The parameter region $`0.002h\text{eV}<M<0.003h\text{eV}`$ which from Fig. 8 is admitted at both the 68.3% and 95.4% levels, (labelled region I), by contrast, corresponds to acceptable values of both $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$. Comparing with Fig. 7, we see that this region has $`0.1\genfrac{}{}{0pt}{}{<}{}q_0\genfrac{}{}{0pt}{}{<}{}0.6`$, corresponding to a Universe with a scalar field still in an early stage of rolling down the potential $`V(\varphi )`$. The label I is thus indicative of the fact that the scalar field is rolling down the potential for the first time (from left to right), while in region II the scalar field is rolling down the potential for the second time (from right to left). In region II $`q_0`$ is positive – however, it corresponds to parameter values for which there would have been a cosmological acceleration at modest redshifts in the past, e.g., at $`z0.2`$, well within the range of the current supernovae dataset.

Since conclusions regarding statistically preferred regions of the parameter space can change if evolution of the sources occurs, we think it is important that this possibility is also examined, as we will now do.

### B Models with evolution

In view of the recent results of Refs. , supernovae Ia may exhibit some evolutionary behaviour, at least as far as their rise times are concerned. A possible evolution in the peak luminosity is therefore a possibility which must be seriously investigated.

In the absence of a detailed physical model to explain precisely how the source peak luminosities vary with redshift, one approach is to assume some particular empirical form for the source evolution, and to examine the consequences. Such an analysis has been recently preformed by Drell, Loredo and Wasserman  in the case of Friedmann-Lemaître models with constant vacuum energy. We will undertake an equivalent analysis for the case of PNGB quintessential cosmologies.

Following Drell, Loredo and Wasserman  we will assume that the intrinsic luminosities of SNe Ia scale as a power of $`1+z`$ as a result of evolution. This model introduces a continuous magnitude shift of the form $`\beta \mathrm{ln}(1+z)`$ to the SNe Ia sample. Eq. (55) then becomes

$$\chi ^2(M,f)=\underset{i=1}{\overset{60}{}}\frac{(m_{B,i}^{\text{eff}}g_i\nu \beta \mathrm{ln}(1+z_i))^2}{\sigma _i^2}.$$
(63)
The parameter $`\beta `$ will be assumed to have a Gaussian prior distribution with mean $`\beta _0`$ and standard deviation $`b`$. Physically the parameter $`\beta _0`$ represents a redshift-dependent evolution of the peak luminosity of the supernovae sources, which might be expected to arise as a result of the chemical evolution of the environment of the supernovae progenitors as abundances of heavier elements increase with cosmic time. Ultimately, one should hope to account for this evolution with astrophysical modelling of the supernovae explosions . The parameter $`b`$ would then account for a local variability in the supernovae environments between regions of individual galaxies at the same redshift which are richer or poorer in metals, or with progenitor populations of different ages and masses etc.

We now have two parameters to marginalize over, $`\nu `$ and $`\beta `$. As in the case of models with no evolution, we will marginalize over $`\nu `$ using a flat prior. We use a Gaussian prior for $`\beta `$ with mean $`\beta _0`$ and standard deviation $`b`$, so that

$$p(\beta ^{})=\frac{1}{b\sqrt{2\pi }}\mathrm{e}^{\beta ^2/2b^2}$$
(64)
where

$$\beta ^{}=\beta \beta _0.$$
(65)
The marginal likelihood is calculated by multiplying the prior (64) by the likelihood resulting from Eq. (59), and integrating over $`\beta `$. The resulting likelihood is

$`(M,f)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }\nu }}{\displaystyle 𝑑\beta ^{}𝑑\nu p(\beta ^{})\mathrm{e}^{\chi ^2/2}}`$ (66)
$`=`$ $`{\displaystyle \frac{s\overline{\sigma }\sqrt{2\pi }}{\mathrm{\Delta }\nu b}}\mathrm{e}^{Q/2}`$ (67)
where

$`Q(M,f)`$ $`=`$ $`{\displaystyle \frac{\widehat{\beta }^2}{\overline{\sigma }^2}}+{\displaystyle \underset{i=1}{\overset{60}{}}}{\displaystyle \frac{(h_is^2)^2}{\sigma _i^2}},`$ (68)
$`h_i(z_i;M,f)`$ $`=`$ $`m_{B,i}^{\text{eff}}g_i\beta _0\mathrm{ln}(1+z_i),`$ (69)
$`(M,f)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{60}{}}}{\displaystyle \frac{h_i}{\sigma _i^2}}.`$ (70)
The conditional best-fit of $`\beta ^{}`$ is given by

$`\widehat{\beta }(M,f)`$ $`=`$ $`\overline{\sigma }^2{\displaystyle \underset{i=1}{\overset{60}{}}}{\displaystyle \frac{h_i[\mathrm{ln}(1+z_i)s^2𝒢]}{\sigma _i^2}},`$ (71)
$`𝒢`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{60}{}}}{\displaystyle \frac{\mathrm{ln}(1+z_i)}{\sigma _i^2}},`$ (72)
and the $`\beta ^{}`$ uncertainity, $`\overline{\sigma }`$, is given by

$$\frac{1}{\overline{\sigma }^2}=\frac{1}{b^2}s^2𝒢^2+\underset{i=1}{\overset{60}{}}\frac{\mathrm{ln}^2(1+z_i)}{\sigma _i^2}.$$
(73)
Although $`\overline{\sigma }`$ is independent of $`M`$ and $`f`$, the marginal likelihood is no longer proportional to the profile likelihood because $`Q`$ is now given by (68) rather than by the chi-square type statistic (60).

We have performed a detailed numerical analysis on the P98 dataset, varying $`\beta _0`$, $`b`$ and $`w_i`$. As a result we find a best-fit value of $`\beta _0\beta ^{}0.414`$, for $`w_i=1.5`$ in the PNGB models. This would correspond to supernovae being intrinsically dimmer by $`0.17`$ magnitudes at a redshift of $`z=0.5`$, which is an effect of the typical order of magnitude being addressed in current attempts to better model the supernova explosions . Furthermore, we find that inclusion of a non-zero variance, $`b^2`$, does not alter the prediction of the best-fit value of $`\beta _0`$, although it naturally does lead to a broadening of the areas of parameter space included at the $`2\sigma `$ level. There is relatively little broadening of the region of parameter values included at the $`1\sigma `$ level in the $`\beta _0=\beta ^{}`$ plane, however, which is no doubt a consequence of the steepness of the $`q_0`$ contours in Fig. 7(a) in the area corresponding to region II.

In Figs. 9 and 10 we display the joint credible regions for $`M`$ and $`f`$, for two slices through the 3-dimensional $`(M,f,\beta _0)`$ parameter space for $`w_i=1.5`$: (a) the best-fit case $`\beta _0=\beta ^{}`$; and (b) $`\beta _0=0`$. Analogously, the best-fit case is shown in Fig. 11 for $`w_i=0.2`$. We see from Fig. 10 that once the likelihood is normalized relative to $`\beta ^{}`$ no regions remain in the $`\beta =0`$ parameter plane which are admitted at the $`1\sigma `$ level when $`b=0`$. Furthermore, even when a non-zero standard deviation, $`b`$, is included, region II of the $`(M,f)`$ parameter plane is favoured at the $`1\sigma `$ level, in contrast to Fig. 8.

We have also undertaken an analysis of the models with $`\beta _0=0`$ but variable $`b`$, similarly to the study of Ref. . In that case, we once again find that region II of the $`(M,f)`$ parameter plane is admitted at the $`1\sigma `$ level if $`b=0.25`$ or $`b=0.5`$. The dependence of the value of

$$\overline{Q}2\mathrm{ln}\left(\mathrm{\Delta }\nu \right)=Q2\mathrm{ln}\left(\frac{\sqrt{2\pi }s\overline{\sigma }}{b}\right)$$
(74)
on the value of $`b`$ is displayed in Fig. 12, for $`\beta _0=0`$ as compared with the best-fit case $`\beta _0=\beta ^{}`$. The quantity $`\overline{Q}`$ is analagous to the chi-square statistic of the maximum likelihood method. We see that the $`\beta _0=0`$ models favour a non-zero value of $`b0.36`$ by a very small margin as compared to the $`b=0`$ case. For $`b=0.36`$ the points of greatest likelihood lie mainly in region II, in contrast to the $`b=0`$ case in Fig. 8.

Varying the inital condition $`w_i`$ does not appear to affect the best-fit value of $`\beta _0`$ significantly. For $`w_i=0.2`$ (cf.,Fig. 11), for example, the best-fit value was $`\beta ^{}=0.435`$, a difference of 5% from the $`w_i=1.5`$ case. Furthermore, the numerical value of the least value of $`\overline{Q}`$ was only 0.2% greater in the $`w_i=0.2`$ case. We have not attempted to find a best-fit value for $`w_i`$.

## V Constraints from Lensing Statistics

Gravitational lensing of distant light sources due to the accumulation of matter along the line of sight provide another relatively sensitive constraint on the cosmological models of interest. For cosmology the situation of most interest is the lensing of high luminosity quasars by intervening galaxies. The abundance of multiply imaged quasars and the observed separation of the images to the source puts constraint on the luminosity–redshift relation and hence the model parameters. Basically, if the volume of space to a given redshift is larger then on average one can expect more lensing events. This leads to a statistical test, which has been used to put bounds on $`\mathrm{\Lambda }`$  and to test properties of some decaying $`\mathrm{\Lambda }`$ or quintessence models .

Gravitational lensing statistics are useful since they provide a tests which potentially provides opposing constraints to those obtained from supernovae magnitude–redshift tests. In particular, in the case of models with a vacuum energy provided by a cosmological constant, the high redshift supernovae have been interpreted as favouring relatively large values of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ \- Perlmutter et al.  give a value of $`\mathrm{\Omega }_\mathrm{\Lambda }=0.72`$ at $`1\sigma `$ level - whereas the gravitational lensing data leads to upper bounds on $`\mathrm{\Omega }_\mathrm{\Lambda }`$: Kochanek  quotes $`\mathrm{\Omega }_\mathrm{\Lambda }<0.66`$ at the $`2\sigma `$ level. A combined likelihood analysis has been performed by various authors .

Gravitational lensing constraints on the PNGB models have been very recently given by Frieman and Waga  for $`w_i=1.5`$. However, Frieman and Waga considered a more retricted range of parameter space, $`M<0.005h\text{eV}`$, since they did not consider the possibility of source evolution in the case of the type Ia supernovae and therefore took values of $`M>0.005h\text{eV}`$ to be ruled out. We wish to extend the range of $`M`$ for the gravitational lensing statistcs to consider parameter values corresponding to region II in the supernovae constraint graphs, Figs. 8, 10, so as to compare the constraints from different tests.

We have thus simply followed the calculation described by Waga and Miceli , who performed a statistical lensing analysis of optical sources described earlier by Kochanek . They used a total of 862 ($`z>1`$) high luminosity quasars plus 5 lenses from seven major optical surveys . (Another alternative not considered here is to analyse data from radio surveys – see, e.g., .) Undertaking a similar analysis for the increased parameter range, we arrive at Fig. 13, which shows the 68.3% and 95.4% joint credible regions for $`M`$ and $`f`$, for two values of $`w_i`$. We refer the reader to Refs.  for details of the calculation. The only regions of parameter space excluded at the $`2\sigma `$ level turn out to be areas of parameter space for which the deceleration is presently negative (cf., Fig. 7), with the scalar field still commencing its first oscillation at the present epoch.

## VI Discussion

Let us now consider the overall implications of the constraints observed above.

Firstly, since empirical models with source evolution do appear to fit the data somewhat better, it would appear that we do have weak evidence for an underlying evolution of the peak luminosity of the type Ia supernovae sources, at least in the context of the PNGB quintessence models. It might be interesting to compare the case of other quintessence models, or the case of a cosmological constant. However, the PNGB model is qualitatively different from such models since its final state corresponds to one in which the ultimate destiny of the universe is to expand at the same rate as a spatially flat Friedmann–Robertson–Walker model, rather than to undergo an accelerated expansion. This is of course precisely why we chose the PNGB models as the basis of our investigation, rather than models in which a late–time accelerated expansion had been built in by hand. If we wish to test the hypothesis that the faintness of the type Ia supernovae is at least partly due to an intrinsic variation of their peak luminosities – which is a very real possibility in view of the results of  – then a quintessence model which possesses a variety of possibilities for the present–day variation of the scale factor is probably the best type of model to investigate.

If only supernovae luminosity distances, (cf., Fig. 9, and gravitational lensing statistics, (cf., Fig. 13), are compared then we see that there is a remarkable concordance between the two tests – region II of Fig. 9 coincides with a region included at even the $`1\sigma `$ level in Fig. 13. This is perhaps not surprising, since in view of Fig. 7 region II corresponds to parameter values for which the present day universe has already undergone almost one complete oscillation of the scalar field about the final critical point $`\mathrm{C}_2`$ of Fig. 1. It is thus already well on the way towards its asymptotic behaviour, which closely resembles that of a standard spatially flat Friedmann-Robertson-Walker model.

Due to the oscillatory behaviour, parameter values in region II correspond to models in which there has been a recent cosmological acceleration, (e.g., at $`z0.2`$), but with a $`q(z)`$ which changes sign three times over the larger range of redshifts, $`0<z<4`$, in the quasar lensing sample, and therefore differing significantly from Friedmann–Lemaître models over this larger redshift range. Extending the SNe Ia sample to include objects at redshifts $`0.1<z<0.4`$ and $`z>0.85`$ in substantial numbers would greatly improve the ability to decide between models in regions I and II.

There is some cause for concern, however, if we consider the favoured values of $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$. In the case of the models with empirical evolution of the supernovae sources, we always found that the best-fit parameter values occured at the $`\mathrm{\Omega }_{\varphi 0}1`$ boundary of the ($`M,f`$) parameter space. The overwhelming evidence of many astronomical observations over the past two decades  would tend to indicate that $`\mathrm{\Omega }_{m0}0.2\pm 0.1`$, indicating that a vacuum energy fraction of $`\mathrm{\Omega }_{\varphi 0}0.7`$$`0.8`$ is desirable, and $`\mathrm{\Omega }_{\varphi 0}\genfrac{}{}{0pt}{}{<}{}0.9`$ in any case. Although parameter values with $`\mathrm{\Omega }_{\varphi 0}<0.8`$ certainly fall within both the 2$`\sigma `$ and 1$`\sigma `$ portions of region II of Fig. 9, for all values of $`b`$, there are potentially serious problems if we wish to simultaneously obtain large values of $`H_0t_0`$. In view of recent estimates of the ages of globular clusters , a lower bound of 12 Gyr for the age of the Universe appears to currently indicated. With $`h0.65`$ this would require $`H_0t_0\genfrac{}{}{0pt}{}{>}{}0.8`$. For $`w_i=1.5`$, parameter values with $`H_0t_0>0.8`$ coincide with values $`\mathrm{\Omega }_{\varphi 0}\genfrac{}{}{0pt}{}{>}{}0.9`$ in region II, which is phenomenologically problematic.

The tension between the values of $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ is somewhat mitigated for lower values of $`w_i`$. For $`w_i=0.2`$, for example, we see from Figs. 11 and 13(b) that the $`\mathrm{\Omega }_{\varphi 0}=0.7`$ and $`H_0t_0`$ meet in Region II, and there is a small region of parameters there with $`0.7\genfrac{}{}{0pt}{}{<}{}\mathrm{\Omega }_{\varphi 0}\genfrac{}{}{0pt}{}{<}{}0.9`$ and $`H_0t_0\genfrac{}{}{0pt}{}{>}{}0.8`$, which is also consistent with the other cosmological tests.

Even if the supernovae sources undergo evolution it is clear that parameter values in region I of 8 which are favoured in the absence of evolution of peak SnIa luminosities, are still included at the $`2\sigma `$ level in the models with evolution, in view of Fig. 9.

Perhaps the most significant aspect of our results is the fact that the mere introduction of an additional dispersion, $`b>0.17`$, in the peak luminosities while leaving their mean value fixed (cf., Fig. 12), gives rise to a change in the best-fit region of parameter space from region I to region II. (See  for further details.) One would imagine that an increased dispersion is likely to be a feature of many models of source evolution, even if evolutionary effects are of secondary importance. Thus even if the empirical models with non-zero $`\beta `$ are somewhat artificial, more sophisticated scenarios could well lead to similar changes in regard to the fitting of cosmological parameters in the PNGB model.

Much tighter bounds on the parameter space of quintessence models, including the present model, will be obtained over the next decade as more supernovae data is collected. What we wish to emphasize, however, is that an effective vacuum energy which is cosmologically significant at the present epoch should not simply be thought of in terms of a “cosmic acceleration”. A dynamical vacuum energy with a varying effective equation of state allows for many possibilities for the evolution of the universe, and overly restrictive assumptions, such as equating quintessence to models with a late period of continuous cosmological acceleration, should be avoided. If detailed astrophysical modeling of type Ia supernovae explosions ultimately shows that the dimness of distant supernova events is largely due to evolutionary effects, it does not spell the end for cosmologies with dynamical scalar fields.

## Acknowledgments

We would like to thank Chris Kochanek for supplying us with the gravitational lensing data, Elisa di Pietro, Don Page and Ioav Waga for helpful discussions about various aspects of the paper, and the Australian Research Council for financial support.


text:
# Printed on Photogalvanic Effects in Heteropolar Nanotubes

\[

## Abstract

We show that an electrical shift current is generated when electrons are photoexcited from the valence to conduction bands on a BN nanotube. This photocurrent follows the light pulse envelope and its symmetry is controlled by the atomic structure of the nanotube. We find that the shift current has an intrinsic quantum mechanical signature in which the chiral index of the tube determines the direction of the current along the tube axis. We identify the discrete lattice effects in the tangent plane of the tube that lead to an azimuthal component of the shift current. The nanotube shift current can lead to ultrafast opto-electronic and opto-mechanical applications.

\]

Recent progress in synthesis of nanometer scale materials has led to the discovery of $`\mathrm{B}_\mathrm{x}\mathrm{C}_\mathrm{y}\mathrm{N}_\mathrm{z}`$ nanotubes. They can be formed in nearly homogeneous ordered structures  or in multiwall hybrid structures with alternating BN and C compositions . These materials differ fundamentally from the structurally similar carbon nanotubes by being non-centrosymmetric (NCS) and polar. This opens the possibility for accessing a new class of photovoltaic effects at the molecular scale.

Photovoltaic effects in NCS materials are often based on an asymmetric generation of hot carriers at momenta $`\pm 𝐤`$ leading to a so called “ballistic photocurrent” . In polar NCS materials, photoexcitation across the bandgap with polarized light also produces the so called electrical “shift current” . Microscopically this originates from a net displacement of charge in the unit cell due to light induced interband transitions ($`𝐉_𝐞`$), intraband relaxation ($`𝐉_𝐬`$) and (radiative) transitions to the original bands ($`𝐉_𝐫`$). In bulk materials, the excitation component $`𝐉_𝐞`$ usually prevails.

We show that unpolarized light can induce a shift current $`𝐉_𝐞`$ in polar NCS nanotubes, with a direction along the tube axis which is determined by the chiral index of the tube. We find that this effect has an essentially quantum mechanical origin, where the sign of the current along the tube axis is controlled by the phase matching of electronic Bloch waves around the tube circumference. The discrete lattice structure controls the azimuthal component of this current and thus produces a net helical current on the tube. These photo-effects can lead to an assortment of new opto-electronic, opto-mechanical and magnetic applications, and interesting extensions to ring structures  and heterojunctions . We should stress that the photoeffects we discuss here are found at nonlinear order in the exciting fields, and are therefore physically distinct from chiral currents tilted with respect to the tube axis, predicted in dc bias driven chiral BN nanotubes  or chiral stretched C nanotubes .

A flat BN sheet has a honeycomb lattice with the B and N occupying alternate sublattices as shown in Figure 1. The physics of the shift current can be understood qualitatively by considering the response of this 2D network of bonds to normally incident polarized light. The valence states are polarized towards the N sites and the conduction states towards the B sites. For vertically polarized ($`y`$-polarized) incident light the response of the system is dominated by the bonds which produce a net $`y`$ polarized current. For horizontally ($`x`$-polarized) incident light, the bonds with a nonzero horizontal component dominate, but these also produce a net $`y`$-polarized current. Note that excitation with an unpolarized incident field cannot produce a net current on this lattice since it has a threefold symmetry. However, this symmetry is removed when the sheet is wrapped into a cylindrical nanotube where the physics should be dominated by excitations with the field polarized along the tube axis. Depending on the wrapped structure of the tube, we anticipate an intrinsic photocurrent which can flow along the tube, around the tube, or in a chiral pattern on the tube surface.

These ideas can be quantified by developing a quantum mechanical model which generalizes our long wavelength theory for a nonpolar carbon nanotube  to the heteropolar lattice. We work in a basis of Bloch orbitals $`\mathrm{\Phi }_{k\alpha }(r)=e^{ikr}_ne^{ik(rR_n\tau _\alpha )}\varphi _{n\alpha }/\sqrt{N}=e^{ikr}U_{k\alpha }`$ where $`\alpha =\pm 1`$ denote the two sublattices occupying sites at $`\tau _\alpha `$ in the $`n`$-th unit cell. We study the states near the conduction and valence band edges at the K and K’ points of the Brillouin zone shown in Figure 1. A long wavelength Hamiltonian for these states is obtained by an expansion in small crystal momenta around these points, and yields in our two component basis ($`\mathrm{}=1`$)

$$H_\lambda (\mathrm{\Delta },k,\delta _\lambda )=\left(\begin{array}{cc}\mathrm{\Delta }& \lambda v_F(ki\delta _\lambda )\\ \lambda v_F(k+i\delta _\lambda )& \mathrm{\Delta }\end{array}\right)$$
(1)
where $`\lambda =\pm 1`$ is an index which labels an expansion around the K or K’ points. The Hamiltonian in Eq. (1) is parameterized by three energies: a symmetry breaking site diagonal potential $`\mathrm{\Delta }`$, the kinetic energy $`v_Fk`$ for motion along the tube, and the kinetic energy $`v_F\delta `$ due to the quantization of the transverse momentum around the circumference of the tube. On a cylindrical tube with lattice constant $`a`$ and primitive lattice translation vectors in its tangent plane $`T_1=a(1,0)`$ and $`T_2=a(1/2,\sqrt{3}/2)`$ (Figure 1) it is conventional to index the lattice structure by two integers $`m`$ and $`n`$ which define a superlattice translation vector $`𝒞_{mn}=mT_1+nT_2`$. The transverse momenta on an $`(m,n)`$ tube are quantized to the values $`\delta _N=\delta _0+2\pi N/|𝒞_{mn}|`$ where $`\delta _0=2\pi \mathrm{sgn}(\nu )/3|𝒞_{mn}|`$, depending on the sign of the chiral index of the tube, $`\nu =\mathrm{mod}(nm,3)`$. In Eq. (1) $`\delta _\lambda =\lambda \delta _N`$.

The eigenvectors of $`H`$ in Eq. (1) are

$$g_v=\left(\begin{array}{c}ue^{i\varphi /2}\\ ve^{i\varphi /2}\end{array}\right)\frac{e^{i\theta _{vk}}}{\sqrt{2}};g_c=\left(\begin{array}{c}ve^{i\varphi /2}\\ ue^{i\varphi /2}\end{array}\right)\frac{e^{i\theta _{ck}}}{\sqrt{2}}$$
(2)
where $`E=\sqrt{v_F^2(k^2+\delta _\lambda ^2)+\mathrm{\Delta }^2}`$, $`\varphi =\mathrm{tan}^1(\delta _\lambda /\lambda k)`$, $`u=\sqrt{(E\mathrm{\Delta })/E}`$, and $`v=\sqrt{(E+\mathrm{\Delta })/E}`$. In the rest of this paper we will set the Fermi velocity $`v_F=1`$ so that energies and momenta are measured in the same units. Eq. (2) explicitly includes a gauge function $`\theta _{mk}`$ since the overall phase of the Bloch function is not fixed, and our calculation will require differentiation with respect to $`k`$. The Bloch eigenfunctions of our problem are expressed as the product of three factors $`\psi _{km}(r)=e^{ikr}_\alpha g_{m\alpha }(k)U_{k\alpha }(r)`$.

We extend the model Hamiltonian in Eq.(1) to include coupling of electrons to the oscillating optical fields $`𝐄(t)=𝐄e^{i\omega t}+c.c.`$ through the dipole operator  $`H_{int}=e𝐄(t)𝐫`$. For a field polarized along the tube direction the matrix elements of the dipole operator between band eigenstates $`m`$ and $`m^{}`$ at crystal momenta $`k`$ and $`k^{}`$ are expressed using a formulation due to Blount

$`r_{m^{}m}(k^{},k)`$ $`=`$ $`\psi _{k^{}m^{}}|r|\psi _{km}=i\psi _{k^{}m^{}}|_k|\psi _{km}`$ (3)
$`+`$ $`i\delta (kk^{}){\displaystyle \underset{\alpha \beta }{}}\delta _{\alpha \beta }g_{m^{}\beta }^{}(k^{})_kg_{m\alpha }(k).`$ (4)
The second term on the r.h.s of (4) forms the connection $`\xi _{m^{}m}=ig_m^{}^{}(k)_kg_m(k)`$ due to the $`k`$ dependence of the eigenstates in (2) .

In general the shift current $`𝐉`$ can be expressed in terms of the band off-diagonal matrix elements of the velocity operator $`𝐯_{\beta \alpha }`$ ($`𝐯=i[H,𝐫]`$), and the band off-diagonal term in the density matrix $`\rho _{\alpha \beta }`$, calculated to second order in the exciting fields $`𝐄`$, as follows

$$𝐉=2\mathrm{e}\underset{mn}{}\frac{dk}{2\pi }𝐯_{mn}(k)\rho _{nm}(k)=𝐉_e+𝐉_s+𝐉_r.$$
(5)
We focus on the component of the excitation current $`\mathrm{J}_\mathrm{e}`$ along the tube direction. After evaluating the sum in (5), we arrive at the transparent result

$$\mathrm{J}_\mathrm{e}=2e\frac{dk}{2\pi }\dot{f}_{cv}(k)R_{cv}(k)e\dot{n}.$$
(6)
Here $`\dot{f}_{cv}(k)`$ is the transition rate at wave vector $`k`$, and the shift vector $`R_{cv}(k)`$ is given by

$$R_{cv}(k)=_k\theta _{vc}(k)+\xi _{cc}(k)\xi _{vv}(k)$$
(7)
with $`\xi _{vc}(k)=|\xi _{vc}(k)|e^{i\theta _{vc}(k)}`$. The shift vector $`R_{cv}(k)`$ is invariant under the gauge transformations $`\mathrm{exp}(i\theta _{mk})`$ and we will evaluate it in a gauge with $`_k\theta _{mk}=0`$.

For an incident electric field polarized along the tube direction, interband excitations are allowed only between band states with the same transverse momentum $`\delta `$. We consider transitions between the lowest two azimuthal subbands ($`N=0`$). Then using the eigenfunctions in Eq.(2) we find

$$\xi _{vc}(k)=\frac{1}{2}uv\left(\frac{\mathrm{\Delta }}{E^2\mathrm{\Delta }^2}_kE+i_k\varphi \right),$$
(8)
which gives the off-diagonal contribution in Eq. (7)

$$_k\theta _{vc}=_k\mathrm{tan}^1\left(\frac{\delta \sqrt{\delta ^2+\mathrm{\Delta }^2+k^2}}{k\mathrm{\Delta }}\right).$$
(9)
The diagonal elements in Eq. (7) are

$`\xi _{mm}(k)=ig_m^{}_kg_m={\displaystyle \frac{\mathrm{\Delta }}{2E}}\left({\displaystyle \frac{\delta }{\delta ^2+k^2}}\right),`$ (10)
where the sign is negative for valence states and positive for conduction states. Combining Eq. (9) and Eq. (10) we obtain the shift vector

$$R_{cv}(k)=\frac{2\delta \mathrm{\Delta }}{(\delta ^2+k^2)\sqrt{\delta ^2+\mathrm{\Delta }^2+k^2}}.$$
(11)
Comparison with Eq. (10) shows that the shift current (Eq.(6)) opposes the direction of the ground state polarization of the tube.

This shift vector is odd in the symmetry breaking potential $`\mathrm{\Delta }`$ and odd in the transverse momentum $`\delta `$. Therefore, two BN tubes with nearly the same wrapped lattice structure but opposite chiral indices $`\nu `$ exhibit reversed ground state polarizations and opposite shift currents. Armchair BN nanotubes with wrapping indices $`(m,m)`$ have chiral index $`\nu =0`$ and do not exhibit a longitudinal shift current (as shown in the middle panel in Figure 2 for the (5,5) tube). Zigzag BN tubes with wrapping indices $`(m,0)`$ can be grouped into three families with positive, negative or zero shift currents distinguished by the sign of the chiral index $`\nu `$ (an example for the (8,0) tube with $`\nu =1`$ is in the left panel of Figure 2). It is important to notice that an isolated flat BN sheet has a threefold symmetry axis perpendicular to the BN plane and therefore it has zero electric polarization in its tangent plane, by symmetry. Thus it has no shift current for excitations at the band edge. The nonzero shift vector in Eq. (11) when $`k0`$ is a remarkable long wavelength quantum mechanical effect that reveals the quantization of the transverse momentum on the wrapped tube.

The shift current $`\mathrm{J}_\mathrm{e}`$ in Eq. (6) is determined by the average shift vector $``$, describing the change of polarity in the excitation, and by the carrier injection rate $`\dot{n}`$. The BN bond has significant ionic character with  2$`\mathrm{\Delta }`$ 4 eV. This gives for excitations near the band edge, $`k0`$, a shift vector with a magnitude $`0.06\mathrm{nm}`$ for a (17,0) tube. For an incident intensity $`\mathrm{S}=100\mathrm{kW}/\mathrm{cm}^2`$ we can also obtain the interband excitation rate; we estimate that $`\dot{n}70\mathrm{nm}^1\mu \mathrm{s}^1`$ produces a net electronic shift current of $`\mathrm{J}_\mathrm{e}`$ 0.67 pA. This value is slighly larger than one obtains for excited bulk polar semiconductors. In a coaxial multiwall tube $`\mathrm{R}_{\mathrm{cv}}`$ can change magnitude and sign on each wall. Although opposing currents on different tube walls reduce the macroscopic current a complete cancellation is very unlikely since the magnitude of $`\mathrm{R}_{\mathrm{cv}}`$ depends on the tube radius. We also expect that the recombination current $`\mathrm{J}_\mathrm{r}`$ will oppose $`\mathrm{J}_\mathrm{e}`$ since both currents rely on optical transitions with longitudinal polarization. However, this should not prevent the observation of $`\mathrm{J}_\mathrm{e}`$ using pulsed laser excitation since the timescale for $`\mathrm{J}_\mathrm{r}`$ is of order nanoseconds.

This formalism can be extended to study chiral shift currents with nonvanishing longitudinal and azimuthal components. Interestingly, the long wavelength Hamiltonian in Eq. (1) gives zero azimuthal current for any chiral index $`\nu `$. As $`Q_\pm =k\pm i\delta 0`$ this occurs because of a cancellation of contributions from the K and K’ points whose long wavelength Hamiltonians have opposite handedness. However, azimuthal currents arise due to discrete lattice effects occurring at higher orders in an expansion in $`Q_\pm `$. To study this, we define two unit vectors at each point in the tube surface, perpendicular to the local outward unit normal $`\widehat{n}`$: $`\widehat{e}_x`$ which points along the tube axis, and $`\widehat{e}_y=\widehat{n}\times \widehat{e}_x`$ which points “counterclockwise” on the tube surface. A chiral angle $`\theta `$ is defined as the angle between $`\widehat{e}_x`$ and the $`T_1`$ lattice direction as shown in Figure 1 (with this convention armchair nanotubes have $`\theta `$ = 0). Then the lowest nonvanishing contribution to an azimuthal ($`y`$-polarized) shift current must have the form $`\mathrm{cos}(6\varphi 3\theta )`$ where $`\varphi =\mathrm{tan}^1(\delta /k)`$. The dependence on $`6\varphi `$ means that azimuthal currents first occur in the continuum theory at sixth order in $`|Q|`$.

We can study these terms by replacing the kinetic energy terms in the Hamiltonian of Eq. (1) by the discrete lattice counterpart $`iv_F(k+i\delta )=iv_FQ_+t_ne^{i(K+Q)\tau _n}`$ where $`t`$ is the hopping amplitude between sites connected by the bond vectors $`\tau _n`$. A systematic expansion of the azimuthal shift current in powers of $`Q_\pm `$ then gives contributions with the correct lattice symmetry and shows that the ratio of the azimuthal and longitudinal in-plane current densities scales as

$$p=\frac{J_y}{J_x}=\frac{(k^2+\delta ^2)^3\tau ^3}{\delta (\mathrm{\Delta }^2+k^2+\delta ^2)}\mathrm{cos}(6\varphi 3\theta ).$$
(12)
The azimuthal current for any zigzag wrapping is zero because of the angular factor, while its longitudinal current, shown in Figure 2, is determined by the chiral index $`\nu `$. For armchair tubes and all non-zigzag tubes with chiral index $`\nu =0`$, we have $`\delta =0`$, so the shift current is purely azimuthal with a direction determined by the sign of the potential $`\mathrm{\Delta }`$. These energy dependent currents vanish at the band edge as $`(E^2\mathrm{\Delta }^2)^3`$. Non-zigzag tubes with chiral index $`\nu 0`$ have a nonzero $`\delta `$ and thus the azimuthal current is nonzero with $`p(\delta /\mathrm{\Delta })^2(\delta \tau )^3`$. There the shift current circulates in a helical pattern on the tube surface as shown in the right panel of Figure 2, with a chirality determined by the index $`\nu `$. The pitch is also energy dependent and increases according to $`p(E)=p(0)+a(E^2\mathrm{\Delta }^2)`$. Thus, as shown in Figure 2, the zigzag nanotube generate a photocurrent like a wire, the armchair tubes like a coil, and the chiral tubes can exhibit a helical current on the tube surface.

Alternatively, the shift current can be studied in a semiclassical model, by considering an expansion in $`t/\mathrm{\Delta }`$, where $`t`$ is the nearest neighbor hopping amplitude, and $`\mathrm{\Delta }`$ is the site diagonal potential. In the limit $`t/\mathrm{\Delta }<<1`$ the tube can be regarded as a network of independent bonds whose excitations are superimposed. Here the local symmetry of the zigzag tube requires zero net azimuthal current while the armchair wrapping has a nonzero azimuthal current with its direction determined by the the sign of $`\mathrm{\Delta }`$. However, it is difficult to correctly describe the physics of the longitudinal current using this model. This is because the longitudinal current is controlled by the quantization of the transverse momentum $`\delta `$. The first terms in the semiclassical expansion that are sensitive to $`\delta `$ occur at order $`(t/\mathrm{\Delta })^N`$ where $`N`$ is the number of bonds ($`40`$) around the tube circumference. This reflects the fact that the longitudinal current is fundamentally a nonclassical quantity for this system.

The longitudinal shift current reduces the ground state electric polarization of the nanotube. This can lead to a fast readjustment of atomic positions on the tube walls. Generation of voltages by a mechanical elongation  of NCS nanotubes is a related effect. In both situations the response can be stronger in buckled nanotubes  or if different types of atoms occur in different layers, as in MoS2, WS2  or GaSe nanotubes . Recently, mechanical motions, strikingly similar to those found in some biomolecular systems , were observed in irradiated nanotube bundles . Interesting mechanical response to an applied dc bias, is observed in carbon nanotubes immersed in a solution of NaCl . Here, the motion is caused by attraction of ions of a given polarity to the nanotube. The speed but short duration of the shift current-induced mechanical response complements the slow but steady state effect of a dc bias.

Experimental observation of the shift current in BN nanotubes would provide a striking illustration of fundamental concepts in the modern quantum theory of polarization. The effect is also promising for new applications since it couples the microscopic physics on short length scales, which can be tuned by local fields and mechanical loads, to the long distance properties of the system (polarization, photovoltage, etc.) Finally, there is a rapidly growing family of related submicron one dimensional materials to which these ideas can be applied.

Acknowledgments
PK acknowledges J.E. Sipe and M. Shapiro for financial support. This work was also supported by the DOE under Grant number DE-FG02-84ER45118 (EJM) by the NSF under Grant number DMR 98-02560 (EJM) and by the ONR under Grant Number N00014-99-1-0252 (DT).


text:
# Untitled Document

% auto-ignore

The paper is withdrawn by the authors due to flaw in the proof of lemma 3.