text:
# Orphan GRB radio afterglows: Candidates and constraints on beaming

## 1 Introduction

Our understanding of gamma-ray bursts (GRBs) has been revolutionized by the discovery of x-ray afterglows by BeppoSAX (e.g. Costa et al. (1997)), and the following detections of optical transients (e.g. van Paradijs et al. (1997)). These efforts led eventually to a confirmation of the cosmological nature of GRBs, through both direct redshift measurements (e.g. Metzger et al. (1997)), and imaging of the host galaxies (e.g. Sahu et al. 1997). The enormous power released in the explosion implies extremely large compactness of the source and, therefore, most GRB models involve compact or collapsed objects (e.g. Goodman (1986); Paczyński (1986); Eichler et al. (1989); Woosley (1993); Levinson & Eichler (1993); Paczyński (1998); see Mészáros 1999 for a review). In spite of the impressive successes of expanding relativistic “fireball” models (e.g.Paczyński & Rhoads (1993); Katz (1994); Mészáros & Rees (1997); Vietri (1997); Waxman 1997a ; Sari, Piran & Narayan (1998); see Mészáros (2002) for review), the precise nature of GRB progenitors remains unknown.

The total energy emitted in the GRB explosion and the rate at which such explosions occur in the Universe are important keys to the understanding of the progenitors. The determination of these factors is complicated, however, by virtue of relativistic beaming. During the phase of $`\gamma `$-ray emission the fireball expands with a large Lorentz factor, $`\mathrm{\Gamma }10^{2.5}`$, so that a distant observer receives radiation from a conical section of the fireball of opening angle $`10^{2.5}`$ around the line of sight. Thus, estimates of total energy and rate based on $`\gamma `$-ray observations are highly uncertain. The first evidence that the fireball may be jetted was provided by radio observations of GRB970508 (Waxman, Kulkarni & Frail (1998)). These observations imply a jet of relatively wide opening angle, which expands sideways and approaches sub-relativistic, spherical expansion $`1`$ yr following the GRB (Frail, Waxman & Kulkarni (2000)). The analysis of radio observations during the sub-relativistic phase allowed to determine the total GRB energy to be $`10^{51}`$ erg.

Jetted GRBs have been widely invoked to explain the optical light curves of the afterglow emission, most notably so for the source GRB990123 (Stanek et al. (1999); Harrison et al. (1999)). Frail et al. (2001) have analyzed a sample of GRB afterglows with known redshifts. They find that most bursts are jetted, with jet opening angle $`\theta 0.1`$ and average “beaming factor,” defined as $`f_b^1(\theta ^2/2)^1`$, of $`f_b^1500`$. They also find that the total $`\gamma `$-ray energy release, when corrected for beaming as inferred from the afterglow light curves, is narrowly clustered around $`0.5\times 10^{51}`$ erg. This result, although somewhat model dependent, also suggests that the conversion efficiency of fireball kinetic energy to radiation is high, in agreement with the conclusion of Freedman & Waxman (2001), who derived total fireball energy using early X-ray afterglow data. The energy estimates of Frail et al. are in agreement with those derived by Freedman & Waxman, and also with the total energy derived for GRB970508 (Frail, Waxman & Kulkarni (2000)).

Using radio surveys to constrain GRB beaming, rate and energetics through the search for radio emission from GRBs that were not necessarily detected in $`\gamma `$-rays, has been proposed by Perna & Loeb (1998) , Woods & Loeb (1999) and by Paczyński (2001). Perna & Loeb have suggested to search for such ”orphan radio afterglow” emission from GRB jets pointing toward us, while Woods & Loeb have suggested to look for $`10^3`$ yr old, non-relativistic remnants. Paczyński suggested, based on the radio observations of GRB970508, to search for emission from $`1`$ yr old nearby GRB remnants, which have undergone the transition from a relativistic jet to spherical sub-relativistic expansion, during which emission is isotropic. Given the hints for an association of GRBs with supernovae (Galama et al. (1998); Bloom et al. (1999); Reichart (1999)), it was suggested by Paczyński to search for strong radio emitters among nearby supernovae.

In this paper we use the model, proposed in (Frail, Waxman & Kulkarni (2000); Livio & Waxman (2000)) for the long term radio afterglow emission of GRBs, to calculate the expected number of orphan afterglows in radio surveys. We find that the detectable orphan population is dominated by $`1`$ yr old GRB remnants which have just undergone the transition to sub-relativistic, spherical expansion. We determine the dependence of the expected number of sources on uncertain model parameters, derive a lower limit to the expected number of sources, and demonstrate that independent constraints on the beaming factor of early GRB emission can be imposed using radio surveys. We searched for point-like radio transients, by comparing the FIRST and NVSS wide field catalogues, taken a few years apart, looking for sources present in one catalogue and absent in the other. Our approach is different than that proposed by Paczyński, namely we look directly for radio transients without relying on an association with nearby supernovae, since the link between GRBs and luminous supernovae is not well established: The only clear case, 1998bw, had an atypical (very weak) GRB, while there are three low redshift, $`z0.45`$, GRBs111see http://www.aip.de/ jcg/grbgen.html with no prominent supernova reported.

In § 2 we provide a brief description of the afterglow model and its main results used in later sections. In § 3 we calculate the number of radio afterglows expected to be detected above a certain flux. In § 4 we present a search for point-like radio transients using the FIRST and NVSS catalogs, and discuss possible contamination by radio supernovae and AGN. The implications of our results are discussed in § 5.

## 2 An outline of the afterglow model

The basic model assumes an ultra-relativistic, jet like ejecta propagating into an ambient medium of density $`n=1n_0\mathrm{cm}^3`$, slowing down and expanding sideways, ultimately becoming non-relativistic. As long as the jet Lorentz factor $`\mathrm{\Gamma }`$ is larger than the inverse of the jet’s opening angle $`\theta `$, it behaves as if it were a conical section of a spherical fireball. Once $`\mathrm{\Gamma }`$ drops below $`\theta ^1`$, the jet expands sideways, and its behavior deviates from that exhibited in the conical phase (Rhoads (1997, 1999)). After a transition stage, in which the jet expands sideways, the flow approaches spherical symmetry and becomes sub-relativistic (Waxman, Kulkarni & Frail (1998); Frail, Waxman & Kulkarni (2000)). The transition from a jet to a spherical sub-relativistic evolution takes place over a time scale,

$$t_{SNT}6\times 10^6(E_{51}/n_0)^{1/3}\mathrm{s}$$
(1)
(Frail, Waxman & Kulkarni (2000); Livio & Waxman (2000)). Here, $`E_{51}`$ is the total energy carried by the jet, in units of $`10^{51}`$ erg, corrected for beaming, i.e. it is related to the isotropic equivalent energy $`E_{\mathrm{iso}}`$ by $`E=(\theta ^2/2)E_{\mathrm{iso}}=f_bE_{\mathrm{iso}}`$, where $`\theta `$ is the jet opening angle. The subscript SNT refers to the Sedov-von Neumann-Taylor self-similar solution that describes the sub-relativistic flow.

After the transition, the power is dominated by synchrotron emission from a sub-relativistic fireball (e.g. Frail, Waxman & Kulkarni 2000). Denoting by $`\xi _e`$ ($`\xi _B`$) the fraction of thermal energy behind the shock that is carried by electrons (magnetic field), and assuming that the electrons are accelerated to a power-law energy distribution, $`dn_e/d\gamma _e\gamma _e^p`$ with $`p=2`$, it can be shown (Livio & Waxman (2000)) that the flux at frequencies above the synchrotron peak frequency at $`t=t_{\mathrm{SNT}}`$,

$$\nu _{}1\left(\frac{1+z}{2}\right)^1\left(\frac{\xi _e}{0.3}\right)^2\left(\frac{\xi _B}{0.3}\right)^{1/2}n_0^{1/2}\mathrm{GHz},$$
(2)
is given by (see also Appendix of Frail, Waxman & Kulkarni 2000)

$$f_\nu 2h^2\left(\frac{1+z}{2}\right)^{1/2}\left(\frac{\xi _e}{0.3}\right)\left(\frac{\xi _B}{0.3}\right)^{3/4}\left(\frac{d_L}{R_0}\right)^2n_0^{3/4}E_{51}\nu _9^{1/2}\left[\frac{t}{t_{\mathrm{SNT}}(1+z)}\right]^{9/10}\mathrm{mJy}.$$
(3)
Here, $`d_L`$ is the luminosity distance, $`R_0=c/H_0=10^{28}h^1`$ cm, and $`\nu _9`$ is the observed frequency measured in GHz. Analysis of GRB970508 radio data implies (Frail, Waxman & Kulkarni (2000)) $`n_01`$, $`\xi _e\xi _B1/3`$.

## 3 Expected number of orphan afterglow sources in radio surveys

Consider a survey with a flux threshold denoted by $`f_{\nu min}`$, and denote by $`\tau _m`$ the time, in units of $`t_{SNT}(1+z)`$, at which the radio flux of a GRB at a given luminosity distance, $`d_L`$, drops below detection limit. From eq. (3) we then obtain

$$\tau _m=\eta (1+z)^{5/9}(d_L/R_0)^{20/9},$$
(4)
where

$$\eta =h^{20/9}\left(\frac{f_{\nu min}}{1mJy}\right)^{10/9}\left(\frac{\xi _e}{0.3}\right)^{10/9}\left(\frac{\xi _B}{0.3}\right)^{10/12}\left(\frac{2}{\nu _9}\right)^{5/9}n_0^{10/12}E_{51}^{10/9}.$$
(5)
We shall assume, in what follows, that the total energy output and beaming factor of GRBs are standard and, moreover, that the fraction of total energy release that goes into driving the external blast wave generating the afterglow emission, is the same for all GRBs. We make this assumption in order to simplify our analysis, and since the distribution of the above parameters is not well known. It is important to emphasize at this point that, as mentioned above, the analysis presented in Frail et al. (2001) indicates that the distribution of GRB true (i.e. corrected for beaming) gamma-ray energies is narrowly clustered around $`E_{\gamma ,51}0.5`$, with a corresponding average beaming factor of $`<f_b^1>=<(\theta ^2/2)^1>500`$. These results also imply high, order unity, radiative efficiency, consistent with the result obtained by Freedman & Waxman (2001) using a different method, and implying that the true total fireball energy is a few times $`E_\gamma `$, i.e. narrowly peaked around $`E_{51}1`$. Thus, we expect corrections to our results, due to deviations of parameters from their standard values, to be small. Moreover, our results, given as explicit functions of the standard parameter values, can be easily generalized to the case where parameter values are drawn from a general distribution, by replacing these values with their corresponding averages.

Let $`\tau _i`$ be the time at which the radio source just becomes isotropic. Then for a given flux threshold, the maximum luminosity distance below which a substantial number of sources can be detected is obtained from eq. (4) for $`\tau _m=\tau _i`$:

$$d_{Lmax}(z_{max})=R_0(\eta /\tau _i)^{9/20}(1+z_{max})^{1/4}.$$
(6)
For a given choice of GRB parameters, a flux threshold $`f_{min}`$ at the observed frequency $`\nu _9`$, and value of $`h`$, the parameters $`\eta `$ and $`\tau _i`$ are fixed. Equation (6) can then be numerically solved for the assumed cosmology to yield the maximum redshift $`z_{max}`$ encompassed by the survey.

We denote by $`\dot{n}`$ the observed number of GRB events per year per Gpc3 at $`z=0`$ (which, assuming that the GRB rate follows the redshift evolution of the star formation, is estimated to be $`\dot{n}0.5`$ (Schmidt (2001))), and by $`\mathrm{\Phi }(z)`$ the redshift evolution factor. The expected number of radio afterglows between $`z`$ and $`z+dz`$ can then be expressed as,

$$\frac{dN_R}{dz}=N_0d_L^2\frac{dl}{dz}(\tau _m\tau _i)\mathrm{\Phi }(z),$$
(7)
where $`dl=R_0^1cdt`$ is the light distance in units of $`R_0`$, and

$$N_0=f_b^14\pi R_0^3\dot{n}t_{SNT}3.5\times 10^4(500f_b)^1h^3\dot{n}(E_{51}/n_0)^{1/3}.$$
(8)
The total number expected, $`N_R`$, can be obtained by integrating eq (7) from $`z=0`$ to $`z=z_{max}`$, with $`z_{max}(\eta /\tau _i)`$ determined from eq. (6).

As an example, consider a survey with a flux threshold of 5 mJy. Assuming $`\tau _i=3`$, $`h=0.75`$ and the above choice for the remaining parameters, we find $`d_{Lmax}0.2R_0=0.8`$ Gpc, and a corresponding redshit $`z_{max}0.2`$. At such small redshifts cosmological effects can be neglected to a reasonable approximation, and the luminosity distance can be taken to be equal to the proper distance. In this case eq. (7) can be integrated analytically. One then finds, using $`n=10^1n_1\mathrm{cm}^3`$

$`N_R`$ $`N_0\eta ^{27/20}\tau _i^{7/20}`$ (11)
$`18(500f_b)^1(\dot{n}/0.5)\left({\displaystyle \frac{f_{\nu min}}{5mJy}}\right)^{3/2}\left({\displaystyle \frac{\xi _e}{0.3}}\right)^{3/2}\left({\displaystyle \frac{\xi _B}{0.03}}\right)^{9/8}n_1^{19/24}E_{51}^{11/6}\nu _9^{3/4}(\tau _i/3)^{7/20}`$
$`10^4f_b^{5/6}(\dot{n}/0.5)\left({\displaystyle \frac{f_{\nu min}}{5mJy}}\right)^{3/2}\left({\displaystyle \frac{\xi _e}{0.3}}\right)^{3/2}\left({\displaystyle \frac{\xi _B}{0.03}}\right)^{9/8}n_1^{19/24}E_{\mathrm{iso},54}^{11/6}\nu _9^{3/4}(\tau _i/3)^{7/20}.`$
As seen, $`N_R`$ depends only weakly on $`\tau _i`$. Consequently, any uncertainty in the time at which the radio emission becomes roughly isotropic should not affect our result considerably. We have also given the result in terms of the isotropic equivalent energy, $`E_{\mathrm{iso}}`$, for which we have chosen a normalization of $`10^{54}`$ erg, since for the sample of GRBs with known redshifts given in Frail et al. 2001, one finds $`<E_{\gamma ,\mathrm{iso}}>=3\times 10^{53}`$ erg and $`<E_{\gamma ,\mathrm{iso}}^{11/6}>^{6/11}=5\times 10^{53}`$ erg. Since the true energy inferred from the observations, which is constrained by the isotropic equivalent energy $`E_{\gamma ,\mathrm{iso}}`$ measured in GRBs with known redshifts, is inversely proportional to the beaming factor, $`N_Rf_b^{5/6}`$. This means that the number of orphans expected to be detected in a flux limited survey decreases with increasing beaming factor, unless the survey is sufficiently deep to detect all sources. This is, of course, due to the fact that larger beaming factor implies, for given isotropic equivalent energy, a lower true energy, and hence a lower radio flux at the transition to spherical expansion, as demonstrated by Eq. (3).

The number of orphan radio GRBs, Eq. (11), depends on the model parameters $`E`$, $`n`$, $`\xi _e`$ and $`\xi _B`$. Afterglow observations strongly suggest that $`\xi _e`$ is close to equipartition, i.e. $`\xi _e1/3`$ (Freedman & Waxman 2001), and that $`E`$ is narrowly peaked around $`E_{51}=1`$ (Frail et al. 2001). The main uncertainty is thus in the parameters $`n`$ and $`\xi _B`$. The peak flux of a GRB afterglow, for an observer lying along the jet axis, is proportional to $`E_{\mathrm{iso}}\sqrt{n\xi _B}`$, and for typical luminosity distance of $`3\times 10^{28}`$ cm it is $`10\sqrt{\xi _Bn_0/10^3}E_{\mathrm{iso},54}`$ mJy (Waxman 1997b ; Gruzinov & Waxman (1999); Wijers & Galama (1998)). Observed afterglow fluxes generally imply $`\xi _Bn_010^3`$ for $`E_{\mathrm{iso},54}0.1`$, and values $`\xi _Bn_01`$ are obtained in several cases (e.g. Waxman 1997b ; Wijers & Galama (1998); Panaitescu & Kumar (2001)). Our parameter choice in Eq. (11), $`(\xi _B/0.3)=0.1`$ and $`n_0=0.1`$, is therefore conservative. Even for such a conservative choice of parameters we expect more than a dozen radio afterglows to be detected by an all sky survey like the NVSS, provided the beaming factor is not much larger than that anticipated from other considerations. Note, that our choice $`\dot{n}=0.5`$ is also conservative: If the GRB rate does not evolve with redshift, rather than following the star formation rate, the $`z=0`$ rate is about 8 times larger, $`\dot{n}4`$, and the expected number of radio GRB remnants is correspondingly larger by the same factor. Finally, the brightest source should have a flux of

$$f_{\nu max}30(500f_b)^{2/3}(\dot{n}/0.5)^{2/3}\left(\frac{\xi _e}{0.3}\right)\left(\frac{\xi _B}{0.03}\right)^{3/4}n_1^{19/36}E_{51}^{11/9}\nu _9^{1/2}(\tau _i/3)^{7/30}\mathrm{mJy}.$$
(12)
Since the flux in this phase declines roughly as $`(t/t_{SNT})^1`$, such a source should remain above 1 mJy for at least several years. This motivates follow-up observations of the brightest transients found in a survey.

The derivation of eq. (11) does not take into account the contribution from sources having a lifetime shorter than $`\tau _it_{SNT}`$. These sources are expected to be beamed, and so the fraction that can be observed (those directed along our sight line) should be of course correspondingly smaller. On the other hand, the radio flux emitted at times earlier than $`\tau _it_{SNT}`$ is larger and, consequently, these beamed sources can be seen out to a larger distance. We may estimate the number of beamed radio sources anticipated to be detected by the survey using the following argument. The ratio of the number of sources having an opening angle $`\theta `$ that are pointing in our direction, to the number of sources that had just become isotropic, is given by $`(2\pi \theta ^2/4\pi )[T(\theta )/t_{SNT}]`$ for a two sided jet, where $`T(\theta )`$ is the observed expansion time to opening angle $`\theta `$: $`T(\theta )\theta ^2`$ (e.g. Rhoads (1999)). Noting that the radio flux measured by an observer within the beam drops as $`t^{1/3}`$ during the expansion of the jet (e.g. Rhoads (1999)), the flux of a beamed source is higher than that of a source that just became spherical (lying at the same distance) by a factor $`\theta ^{2/3}`$. This implies that beamed sources can be seen out to a distance larger by a factor $`\theta ^{1/3}`$ compared to those that just became spherical. Thus, the ratio of the number of beamed sources to those that just became spherical, observed above a certain flux, is $`(\theta ^{1/3})^3(\pi \theta ^2/4\pi )[T(\theta )/t_{SNT}]\theta ^3`$. Thus, the number of detectable sources is dominated by sources which have already expanded to spherical symmetry.

For deeper surveys cosmological effects cannot be ignored and must be taken into account. We have solved eqs. (6) and (7) numerically for different choices of cosmological parameters. For the range of parameters presently favorable we find that $`N_R`$ depends only weakly on the assumed cosmology. In the following we present the results obtained for $`h=0.75`$, $`\mathrm{\Omega }_M=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. We also made the popular assumption that the comoving GRB density distribution traces the star formation rate, and invoked a strong redshift evolution for the later, as suggested by Madau & Pozetti (2000). Specifically, we take $`\mathrm{\Phi }(z)=(1+z)^3`$ in the following calculations.

The solution of eq. (6) is depicted in fig 1a, where the dependence of the maximum redshift $`z_{max}`$ on $`\eta /\tau _i`$ is plotted. Fig. 1b shows the dependence of $`N_R/N_0`$ on $`\eta `$ for different values of $`\tau _i`$.

For the cosmology adopted above, the radio flux of a source at a redshift of $`z=2`$ is given by

$$f_\nu 50\left(\frac{\xi _e}{0.3}\right)\left(\frac{\xi _B}{0.3}\right)^{3/4}n_0^{3/4}E_{51}\nu _9^{1/2}\left(\frac{t}{t_{SNT}}\right)^{9/10}\mu \mathrm{Jy}$$
(13)
and is smaller by a factor of 4 at $`z=3`$. Thus for our above (conservative) choice of afterglow parameters a flux limit of $`1\mu Jy`$ is required to see all the orphans if the GRB population is dominated by those at $`z2`$.

## 4 A Search for NVSS/FIRST Radio Transients

### 4.1 Analysis

We conducted a search for point-like radio transients in the FIRST radio catalog (White et al. 1997). We searched for compact FIRST radio sources that do not appear in the NVSS radio catalog (Condon et al. 1998). The sample we obtain this way is incomplete, since it rejects afterglow sources that were detected by the NVSS, and were bright enough to remain above 6 mJy during the FIRST observation. We examine this incompleteness factor further below and show that it is smaller than 5 percent. As the FIRST and NVSS surveys were obtained in different VLA configurations, the comparison was done carefully. We use this comparison to put an upper limit on the number of GRB radio afterglows.

The FIRST is a 20cm radio survey conducted using the NRAO VLA in B-configuration, with $`5^{\prime \prime }`$ beam-size. The final maps have a typical noise RMS of $`0.15`$ mJy. This value allows $`5\sigma `$ detection of $`1`$ mJy sources. The survey’s resolution enables an astrometric accuracy better than $`1\mathrm{"}`$ ($`90\%`$ CL). The FIRST radio survey has begun in 1994, and when completed it will cover $`10,000`$ deg2, mostly in the north Galactic cap. In the search described below, we used the FIRST 2001, Oct 15 version, covering $`8,565`$ deg2. This catalog has $`771,076`$ sources, hence, a density of about $`90`$ sources per deg2.

The NVSS 20cm radio survey (Condon et al. 1998) was conducted between 1993 to 1996 and covers all the sky north of declination $`40`$ deg. It utilized the VLA D- and DnC-configurations, with a circular beam-size of $`45^{\prime \prime }`$. This is significantly larger than the median angular size ($`10^{\prime \prime }`$) of faint extragalactic sources. With this beam-size, the faintest NVSS sources have an astrometric accuracy better than $`7\mathrm{"}`$ ($`68\%`$ CL). The images have RMS brightness fluctuations of $`0.45`$ mJy beam-1, that allows for $`5\sigma `$ detection of $`2.5`$ mJy sources. The NVSS $`99\%`$ completeness level is about $`3.4`$ mJy. With density of $`\rho 60`$ sources per deg-2, it has total of $`2\times 10^6`$ radio sources.

Radio sources could be intrinsically variable. Therefore, using the NVSS $`99\%`$ completeness of $`3.4`$ mJy as our threshold, could introduce a large number of variable sources (e.g., that were detected in the FIRST survey, but were just under detection in the NVSS) into our transients list. In order to compare the two catalogs, we need to choose a flux limit that rejects most of the variable objects. We plot in Figure 2 all FIRST point-like sources (above $`4`$ mJy) with NVSS counterparts. The solid line shows the location of sources for which the FIRST flux equals the NVSS flux, and the horizontal dashed line shows the NVSS, $`99\%`$ completeness flux level (i.e., $`3.4`$ mJy). Based on this plot we chose a threshold of $`6`$ mJy, for which the scatter around the solid line does not produce large number of false transients sources.

In the following we assumed that all the FIRST images, were taken after the corresponding NVSS images. This assumption will be examined more carefully later. Unless a GRB afterglow happens in a radio galaxy, it is expected to: (i) be a point source; and (ii) be absent from the NVSS. Since radio galaxies are a tiny fraction of the overall galaxy population, these assumptions hold for most of the afterglow population.

We selected all the FIRST point-like sources with peak flux above $`6`$ mJy ($`21.5\%`$). From this subsample, we selected all the unresolved objects ($`4.3\%`$ of the $`6`$ mJy subsample). However, some of the resolved FIRST objects could actually be point sources, classified as resolved due to measurements errors. Adopting the FIRST size error criterion (White et al. 1997), for $`9.2\%`$ of the FIRST sources above $`6`$ mJy it is not possible to reject the hypothesis that they are point-like at the $`99.7\%`$ CL. Therefore, the number of sources above $`6`$ mJy for which the point-source hypothesis could not be rejected at the $`3\sigma `$ level is $`2.1`$ ($`=9.2/4.3`$) times larger than the number of point sources we used in our sample. We ended with a catalog of $`7181`$ FIRST point-like sources with peak flux above $`6`$ mJy. For these sources, we searched the NVSS for counterpart within $`30\mathrm{"}`$ of the FIRST position. Although the NVSS astrometry is better than $`7\mathrm{"}`$, the use of such a high threshold was necessary because in cases of two (or more) nearby point-like radio sources, the NVSS could detect only one elongated source, with its center offset by more than $`7\mathrm{"}`$ from one of the FIRST point source positions. The search yields $`110`$ point-like FIRST sources with peak flux above $`6`$ mJy for which there is no NVSS counterpart.

We examined the POSS-II $`E`$ images and the FIRST and NVSS radio maps of all these transient candidates. We found that: $`2\%`$ result from holes in the NVSS maps; $`14\%`$ are artifacts due to un-cleaned aliases in the radio maps; $`61\%`$ are multiple sources for which the NVSS found an elongated source centered more than $`30\mathrm{"}`$ from one of the FIRST radio positions; and for the rest ($`23\%`$, $`25`$ sources) we could not reject the hypothesis that these are radio transients or variable radio sources of some kind (e.g., AGNs, SNe, GRBs, Pulsars). A large fraction of these ’sources’ are possibly: radio artifacts; NVSS non-detections due to locally high background (e.g., in the vicinity of bright radio source); or physical phenomena other than GRBs (e.g., AGNs, SNe). If we change the flux limit from $`6`$ to $`10`$ mJy, the number of candidates decreases from $`21`$ to $`6`$.

The following table lists the $`25`$ candidates along with their basic properties.

For each of these objects, we conducted a search for known counterparts within a radius of $`5^{}`$ in the NED database and within $`1^{}`$ in the HEASARC database. Our findings are listed in the table.

$`2575`$ deg2 of the FIRST survey were taken before the completion of the NVSS survey. Therefore, we reject all the sources in this area (i.e., $`22.2^{}<\delta <42.5^{}`$). Omitting this strip (PSR J1652+2651 is in this strip), the area covered by our search is therefore $`5990`$ deg2, and the number of sources is: $`16`$ and $`3`$, for the $`6`$ mJy and $`10`$ mJy flux thresholds respectively. The upper $`95\%`$ Poisson CL on these numbers are: $`26`$ and $`9`$ respectively. We consider these numbers to be conservative upper-limits on the number of possible GRB afterglows.

Finally, we have to correct our upper limits for the various completeness factors involved: (i) Missing point sources due to FIRST measurements errors. As we showed above, a very conservative estimate for this number is $`C_{missing}=9.2/4.33=2.1`$; (ii) Area completeness (e.g., the fraction of the sky covered by our survey) $`C_{area}=41252/5990=6.8868`$; (iii) Given The crude NVSS resolution, the probability that a transient source will be detected up to $`r=90^{\prime \prime }`$ (twice the NVSS FWHM), from a known un-related NVSS source, and therefore could be confused with it, is $`P(<r)=1\mathrm{exp}(\pi \rho r^2)`$ $`=11\%`$. Therefore, the source confusion completeness factor is: $`C_{conf}=1.11`$; (iv) NVSS completeness - For about $`2\%`$ of the FIRST sources there are no corresponding NVSS images (holes in the catalog), hence we set $`C_{comp}=1.02`$. (v) Incompleteness due to our search criteria - afterglow sources detected by the NVSS that were bright enough to remain above 6 mJy during the FIRST observation are rejected by our search. The number of such sources depends on the average time separation between the NVSS and FIRST observations. Denoting by $`t_{sep}`$ the time interval between two consecutive observations of a typical afterglow source, eq. (3) implies that the flux during the first observation was larger than during the second observation by a factor of $`(t_{sep}/t_{SNT}+1)^{9/10}`$, assuming the source was near maximum flux during the first detection (i.e., its age was roughly $`t_{SNT}`$). Consequently, the fraction of bright afterglow sources in the NVSS catalog that may be detected above 6 mJy by the FIRST at time $`t_{sep}`$ after 1996, must be smaller than $`\delta (t_{sep})=(t_{sep}/t_{SNT}+1)^{27/20}`$. The smallest time separation between the catalogs in the strips we consider is about $`t_{min}=6`$ months, and the largest is about $`t_{max}=4.5`$ years. Taking the FIRST detections to be uniformly distributed in time, integrating the fraction $`\delta (t_{sep})`$ from $`t_{min}`$ to $`t_{max}`$ yields a total fraction of less than 5% that our search rejected. Thus, we set $`C_{search}=1.05`$. Taking together all these factors, we can put a conservative $`95\%`$ CL upper limit of $`26\times C_{area}\times C_{missing}\times C_{conf}\times C_{comp}\times C_{search}=447`$ for the all-sky number of GRB afterglows with peak specific-flux above $`6`$ mJy (or $`154`$ for the $`10`$ mJy threshold).

### 4.2 Radio Supernovae

The characteristic behavior of radio supernovae (RSNe) (see Weiler et al. (2002) for a recent review) appears to be similar to the expected behavior of radio afterglows (RGRBs). In particular, the radio light curve of type II RSNe rises on a characteristic time scale of $`100`$ days (Weiler et al. (1998)), comparable to $`t_{SNT}`$ (see eq. ), and has an overall similar shape to that of RGRBs following the peak. This is not at all surprising since both are produced by the same mechanism, namely a blast wave, albeit under somewhat different conditions. This similarity would render any attempt to distinguish RGRBs from RSNe difficult and, therefore, motivates a careful examination of potential differences between these two classes of objects.

We identify two features that may allow one to distinguish RGRBs from RSNe. Firstly, the radio flux is anticipated to decline smoothly with a $`t^1`$ dependence during the sub-relativistic phase in RGRBs, while the radio light curves seen in RSNe are often much more complex. However, some RSNe do show a smooth decline, and those may be confused with RGRBs. Secondly, the absolute peak radio luminosity of RSNe may differ considerably from the radio luminosity of the brightest RGRBs. For a typical GRB the energy injected into the blast wave is $`E10^{51}`$ erg, and with the above choice of parameters: $`n=0.1\mathrm{cm}^3`$, $`\tau _i=3`$, $`\xi _B=0.03`$ and $`\xi _e=0.3`$, the RGRB luminosity density at 5 GHz is $`7\times 10^{28}\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1`$ at a time of $`0.5`$ yr after the explosion. This luminosity is $`100`$ times larger than the typical peak luminosity (of detected) RSNe, of order $`10^{27}\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1`$ (Weiler et al. (1998)), and several times brighter than the brightest known RSNe, SN1988z, which had a peak flux density of $`2\times 10^{28}\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1`$. This suggests that in a flux limited survey RSNe and RGRBs should have a distinct redshift distribution, with RSNe having preferentially smaller redshifts.

We estimate below the number of RSNe expected to be detected in a search of the type described in this paper for RGRBs, and compare their expected redshift distribution to that calculated for RGRBs. In order to make an accurate prediction for RSNe, their (peak) luminosity function and the shape of the corresponding light curves must be known. Unfortunately, due to the small sample of RSNe presently available (12 type II and 4 type I), and the lack of a supernova model that is capable of predicting the distribution of RSN properties, the determination of the peak luminosity function is highly uncertain. We therefore provide below a crude estimate based on the Weiler et al. (1998) compilation of RSNe radio data. Based on Figure 3 of Weiler et al., we adopt a RSNe peak luminosity function with equal number of RSNe per logarithmic peak luminosity interval over the 5 GHz peak luminosity range $`L_1=10^{23}\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1`$ to $`L_2=10^{28}\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1`$. In addition, we invoke a correlation between 5 GHz peak luminosity $`L_p`$ and peak time $`T_p`$ of the form, $`L_p=10^{28}(T_p/10^3\mathrm{day})^{1.4}\mathrm{erg}/\mathrm{sHz}`$, as demonstrated in Weiler et al. (1998). These assumptions hold, based on the available data, for type II RSNe. The properties of Type I RSNe are less certain. However, since type II RSNe dominate the (detected) RSNe population, and since the type I RSNe peak luminosities are close to the typical type II peak luminosities, we expect our estimate below of the type II RSNe “background” rate to provide a (crude) estimate for the total RSNe background.

Assuming that the RSNe flux drops as $`T_p/t`$ for $`t>T_p`$, the number of 5 GHz sources on the sky above a certain flux density $`f_{\nu ,\mathrm{min}}`$ at any given time can be readily computed:

$$N_{RSNe}=\frac{112\pi }{93\mathrm{ln}(L_2/L_1)}R_{RSNe}^3\dot{n}_{RSNe}T_2.$$
(14)
Here, $`\dot{n}_{RSNe}`$ is the RSNe rate per unit volume, $`T_2`$ is the peak time corresponding to peak luminosity $`L_2`$, and $`R_{RSNe}`$ is the maximum distance at which a RSN can be observed,

$$R_{RSNe}\left(\frac{L_2}{4\pi f_{\nu ,\mathrm{min}}}\right)^{1/2}=41\left(\frac{L_2/10^{28}\mathrm{erg}/\mathrm{sHz}}{f_{\nu ,\mathrm{min}}/5\mathrm{m}\mathrm{J}\mathrm{y}}\right)^{1/2}\mathrm{Mpc}.$$
(15)
Since the typical distances of detectable RSNe are $`10`$ Mpc, we have assumed Euclidean geometry for the calculation above. It is also straight forward to calculate the average distance of detected RSNe,

$$<R>_{RSNe}=\frac{93}{304}R_{RSNe}=13\left(\frac{L_2/10^{28}\mathrm{erg}/\mathrm{sHz}}{f_{\nu ,\mathrm{min}}/5\mathrm{m}\mathrm{J}\mathrm{y}}\right)^{1/2}\mathrm{Mpc}.$$
(16)
These distances should be compared with those obtained for RGRBs: for the conservative GRB parameters quoted above, the maximum and average distances of RGRBs detectable above 5 mJy at 5 GHz are $`R_{GRB}=140`$ Mpc and $`<R>_{GRB}=R_{GRB}/3=46`$ Mpc respectively.

Using the above equations, and taking $`\dot{n}_{RSNe}=1.2\times 10^5\mathrm{Gpc}^3\mathrm{yr}^1`$, corresponding to the total rate of core-collapse SNe (Filippenko (2001)), we find $`N_{RSNe}=3(f_{\nu ,\mathrm{min}}/5\mathrm{m}\mathrm{J}\mathrm{y})^{3/2}`$. Comparing this result with the estimates for RGRBs given in §2, we find that the number of RGRBs detected in a survey of the type described here should be at least comparable, and most likely significantly larger, than the number of detected RSNe. The characteristic distances of detectable RGRBs should be $`4`$ times larger than those expected for RSNe.

The anticipated proximity of RSNe in the FIRST survey implies that the host galaxy should be easily detectable. In fact, even very faint galaxies should be detectable above the POSS-II limit out to a distance of 100 Mpc. This means that all the radio transients listed in table 1 for which no optical counterpart is detected are unlikely to be RSNe, as even dwarf hosts should be clearly evident.

### 4.3 Radio loud AGNs

Radio loud AGNs are known to exhibit large amplitude variations over a broad range of timescale and frequencies, with flare durations ranging from minutes to years. The characteristic durations of radio outbursts (with the exception of ID variability) are between weeks to years (Valtaoja et al. (1999); Venturi et al. (2001)), with amplitudes that occasionally exceed a factor of 10. Thus, it is conceivable that some of the transients in our selected sample are radio loud AGNs. Indeed, our initial list of 110 transients contains a number of such sources.

Most of the candidates in table 1 have no optical counterpart brighter than $`M_R=21`$. What is the likeliness that some or most of them are in fact AGNs? One way to characterize the radio loudness is in terms of the radio-to-optical luminosity ratio, $`R`$, which is commonly defined using radio and optical luminosities at some fiducial rest-frame wavelengths. Existing samples of AGNs exhibit a distribution in radio loudness that depends on selection criteria. The distribution of $`R`$ in the FBQS (White et al. (2000)) for example peaks at $`R10`$, with very few sources (0.3 %) having $`R>10^4`$. The distribution of radio selected quasars has a median at $`R10^3`$ (Wadadekar & Kembhavi, (1999)). It is not clear at present whether this is a result of a selection effect, and whether there exists a population of extremely radio loud AGNs. However, based on our current knowledge, we would naively expect that deep enough observations at optical wavelengths would be able to ultimately reveal the underline AGN. future X-ray detections of some of our candidates will also imply that these are likely AGNs. A precise determination of the minimum value of $`R`$ for the objects in table 1 requires a knowledge of their redshifts, which are not available. If $`R`$ is defined as in Wadadekar & Kembhavi (1999) then a rough estimate yields $`R>300`$ for those objects. In addition, they should have no extended structure above $`1mJy`$ on scales larger than a few arcseconds.

Another potential difficulty for future surveys may arise due to the presence of faint, BL Lac type sources. If a population of radio BL Lacs exist for which the continuum flux associated with the nuclear activity drops below that of the host galaxy during quiescent periods, then a radio transient associated with such an object may be confused with a RSN or orphan afterglow if optical observations were made when the AGN was in a low state. An example is the source SDSS J124602.54 which has been proposed as a candidate optical orphan afterglow, and later confirmed to be a BL Lac (Gal-Yam et al. (2002)). In order to select out such sources in future surveys, monitoring in the radio and/or optical bands may be necessary.

## 5 Discussion

We have calculated the number of orphan radio afterglows expected to be detected in a flux limited survey. For the GRB parameters derived by Frail et al. (2001); total jet energy $`E10^{51}`$ erg and average beaming factor $`<f_b^1>(\theta ^2/2)^1500`$, the number of sources expected to be present over all sky at any given time with flux exceeding $`f_{\nu \mathrm{min}}`$ is \[see Eq. (11)\]

$$N_R20\left(\frac{\dot{n}}{0.5\mathrm{Gpc}^3\mathrm{yr}^1}\right)\left(\frac{f_{\nu min}}{5\mathrm{m}\mathrm{J}\mathrm{y}}\right)^{3/2}\left(\frac{\xi _B}{0.03}\right)^{9/8}n_1^{19/24}.$$
(17)
Here $`\dot{n}`$ is the GRB rate at redshift $`z=0`$, $`n=0.1n_1\mathrm{cm}^3`$ is the number density of the ambient medium into which the blast wave expands, and $`\xi _B`$ is the energy density of the magnetic field behind the shock in units of the equipartition value. $`\dot{n}`$ is normalized in Eq. (17) to the value obtained under the assumption that the GRB rate follows the redshift evolution of star formation rate (Schmidt (2001)). If the GRB rate does not evolve with redshift, $`\dot{n}`$ is a factor of 8 larger than the value used in Eq. (17). Since afterglow observations generally imply $`\xi _Bn10^3\mathrm{cm}^3`$, and values $`\xi _Bn1\mathrm{c}\mathrm{m}^3`$ are obtained in several cases \[see discussion in the paragraphs following Eq. (11)\], our choice of $`n`$ and $`\xi _B`$ values are conservative. Even for such a conservative choice of parameters we expect $`20`$ radio afterglows to be detected by an all sky survey with $`f_{\nu \mathrm{min}}5`$ mJy. These sources are detected out to a distance

$$d_{\mathrm{max}}200n_1^{3/8}\left(\frac{\xi _B}{0.03}\right)^{3/8}\left(\frac{f_{\nu min}}{5\mathrm{m}\mathrm{J}\mathrm{y}}\right)^{1/2}\mathrm{Mpc}.$$
(18)
We conducted a search for point-like radio transients with flux densities larger than 6 mJy using the FIRST and NVSS surveys, and discovered 25 radio afterglow candidates, listed in Table 1. The two main types of sources that may produce radio transients, which may be confused in our analysis with a radio afterglow, are radio supernovae (RSNe) and radio loud AGNs. While we have shown that our candidate sources are unlikely to be RSNe (see sec. §4.2), the possibility that most of them are radio loud AGNs cannot be ruled out (see sec. §4.3). Thus, our detected sources allow to set an upper limit to the number of orphan radio afterglows. Correcting for various completeness factors (see sec. §4.1), we obtained an upper limit of 447 all sky radio afterglows above 6 mJy, and 154 above 10 mJy.

We have shown that for a given isotropic equivalent burst energy, $`E_{\mathrm{iso}}=f_b^1E`$, the number of afterglows detected in a flux limited survey is smaller for larger beaming factor $`f_b^1`$, i.e. for smaller jet opening angle $`\theta `$ \[see Eq. (11)\]. The upper limit we derived on the number of radio afterglow sources therefore implies a lower limit for the beaming factor,

$$f_b^110\left(\frac{\dot{n}}{0.5\mathrm{Gpc}^3\mathrm{yr}^1}\right)^{6/5}\left(\frac{\xi _B}{0.03}\right)^{27/20}n_1^{19/20}.$$
(19)
Here we have used Eq. (11) with $`<E_{\mathrm{iso}}^{11/6}>^{6/11}=5\times 10^{53}`$ erg, based on GRBs with known redshifts (Frail et al. 2001). This result is consistent with the value $`<f_b^1>=500`$ inferred by Frail et al. (2001). Better determination of the distribution of $`n`$ and $`\xi _B`$, based on afterglow observations of identified GRBs, will improve the constraint on beaming. Alternatively, if the beaming factor is determined independently using other methods, then eq. (19) can be used to constrain the afterglow parameters. For instance, invoking $`<f_b^1>=500`$ implies

$$\left(\frac{\xi _B}{0.03}\right)^{27/19}n<5\left(\frac{\dot{n}}{0.5\mathrm{Gpc}^3\mathrm{yr}^1}\right)^{24/19}\mathrm{cm}^3.$$
(20)
In reality the afterglow parameters in a sample of radio GRBs should have a spread, and so the letter condition may be taken as a constraint on the mean.

Additional observational efforts aimed at identifying the candidates in table 1 are strongly motivated by two arguments. First, some of the candidates may turn out to be orphan afterglows. The fact that most candidates in table 1 have no optical counterparts brighter than the POSS-II limit, implies that if some of them are indeed GRB radio orphans then their host galaxies are not as bright as would have been expected under the assumption that GRBs are associated with star-forming regions, given the distances estimated in Eq. (18) for the conservative values of $`n`$ and $`\xi _B`$. Nevertheless, the observed number of orphan candidates is consistent with larger values of these parameters, as shown in Eq. (20), for which the distances are larger and the absence of optical counterparts is not surprising. Second, further observations are also motivated by the fact that rejection of sources from our list in Table 1 as possible candidates would impose a more stringent constraint for the beaming factor. Rejection of all the sources in Table 1 would imply an all sky upper limit of 49 sources, and a corresponding beaming factor $`f_b^1>100`$.

Future radio surveys may go deeper than the FIRST and would be of comparable or better quality. The Allen Telescope Array (formerly the ”1-hectar radio telescope”) will allow a full coverage of the northern sky every night with $`1`$ mJy sensitivity222see http://astron.berkeley.edu/ral/, and the square-kilometer-array may provide even better sensitivity333http://www.nfra.nl/skai/. A search for point-like radio transients similar to that outlined in §4.1, would probably yield a large number of candidates. The analysis presented in §3 predicts that the number of all sky radio GRBs above 1 mJy is likely to lie in the range between a few hundreds and a few thousands, depending on afterglow parameters \[see Eq. (17)\], provided the average beaming factor of GRBs is not significantly larger than currently estimated from other considerations. Since the detectable orphan population is dominated by $`1`$ yr old GRB remnants, which have just undergone the transition to sub-relativistic spherical expansion, these sources should be transient, exhibiting $`1/t`$ flux decline on time scale of months. Expansion at the speed of light over $`1`$ yr implies that the angular size of the nearby remnants, lying at a distance $`100`$ Mpc, is $`1`$ mas. The nearby remnants may therefore be resolved using the VLBA and, moreover, since they have just undergone the transition to spherical expansion, their structure may show significant anisotropy (see Paczyński 2001).

RSNe and radio loud AGNs would probably contribute a large number of transients, and may render the identification of radio afterglows difficult. As shown in §4.2, RSNe and RGRBs can, in principle, be distinguished statistically by their redshift distributions. At the mJy level, the average distance of the radio afterglows is expected to lie between 0.1 to 1 Gpc, depending on afterglow parameters \[see Eq. (18)\], so it should be possible to detect the host galaxies and obtain their redshifts. RSNe should typically lie at significantly shorter distances (see §4.2). Detection of such bimodality in redshift distribution will greatly assist in identifying radio GRBs. Furthermore, the relativistic expansion of the GRB remnants also implies, as explained above, that unlike RSNe they may be resolved with the VLBA, possibly revealing anisotropic structure. Efficient rejection of flaring blazars would require follow-up observations in the radio, optical and X-ray bands. In particular, the orphan radio afterglows should exhibit a $`t^1`$ decay of the radio flux.

Finally, it is conceivable that GRBs emit substantial amounts of energy in the form of gravitational waves and high-energy neutrinos (Waxman & Bahcall (1997); Rachen & Mészáros (1998); Waxman & Bahcall (1999, 2000); Mészáros, & Waxman (2001); see Waxman (2001) for review). Association of orphan radio afterglows with potential signals in gravitational wave and neutrino detectors may help increasing the signal-to-noise ratio. In a model proposed recently (Van Putten (2001); Van Putten & Levinson (2002)), the GRB results from extraction of the rotational energy of a Kerr black hole through its interaction with a magnetized torus. This model predicts that the major fraction of the rotational energy available will be emitted in the form of gravitational waves in the band 0.5 - 1.5 KHz, for a reasonable range of black hole masses. It further proposes that measuring the product of the total energy carried by the gravitational waves and the corresponding frequency, which reflects the compactness of the system, can provide an existence test for Kerr black holes. However, such a measurement requires a knowledge of the distance to the source. Since the gravitational wave emission is isotropic, the chance of associating a LIGO source with a GRB event is small if GRBs have a beaming factor as large as currently believed. Nevertheless, if this model is correct, then LIGO sources should be associated with the orphan radio afterglows that should appear several months after the burst of gravitational radiation. At the distance within which LIGO can detect the predicted signal ($``$ 100 Mpc), the radio source is expected to be very bright (see eq. ) and can be monitored for at least several years. Even though the angular resolution of the gravitational wave telescopes is expected to be quite poor, the chance that such a bright transient would coincide with a LIGO source should still be very small, and so such an association would be statistically significant. Our analysis for instance yielded an upper limit of about $`10^2`$ transients per square degree above 10 mJy. The association of a radio afterglow with a LIGO source would also provide a means to identify the host galaxy, given the superior resolution of the radio telescopes, and hence a distance measurement for the LIGO source.

We thank K. Weiler, A. Loeb B. Paczyński and R.D. Blandford for useful discussions. EW is partially supported by AEC grant and a Minerva grant, and is an incumbent of the Beracha Career development Chair.


format: text

text:
# Optimal Control of One-Qubit Gates

## 1 Introduction

Recent advances in experimental physics allowing for manipulation and measurement of single quantum systems, have stimulated a flurry of investigations on the control of quantum systems, and more or less formal schemes have been advanced. The conditions under which a given quantum system is completely controllable have been explored  and some limits of quantum controllability  have been found. Quantum control theory have several important applications including quantum state engineering , control of chemical reactions , laser cooling of molecular degrees of freedom , quantum register initialization  and the fabrication of robust quantum memories .

A major application of the theory of quantum control is the subject of quantum computation. The physical implementation of a quantum computer is a major challenge and many proposals including ion traps, optical cavities, and quantum dots have been advanced. Promissing practical implementations should be scalable, facing the problem of heat dissipation which gets worse along with the shrinking of the size of the proposed physical system. In this work we address the problem of carrying an initial qubit (more precisely of an initial Bloch vector) to an specified final qubit in a given amount of time using the minimum amount of energy possible. This paper is organized as follows: first we review the controllability of Bloch vector, and then we formulate and solve the problem of optimal control under the criteria of minimum energy, minimum energy derivative and a combination of both.

Two level systems adequately model many physical systems (spin 1/2, photon polarization, atoms in (quase) monochromatic electromagnetic fields, etc), despite its simplicity. A general 1-qubit gate can be represented as a 2-level quantum system, and the most general Hamiltonian for such a system can be written as $`H=h_0I+\stackrel{}{h}\stackrel{}{\sigma }`$, where $`\stackrel{}{\sigma }`$ are Pauli´s matrices, $`h_0`$ determines the zero energy reference, and $`\stackrel{}{h}`$ is a classical vector. Since we can use any $`\stackrel{}{h}`$ we want, and this is the field we use to control our system, here on we call $`\stackrel{}{h}`$ our vector control. We moreover assume $`h_0=0`$.

## 2 Bloch vector Controllability

All of the information of the quantum state of a two level system is completely determined by its density matrix $`\rho `$, or equivalently by its Bloch vector $`\stackrel{}{s}(t)=\frac{1}{2}\mathrm{Tr}(\rho \stackrel{}{\sigma })`$. The dynamics of Bloch vector, given by well known Bloch equation $`\dot{\stackrel{}{s}}(t)=\stackrel{}{b}\times \stackrel{}{s}`$, where $`\stackrel{}{b}=2\stackrel{}{h}/\mathrm{}`$, can be put in the more explicit way $`\dot{\stackrel{}{s}}=(b_x𝒥_x+b_y𝒥_y+b_z𝒥_z)\stackrel{}{s}`$, where the $`𝒥`$s are the rotation generators. In this case we can formulate the problem of taking an initial state $`\rho _i`$ to a final state $`\rho _f`$ in a specified amount of time $`T=t_ft_i`$. After rescaling, we take $`t_i=0,t_f=1`$. Remember that we assume that all three components of $`\stackrel{}{h}`$ (or of $`\stackrel{}{b}`$) are control fields. Thus, since we have all three rotation generators, the system is completely controllable in the sense that every rotation can be reached from the identity . In other words, any final vector can be reached from any other initial vector (of the same length, or the same degree of mixture) in any finite time, provided there are no constraints on the size of the control fields. Moreover, the motion equation can be inverted  to give $`\stackrel{}{b}=\stackrel{}{s}\times \dot{\stackrel{}{s}}s^2f(t)\stackrel{}{s}`$, where $`f(t)`$ is an arbitrary function, which shows that not only the control problem is solvable, but that it is solvable even when a path to realize such control is given. Moreover, even in this event the solution is not unique.

## 3 Optimal Control

As mentioned before, the larger the fields used for control the greater the amount of heat to be dissipated. Then, proved the complete controllability of the one qubit gate, it is meaningful asking which is the control vector field required to perform an arbitrary rotation operation with the minimum amount of energy spent. We address this problem an optimal “classical” control problem for the Bloch vector, that is we shall extremize the cost functional

$$S=_{t_0}^{t_f}𝑑t\left\{\frac{1}{2}\stackrel{}{b}\stackrel{}{b}+\stackrel{}{\lambda }\left(\dot{\stackrel{}{s}}\stackrel{}{b}\times \stackrel{}{s}\right)\right\},$$
(1)
where $`\stackrel{}{b}=2\stackrel{}{h}/\mathrm{}`$, and $`\lambda `$ is a (vector) Lagrange multiplier, and the momentum corresponding to $`\stackrel{}{s}`$. Notice that the form of the cost functional above is not arbitrary: were the vector control field a magnetic field acting on a spin half system, or an electromagnetic field for a charged two level system, the (electro)magnetic energy would have had the assumed form. For the sake of simplicity we scale the variables to get the following cost functional

$$S=_0^1𝑑t\left\{\frac{1}{2a}\stackrel{}{b}\stackrel{}{b}+\stackrel{}{\lambda }\left(\dot{\stackrel{}{s}}\stackrel{}{b}\times \stackrel{}{s}\right)\right\},$$
(2)
where the constant $`a`$ have been introduced to have $`\stackrel{}{s}`$ and $`\stackrel{}{\lambda }`$ in the same units.

Not all of the resulting Euler–Lagrange equations

$$\stackrel{}{b}=a\stackrel{}{s}\times \stackrel{}{\lambda },\dot{\stackrel{}{s}}=\stackrel{}{b}\times \stackrel{}{s},\dot{\stackrel{}{\lambda }}=\stackrel{}{b}\times \stackrel{}{\lambda },$$
(3)
are true dynamical equations: the first equation is a constraint. A few simple calculations show that the system of equations (3) possesses the following constants of motion $`s^2=\stackrel{}{s}\stackrel{}{s}`$, $`\lambda ^2=\stackrel{}{\lambda }\stackrel{}{\lambda }`$ and $`\nu =\stackrel{}{\lambda }\stackrel{}{s}=\lambda s\mathrm{cos}(\theta )`$, where $`\theta `$ is the angle between $`\stackrel{}{\lambda }`$ and $`\stackrel{}{s}`$. The equations of motion can be put in the equivalent form

$$\frac{d}{d[at]}\left(\begin{array}{c}\stackrel{}{s}\\ \stackrel{}{\lambda }\end{array}\right)=\left(\begin{array}{cc}\stackrel{}{\lambda }\stackrel{}{s}& s^2\\ \lambda ^2& \stackrel{}{\lambda }\stackrel{}{s}\end{array}\right)\left(\begin{array}{c}\stackrel{}{s}\\ \stackrel{}{\lambda }\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}(\theta )& 1\\ 1& \mathrm{cos}(\theta )\end{array}\right)\left(\begin{array}{c}\stackrel{}{s}\\ \stackrel{}{\lambda }\end{array}\right),$$
(4)
whose solution reads

$$\left(\begin{array}{c}\stackrel{}{s}\\ \stackrel{}{\lambda }\end{array}\right)=\frac{1}{\mathrm{sin}(\theta )}\left(\begin{array}{cc}\mathrm{sin}(\theta at\mathrm{sin}\theta )& \mathrm{sin}(at\mathrm{sin}\theta )\\ \mathrm{sin}(at\mathrm{sin}\theta )& \mathrm{sin}(\theta +at\mathrm{sin}\theta )\end{array}\right)\left(\begin{array}{c}\stackrel{}{s}(0)\\ \stackrel{}{\lambda }(0)\end{array}\right).$$
(5)
We have assumed that $`\stackrel{}{s}`$ and $`\stackrel{}{\lambda }`$ are unitary vectors. Notice that we have obtained an orthogonal transformation. In fact this could have been anticipated by showing that $`\stackrel{}{b}`$ is constant both in norm and direction. Let us see that the solution is valid for any value of $`\theta `$, so we have freedom to choose $`\theta `$ to our best convenience. For example, we can choose $`\stackrel{}{\lambda }(0)=\stackrel{}{s}(1)`$. It follows $`\mathrm{cos}\theta =\stackrel{}{s}(0)\stackrel{}{s}(1)`$, and $`\mathrm{sin}(a\mathrm{sin}(\theta ))=\mathrm{sin}(\theta )`$. This choice leads to the solution

$$\stackrel{}{s}(t)=\frac{\mathrm{sin}(b(1t))}{\mathrm{sin}(b)}\stackrel{}{s}_i+\frac{\mathrm{sin}(bt)}{\mathrm{sin}(b)}\stackrel{}{s}_f,$$
(6)
where

$$\stackrel{}{b}(t)=\stackrel{}{b}(0)=\frac{\mathrm{ArcSin}\left(\mathrm{sin}(\theta )\right)}{\mathrm{sin}(\theta )}\stackrel{}{s}_i\times \stackrel{}{s}_f=(\theta +2\pi n)\stackrel{}{s}_{}=b_1(t;n)\stackrel{}{s}_{}.$$
(7)
Equations (6) and (7) constitute the solution to the proposed problem, and agree with the result of reference . Let us see that in the scaled variables the magnitude of the field is essentially equal to the value of the angle between the initial and the final Bloch vectors. The solution corresponds to a rotation around $`()\stackrel{}{s}_{}`$ at constant speed, from the initial Bloch vector to the final one. We see then that the solution is not unique: due to the multivaluedness of the function ArcSin, there is an infinite number of solutions. Each solution corresponds to a local minimum, the global minimum being the solution whose magnitude is equal to the angle between the initial and the final Bloch vector. The other solutions correspond to an integer number of turns followed by this angle. Finally, the other solutions correspond to reaching the final Bloch vector from the initial one but rotating the other way. Choices different to $`\stackrel{}{\lambda }(0)=\stackrel{}{s}(1)`$ are also meaningful, but should lie on the plane which contains the initial and final Bloch vectors. Let us see that the solution found breaks down when the initial and final vector are antipodal. In this case there are an infinite number of solutions all of them spending the same amount of energy. In this case the choice of an initial Lagrange multiplier vector which is not (anti)parallel to the initial Bloch vector leads to a particular solution.

Now, we remark that the minimum energy control corresponds to minimum path on the Bloch sphere. The length $`l`$ transversed by the tip of the Bloch vector is given by

$$l=_0^1\left|\frac{d\stackrel{}{s}}{dt}\right|𝑑t=_0^1\left|\stackrel{}{b}\times \stackrel{}{s}\right|𝑑t=s_0^1b(t)|\mathrm{sin}(\theta (t))|𝑑ts_0^1b(t)𝑑t.$$
(8)
This means that, for a given magnitude of the control field, the transversed length is maximized when the control vector field makes a right angle with the Bloch vector. In other words, if we restrict to control fields that are perpendicular to Bloch vector, and if Bloch vector transverses a longer path, it is necessary to have a larger average control field. This remark allows us to restrict ourselves to consider the vector control to point out in the direction of $`\stackrel{}{s}_i\times \stackrel{}{s}_f`$. Making the Ansatz

$$\stackrel{}{b}(t)=b(t)\stackrel{}{s}_{},$$
(9)
after some algebra we obtain

$$\stackrel{}{s}(t)=\frac{\mathrm{sin}(\theta _0^tb(t^{})𝑑t^{})}{\mathrm{sin}(\theta )}\stackrel{}{s}_i+\frac{\mathrm{sin}(_0^tb(t^{})𝑑t^{})}{\mathrm{sin}(\theta )}\stackrel{}{s}_f,$$
(10)
were $`b`$ should satisfy the equality $`_0^1b(t)𝑑t=\theta `$. If we write $`b(t)`$ as $`\theta +\delta (t)`$ we see that the average value of $`\delta (t)`$ over the unitary interval is zero and that $`_0^1b^2(t)𝑑t=\theta ^2+_0^1\delta ^2(t)𝑑t`$. Thus, the minimum is attained when the control vector field is constant from the initial till the final time. It is possible to give an alternative argument which shows that the solution of minimum fluence is the same as the shortest geodesic. Inverting the Bloch equation we obtain the control field $`\stackrel{}{b}=\stackrel{}{s}\times \dot{\stackrel{}{s}}+f\stackrel{}{s}`$, and the energy spent in the control $`_0^1b^2(t)𝑑t=_0^1\left((\dot{\stackrel{}{s}})^2(\stackrel{}{s}\dot{\stackrel{}{s}})^2+f^2\right)𝑑t`$. Since the second subintegral term is identically zero, and the function $`f(t)`$ should be zero for the extrema, we see that the fluence minimization and geodesic minimization (see Eq. 8) are almost the same, and reach their extrema together. Had we chosen the squared length instead of the length, both expressions would have been identical.

Since the times to perform quantum computation are generally short in low dimensional condensed matter systems, which are the most promising candidates, one also should analize possible limitations set by the rate at which control fields can be set. In particular, it is worth paying attention that the solution (7) is a discontinuous one, zero before the initial time, constant between the initial and final times, and zero again from the final time on. Had we used the square of the time derivative of the control field instead of the square of the field itself, defining the cost functional

$$S=_0^1𝑑t\left\{\frac{1}{2\mathrm{\Omega }^2}\frac{d\stackrel{}{b}}{dt}\frac{d\stackrel{}{b}}{dt}+\stackrel{}{\lambda }\left(\dot{\stackrel{}{s}}\stackrel{}{b}\times \stackrel{}{s}\right)\right\},$$
(11)
the solutions obtained above would have been also solutions of the new problem. In this case a whole set of new solutions arise, which are of constant magnitude but whose direction changes with time. It is easy to construct such a kind of solutions. Let $`\{\stackrel{}{b}(t),\stackrel{}{s}(t)\}`$ a solution of the new problem, but with a time dependent $`b^2(t)`$, then $`\{\stackrel{~}{\stackrel{}{b}}(t)=\stackrel{}{b}(t)+f(t)\stackrel{}{s},\stackrel{}{s}\}`$, is also a solution, no matter how the function $`f`$ is chosen. In particular, we can adjust $`f`$ to obtain $`\stackrel{~}{b}^2`$ a constant. For instance, if we set $`\stackrel{}{s}(t)=\mathrm{cos}(\varphi (t))\stackrel{}{s}_0(t)+\mathrm{sin}(\varphi (t))\stackrel{}{s}_{}`$, where $`\stackrel{}{s}_0`$ is the solution for the problem of minimum fluence, and $`\mathrm{cos}(\varphi (t))`$ a function with value $`1`$ both at $`t=0`$ and at $`t=1`$, we have the control field $`\stackrel{}{b}(t)=\theta \mathrm{cos}^2(\varphi )\stackrel{}{s}_{}\dot{\varphi }\stackrel{}{s}_\tau \theta \mathrm{sin}(\varphi )\mathrm{cos}(\varphi )\stackrel{}{s}_0`$ where $`\stackrel{}{s}_\tau =\stackrel{}{s}_{}\times \stackrel{}{s}_0`$ is a unitary vector needed to define a time dependent right triade $`\{\stackrel{}{s}_0,\stackrel{}{s}_\tau ,\stackrel{}{s}_{}\}`$. One can choose $`f`$ as $`f=\pm \sqrt{B^2\theta ^2\mathrm{cos}^2(\varphi )\dot{\varphi }^2}`$, with $`B^2`$ the maximum value of $`b^2`$, so at instants where the maximum is attained, $`f`$ vanishes. For the sake of definiteness we use $`\varphi (t)=\theta \mu t(1t)`$ which yields $`f^2=\theta ^2(\mathrm{sin}^2(\varphi (t))+\mu ^2t(2t))`$, and produces a new constant norm vector control $`\stackrel{~}{\stackrel{}{b}}`$ with magnitude $`\theta \sqrt{1+\mu ^2}`$. This solution, of course, spends more energy than the found before to perform the same control.

Solutions with vanishing magnitude at the initial and final instants of time also exist. Observe that in this case the equation for $`\stackrel{}{b}`$ is

$$\ddot{\stackrel{}{b}}=\mathrm{\Omega }^2\stackrel{}{s}\times \stackrel{}{\lambda }$$
(12)
which is a true dynamical equation. This allows for some extra flexibility: now we can add initial and final conditions on the value of the control field. From the point of view of the energy injected to the system, the most physically sensible conditions are those of vanishing control field both at the initial and the final times. We observe that in this case the solution should follow the shortest geodesic between the initial and the final Bloch vectors. In fact, since the control field begins and ends with a vanishing value, it should grow and decrease as slowly as possible but fast enough as to reach the maximum value necessary to have an average magnitude of at least $`\theta `$. If it had grown to a greater value it would have grown at a larger pace, so it would have been not the minimum solution sought.

Differentiating the equation (12) we obtain $`\frac{d^3\stackrel{}{b}}{dt^3}=\stackrel{}{b}\times \frac{d^2\stackrel{}{b}}{dt^2}`$, which immediately tells us that the second derivative of the vector control has constant norm. The first integral of this third order equation, $`\frac{d^2\stackrel{}{b}}{dt^2}=\stackrel{}{b}\times \frac{d\stackrel{}{b}}{dt}\mathrm{\Omega }^2\stackrel{}{s}_i\times \stackrel{}{s}_f`$, where we have assumed $`\lambda (0)=\stackrel{}{s}_f`$, can be solved under the assumption that the control vector is a second degree polynomial in $`t`$, with the result $`\stackrel{}{b}(t)=\frac{\mathrm{\Omega }^2|\mathrm{sin}(\theta )|}{2}t(1t)\stackrel{}{s}_{}=b_2(t)\stackrel{}{s}_{}`$ . Numerical solution of these equations, without the ansatz made above also lead to the same solutions. Before expressing $`\stackrel{}{b}`$ completely in terms of $`\stackrel{}{s}_i`$ and $`\stackrel{}{s}_f`$, we proceed to discuss the more general physical criterion in which one is interested on energy saving but with a limited rate of change of the vector control, through the cost functional

$$S=_0^1𝑑t\left\{\frac{1}{2a}\left(\stackrel{}{b}\stackrel{}{b}+\frac{1}{\omega ^2}\frac{d\stackrel{}{b}}{dt}\frac{d\stackrel{}{b}}{dt}\right)+\stackrel{}{\lambda }\left(\dot{\stackrel{}{s}}\stackrel{}{b}\times \stackrel{}{s}\right)\right\}.$$
(13)
The experience gained with the previous examples shows that the solution control fied should point (anti)parallel to $`\stackrel{}{s}_{}`$. Some algebra leads to the solution

$$\stackrel{}{b}(t)=a|\mathrm{sin}(\theta )|\left(1\frac{\mathrm{cosh}\omega (t\frac{1}{2})}{\mathrm{cosh}(\frac{\omega }{2})}\right)\stackrel{}{s}_{}=b_3(t)\stackrel{}{s}_{}.$$

We notice that all of the solutions so far found have the form of equation (9) and therefore have the solution (10). We only have to take care of the final value of $`\stackrel{}{s}`$. This leads to the following more explicit forms for $`b(t)`$

$$b_2(t;n)=6(\theta +2\pi n)t(1t),b_3(t;n)=\frac{\theta +2\pi n}{1\frac{\mathrm{tanh}(\omega /2)}{\omega /2}}\left(1\frac{\mathrm{cosh}(\omega (t\frac{1}{2}))}{\mathrm{cosh}(\frac{\omega }{2})}\right).$$
(14)
For the second case considered, the intuitive choice, $`b=\pi \theta \mathrm{sin}(\pi t)/2`$ produces a value of the cost functional only 1.5% above that of the optimal solution. Finally, for control fields of the form of equation (9) the cost functional can be written in purely geometric terms. If we set $`\dot{\varphi }(t)=b(t)`$, $`\varphi (0)=0`$, then $`S`$ can be expressed as

$$S=\frac{1}{2a}_0^{\overline{b}}b(\varphi )\left(1+\left(\frac{1}{\omega }\frac{db}{d\varphi }\right)^2\right)𝑑\varphi ,$$
(15)
where $`\overline{b}`$ is the average magnitude of the control field, and $`\varphi `$ the accumulated angle (or the arc length) transversed by the Bloch vector. Equation (15), just like equation (13), contains the other two cases: the first in the limit $`1/\omega 0`$, and the second in the limit $`1/a0`$ but with $`a\omega ^2=\mathrm{\Omega }^2`$ fixed. Of course, $`b_1(t;n)=lim_{1/\omega 0}b_3(t;n)`$ and $`b_2(t;n)=lim_{1/a0,a\omega ^2=\mathrm{\Omega }^2}b_3(t;n)`$.

We have formulated and solved in an analytic way, the problem of rotation of the Bloch vector (which characterizes completely the state of a two level system) from a prescribed initial vector to a prescribed final vector, in a given amount of time, using an optimal control scheme which minimizes the energy spent by the control fields, or the magnitude of the rate of change of the control fields or a linear combination of both. We have found control fields perpendicular to both the initial and final Bloch vectors, and multiple local minima corresponding to arrival from the initial to the final Bloch vector in one or other senses or after one or more complete turns.

This work was partly funded by DIB-UN (División de Investigaciones, Sede Bogotá, Universidad Nacional).


text:
# Semiflexible Polymer in a Strip

## I Introduction

Polymers are often modeled as self- and mutually avoiding walks placed on a lattice, and much has been learned about their thermodynamic properties through such models . Grand-canonical models of this kind, where the number of monomers incorporated into the chain is allowed to fluctuate controlled by an activity $`z=\mathrm{exp}(\mu /k_BT)`$, display a phase transition at some value of the activity (for infinite chains, that is, in the polymer limit). This transition is discontinuous in the one dimensional case $`d=1`$ , and continuous for $`d>2`$. Rather precise estimates of the critical value of $`z`$ were found in two dimensions through transfer matrix calculations  and series expansions . Also, exact values for the critical exponents are available in this case .

More recently, properties of such models on lattices limited by walls have attracted much interest , following studies of magnetic models in the same situation . The short range interaction between the wall and the polymer may be introduced associating a Boltzmann factor $`\omega =\mathrm{exp}(ϵ/k_BT)`$ to each monomer placed on the wall, so that $`ϵ<0`$ corresponds to attracting walls while repulsive walls are described by $`ϵ>0`$. The grand-canonical partition function for a model of a single chain is

$$Y(x,\omega )=z^N\omega ^{N_w},$$
(1)
where $`N`$ is the number of monomers in the chain, the sum is over all configurations of the chain with the initial monomer placed on the wall, and $`N_wN`$ stands for the number of monomers located on the wall. Such a model shows interesting features, and even in the limit where the self avoidance constraint is neglected (the so-called ideal chains) one finds that for sufficiently large values of $`\omega >\omega _0`$ the surface polymerization transition will occur at a lower value of the activity $`z`$ than the one in the bulk . The point $`(z_0,\omega _0)`$ in the phase diagram where the bulk (also called ordinary) transition line meets the surface transition line is called the adsorption transition point. In two dimensions, such models have been studied through transfer matrix calculations  and series expansions . Additional walls may be added, confining the polymer inside a strip, slab or pore , and the force applied on the walls in such situations is of interest even from the point of view of applications of polymers as adhesives . The model of ideal chains confined in a slab has been studied in the past , and it was found that the force on the walls is attractive if $`\omega `$ exceeds the adsorption value $`\omega _0`$. In the case of self-avoiding chains confined in a strip on the square lattice, transfer matrix calculations show that attractive forces appear for $`\omega `$ below $`\omega _0`$ . Also, the profile of the monomer density inside the strip was obtained , and, unlike what happens for ideal chains, for a self-avoiding chain profiles which are neither convex nor concave are found for a range of values of $`\omega `$. Finally, similar techniques have also been applied to shed light on the scaling behavior of such models for strips with periodic boundary conditions (cylinders) as the width of the strips becomes large .

Another generalization of the original polymer model is to introduce an energy associated to the operation of bending the chain. On hypercubic lattices, the elementary bends will always be at right angles, and an energy $`ϵ_b`$ may be associated to each of them. This semiflexible polymer problem (also called persistent or biased walks), has been studied some time ago . Recently, this model has attracted renewed interest , since it may describe some relevant aspects in the protein folding problem . The thermodynamic properties of the model have been studied on the Bethe lattice  and, equivalently, in the Bethe approximation . The end-to-end distance of semi-flexible chains on Bethe and Husimi lattices was obtained . In this paper we study the thermodynamic behavior of a semiflexible polymer confined inside a strip. The partition function of the model may be written as

$$Y(x,\omega )=z^N\omega ^{N_w}\omega _b^{N_b},$$
(2)
where $`N_b`$ is the number of elementary bends in the configuration and $`\omega _b=\mathrm{exp}(ϵ_b/k_BT)`$ is the Boltzmann factor associated to each elementary bend and the sum is over all configurations of the chain. We define a transfer matrix for the model and obtain the grand-canonical partition function in the thermodynamic limit, determined by the largest eigenvalue of this matrix. In order to be able to obtain the distribution of monomers in the transverse section of the strip, we define a position dependent activity for the monomers. Also, the force applied by the polymer on the walls is calculated, as a function of the width of the strip, $`\omega `$ and $`\omega _b`$. All thermodynamic properties are calculated at the polymerization transition, that is, at the value of the activity of a monomer for which the number of monomers incorporated into the polymer diverges.

In section II the model is defined in detail and its solution is presented. The results we find for the thermodynamic behavior of the model are shown in section III as well as final conclusions and discussions.

## II Definition of the model and its solution

The self avoiding chain is constrained inside a strip of width $`m`$ defined on the square lattice in the $`(x,y)`$ plane, so that $`0xm`$. The chain runs through the whole strip, from $`y\mathrm{}`$ to $`y+\mathrm{}`$. We may define a transfer matrix for the problem following a prescription proposed by Derrida, in a way to take into account the self-avoidance constraint exactly . The connectivity properties of all vertical bonds of the chain arriving at a line $`y_0`$ from below are specified through the indication of

1. The (unique) bond connected to the initial monomer of the chain (placed in $`y\mathrm{}`$) through a path lying entirely below the line $`y_0`$ (passing though sites with $`y<y_0`$);
2. The pairs of bonds connected to each other through a path lying entirely below the line $`y_0`$.

In figure 1 the five configurations for the case $`m=2`$ are depicted, and a portion of the chain placed inside the strip is shown with each configuration indicated. Configurations 1 and 3, as well as 4 and 5, are related to each other by reflection symmetry. We define column dependent activities $`z_i`$, $`i=1,2,\mathrm{}n_a(m)`$ according to this reflection symmetry, as indicated also in figure 1. The number of activities $`n_a(m)`$ is equal to $`m/2+1`$ for even values of $`m`$ and $`(m+1)/2`$ for odd values of $`m`$.

For a fixed connectivity configuration of a set of $`m+1`$ vertical bonds arriving at the line $`y_0`$, the possible configurations of the set of vertical bonds arriving at $`y_0+1`$ may be obtained, as well as the contribution to the partition function from the sites comprised between both sets of vertical bonds. This contribution will be given by

$`\omega _b^{N_{b,y_0}}{\displaystyle \underset{i=1}{\overset{n_a(m)}{}}}z_i^{N_{i,y_0}}`$

where $`N_{i,y_0}`$ is the number of monomers with activity $`z_i`$ in line $`y_0`$ and $`N_{b,y_0}`$ is the number of elementary bends in this line. These contributions, shown for a particular example in figure 2, define a line of the transfer matrix. For $`m=2`$ the transfer matrix will be given by

$$𝐓=\left(\begin{array}{ccccc}z_2& z_1z_2\omega _b^2& z_1z_2^2\omega _b^2& z_1z_2^2\omega _b^2& 0\\ z_1z_2\omega _b^2& z_1& z_1z_2\omega _b^2& 0& 0\\ z_1z_2^2\omega _b^2& z_1z_2\omega _b^2& z_2& 0& z_1z_2^2\omega _b^2\\ 0& 0& z_1z_2^2\omega _b^2& z_1z_2^2& 0\\ z_1z_2^2\omega _b^2& 0& 0& 0& z_1z_2^2\end{array}\right).$$
(3)
The grand canonical partition function of the model,considering periodic boundary conditions in the $`y`$ direction, will be given by

$$\mathrm{\Xi }=(Tr)𝐓^{N_y},$$
(4)
where $`N_y`$ is the total length of the strip in the $`y`$ direction, so that the total number of sites is given by $`N_s=N_y(m+1)`$. The number of monomers with activity $`z_i`$ is

$$N_i=\frac{z_i}{\mathrm{\Xi }}\frac{\mathrm{\Xi }}{z_i},$$
(5)
and the total number of monomers in the chain will be

$$N=\underset{i=1}{\overset{n_a(m)}{}}N_i.$$
(6)
The fraction of monomers placed at column $`x`$ is

$$\rho (x)=\frac{N_{x+1}}{(2\delta _{x,m/2})N},$$
(7)
where $`x=0,1,\mathrm{},n_a(m)1`$ and $`\rho (x)=\rho (mx)`$. The Kronecker delta in the denominator contributes only for even values of $`m`$.

In the thermodynamic limit $`N_y\mathrm{}`$ the partition function 4 is dominated by the largest eigenvalue of the transfer matrix $`\lambda _1`$, so that $`\mathrm{\Xi }\lambda _1^{N_y}`$ and in this limit

$$N_i=N_y\frac{z_i}{\lambda _1}\frac{\lambda _1}{z_i}.$$
(8)
The first order polymerization transition in the finite strip will take place when the thermodynamic potential

$$\psi =k_BT\mathrm{ln}(\mathrm{\Xi })$$
(9)
is equal to the one for the empty lattice $`\psi _0=0`$. Thus, since in the thermodynamic limit we have

$$\psi /N_s=k_BT(m+1)\mathrm{ln}(\lambda _1)$$
(10)
the polymerized phase will coexist with the non-polymerized phase for $`\lambda _1=1`$. Therefore, all thermodynamic quantities below will be calculated for $`z_i=z,i=1,\mathrm{},n_a(m)1`$ and $`z(n_a(m))=\omega z`$, where the activity $`z`$ is then fixed at the coexistence value $`z_c`$ so that, for a given value of $`\omega _b`$, we have $`\lambda _1=1`$.

Finally, the force applied on the walls is given by

$$F=\frac{1}{a}\left(\frac{\mathrm{\Xi }}{m}\right)_{z,\omega ,\omega _b},$$
(11)
where $`a`$ is the lattice parameter and positive values correspond to attractive forces. An adimensional force per monomer at the coexistence may be then defined as

$$f=\frac{Fa}{k_BTN}=\frac{\left(\frac{\mathrm{\Xi }}{m}\right)_{z,\omega ,\omega _b}}{z\left(\frac{\mathrm{\Xi }}{z}\right)_{m,\omega ,\omega _b}}=\frac{1}{z_c}\left(\frac{z_c}{m}\right)_{\omega ,\omega _b}.$$
(12)
Since our results correspond to integer values of $`m`$, the force $`f`$ was estimated making the discrete approximation

$`f(m+1/2,\omega ,\omega _b)`$ $``$ $`{\displaystyle \frac{2}{z_c(m+1,\omega ,\omega _b)+z_c(m,\omega ,\omega _b)}}`$ (14)
$`\times [z_c(m+1,\omega ,\omega _b)z_c(m,\omega ,\omega _b)].`$
A simple calculation may be performed in the limit of rigid rods $`\omega _b=0`$, where bends are not allowed. In this limit the transfer matrix is diagonal of size $`(m+1)\times (m+1)`$. For $`\omega >1`$ all monomers are on the walls, for $`\omega =1`$ they are uniformly distributed, whereas for $`\omega <1`$ an uniform distribution in the sites away from the walls is found. As expected, in this limit the force on the walls vanishes. In the numerical results below we consider $`\omega _b1`$.

## III Numerical results and conclusion

The transfer matrices for the model were obtained exactly for strips of widths ranging from 3 to 9. After using the reflection symmetry, the sizes of the matrices are, respectively, equal to 6, 16, 38, 100, 256, 681, and 1805. We obtained the density profile at the coexistence condition for values of $`\omega `$ and $`\omega _b`$ mostly between 1 and 3.

For neutral walls ($`\omega =1`$) the density profile is always concave, with a higher density in the center. This may be understood since the region away from the walls is favored entropically. As $`\omega `$ becomes larger, monomers on the walls are energetically favored, so for sufficiently large values of $`\omega `$ a convex density profile is expected. For ideal flexible chains ($`\omega _b=1`$), at the adsorption value $`\omega =4/3`$ the density profile is flat for all sites which are not on the walls . Such a flat transition profile is not observed for self-avoiding chains, where convex profiles are separated from concave ones by an interval of values of $`\omega `$, located well below the adsorption transition value, where the profile is neither convex nor concave. This interval is quite narrow for flexible chains . In figure 3 the densities are plotted as functions of $`\omega `$ for two values of $`\omega _b`$. It is clear that as $`\omega _b`$ is increased, favoring bends, the interval of values of $`\omega `$ with a profile without well defined convexity grows.

The values of $`\omega `$ below which the density profile is concave ($`\omega _1`$)and those above which the profile is convex ($`\omega _2`$) are plotted in figure 4 as functions of $`\omega _b`$ for two values of $`m`$. As is apparent the range with no defined convexity grows with $`\omega _b`$ in a nearly linear way. As a general rule, we found that as $`\omega _b`$ is increased from 1, the first pair of densities to cross, destroying concavity, is always correspondent to the columns which are first and second neighbors to the walls. Also, the last crossing, which turns the profile convex, is between densities corresponding to the columns at the wall and the neighbor columns. In the limit $`\omega _b\mathrm{}`$ we found that $`\omega _2\mathrm{}`$.

Finally, we calculated the force per monomer on the walls for pair of widths $`m`$, $`m+1`$, using expression 14. As already found for flexible self-avoiding chains , attractive forces appear for sufficiently large values of $`\omega `$. Figure 5 shows results for the force as a function of $`\omega `$ for some values of $`\omega _b`$. The origin of attractive forces in the system are portions of the chain limited by monomers adsorbed on opposite walls. Thus, as expected, the curves $`f\times \omega `$ display a maximum, because the force vanishes as $`\omega \mathrm{}`$, since in this limit such “bridge” segments are absent.

In figure 6 the maximum force per monomer (with respect to $`\omega `$) is plotted as a function of $`\omega _b`$. It is apparent that the maximum of these curves is located at $`\omega _b>1`$, and thus we conclude that polymers for which bends are somewhat favored give rise, in general, to larger attractive forces than flexible ones. As $`\omega _b`$ is increased, the maximum attractive force occurs at higher values of $`\omega `$

The force per monomer as a function of the width of the strip $`m`$, for $`\omega `$ above a threshold, shows a stable equilibrium point at very low separation $`m_1`$ and an unstable equilibrium point at a larger separation $`m_2`$ . The force is attractive in the interval $`m_1<x<m_2`$.

The stable equilibrium point $`m_1(\omega ,\omega _b)`$ is found in the range $`0.5m_11.5`$, showing little variation as a function of the Boltzmann factors, as long as $`\omega `$ is larger than the threshold value. The unstable equilibrium point $`m_2(\omega ,\omega _b)`$ shows a rather strong monotonic dependence on both variables, being an increasing function of $`\omega `$ and a decreasing function of $`\omega _b`$.

In conclusion, we may summarize the behavior of semiflexible chains confined inside strips observing that as the presence of elementary bends in the chains is favored, the range of values of $`\omega `$ for which the density profile has no defined convexity grows. Also, the largest attractive forces are found for $`\omega _b>1`$, that is, when the bending of the chains is favored. It should be noticed that the nonmonotonic behavior observed for the tension as a function of the width of the strip, with two equilibrium points, happens in a regime where the finite size scaling behavior has not yet been reached .

## ACKNOWLEDGMENTS

Partial financial support by the brazilian agencies CNPq and FAPERJ is gratefully acknowledged.