text:
# GRB Formation Rates inferred from the Spectral Peak Energy – Luminosity Relation

## 1 Introduction

Many ground based telescopes have tried to detect optical afterglows of GRBs and to measure their redshifts by using the spectral absorption and/or emission lines of the interstellar matter in the host galaxy. However, the number of GRBs with measured redshift is only a fraction of all GRBs detected with BATSE, BeppoSAX, HETE-II and INTEGRAL satellites; we have still only about 30 GRBs with the known redshifts. The most of them occur at the cosmological distance, and the current record holder is GRB 000131 at $`z=4.5`$ (Andersen et al., 2000). According to the brightness distribution of GRBs with the known redshifts, the above satellites should have already detected much more distant GRBs, such as at $`z20`$ (Band, 2003). If we can establish a method for estimating the intrinsic brightness from the characteristics of the prompt gamma-ray emission, we can use the brightness of the GRB as a lighthouse to determine the unknown redshifts of majority of GRBs, which enables us to explore the early universe out to $`z20`$.

Using the geometrical corrections of the collimated jets, Frail et al. (2001) and Bloom et al. (2003) revealed that the bolometric energies released in the prompt emission are tightly clustering around the standard energy of $`1\times 10^{51}\mathrm{erg}`$. Thus, the explosion energy of GRBs can be used as a standard candle as the supernovae. However, the apparent brightness of GRBs strongly depends on the jet opening angle and the viewing angle. To use the GRB as the standard candle, we need to correct such effects.

Several authors tried to establish a method for estimating the isotropic luminosity from the observed GRB properties. Using the lightcurves of the prompt gamma-ray emissions, some pioneering works have been done; the variability–luminosity relation reported by Fenimore & Ramirez-Ruiz (2000), which indicates that the variable GRBs are much brighter than the smooth ones. The spectral time-lag, which is the interval of the peak arrival times between two different energy bands, also correlates with the isotropic luminosity (Norris et al., 2000). These properties based on the time-series data might be due to the effect of the viewing angle from the GRB jet (e.g., Ioka & Nakamura, 2001; Norris, 2002; Murakami et al., 2003).

On the other hand, based on the spectral analyses with the K-correction (Bloom et al., 2001), Amati et al. (2002) found the correlation between the isotropic-equivalent energy radiated in GRBs and the peak energies $`E_p`$, which is the energy at the peak of $`\nu F_\nu `$ spectrum. Atteia (2003) suggested a possibility of the empirical redshift indicator, which is based on the $`E_p`$ and the arrival number of photons.

Applying these luminosity indicators to the GRBs with the unknown redshifts, their redshifts can be estimated from the apparent gamma-ray brightness. As a natural application of the obtained redshift distribution, the GRB formation rates are discussed by several authors (e.g., Fenimore & Ramirez-Ruiz, 2000; Norris et al., 2000; Schaefer et al., 2001; Lloyd-Ronning et al., 2002; Murakami et al., 2003). Especially, using mathematically rigid method (Efron & Petrosian, 1992; Petrosian, 1993; Maloney & Petrosian, 1999), Lloyd-Ronning et al. (2002) have estimated the GRB formation rates from the variability–luminosity relation. These works give basically the same results; the GRB formation rate rapidly increases with the redshift at $`0z2`$, and it keeps on rising up to higher redshift ($`z12`$). The GRB formation rate did not decrease with $`z`$ in contrast with the star formation rates (SFRs) measured in UV, optical and infrared band (e.g., Madau et al., 1996; Barger et al., 2000; Stanway et al., 2003).

In this paper, we establish a new calibration formula of the redshift based on the $`E_p`$–luminosity relation of the brightest peak of the prompt gamma-ray emission from 9 GRBs with the known redshifts. Importantly, the uncertainty of our formula is much less than those of the previous works (lag, variability) as shown in section 3. Applying the obtained calibration, in section 4, we estimate the redshifts of 684 GRBs with the unknown redshifts. We then demonstrate the GRB formation rate out to $`z12`$ and the luminosity evolution using the non-parametric method (Efron & Petrosian, 1992; Petrosian, 1993; Maloney & Petrosian, 1999; Lloyd-Ronning et al., 2002). We emphasize that the GRB formation rate derived by us is based on the spectral analysis for the first time, and its uncertainty is well controlled in a small level. Throughout the paper, we assume the flat-isotropic universe with $`\mathrm{\Omega }_m=0.32`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.68`$ and $`H_0=72\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (Bennett et al., 2003; Spergel et al., 2003).

## 2 Data Analysis

We used 9 GRBs (970508, 970828, 971214, 980703, 990123, 990506, 990510, 991216, and 000131) in the BATSE archive with the known redshifts, and focused on the brightest peak in each GRB. We performed the spectral analysis with the standard data reduction for each GRB. We subtracted the background spectrum, which was derived from the average spectrum before and after the GRB in the same data set.

We adopted the spectral model of smoothly broken power-law (Band et al., 1993). The model function is described below.

$`N(E)=\{\begin{array}{cc}A\left(\frac{E}{100\mathrm{keV}}\right)^\alpha \mathrm{exp}(\frac{E}{E_0}),\hfill & \\ A\left(\frac{E}{100\mathrm{keV}}\right)^\beta \left(\frac{(\alpha \beta )E_0}{100\mathrm{keV}}\right)^{\alpha \beta }\mathrm{exp}(\beta \alpha ),\hfill & \end{array}`$ (3)
for $`E(\alpha \beta )E_0`$ and $`E(\alpha \beta )E_0`$, respectively. Here, $`N(E)`$ is in units of $`\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1\mathrm{keV}^1`$, and $`E_0`$ is the energy at the spectral break. $`\alpha `$ and $`\beta `$ is the low- and high-energy power-law index, respectively. For the case of $`\beta <2`$ and $`\alpha >2`$, the peak energy can be derived as $`E_p=(2+\alpha )E_0`$, which corresponds to the energy at the maximum flux in $`\nu F_\nu `$ spectra. The isotropic luminosity can be calculated with the observed flux as $`L=4\pi d_\mathrm{L}^2F_\gamma `$, where $`d_\mathrm{L}`$ and $`F_\gamma `$ are the luminosity distance and the observed energy flux, respectively.

## 3 $`E_p`$–Luminosity relation

In figure 1, we show the observed isotropic luminosity in units of $`10^{52}\mathrm{erg}\mathrm{s}^1`$ as a function of the peak energy, $`E_p(1+z)`$, in the rest frame of each GRB. For one of our sample (GRB 980703), a lower limit of $`E_p(1+z)`$ is set because the spectral index $`\beta >2`$. The results of BeppoSAX reported by Amati et al. (2002) are also plotted in the same figure. Here, we converted their peak fluxes into the same energy range of 30 – 10000 keV in our analysis with their spectral parameters. There is a good positive correlation between the $`E_p(1+z)`$ and the $`L`$. The linear correlation coefficient including the weighting factors is 0.957 for 15 degree of freedom (17 samples with firm redshift estimates) for the $`\mathrm{log}[E_p(1+z)]`$ and the $`\mathrm{log}[L]`$. The chance probability is $`1.85\times 10^9`$. When we adopt the power-law model as the $`E_p`$–luminosity relation, the best-fit function is

$`{\displaystyle \frac{L}{10^{52}\mathrm{ergs}}}=(4.29\pm 0.15)\times 10^5\left[{\displaystyle \frac{E_p(1+z)}{1\mathrm{keV}}}\right]^{1.94\pm 0.19}`$ (4)
where the uncertainties are 1 $`\sigma `$ error. This relation agrees well with the standard synchrotron model, $`LE_p^2`$ (e.g., Zhang & Mészáros, 2002; Lloyd et al., 2000).

## 4 Redshift Estimation and GRB Formation Rate

The $`E_p`$–luminosity relation derived in the previous section seems to be much better estimator of the isotropic luminosity compared with the spectral time-lag and the variability of GRBs (Norris et al., 2000; Fenimore & Ramirez-Ruiz, 2000; Schaefer et al., 2001) since the chance probability is extremely low. In this section, using the $`E_p`$–luminosity relation, we try to estimate the isotropic luminosities and the redshifts of the BATSE GRBs with the unknown redshifts.

We first selected about 1000 GRBs in a class of the long duration of $`T_{90}>2\mathrm{sec}`$ detected by BATSE, and extracted the brightest pulse in each GRB. We performed the spectral analysis for these peaks with the same method described in section 2. Once we obtained observed energy-flux $`F_\gamma `$ and $`E_p`$ at the observer’s rest frame, we can estimate the redshift with equation (4) and the luminosity distance as a function of the redshift. We could not determine $`E_p`$ for $`5`$ % of samples because of $`\beta >2`$ within 1 $`\sigma `$ uncertainty, and excluded them. After setting the flux limit of $`F_{\mathrm{limit}}=1\times 10^7\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, which is based on the dimmest one (GRB 970508), on the data set of $`(z,L)`$ plane, 684 samples are remained. In figure 2, we show the sample distribution in $`(z,L)`$ plane with the truncation by the flux limit.

The normalized cumulative luminosity function at each redshift bin is also shown in figure 2 . The shape of these luminosity functions are similar to each other, but the break-luminosity seems to increase toward the higher redshift. This indicates that a luminosity evolution is hidden in the $`(z,L)`$ plane in figure 2. Therefore, we have to remove the effect of the luminosity evolution from the data set before discussing the real luminosity function and the GRB formation rate. This is because the univariate distribution of the redshift and the luminosity can be estimated only when they are independent of each other.

We used the non-parametric method (Lynden-Bell, 1971; Efron & Petrosian, 1992; Petrosian, 1993; Maloney & Petrosian, 1999), which were used for the Quasar samples and first applied to the GRB samples by Lloyd-Ronning et al. (2002). The total luminosity function $`\mathrm{\Phi }(L,z)`$ can be written $`\mathrm{\Phi }(L,z)=\rho (z)\varphi (L/g_k(z),\alpha _s)/g_k(z)`$ without the loss of generality. Here, each function means the luminosity evolution $`g_k(z)`$, the density evolution $`\rho (z)`$ and the local luminosity function $`\varphi (L/g_k(z),\alpha _s)`$, respectively. Although the parameter $`\alpha _s`$ represents the shape of the luminosity function, we will ignored the effect of this parameter because the shape of the luminosity function is approximately same as shown in figure 2. We do not mention more about the method, but the details are found in Maloney & Petrosian (1999).

Assuming the functional form of the luminosity evolution as $`g_k(z)=(1+z)^k`$, which is also used in Maloney & Petrosian (1999) and Lloyd-Ronning et al. (2002), we convert the data set of $`(z,L)`$ to $`(z,L^{})`$, where $`L^{}=L/g_k(z)`$. We search the best value of $`k`$ giving the independent set of $`(z,L^{})`$ within the significance of $`1\sigma `$ error, and it is found to be $`g_k(z)=(1+z)^{1.85\pm 0.08}`$. Once we obtain the independent data set of the $`(z,L^{})`$, we can generate the cumulative luminosity function $`\psi (L^{})`$ with a simple formula (see equation (14) of Maloney & Petrosian, 1999). In figure 3, we show the cumulative luminosity function calculated by the formula. This function is approximately described by $`\psi (L^{})L^{0.3}`$ for $`0.1<L_{52}^{}<1`$ and $`\psi (L^{})L^{1.2}`$ for $`1<L_{52}^{}<50`$, where $`L_{52}^{}=L^{}/10^{52}\mathrm{ergs}`$.

We can also obtain the cumulative number distribution $`\psi (z)`$ as a function of $`z`$. The differential form of the function is useful for the purpose of comparison with the star formation rates in other wave bands. We convert $`\psi (z)`$ into the differential form with the following equation.

$`\rho (z)={\displaystyle \frac{d\psi (z)}{dz}}(1+z)\left({\displaystyle \frac{dV(z)}{dz}}\right)^1,`$ (5)
where the additional factor of $`(1+z)`$ comes from the cosmological time dilation, and $`dV(z)/dz`$ is a differential comoving volume. In figure 4, we show the relative comoving GRB rate $`\rho (z)`$ in unit proper volume.

## 5 Discussion

We have investigated the spectral property of the brightest peak of each GRB with the known redshifts, and have found a fine correlation between the peak energy $`E_p(1+z)`$ and the isotropic luminosity. While this correlation against a smaller sample has been pointed out by some authors (e.g., Amati et al., 2002; Atteia, 2003; Schaefer, 2003a, b), we have succeeded in combining the results of BeppoSAX and BATSE to describe equation (4).

Using the new $`E_p`$–luminosity relation, we have estimated the redshifts of 684 GRBs with the unknown redshifts. We found the existence of the luminosity evolution $`g_k(z)=(1+z)^{1.85\pm 0.08}`$ for GRB samples as shown in figure 2. The null-hypothesis of the luminosity evolution is rejected about $`9\sigma `$ significance. Lloyd-Ronning et al. (2002) also suggest the presence of the luminosity evolution as $`g_k(z)=(1+z)^{1.4\pm 0.5}`$. These two values are consistent with each other. The luminosity evolutions are found in other objects. For example, Caditz & Petrosian (1990) and Maloney & Petrosian (1999) estimated the luminosity evolution of the QSO samples as $`g_k(z)=(1+z)^3`$ and $`(1+z)^{2.58}`$, respectively. Based on the observation of Subaru Deep Field and a photometric redshift estimation for $`K^{}`$-band selected galaxy samples, Kashikawa et al. (2003) found the strong luminosity evolution in the rest UV band of their galaxies.

The form of the cumulative luminosity function is independent of $`z`$ except for the break luminosity, which changes with $`z`$. We propose that figure 3 might include important information on the jets responsible for the prompt $`\gamma `$-ray emissions and a distribution of their opening angles. Consider an extremely simple model for a uniform jet with an opening half-angle $`\theta _j`$ and a constant geometrically-corrected luminosity $`L_0`$, which is viewed from an angle of $`\theta _v`$. Then, in a crude approximation the luminosity $`L`$ is given by

$`L=\{\begin{array}{cc}2L_0\theta _j^2\hfill & \mathrm{for}\theta _v<\theta _j\hfill \\ 2L_0(\theta _j^6/\theta _v^8)\hfill & \mathrm{for}\theta _v>\theta _j.\hfill \end{array}`$ (8)
For the case of $`\theta _v>\theta _j`$, $`L`$ is proportional to $`\delta ^4`$, where $`\delta =[\gamma (1\beta \mathrm{cos}\theta _v)]^1\theta _v^2`$, so that the luminosity has the dependence of $`\theta _v^8`$ (Ioka & Nakamura, 2001; Yamazaki et al., 2002, 2003). The dependence of $`\theta _j^6`$ is determined in order that two functions in equation (8) are continuously connected at $`\theta _v=\theta _j`$. We also consider the distribution of $`\theta _j`$ in the form $`f(\theta _j)d\theta _j\theta _j^qd\theta _j`$ when $`\theta _{\mathrm{min}}<\theta _j<\theta _{\mathrm{max}}`$. Then, easy calculations show that in the case of $`q<5/2`$, we have

$`N(>L)\{\begin{array}{cc}L^{1/4}\hfill & \mathrm{for}L<2L_0\theta _{\mathrm{max}}^2\hfill \\ L^{(q3)/2}\hfill & \mathrm{for}2L_0\theta _{\mathrm{max}}^2<L<2L_0\theta _{\mathrm{min}}^2.\hfill \end{array}.`$ (11)
This is a broken power low with the break luminosity $`2L_0\theta _{\mathrm{max}}^2`$. Then if $`\theta _{\mathrm{max}}^2L_0g_k(z)`$ with $`q=0.6`$, we can roughly reproduce figure 3. This suggests that either the maximum opening half angle of the jet decreases or $`L_0`$ increases as a function of the redshift.

This work is the first study to generate the GRB formation rate with $`E_p`$–luminosity relation. The result indicates that the GRB formation rate always increases toward $`z12`$. This is consistent with the previous works using the GRB variability (Fenimore & Ramirez-Ruiz, 2000; Lloyd-Ronning et al., 2002) and the spectral time-lag (Norris et al., 2000; Schaefer et al., 2001; Murakami et al., 2003). Quantitatively, GRB formation rate is proportional to $`(1+z)^{56}`$ for $`z1`$ and to $`(1+z)^1`$ for $`z1`$. On the other hand, the SFR measured in UV, optical, and infrared is proportional to $`(1+z)^{23}`$ for $`z1`$ (e.g., Lilly et al., 1996; Cowie et al., 1999; Glazebrook et al., 2003) and to $`(1+z)^{10}`$ for $`z1`$ (e.g., Madau et al., 1996; Barger et al., 2000; Stanway et al., 2003; Kashikawa et al., 2003). Therefore, we find that the ratio of the GRB formation rate to the SFR evolves along redshift, say $`(1+z)^{23}`$.

Recently, it has been strongly suggested that the long duration GRBs arise from the collapse of a massive star (Galama et al., 1998; Hjorth et al., 2003; Price et al., 2000; Uemura et al., 2003; Stanek et al., 2003). Hence our result implies that either the formation rate of the massive star or the fraction of GRB progenitor in massive stars at the high redshift should be significantly greater than the present value. However, if the SFR rapidly increases along redshift as suggested by Lanzetta et al. (2003), the fraction of GRB progenitors does not change so much.

The existence of the luminosity evolution of GRBs, i.e., $`g_k(z)=(1+z)^{1.85}`$, may suggest the evolution of GRB progenitor itself (e.g., mass) or the jet evolution. Although the jet opening angle evolution was suggested (Lloyd-Ronning et al., 2002), in our extremely simple model, either the maximum jet opening angle decreases or the jet total energy increases. In the former case the GRB formation rate shown in figure 4 may be an underestimate since the chance probability to observe the high redshift GRB will decrease. If so, the evolution of the ratio of the GRB formation rate to the SFR becomes more rapid. On the other hand, in the latter case, GRB formation rate of figure 4 gives a reasonable estimate.

We thank HEASOFT help desk, especially Bryan K. Irby, for the support to the BATSE analysis tools. This research is supported by Grants-in-Aid from the Japanese Ministry of Education, Culture, Sports, Science and Technology (15740149, 14204024, 14047212 and 12640302). R.Y., A.K.I. and K.I. are supported by JSPS Research Fellowship for Young Scientists.


format: text

text:
# Lemma 1

The Interleaved Multichromatic Number of a Graph

Valmir C. Barbosa

Universidade Federal do Rio de Janeiro

Programa de Engenharia de Sistemas e Computação, COPPE

Caixa Postal 68511

21941-972 Rio de Janeiro - RJ, Brazil

valmir@cos.ufrj.br

Abstract

For $`k1`$, we consider interleaved $`k`$-tuple colorings of the nodes of a graph, that is, assignments of $`k`$ distinct natural numbers to each node in such a way that nodes that are connected by an edge receive numbers that are strictly alternating between them with respect to the relation $`<`$. If it takes at least $`\chi _{\mathrm{int}}^k(G)`$ distinct numbers to provide graph $`G`$ with such a coloring, then the interleaved multichromatic number of $`G`$ is $`\chi _{\mathrm{int}}^{}(G)=inf_{k1}\chi _{\mathrm{int}}^k(G)/k`$ and is known to be given by a function of the simple cycles of $`G`$ under acyclic orientations if $`G`$ is connected.1 V. C. Barbosa and E. Gafni, ACM Trans. on Programming Languages and Systems 11 (1989), 562–584. This paper contains a new proof of this result. Unlike the original proof, the new proof makes no assumptions on the connectedness of $`G`$, nor does it resort to the possible applications of interleaved $`k`$-tuple colorings and their properties.

Keywords: Chromatic number, multichromatic number, interleaved multichromatic number, acyclic orientations.

1. Introduction

Let $`G=(N,E)`$ be an undirected graph with $`|E|1`$. For $`k1`$, a $`k`$-tuple coloring of $`G`$ is an assignment of $`k`$ distinct natural numbers (colors) to each of the nodes of $`G`$ in such a way that no two nodes connected by an edge are assigned an identical color. The least number of colors with which $`G`$ can be $`k`$-tuple colored is its $`k`$-chromatic number, denoted by $`\chi ^k(G)`$. For $`k=1`$, the $`1`$-chromatic number of $`G`$ is the graph’s chromatic number, denoted by $`\chi (G)`$.

There is also a notion of “efficiency” that goes along $`k`$-tuple colorings of $`G`$, which is the notion that, for two distinct positive integers $`k_1`$ and $`k_2`$, it is more “efficient” to $`k_1`$-tuple color $`G`$ than to $`k_2`$-tuple color it if and only if $`\chi ^{k_1}(G)/k_1<\chi ^{k_2}(G)/k_2`$. This notion is formalized by the definition of the multichromatic number of $`G`$ , denoted by $`\chi ^{}(G)`$ and given by

$$\chi ^{}(G)=\underset{k1}{inf}\frac{\chi ^k(G)}{k}.$$
$`(1)`$
Because the infimum in (1) can be shown to be always attained , a $`k`$ for which the ratio in (1) is minimum gives a most “efficient” $`k`$-tuple coloring of $`G`$.

One special case of $`k`$-tuple colorings that has great practical interest is the case of interleaved $`k`$-tuple colorings, defined as follows. For $`iN`$, let $`c_i^1,\mathrm{},c_i^k`$ be the $`k`$ colors assigned to node $`i`$. This $`k`$-tuple coloring is interleaved if and only if either $`c_i^1<c_j^1<\mathrm{}<c_i^k<c_j^k`$ or $`c_j^1<c_i^1<\mathrm{}<c_j^k<c_i^k`$ for all $`(i,j)E`$. Every $`1`$-tuple coloring is interleaved, and the counterparts for interleaved colorings of the graph’s $`k`$-chromatic and multichromatic numbers are, respectively, its interleaved $`k`$-chromatic number and interleaved multichromatic number. These are denoted, respectively, by $`\chi _{\mathrm{int}}^k(G)`$ and $`\chi _{\mathrm{int}}^{}(G)`$. Following (1), the latter is defined as

$$\chi _{\mathrm{int}}^{}(G)=\underset{k1}{inf}\frac{\chi _{\mathrm{int}}^k(G)}{k}.$$
$`(2)`$
Clearly, $`\chi ^{}(G)\chi _{\mathrm{int}}^{}(G)\chi (G)`$.

For connected $`G`$, the infimum in (2) has also been shown to be always attained and given as follows . Let $`\mathrm{K}`$ be the set of all simple cycles of $`G`$. For $`\kappa \mathrm{K}`$, let $`\kappa ^+`$ and $`\kappa ^{}`$ be the two possible traversal directions of $`\kappa `$. If $`\omega `$ is an acyclic orientation of $`G`$ (that is, an assignment of directions to edges that forms no directed cycles), then we denote by $`m(\kappa ^+,\omega )`$ the number of edges in $`\kappa `$ oriented by $`\omega `$ in the $`\kappa ^+`$ direction. Likewise for $`m(\kappa ^{},\omega )`$. What has been shown for nonempty $`\mathrm{K}`$ is that

$$\frac{1}{\chi _{\mathrm{int}}^{}(G)}=\underset{\omega \mathrm{\Omega }(G)}{\mathrm{max}}\underset{\kappa \mathrm{K}}{\mathrm{min}}\frac{\mathrm{min}\{m(\kappa ^+,\omega ),m(\kappa ^{},\omega )\}}{|\kappa |},$$
$`(3)`$
where $`\mathrm{\Omega }(G)`$ is the set of all acyclic orientations of $`G`$ and $`|\kappa |`$ is the size of $`\kappa `$, that is, $`|\kappa |=m(\kappa ^+,\omega )+m(\kappa ^{},\omega )`$. If $`\mathrm{K}=\mathrm{}`$ (i.e., $`G`$ is a tree), then $`\chi _{\mathrm{int}}^{}(G)=2`$.

In , $`1/\chi _{\mathrm{int}}^{}(G)`$ is analyzed as a measure of concurrency (or parallelism) for certain distributed computations over systems that can be represented by connected graphs, and (3) is obtained through proofs that make heavy use of both the graph’s connectedness and properties of those distributed computations. What we do in this paper is to provide an alternative proof for (3), one that does not make assumptions on the connectedness of $`G`$ and is built from “first principles,” that is, independently of the possible applications of the concept of “efficient” interleaved $`k`$-tuple colorings. This is done in Section 2, and in Section 3 we give concluding remarks and argue that (3) provides both “efficient” interleaved $`k`$-tuple colorings and optimal $`1`$-tuple colorings with a common underlying relationship to the set $`\mathrm{\Omega }(G)`$ of the graph’s acyclic orientations.

2. A new proof

Our proof is based on the lexicographic product  of $`G`$ and the complete graph on $`k`$ nodes. This lexicographic product is another undirected graph, denoted by $`G^k`$ and having $`k`$ nodes for each node in $`G`$. For $`iN`$, we denote the nodes of $`G^k`$ by $`i^1,\mathrm{},i^k`$. The edges of $`G^k`$ are deployed in such a way that, for fixed $`i`$, an edge exists connecting every pair of nonidentical nodes of $`i^1,\mathrm{},i^k`$, and for $`(i,j)E`$, an edge exists between each one of $`i^1,\mathrm{},i^k`$ and each one of $`j^1,\mathrm{},j^k`$. Clearly, every $`k`$-tuple coloring of $`G`$ corresponds to a $`1`$-tuple coloring of $`G^k`$, and conversely.

Now consider the set $`\mathrm{\Omega }(G^k)`$ of the acyclic orientations of $`G^k`$. We say that $`\phi \mathrm{\Omega }(G^k)`$ is a layered orientation of $`G^k`$ if and only if every edge $`(i^{k_1},j^{k_2})`$ for which $`k_1>k_2`$ is oriented from $`i^{k_1}`$ to $`j^{k_2}`$, and furthermore every edge $`(i^{k_1},j^{k_1})`$ gets the same orientation for $`k_1=1,\mathrm{},k`$. If we visualize $`G^k`$ as $`k`$ instances of $`G`$ stacked upon one another in such a way that the instance at the bottom contains superscript-$`1`$ nodes, the one right above it superscript-$`2`$ nodes, and so on, then what a layered acyclic orientation of $`G^k`$ does is to provide the $`k`$ instances of $`G`$ with identical acyclic orientations from $`\mathrm{\Omega }(G)`$ and orient all other edges (the ones that connect nodes from different instances of $`G`$) from top to bottom. The set of all layered acyclic orientations of $`G^k`$ is denoted by $`\mathrm{\Omega }^L(G^k)`$.

In $`G^k`$, we refer to each of the aforementioned instances of $`G`$ as a layer, and number each layer with the superscript of the nodes that it contains. Given a layered acyclic orientation, every maximal directed path starts in layer $`k`$ (because that is where all sources are, these being nodes whose adjacent edges are all oriented outward) and ends in layer $`1`$ (which contains all sinks, that is, nodes whose adjacent edges are all oriented inward).

For a generic undirected graph $`H`$ and for an acyclic orientation $`\omega `$ of $`H`$, let $`l_\omega `$ be the number of nodes in the longest directed path in $`H`$ according to $`\omega `$. Also, let a $`1`$-tuple coloring of $`H`$ be called monotonic with respect to an acyclic orientation $`\omega `$ of $`H`$ if and only if the colors it assigns decrease along any directed path according to $`\omega `$. If $`\mathrm{\Omega }(H)`$ is the set of acyclic orientations of $`H`$, then we have the following, for whose proof the reader is referred to .

###### Lemma 1

To every $`\omega \mathrm{\Omega }(H)`$ there corresponds a monotonic $`1`$-tuple coloring of $`H`$ with at most $`l_\omega `$ colors. Conversely, to every $`1`$-tuple coloring of $`H`$ by a total of $`c`$ colors there corresponds an orientation $`\omega \mathrm{\Omega }(H)`$ with respect to which the coloring is monotonic and such that $`l_\omega c`$.

By Lemma 1, the chromatic number of graph $`H`$ is given by the number of nodes on the longest directed path that is shortest among all possible acyclic orientations of $`H`$ . In the case of $`G^k`$, and considering that $`\chi ^k(G)=\chi (G^k)`$, we have $`\chi ^k(G)=\mathrm{min}_{\phi \mathrm{\Omega }(G^k)}l_\phi `$. Lemma 1 is the basis of our first supporting result, which states that it suffices to look at the orientations in $`\mathrm{\Omega }^L(G^k)\mathrm{\Omega }(G^k)`$ to obtain $`\chi _{\mathrm{int}}^k(G)`$.

###### Lemma 2

For $`k1`$, $`\chi _{\mathrm{int}}^k(G)=\mathrm{min}_{\phi \mathrm{\Omega }^L(G^k)}l_\phi `$.

Proof: Let $`\phi `$ be a layered acyclic orientation of $`G^k`$. By Lemma 1, $`G^k`$ can be colored monotonically by a $`1`$-tuple coloring that employs a total of at most $`l_\phi `$ colors. If the instance of $`(i,j)E`$ within each layer is oriented by $`\phi `$ from the corresponding instance of $`i`$ to the corresponding instance of $`j`$, then by definition of a layered orientation a directed path exists from $`i^k`$ to $`j^1`$ that alternates instances of $`i`$ and $`j`$ along its way (that is, a path like $`i^kj^ki^{k1}j^{k1}\mathrm{}i^1j^1`$). It then follows from the monotonicity of $`1`$-tuple colorings of $`G^k`$ that such a coloring corresponds to an interleaved $`k`$-tuple coloring of $`G`$, thus

$$\chi _{\mathrm{int}}^k(G)\underset{\phi \mathrm{\Omega }^L(G^k)}{\mathrm{min}}l_\phi .$$

Conversely, consider an interleaved $`k`$-tuple coloring of $`G`$ by a total of $`c`$ colors. This coloring corresponds to possibly several $`1`$-tuple colorings of $`G^k`$ by the same $`c`$ colors. Of these, and for $`iN`$, consider the $`1`$-tuple coloring that assigns the highest of the $`k`$ colors of $`i`$ to $`i^k`$, then the next highest to $`i^{k1}`$, and so on. By Lemma 1, an orientation $`\phi `$ of $`G^k`$ exists with respect to which this $`1`$-tuple coloring of $`G^k`$ is monotonic and as such orients all edges of $`G^k`$ from the node with higher color to the one with lower. For an arbitrary edge $`(i,j)`$ of $`G`$, this, together with the fact that the assumed $`k`$-tuple coloring of $`G`$ is interleaved, implies that $`\phi `$ is layered. In addition, Lemma 1 also ensures that $`l_\phi c`$, hence

$$\chi _{\mathrm{int}}^k(G)\underset{\phi \mathrm{\Omega }^L(G^k)}{\mathrm{min}}l_\phi .$$

Our next supporting result is a statement on the morphology of longest directed paths according to layered acyclic orientations of $`G^k`$. What it states is that such a directed path’s every edge is either confined to a layer or joins two adjacent layers. In the latter case, it also states that, if such an edge exists from $`i^{k_1}`$ to $`j^{k_11}`$, then a directed edge also exists from $`j^{k_1}`$ to $`i^{k_1}`$.

###### Lemma 3

For $`kk_1k_21`$ and $`\phi \mathrm{\Omega }^L(G^k)`$, let $`p`$ be a longest directed path according to $`\phi `$, and let $`i^{k_1}`$ and $`j^{k_2}`$ be such that a directed edge exists on $`p`$ from $`i^{k_1}`$ to $`j^{k_2}`$. Then, for $`\mathrm{}=1,\mathrm{},k`$, either $`k_1=k_2`$ and a directed edge exists from $`i^{\mathrm{}}`$ to $`j^{\mathrm{}}`$, or $`k_1=k_2+1`$ with $`ij`$ and a directed edge exists from $`j^{\mathrm{}}`$ to $`i^{\mathrm{}}`$.

Proof: Note, first, that $`k_1k_2+1`$, since $`k_1>k_2+1`$ contradicts the hypothesis that $`p`$ is longest. If $`k_1=k_2`$, then the lemma follows directly from the definition of a layered acyclic orientation of $`G^k`$. We then assume that $`k_1=k_2+1`$.

What remains to be shown is that two specific scenarios cannot happen if $`\phi `$ is layered. The first one is that $`i=j`$, that is, $`p`$ contains an edge directed from $`i^{k_1}`$ to $`i^{k_11}`$. That this cannot happen follows from the fact that replacing this edge on $`p`$ with either the directed path $`i^{k_1}r^{k_1}i^{k_11}`$ or the directed path $`i^{k_1}r^{k_11}i^{k_11}`$ for some $`(i,r)E`$ yields a directed path longer than $`p`$. Because $`|E|1`$, at least one of the two possibilities is certain to occur, so we have a contradiction.

The second scenario corresponds to $`(i,j)E`$ when every layer $`\mathrm{}`$ contains an edge directed from $`i^{\mathrm{}}`$ to $`j^{\mathrm{}}`$. In this case, an edge exists directed from $`i^{k_1}`$ to $`j^{k_11}`$ on $`p`$. Once again, this edge can be replaced with the directed path $`i^{k_1}j^{k_1}i^{k_11}j^{k_11}`$ to yield a directed path longer than $`p`$, also a contradiction.

By Lemma 3, a longest directed path $`p`$ in $`G^k`$ according to a layered acyclic orientation $`\phi `$ looks like the following. It starts at a source in layer $`k`$, winds its way through the $`k`$ layers, and ends at a sink in layer $`1`$. In addition, for each of the $`k1`$ edges that lead from one layer to the next, say for example an edge from $`i^{k_1}`$ to $`j^{k_11}`$ for $`1<k_1k`$, an edge exists in layer $`\mathrm{}`$ directed from $`j^{\mathrm{}}`$ to $`i^{\mathrm{}}`$ for $`\mathrm{}=1,\mathrm{},k`$.

If we let $`\omega \mathrm{\Omega }(G)`$ be the acyclic orientation of $`G`$ that corresponds to the orientation of a layer by $`\phi `$ (the same for all layers), then $`p`$ can be seen to leave a trace on $`G`$ oriented by $`\omega `$ as well. This trace is an undirected path in $`G`$ (not necessarily a simple one) that can be traversed by following $`\omega `$ on the edges that correspond to $`p`$ being confined to a layer, and by going against $`\omega `$ on the $`k1`$ edges that correspond to the transition of $`p`$ from one layer to the next. Of course, leaving such a trace on $`G`$ oriented by $`\omega `$ is no prerogative of longest directed paths, but rather applies to any other directed path whose edges are either confined to a same layer or join adjacent layers as described in Lemma 3.

Next we give our main result, in which (3) is established. Recall that $`\mathrm{K}`$ stands for the set of all simple cycles of $`G`$.

###### Theorem 4

If $`G`$ is a forest, then $`\chi _{\mathrm{int}}^{}(G)=2`$. Otherwise,

$$\chi _{\mathrm{int}}^{}(G)=\underset{\omega \mathrm{\Omega }(G)}{\mathrm{min}}\underset{\kappa \mathrm{K}}{\mathrm{max}}\frac{|\kappa |}{\mathrm{min}\{m(\kappa ^+,\omega ),m(\kappa ^{},\omega )\}}.$$

Proof: By definition of $`\chi _{\mathrm{int}}^{}(G)`$ (cf. (2)) and by Lemma 2, we have

$$\begin{array}{cc}\hfill \chi _{\mathrm{int}}^{}(G)& =\underset{k1}{inf}\frac{\chi _{\mathrm{int}}^k(G)}{k}\hfill \\ & =\underset{k1}{inf}\frac{1}{k}\underset{\phi \mathrm{\Omega }^L(G^k)}{\mathrm{min}}l_\phi .\hfill \end{array}$$

Also, the definition of layered acyclic orientations of $`G^k`$ implies that there exists a one-to-one correspondence between such orientations and the acyclic orientations of $`G`$. If, for $`k1`$ and $`\omega \mathrm{\Omega }(G)`$, we let $`\phi _\omega ^k`$ be the corresponding orientation in $`\mathrm{\Omega }^L(G^k)`$, then we can further write

$$\begin{array}{cc}\hfill \chi _{\mathrm{int}}^{}(G)& =\underset{k1}{inf}\frac{1}{k}\underset{\omega \mathrm{\Omega }(G)}{\mathrm{min}}l_{\phi _\omega ^k}\hfill \\ & =\underset{\omega \mathrm{\Omega }(G)}{\mathrm{min}}\underset{k1}{inf}\frac{l_{\phi _\omega ^k}}{k}.\hfill \end{array}$$

If $`G`$ is a forest, then the theorem follows trivially (this is also a consequence of the facts that $`\chi ^{}(G)=\chi (G)=2`$ if $`G`$ is a forest and that $`\chi ^{}(G)\chi _{\mathrm{int}}^{}(G)\chi (G)`$). If not, then it suffices to argue that

$$\underset{\kappa \mathrm{K}}{\mathrm{max}}\frac{|\kappa |}{\mathrm{min}\{m(\kappa ^+,\omega ),m(\kappa ^{},\omega )\}}=\underset{k1}{inf}\frac{l_{\phi _\omega ^k}}{k}$$

for all $`\omega \mathrm{\Omega }(G)`$.

Let $`\widehat{\mathrm{K}}\mathrm{K}`$ be the set of all simple cycles $`\kappa `$ of $`G`$ for which $`|\kappa |/\mathrm{min}\{m(\kappa ^+,\omega ),m(\kappa ^{},\omega )\}`$ is maximum over $`\mathrm{K}`$, and for all $`\kappa \widehat{\mathrm{K}}`$ assume that $`m(\kappa ^+,\omega )m(\kappa ^{},\omega )`$. For $`k1`$, every $`\kappa \widehat{\mathrm{K}}`$ can be seen to give rise to directed paths in $`G^k`$ whose edges either lead to other nodes in the same layer or else descend to the next layer as described in Lemma 3. As we discussed following the presentation of that lemma, each such path leaves a trace on $`G`$ oriented by $`\omega `$. One possibility is the path that we denote by $`p_k(\kappa ,\omega )`$ and whose trace has the following characteristics. It starts at the origin of the longest segment of $`\kappa `$ oriented by $`\omega `$ in the $`\kappa ^+`$ direction, then winds its way around $`\kappa `$ a number of times given by $`t=(k1)/m(\kappa ^{},\omega )`$, and finally goes around $`\kappa `$ one last time until $`(k1)modm(\kappa ^{},\omega )`$ additional edges have been traversed against $`\omega `$ and one last segment has been traversed in agreement with $`\omega `$. Along this last traversal around $`\kappa `$, let $`ϵ^+`$ denote the number of edges traversed in agreement with $`\omega `$ and $`ϵ^{}`$ of those traversed against $`\omega `$. Evidently, $`ϵ^{}m(\kappa ^{},\omega )1`$.

Note that this last traversal around $`\kappa `$ starts, like $`p_k(\kappa ,\omega )`$, at the origin of the longest segment of $`\kappa `$ oriented by $`\omega `$ in the $`\kappa ^+`$ direction. If, in addition, it is chosen to be longest among all possibilities, then we have

$$ϵ^+(1+ϵ^{})\frac{m(\kappa ^+,\omega )}{m(\kappa ^{},\omega )}.$$

Equivalently,

$$ϵ^+m(\kappa ^{},\omega )(1+ϵ^{})m(\kappa ^+,\omega ),$$

or yet,

$$\begin{array}{c}\text{ }\left(1+tm(\kappa ^+,\omega )+tm(\kappa ^{},\omega )+ϵ^++ϵ^{}\right)m(\kappa ^{},\omega )\text{ }\hfill \\ \text{ }\left(1+tm(\kappa ^{},\omega )+ϵ^{}\right)\left(m(\kappa ^+,\omega )+m(\kappa ^{},\omega )\right),\text{ }\hfill \end{array}$$

which yields

$$\begin{array}{cc}\hfill \frac{1+tm(\kappa ^+,\omega )+tm(\kappa ^{},\omega )+ϵ^++ϵ^{}}{1+tm(\kappa ^{},\omega )+ϵ^{}}& \frac{m(\kappa ^+,\omega )+m(\kappa ^{},\omega )}{m(\kappa ^{},\omega )}\hfill \\ & =\frac{|\kappa |}{m(\kappa ^{},\omega )}.\hfill \end{array}$$

If we let $`l_k(\kappa ,\omega )`$ denote the number of nodes of $`p_k(\kappa ,\omega )`$, then the left-hand side of the latter inequality is equal to $`l_k(\kappa ,\omega )/k`$. Considering further that $`l_{\phi _\omega ^k}l_k(\kappa ,\omega )`$, we obtain

$$\frac{l_{\phi _\omega ^k}}{k}\frac{|\kappa |}{m(\kappa ^{},\omega )}.$$

Also, for sufficiently large $`k`$, a longest directed path in $`G^k`$ (comprising $`l_{\phi _\omega ^k}`$ nodes, by definition) leaves a trace on $`G`$ oriented by $`\omega `$ which, after a preamble that comprises a number of nodes that is bounded as $`k`$ grows, coincides with one of the members, say $`\widehat{\kappa }`$, of $`\widehat{\mathrm{K}}`$. This is ensured by the fact that, for every $`m(\widehat{\kappa }^{},\omega )`$ times that the path has to change layers, the greatest contribution to its number of nodes is $`|\widehat{\kappa }|`$, by definition of $`\widehat{\mathrm{K}}`$.

So we can write, by Lemma 3,

$$\frac{l_{\phi _\omega ^k}}{k}=\frac{tm(\widehat{\kappa }^+,\omega )+tm(\widehat{\kappa }^{},\omega )+o(k)}{tm(\widehat{\kappa }^{},\omega )+o(k)},$$

following the standard use of $`o(k)`$ to signify that $`lim_k\mathrm{}o(k)/k=0`$, for some $`t`$ that grows without bounds along with $`k`$. Clearly, then,

$$\underset{k\mathrm{}}{lim}\frac{l_{\phi _\omega ^k}}{k}=\frac{|\widehat{\kappa }|}{m(\widehat{\kappa }^{},\omega )},$$

which concludes the proof.

3. Discussion

This paper contains a new proof of (3). This new proof, in contrast to the original proof that appears in , does not make assumptions on the connectedness of $`G`$, nor does it resort to the specifics of possible applications of the concept of interleaved $`k`$-tuple colorings. Unlike the proof of , however, the proof in this paper does not establish that the infimum in (2) is always achieved, which is indeed the case. For this, the original proof continues to be the source, as establishing that property is achieved independently of the connectedness of $`G`$.

As expressed in Theorem 4, $`\chi _{\mathrm{int}}^{}(G)`$ corresponds to the minimum, taken over all of the acyclic orientations of $`G`$, of a function of certain undirected paths of $`G`$ under those acyclic orientations. Specifically, if for $`\omega \mathrm{\Omega }(G)`$ we let $`P(\omega )`$ comprise every path $`p`$ in $`G`$ such that $`p`$ contains all nodes in a simple cycle and starts at a node that is a source according to $`\omega `$ if only the edges on the simple cycle are considered, then $`\chi _{\mathrm{int}}^{}(G)`$ can be written as

$$\chi _{\mathrm{int}}^{}(G)=\underset{\omega \mathrm{\Omega }(G)}{\mathrm{min}}\underset{pP(\omega )}{\mathrm{max}}\frac{l(p)}{1+c(p,\omega )},$$
$`(4)`$
where $`l(p)`$ is the number of nodes of $`p`$ and $`c(p,\omega )`$ is the number of edges on $`p`$ that are oriented by $`\omega `$ contrary to the traversal of $`p`$ from its starting node.

Interestingly, the expression in (4) is the same as the aforementioned consequence of Lemma 1 that

$$\chi (G)=\underset{\omega \mathrm{\Omega }(G)}{\mathrm{min}}\underset{pP(\omega )}{\mathrm{max}}l(p),$$
$`(5)`$
if only we now let $`P(\omega )`$ stand for the set of directed paths of $`G`$ according to $`\omega `$, since in this case $`c(p,\omega )=0`$ for all $`\omega \mathrm{\Omega }(G)`$ and all $`pP(\omega )`$.

Since by (4) and (5), respectively, it is possible to assess a graph’s interleaved multichromatic number and chromatic number by looking at functions of the graph’s paths as oriented by the acyclic orientations of $`G`$, one interesting question is whether $`\chi ^{}(G)`$, the multichromatic number of $`G`$, is also amenable to a similar characterization. Unlike the other two cases, however, there does not appear to be a result like Lemma 2, which establishes the fundamental one-to-one correspondence between the acyclic orientations of $`G^k`$ that matter and those of $`G`$. Finding such a characterization for $`\chi ^{}(G)`$ remains then an open question.

Acknowledgments

The author acknowledges partial support from CNPq, CAPES, the PRONEX initiative of Brazil’s MCT under contract 41.96.0857.00, and a FAPERJ BBP grant.

References

1. V. C. Barbosa and E. Gafni, “Concurrency in heavily loaded neighborhood-constrained systems,” ACM Trans. on Programming Languages and Systems 11 (1989), 562–584.

2. F. H. Clarke and R. E. Jamison, “Multicolorings, measures and games on graphs,” Discrete Mathematics 14 (1976), 241–245.

3. R. W. Deming, “Acyclic orientations of a graph and chromatic and independence numbers,” J. of Combinatorial Theory B 26 (1979), 101–110.

4. D. Geller and S. Stahl, “The chromatic number and other functions of the lexicographic product,” J. of Combinatorial Theory B 19 (1975), 87–95.

5. S. Stahl, “$`n`$-tuple colorings and associated graphs,” J. of Combinatorial Theory B 20 (1976), 185–203.


text:
# Multi-frequency & Multi-epoch VLBI study of Cygnus A

## 1. Introduction

Cygnus A is one of the most powerful radio galaxies at a redshift of ($`z=0.057`$). It is the archetypical FR II radio galaxy. In the radio bands, Cygnus A is characterized by two strong lobes separated by $``$ 2 in the sky. Two highly collimated jets connect the lobes with the core (Perley et al. 1984; Carilli et al. 1991). On kiloparsec scales, the jet is oriented along P.A. $`285^{}`$ and the fainter counter-jet along P.A. $`107^{}`$. Due to the large inclination of the jet with respect to the observer, and the correspondingly reduced relativistic effects, Cygnus A is an ideal candidate for detailed studies of its jet physics, which is thought to be similar to those of luminous quasars (e.g. Barthel 1989).

We present and discuss results from multi-frequency VLBI observations at 15, 22 and 43 GHz and from the first multi-frequency phase referenced observations of Cygnus A done with the VLBA at 15 and 22 GHz. To complement our data at 15 GHz we used ten epochs from the VLBA 2 cm Survey (Kellermann et al. 1998; Zensus et al. 2002) to analyse the source kinematics with improved accuracy.

## 2. Observations and Data Reduction

We observed Cygnus A in 1996 with the VLBA+Effelsberg at 15, 22, and 43 GHz, in 2002 at 5 and 15 GHz and in 2003 with the VLBA only at 15 and 22 GHz in phase-referencing mode. All observations were done in dual circular polarization and were correlated in Socorro. The data were reduced in the standard manner using Aips. The imaging of the source implying phase and amplitude self-calibration was done using Difmap.

After imaging we fitted circular Gaussian components to the self-calibrated data in order to parameterise the source structure. Conservatively, we assume errors of $`10`$ % in the flux density arising from the uncertainties of the amplitude calibration and from the formal errors of the model fits. An estimate for the position error is given by $`\mathrm{\Delta }r=\frac{\sigma \mathrm{\Theta }}{2S_\mathrm{P}}`$ (Fomalont 1989), where $`\sigma `$ is the residual noise of the map after the subtraction of the model, $`\mathrm{\Theta }`$ the width of the component, and $`S_\mathrm{P}`$ the peak flux density. This formula tends to underestimate the error if the peak flux density is very high or the width of the component is small. In the case of a small FWHM we used the beam size instead.

## 3. Results and Discussion

To investigate the spectral properties and the kinematics of Cygnus A on parsec scales we cross-identified individual model components along the jet using their relative separation from each other, their flux density and size. Since the observations from 2003 were phase-referenced, we compared the alignment between the two frequencies from the modelfits with those from the phase-referencing and found similar results. The most likely position for the central engine is located between component C3 and J10 (see Fig. 1). This position shows a slightly inverted spectrum and turned out to be stable in our kinematical study. It seems to be the same position, which Krichbaum et al. (1998) already used in their study.

### 3.1. Spectral analysis

Figure 1 (left) shows the hybrid maps at 15 and 22 GHz with the positions of the modelfit components denoted by dotted lines. The images are registered using phase-referencing to the quasar 2005+403 at a distance of $`1.5^{}`$. On the right hand side the spectral indices between 15, 22 and 43 GHz are plotted versus the core distance, assuming N is the core.

Our analysis reveals that the core spectrum is slightly inverted with a spectral index of $`\alpha _{15/22}\alpha _{22/43}0.5`$. The jet shows a steep spectrum with spectral indices of $`\alpha _{15/22}0.5`$ to $`\alpha _{15/22}1.2`$ and $`\alpha _{22/43}1.5`$. The counter-jet has a highly inverted spectrum with $`\alpha _{15/22}`$ up to 2.5 in the inner part ($`r2`$ mas) and a flat spectrum between 15 and 22 GHz further out. Between 22 and 43 GHz the counter-jet spectrum is inverted with a spectral index of $`\alpha _{22/43}1`$. The inverted spectrum on the counter-jet side is likely due to free-free absorption by a foreground absorber. This is supported by the spectral behavior of the jet to counter-jet flux density ratio (Krichbaum et al. 1998; Bach et al. 2002). UV spectroscopy (Antonucci et al. 1994) and optical spectro-polarimetry (Ogle et al. 1997) also show evidence for a hidden broad line region. According to the unified scheme, this is strong evidence for an obscuring torus around the central engine and is consistent with the results from 21 cm absorption line VLBI (Blanco & Conway 1996).

### 3.2. Kinematics and Geometry

Figure 2 shows the component positions against time after a careful identification and linear fits to those, yielding apparent velocities $`\beta _{\mathrm{app}}`$ as follows. On the jet side the components start with 0.1 mas yr-1 near the core and accelerate to 0.25 mas yr-1 at larger separations. In Cygnus A 1 mas corresponds to 0.8 $`h^1`$ pc so that an apparent motion of 0.1 mas yr-1 corresponds to $`\beta _{\mathrm{app}}=0.26h^1`$.

The situation on the counter-jet side is less clear than on the jet side. The components are weaker and more extended than on the jet side. The innermost component C3 shows a marginal motion of $`0.010\pm 0.007`$ mas yr-1 and the next component C2 seems to travel inwards with $`0.02\pm 0.01`$ mas yr-1. In the epochs with the highest quality data, component C2 consists of two nearly equally bright components. The apparent inward motion thus might be due to blending effects of these sub-components. A summary of the component speeds is given in Table 2.

From the right panel in Fig. 2 it seems that the jet components accelerate as they travel down the jet and reach there maximum speed of $`0.5c`$ at a distance of $`r2`$ mas. Using $`\beta _{\mathrm{app}}=\frac{\beta \mathrm{sin}\theta }{1\beta \mathrm{cos}\theta }`$, which is maximized at $`\mathrm{cot}\theta =\beta _{\mathrm{app}}`$ and with $`\beta _{\mathrm{app}}=0.5h^1`$ one can calculate the lower limit of the intrinsic velocity $`\beta _{\mathrm{min}}`$ of $`0.45h^1c`$ and the corresponding angle to the line of sight of $`\theta =64^{}`$. Indeed, the analysis of the jet to counter-jet ratio also favors a large inclination of $`80^{}\pm 8^{}`$ (Krichbaum et al. 1998; Bach et al. 2002).

Due to the small relativistic effects at such conditions, the angle to the line of sight needs to change by more than $`25^{}`$ to explain the observed velocities by geometrical jet curvature. Since the jet appears to be very straight from parsec to kiloparsec scales, it is more likely that we observe a true acceleration, possibly caused by the collimation of the jet in the inner 2 parsec. Alternatively, we might observe a highly stratified jet with the different velocities belonging to different layers in the jet. This idea can also explain the absence of detectable motions on the counter-jet side. Due to the fact that we see the counter-jet from its ‘back’ the emission of the faster components is beamed away from us and we observe only the slower velocities of the outer sheath of the counter-jet.

Bach et al. (2002) detected apparent motion of $`\beta _{\mathrm{app}}=0.15\pm 0.12h^1`$ on the counter-jet side at a distance of $`r8`$ mas from the core at 5 GHz, but this optically thin component is no more visible at 15 GHz. From the assumption of a simple Königl jet (Königl 1981) and an inclination of 65 – 80 we would expect a proper motion of $`\beta _{\mathrm{app}}=0.30.4h^1`$ on the counter-jet side. That we do not see these ‘high’ velocities on the counter-jet points to a more complicated jet structure.

## 4. Conclusion

We carried out the first phase-referencing observations of Cygnus A with the VLBA at 15 and 22 GHz. The analysis of this data and a multi-frequency observation from 1996 at 15, 22 and 43 GHz revealed the spectral properties of the innermost jet structure in Cygnus A. We found a slightly inverted core and a steep jet spectrum. On the counter-jet side the inner part ($`r2`$ mas) shows a highly inverted spectrum with an $`\alpha _{15/22}`$ of up to 2.5 and $`\alpha _{22/43}1`$. In the outer regions the counter-jet shows a flat spectrum between 15 and 22 GHz and is still inverted between 22 and 43 GHz. Together with the frequency dependence of the jet to counter-jet ratio (Krichbaum et al. 1998; Bach et al. 2002) there is strong evidence that this is due to free-free absorption by an obscuring torus.

The apparent acceleration in the jet and the absence of detectable motions on the counter-jet side might reflect a more complicated jet structure than that of the simple Königl jet with a well defined jet flow and questions the assumption that the jet and counter-jet are intrinsically the same. The observations could be explained by a stratification of the jet were we observe different velocities sheaths depending on the optical depth and the orientation of the jet.

#### Acknowledgments.

We thank the group of the VLBA 2cm Survey for providing their data. This work made use of the VLBA, which is an instrument of the National Radio Astronomy Observatory, a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc. and of the 100 m telescope at Effelsberg, which is operated by the Max-Planck-Institut für Radioastronomie in Bonn.

## References

Antonucci, R., Hurt, T., & Kinney, A. 1994, Nature, 371, 313

Bach, U., et al. 2002, 6th EVN Symp., Bonn, ed. Ros E., (Bonn, MPIfR), 155

Barthel P.D. 1989, ApJ 336, 606

Blanco, P.R. & Conway, J. 1996, in Cygnus A – Study of a Radio Galaxy, ed. Carilli, C.L. & Harris, D.E., Cambridge University Press, 69

Carilli C.L., Bartel N., & Linfield R.P. 1991, AJ, 102, 1691

Fomalont, E.B. 1989 in Synthesis Imaging in Radio Astronomy, ASP Conf. Ser., Vol 6 Socorro, ed. Perley, R.A., Schwab, F.R. & Bridle, A.H., (San Francisco: ASP), 215

Kellermann, K.I., et al. 1998, AJ, 115, 1295

Königl, A. 1981, ApJ, 243, 700

Krichbaum, T.P., et al. 1998, A&A, 329, 873

Ogle, P.M., et al. 1997, ApJ, 482, L37

Perley R.A., Dreher J.W. & Cowan J.J. 1984, ApJ, 285, L35.

Zensus, J.A., et al. 2002, AJ, 124, 662


text:
# Fractionalization of a flux quantum in a one-dimensional parallel Josephson junction array with alternating 𝜋 junctions

## Abstract

We study numerically and analytically the properties of a one-dimensional array of parallel Josephson junctions in which every alternate junction is a $`\pi `$ junction. In the ground state of the array, each cell contains spontaneous magnetic flux $`\mathrm{\Phi }\mathrm{\Phi }_0/2`$ which shows antiferromagnetic ordering along the array. We find that an externally introduced $`2\pi `$-fluxon $`\mathrm{\Phi }_0`$ in such an array is unstable and fractionalizes into two $`\pi `$ fluxons of magnitude $`\frac{1}{2}\mathrm{\Phi }_0`$. We attribute this fractionalization to the degeneracy of the ground state of the array. The magnitude of the flux in the fractional fluxons can be controlled by changing the critical current of the $`\pi `$ junctions relative to the 0 junctions. In the presence of an external current, the fluxon lattice in the antiferromagnetic ground state can be depinned. We also observe a novel resonant structure in the $`V`$-$`I`$ characteristics above the depinning current due to the interaction between the fluxon lattice and the array.

One of the exciting developments in the field of Josephson devices is the fabrication of the three terminal controllable Josephson junction morpurgo . The supercurrent through a Josephson junction is given by $`I=I_c\mathrm{sin}(\mathrm{\Delta }\varphi )`$ where $`\mathrm{\Delta }\varphi `$ is the gauge-invariant phase difference between the superconductors, and the critical current $`I_c`$ depends upon the junction geometry, normal state resistance $`R_n`$ and the temperature $`T`$. Morpurgo et al.morpurgo  showed that the supercurrent through a superconductor-metal-superconductor (SNS) junction changes on passing a control current through the normal metal. For such a junction, the supercurrent $`I\mathrm{sin}(\mathrm{\Delta }\varphi +\chi )`$ where the additional phase difference $`\chi `$ is dependent on the current through the normal metal. Further theoretical work showed that in the diffusive limit of the junction, the additional phase factor $`\chi `$ can be made $`\pi `$, thus reversing the direction of the supercurrent with respect to the phase difference $`\mathrm{\Delta }\varphi `$ volkov ; wilhelm . Josephson junction with $`\chi =\pi `$ is referred as the $`\pi `$ junction (we use the term 0 junction for the Josephson junction for which $`\chi =0`$). The $`\pi `$ junction has now been realized in several experiments baselmans ; shaikh ; comment1 . The fabrication of such tunable junctions has opened immense possibilities for new applications, as demonstrated recently by the development of controllable $`\pi `$-SQUID baselmans1 .

The next natural step in this field would be to consider Josephson junction array (JJA) containing $`\pi `$ junctions. Theoretically, Kusmartsev kusum  considered a loop containing an odd number of $`\pi `$ junctions and showed that the loop contains spontaneous magnetic flux in the ground state. In the continuum limit, the long Josephson junction with alternating critical current density have been studied which shows self-generated magnetic flux mints ; gold . Recent studies of JJAs with $`\pi `$ junctions mahesh ; kulkarni  have shown some novel features arising out of the interplay between 0 and $`\pi `$ junctions. Moreover, JJA is a unique system which provides experimental realizations of several interesting physical phenomena, some examples of which are field induced superconductor to insulator transition zant , Aharonov-Casher effect elion , coherent emission of radiation barb . One is then led to ask as what new physical phenomenon exist in the 1D JJA containing $`\pi `$ junctions.

In this paper, we study numerically and analytically a new class of 1D JJA: an array of parallel Josephson junctions in which every alternate junction is a $`\pi `$ junction. The ground state contains spontaneous magnetic flux in each cell and are ordered antiferromagnetically along the array. We find that a quantum of flux (fluxon) with a $`2\pi `$ kink in the phase is unstable in such an array, and fractionalizes into two spatially separated $`\pi `$ kink fluxons. We also calculate the $`V`$-$`I`$ characteristics which shows a novel structure due to resonant interaction between the moving antiferromagnetic fluxon lattice and the linear waves emitted by the array.

Consider a 1D array of parallel Josephson junctions containing alternate $`\pi `$\- and 0 junctions (see inset Fig. 1). The Hamiltonian for this system is

$$\frac{H}{E_J}=\underset{i=0}{\overset{N1}{}}\frac{1}{2}\left(\frac{\varphi _i}{t}\right)^2+\frac{\lambda _J^2}{2}(\varphi _{i+1}\varphi _i)^2+[1(1)^i\mathrm{cos}\varphi _i],$$
(1)
where $`\varphi _i`$ is the gauge-invariant phase difference across the $`i`$-th junction. The periodic boundary condition is imposed at the two ends of the array such that $`\varphi _0=\varphi _N`$ ($`N`$ is assumed to be even). In Eq.(1), the first term represents the charging energy and the second term is the energy of the induced magnetic field due to finite self-inductance of the cell (the effect of mutual inductance between the cells is neglected). The last term represents the energy associated with the Josephson currents. The prefactor for the $`\mathrm{cos}\varphi _i`$ term alternates in sign for odd ($`\pi `$-) and even (0-) junctions. The Josephson coupling energy $`E_J=I_c\mathrm{\Phi }_0/2\pi `$ where $`I_c`$ is the critical current of a single junction. The time $`t`$ is in the units of inverse plasma frequency $`\omega _P^1=\sqrt{\frac{\mathrm{\Phi }_0C}{2\pi I_c}}`$, where $`C`$ is the averaged capacitance per unit area of the junction. The effective Josephson penetration depth is given by $`\lambda _J=(\frac{\mathrm{\Phi }_0}{2\pi L_0I_c})^{1/2}`$ where $`L_0`$ is the self inductance of a single cell. The $`\lambda _J`$ determines the screening strength of the array, and is related to the SQUID parameter $`\beta _L=\lambda _J^2`$.

From Eq.(1), the equation of motion for $`\varphi _i`$ is

$$\frac{d^2\varphi _i}{dt^2}+\alpha \frac{d\varphi _i}{dt}+(1)^i\mathrm{sin}\varphi _i+\gamma =\lambda _J^2(\varphi _{i+1}+\varphi _{i1}2\varphi _i),$$
(2)
where a dissipative term $`\alpha \frac{d\varphi _i}{dt}`$ is also added barone . The coefficient $`\alpha =\beta _c^{1/2}`$, where $`\beta _c=2\pi R_n^2I_cC/\mathrm{\Phi }_0`$ is the McCumber parameter. The parameter $`\gamma =I_{ext}/I_c`$ represents the external current through the junction. For the numerical simulation, Eq.(2) is integrated using the fourth order predictor-corrector method. The consistency of the steady state solutions was checked using different initial configurations of $`\varphi _i`$’s. The magnetic flux in the $`i^{th}`$ cell is defined as $`2\pi \frac{\mathrm{\Phi }_i}{\mathrm{\Phi }_0}=(\varphi _{i+1}\varphi _i)`$. We remark that for the case where all junctions are 0 junctions (henceforth referred to as the 0-JJA), Eq.(2) is the discrete perturbed sine-Gordon equation, and have been studied extensively ustinov1 ; ustinov2 . First, we consider the results from the numerical simulation.

Figure 1 shows the ground state flux configuration $`2\pi \mathrm{\Phi }_i/\mathrm{\Phi }_0`$ for $`N=100`$ and $`\lambda _J=1.0`$. The self-induced magnetic flux $`\mathrm{\Phi }_i`$ changes sign across neighboring cells with $`|\mathrm{\Phi }_i|=\mathrm{\Phi }`$ remaining constant. Such a configuration of $`\mathrm{\Phi }_i`$ is reminiscent of the ground state in 1D classical Ising model with antiferromagnetic (AF) coupling, as depicted schematically in the inset of Fig. 1. Therefore, we call this array the antiferromagnetic JJA (AFJJA). The AF ordering of $`\mathrm{\Phi }_i`$ (and hence, $`\varphi _i`$) implies that the self-induced screening currents in neighboring cells are oppositely oriented. The magnitude of the flux $`\mathrm{\Phi }`$ in a cell depends on the screening strength $`\lambda _J`$ as shown in the Fig. 1(inset). With increasing $`\lambda _J`$, the magnetic flux in the neighboring cells tends to overlap, and $`\mathrm{\Phi }0`$ as $`\lambda _J\mathrm{}`$. In the strong screening limit, $`\lambda _J0`$ and $`\mathrm{\Phi }\frac{\mathrm{\Phi }_0}{2}`$.

Next, we consider the consequence of introducing an external fluxon in AFJJAcomment3 . In the 0-JJA, a fluxon corresponds to a $`2\pi `$-kink in the phase profile $`\varphi (x)`$ and the magnetic field ($`\frac{\mathrm{\Delta }\varphi }{\mathrm{\Delta }x}`$) is spatially localized on the length scale $`\lambda _J`$. Figure 2(a) shows the steady state profiles of $`2\pi \mathrm{\Phi }_i/\mathrm{\Phi }_0`$ and $`\varphi _i`$ in AFJJA in the presence of a $`2\pi `$-fluxon. The self-induced magnetic field of the AF ground state has been subtracted from $`2\pi \mathrm{\Phi }_i/\mathrm{\Phi }_0`$. We find that a $`2\pi `$-fluxon is unstable in the AFJJA and fractionalizes into two spatially separated fluxons, each carrying half the quantum of flux. Also, each of the fractional fluxons is a $`\pi `$-kink in $`\varphi _i`$. The magnetic field around the fractional fluxon decays as $`\mathrm{exp}(x/\lambda _{eff})`$ where $`\lambda _{eff}2\lambda _J`$. This should be compared with the $`2\pi `$-fluxon in the 0-JJA \[Fig. 2(b)\] where $`\lambda _{eff}\lambda _J`$. The increase in $`\lambda _{eff}`$ in AFJJA is a consequence of the magnetic flux in the AF ground state.

It is possible to vary the magnitude of the magnetic flux in each fraction by changing $`i_c^{}=I_c^\pi /I_c^0`$ where $`I_c^\pi `$ and $`I_c^0`$ are the critical currents of the $`\pi `$ and 0 junctions, respectivelycomment2 . Figure 3(a) shows the spatial profile of the fractional fluxons for $`i_c^{}=0.8`$ and $`\lambda _J=1.0`$. The phase change across the fractional fluxons is not $`\pi `$ but is dependent on the value of $`i_c^{}`$. The total phase change across both the fractions is always $`2\pi `$, as required by the flux conservation. Figure 3(b) shows the magnitude of the integrated flux $`2\pi \frac{\mathrm{\Phi }_T}{\mathrm{\Phi }_0}`$ in each fractional fluxons as a function of $`i_c^{}`$ for $`\lambda _J=1.0`$. The fractionalization occurs for $`i_{cl}^{}<i_c^{}<i_{cu}^{}`$, where $`i_{cl}^{}`$ and $`i_{cu}^{}`$ are the two critical values. The slope $`\frac{d\mathrm{\Phi }_T}{di_c^{}}`$ at the critical values $`i_{cl}^{}`$ and $`i_{cu}^{}`$ appears to diverge, suggesting a transition between the fractionalized state and the single fluxon state. In experiments, the magnitude of the flux at the center of the fluxon $`\mathrm{\Phi }_m`$ can be measured more easily. Figure 3(b) shows the behavior of $`\mathrm{\Phi }_m(i_c^{})`$.

In Fig. 4, we show the numerically obtained parameter space $`\lambda _J`$-$`i_c^{}`$. We have assumed that $`i_c^{}`$ can be varied independent of $`\lambda _J`$. The region of fractional fluxons is bounded by $`i_{cl}^{}(\lambda _J)`$ and $`i_{cu}^{}(\lambda _J)`$. It is easy to understand the absence of fractional fluxons in the limit $`i_c^{}0`$ since the array becomes a 0-JJA (with lattice constant twice the original array) which allows only $`2\pi `$-fluxons. In the opposite limit $`i_c^{}\mathrm{}`$ such that $`I_c^00`$ and $`I_c^\pi `$ is finite, there are two $`\pi `$ junctions in each cell and the array can be shown to be equivalent to the 0-JJA, and the fractionalization is again not expected. In obtaining the parameter space in Fig. 4, $`I_c^0`$ is assumed to be finite and fixed which leads to fractionalization for $`\lambda _J<0.7`$ even as $`i_c^{}\mathrm{}`$.

The simulation results discussed above can also be understood analytically. Consider the case of $`i_c^{}=1`$. Define

$$\varphi _{2m+1}=u_m,\mathrm{and}\varphi _{2m}=v_m,$$
(3)
where $`m=0,\frac{N}{2}1`$. Thus, $`u_m`$ and $`v_m`$ are the gauge invariant phase differences across the $`\pi `$ and 0 junctions, respectively. In the absence of any external fluxon, the $`u_m`$ and $`v_m`$ are invariant on translation by the lattice vector along the array. Hence, substituting $`u_m=u`$ and $`v_m=v`$, Eq.(2) becomes

$$\mathrm{sin}u=\mathrm{sin}v=2\lambda _J^2(uv).$$
(4)
There are two trivial solutions of the Eq.(4): $`u=v=0`$ and $`u=v=\pi `$. The non-trivial solution of Eq.(4) is given by

$$u=\pi v,\mathrm{and}\mathrm{sin}v=2\lambda _J^2(\pi 2v).$$
(5)
For a given value of $`\lambda _J`$, the quantities $`v`$ and $`u`$ can be calculated graphically from Eq.(5). It can be easily verified that the non-trivial solution is the ground state for any finite $`\lambda _J`$. The magnetic flux in the cell is given by $`\mathrm{\Phi }=\pm (uv)\frac{\mathrm{\Phi }_0}{2\pi }`$, where the $`+`$ and $``$ sign is for the cell to the left and the right of the $`\pi `$ junction, respectively. Thus, the magnetic flux alternates in sign along the array. The values of $`u`$, $`v`$ and $`|\mathrm{\Phi }|`$ obtained from Eq.(5) are in excellent agreement with the numerically obtained values. In the strong screening limit $`\lambda _J0`$, $`u=\pi `$ and $`v=0`$, and the flux in each cell attains the maximum value $`\mathrm{\Phi }_0/2`$. In the limit $`\lambda _J\mathrm{}`$, the solution of Eq.(5) is $`u=v=\pi /2`$ which is degenerate to the trivial solutions of Eq.(4), and the ground state contains no spontaneous magnetic flux.

To understand the fractionalization of a $`2\pi `$-fluxon, we note that the AF ground state solution of Eq.(4) is two-fold degenerate. If $`\varphi =\{u,v\}`$ obtained from Eq.(5) is one solution, the other solution is obtained by translating $`\varphi `$ by one lattice constant. Thus, the other solution is $`\varphi ^{}=\{u^{},v^{}\}`$, where $`u^{}=v+\pi `$ and $`v^{}=u+\pi `$. The degeneracy of the ground state has an important implication: the elementary excitation for the array is a kink (or domain wall) in $`\varphi _i`$ which interpolates between the two degenerate ground states $`\varphi `$ and $`\varphi ^{}`$. It can be easily verified that the phase change across the ‘kink’ is $`\pi `$ and corresponds to an additional flux $`\mathrm{\Phi }_0/2`$ in the array. In the presence of a $`2\pi `$ fluxon, the energy minimization leads to two $`\pi `$-kinks which are separated by the degenerate ground states. This is shown schematically in the inset of Fig. 2(a). Fractionalization of a $`2\pi `$ fluxon in AFJJA is thus a consequence of the degeneracy of the ground state. A similar phenomenon is observed in polyacetylene su  and certain field theories jackiw ; rajaraman  where fractional topological excitations occur and are related to the ground state degeneracy.

Equations (4) and (5) can be extended for the case $`i_c^{}1`$. We find that the AF ground state is stable and fluxon fractionalization occurs for $`i_c^{}>1`$ when $`\lambda _J<\lambda _{Ju}=\frac{1}{\sqrt{2}}(1\frac{1}{i_c^{}})^{1/2}`$, and for $`i_c^{}<1`$ when $`\lambda _J<\lambda _{Jl}=\frac{1}{\sqrt{2}}(\frac{i_c^{}}{1i_c^{}})^{1/2}`$. This is plotted in Fig. 4 and is in good agreement with the numerically obtained behavior of $`i_{cu}^{}(\lambda _J)`$ and $`i_{cl}^{}(\lambda _J)`$. Further details of the analytical calculations will be given elsewhere.

Next, we study the dynamical properties of the AFJJA. We restrict the analysis to the case $`i_c^{}=1.0`$ and in the absence of any external fluxon. The $`V`$-$`I`$ curve is obtained by sweeping $`\gamma `$ in small steps, and calculating $`V=\alpha \frac{d\varphi }{dt}`$ in the steady state ($`V`$ is in units of $`R_nI_c`$). Recall that for the 0-JJA, all junctions switch from the superconducting state to the normal state simultaneously at $`\gamma =1`$. In the AFJJA, the magnetic flux in the ground state alters this behavior significantly.

Figure 5 shows the $`V`$-$`I`$ curve for $`N=30`$ and $`\lambda _J=1`$. The $`\alpha =0.1`$ for the rest of the discussion belowcomment4 . The important feature of the $`V`$-$`I`$ curve is the appearance of a plateau in $`V`$ above a depinning current $`\gamma _c`$. The transition to the normal state occurs at a higher current $`\gamma _n`$. By analyzing the spatio-temporal dynamics of fluxons on the voltage plateau, we find that the interpenetrating lattice of the flux ($`\mathrm{\Phi }`$) and the antiflux ($`\mathrm{\Phi }`$) move in the opposite directions which appears as a stationary wave of breathing flux-antiflux pairs likharev . This is evident from the 2D space-time plot of $`2\pi \frac{\mathrm{\Phi }_i}{\mathrm{\Phi }_0}`$ which is shown for a section of the array in the inset of Fig. 5. The $`\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}`$ in each cell oscillates between the positive and the negative values, and is in antiphase with the neighboring cells. Thus, at any instant of time, the total flux in the array is zero as expected from the ground state. We also find that for $`\lambda _J>1`$, a linear flux flow regime appears before the plateau, whereas for small $`\lambda _J`$, the $`V`$ shows a sharp step as shown in the inset of Fig. 5. The $`\lambda _J`$ dependence of the width of the voltage plateau $`(\mathrm{\Delta }\gamma )_s=\gamma _n\gamma _c`$ is shown in Fig. 6(a). The $`(\mathrm{\Delta }\gamma )_s`$ is non-monotonic and is maximum for $`\lambda _J2`$. Below a cut-off $`\lambda _J^{}0.45`$, the plateau in $`V`$ disappears and $`\gamma _c=\gamma _n`$.

The origin of the plateau (or step) in $`V`$ can be understood from the fluxon dynamics. In a 0-JJA, the motion of a single $`2\pi `$-fluxon leads to the emission of small amplitude linear waves (plasma waves) due to the discreteness of the array. The resonances between these linear waves and the periodic motion of the fluxon causes a series of plateaus in $`V`$ ustinov1 . Extending this to the case of AFJJA, the plateau in $`V`$-$`I`$ can be attributed to the phase locking between the moving fluxon lattice and the linear waves emitted by the array. The frequency $`\omega _s`$ of the linear waves can be calculated from Eq. (2). In the absence of any external fluxon, the symmetry of the ground state allows only waves with the wavevector $`k=2\pi /(2a)`$ to be coupled resonantly to the moving fluxon lattice (here $`a`$ is the lattice constant). Thus, all 0 and $`\pi `$ junctions oscillate with the same amplitude and phase. Linearizing Eq.(2) for the $`\pi `$ and the 0 junctions using $`u=u_0\mathrm{exp}(i\omega _st)`$ and $`v=v_0\mathrm{exp}(i\omega _st)`$, respectively,

$`\omega _s^2uu`$ $`=`$ $`2\lambda _J^2(vu),`$
$`\omega _s^2v+v`$ $`=`$ $`2\lambda _J^2(uv).`$ (6)
For simplicity, we have used $`\alpha =\gamma =0`$. Adding the above equations gives the relation between $`u`$ and $`v`$, $`u=v(1\omega _s^2)/(1+\omega _s^2)`$. Substituting for $`u`$ in the second equation leads to

$$(1\omega _s^2)(1+\omega _s^2)+4\lambda _J^2\omega _s^2=0.$$
(7)
The above quadratic equation in $`\omega _s^2`$ can be solved to obtain the frequency of the linear waves

$$\omega _s(\lambda _J)=\sqrt{2\lambda _J^2+\sqrt{1+4\lambda _J^4}}.$$
(8)
The $`\omega _s`$ is in units of the plasma frequency $`\omega _P`$. The condition for phase locking of the linear waves with the moving fluxon lattice then becomes $`\omega _sT=2\pi `$, where $`T`$ is the time period corresponding to the motion of the fluxon lattice. From the simulation, the $`T`$ (and hence $`\omega _s`$) on the plateau can be obtained from the time evolution of the magnetic flux $`\mathrm{\Phi }_i(t)`$ in a cell. In Fig. 6(b), the $`\omega _s`$ calculated analytically is compared with that obtained from the simulation for increasing $`\lambda _J`$. A reasonable agreement is observed over a range of $`\lambda _J`$. The discrepancy appears as $`\lambda _J`$ approaches $`\lambda _J^{}`$ below which the plateau in $`V`$ is not observed. The voltage $`V_s`$ on the plateau is given by $`V_s=\alpha \omega _s`$. This is shown in Fig. 6(c) and is in good agreement with the simulation.

Finally, the dependence of $`\gamma _c`$ on the screening strength is shown in Fig. 6(d). A finite $`\gamma _c`$ in JJA is attributed to the pinning potential created at the center of the cell lobb . For AFJJA, the pinning potential or the energy barrier can be defined as $`\mathrm{\Delta }E=E_ME_A`$, where $`E_M`$ is the energy of the array with the fluxon lattice placed on the junctions and $`E_A`$ is the ground state energy. Equating the pinning force $`2\mathrm{\Delta }E`$ to the Lorentz force required to overcome the energy barrier gives the depinning current $`\gamma _c`$ (per junction),

$$\gamma _c=\frac{1}{\pi }\left(2\mathrm{cos}(v)\frac{1}{2}\lambda _J^2(\pi 2v)^2\right).$$
(9)
The phase $`v`$ in the above equation is the solution of the Eq.(5). Fig. 6(d) compares the above expression with the $`\gamma _c`$ obtained from the simulation. The agreement is good at large $`\lambda _J`$, and the deviation appears as $`\lambda _J\lambda _J^{}`$ below which no flux flow is observed.

In conclusion, we have introduced a new class of JJA containing $`\pi `$ junctions. We considered the one-dimensional case in which every alternate junction is a $`\pi `$ junction, and studied its properties numerically and analytically. The ground state of the array contains spontaneous magnetic flux in each cell and are ordered antiferromagnetically along the array. A $`2\pi `$-fluxon in such an array is unstable and fractionalizes into two spatially separated $`\pi `$-fluxons. The fractionalization is related to the ground state degeneracy. The $`V`$-$`I`$ curve shows a voltage plateau due to resonant interaction between the linear waves emitted by the array and the moving antiferromagnetic fluxon lattice present in the ground state.

Acknowledgments : Part of this work was done during the authors stay at the University of California, Davis.

Note : After the completion of the work, we became aware of the experimental fabrication of an array of coupled $`\pi `$-loops using YBaCuO-Nb zigzag structure in Ref.hilg . Some of the results obtained in this paper are applicable for such structures also. We thank the referees for bringing this paper to our attention.


text:
# Integrable spin-boson interaction in the Tavis-Cummings model from a generic boundary twist

## Abstract

We construct models describing interaction between a spin $`s`$ and a single bosonic mode using a quantum inverse scattering procedure. The boundary conditions are generically twisted by generic matrices with both diagonal and off-diagonal entries. The exact solution is obtained by mapping the transfer matrix of the spin-boson system to an auxiliary problem of a spin-$`j`$ coupled to the spin-$`s`$ with general twist of the boundary condition. The corresponding auxiliary transfer matrix is diagonalized by a variation of the method of $`Q`$-matrices of Baxter. The exact solution of our problem is obtained applying certain large-$`j`$ limit to $`su(2)_j`$, transforming it into the bosonic algebra.

Introduction Models representing interactions between a single bosonic mode and spin degrees of freedom find application in many different contexts. In atomic physics they describe atoms interacting with electromagnetic field  and many phenomena like spontaneous emissions in cavity  and Rabi oscillations  in two-level atoms are captured by an important representative of the models mentioned above, the Jaynes-Cummings (J-C) model . Recently the J-C dynamics was intensively studied in the research field of ions in harmonic traps  and then in quantum computation . Finally this kind of models found applications in quasi-2D semiconductors in transverse magnetic field .
Toy model Hamiltonians representing a single bosonic mode interacting with a spin $`s`$ are of the following type

$$H=H_{sph}+\gamma \left(S^+a^{}+S^{}a\right)$$
(1)
where $`H_{sph}=\omega a^{}a+\alpha S^z+\beta \left(S^+a+S^{}a^{}\right)`$; operators $`a^{}`$ and $`a`$ are bosonic operators that commute with spin operators $`S^v`$, $`v=\{z,\pm \}`$. The model (1) with $`\gamma =\beta `$ was originally proposed for quantum optics purposes to describe a dipole-like interaction in atoms-radiation systems and it is known as the Tavis-Cummings (T-C) model; the model reduces to the J-C one for $`s=1/2`$. In solid state physics models of type (1) can describe certain quantum circuits .

There are two cases in which the model can be solved analytically: i)In the limit $`s\mathrm{}`$ there are exact results . They were applied to study the entanglement across the quantum phase transition between normal to super-radiant phase. ii) For “single photon” interactions the model can be simplified employing the Rotating Wave Approximation (RWA) that neglects the so called “counter rotating” terms: $`S^+a^{}`$, $`S^{}a`$. Within the RWA the model (1) can be solved exactly. For generic parameters $`\alpha ,\beta ,\gamma ,\omega `$, and for finite $`s`$ the model, as it stands in Eq. (1), is non-integrable. Merging the model into the main stream of the Quantum Inverse Scattering (QIS) method constitues often a guide to discover unsuspected exactly solvable models with sufficiently generic interaction. According to this method the Hamiltonian is obtained as output of the procedure that remarkably ensures the integrability of the theory. In the simplest cases the QIS method provide integrable Hamiltonians after periodic boundary condition are imposed. The variety of the integrable models can be considerably enriched by considering more general boundary conditions . By this is meant that the monodromy matrix is multiplied by non-trivial matrices that ultimately cause the presence of boundary terms in the Hamiltonian. Of interest in the present paper is the case of constant boundary matrix; this realize the, so called, twisted boundary conditions. For the type of models under consideration the QIS method was employed in Ref. where nonlinear generalizations of $`H_{sph}`$ were studied. These generalizations were obtained by twisting the boundary conditions. The twist matrices were chosen as the same diagonal matrix for both the bosonic and spin degrees of freedom; finally in a certain sense (specified below) they are classical. As a result, although nonlinear, these models contain the standard interaction. Here a more general interaction is obtained applying more general boundary conditions: the twist matrices $`K_B`$ and $`K_S`$ are non-diagonal, different for the boson and the spin, and of “quantum” nature (see (5)) . The Hamiltonian we found is Eq. (9). The QIS method pave the way towards the exact solution of the theory through Bethe Ansatz. In the cases where there exist an obvious “reference” state a direct (algebraic) Bethe Ansatz approach can be applied. For the model we found here, however, there is no simple vacuum state since (9) does not commute with $`S^z+a^{}a`$. A standard route to attack the exact solutions of such kind of spectral problems is to apply the technique that Baxter  invented to obtain the eigenvalues without knowledge of the eigenstates. We obtain the eigenvalues (the calculation of the exact eigenstates will be the object of a future publication) in the following way. We first define an auxiliary problem consisting of two spins with two distinct representations $`j`$ and $`s`$; the boundary conditions are generically twisted; the spin $`j`$ is affected by an “impurity” $`\nu `$. We diagonalize the auxiliary problem by adapting the Baxter method to it. Then the solution of the spin-bosonic problem is obtained performing certain $`j\mathrm{}`$ limit (see (14)) in the results for the auxiliary spin-spin problem. The eigenvalues are given in (23) and the parameters $`\lambda _k`$ are fixed by (24).

The paper is laid out as follows. In the next section we derive the integrable model. In the section III the exact solution is obtained. The section IV is devoted to our conclusions.

Integrability. The starting point of the QIS method is to define quantum Lax matrices $`L(\lambda )`$ and a scattering matrix $`R(\lambda )`$ satisfying the Yang Baxter (YB) equation: $`R(\lambda \mu )L(\lambda )L(\mu )=L(\mu )L(\lambda )R(\lambda \mu ),`$ where $`\lambda `$ is the spectral parameter. For the present case, the Lax operators $`L`$ we consider are

$$L_S(\lambda ):=\left(\begin{array}{cc}\lambda \eta S^z& \eta S^+\\ \eta S^{}& \lambda +\eta S^z\end{array}\right),$$
(2)
$$L_B(\lambda ):=\left(\begin{array}{cc}\lambda \mathrm{\Delta }\eta ^1\eta a^{}a& a^{}\\ a& \eta ^1\end{array}\right),$$
(3)
each satisfying the YB equation with: $`R(\lambda ;\eta )=\eta 1\mathrm{l}1\mathrm{l}+\lambda \mathrm{IP},`$ where $`\eta \text{ }\mathrm{R}`$ and $`\mathrm{IP}`$ is the permutation: $`\mathrm{IP}AB\mathrm{IP}=BA`$. The monodromy matrix is

$$T(\lambda )=K_BL_B(\lambda )K_SL_S(\lambda ),$$
(4)
where $`K_B`$ and $`K_S`$ are $`\text{ }\mathrm{C}`$-number matrices that produce boundary terms (without “internal dynamics”, the matrices $`K`$ not depending on $`\lambda `$). Notice that we have two different boundaries each for the spin and for the boson. In Refs. $`K_BK_S`$ is assumed; the T-C model (without counter rotating terms) $`H_{sph}`$ is obtained for $`\eta 0`$. The matrix $`T(\lambda )`$ fulfills the YB relation: $`R(\lambda \mu )T(\lambda )T(\mu )=T(\mu )T(\lambda )R(\lambda \mu ),`$ due to the fact that $`[R,K_BK_B]=[R,K_SK_S]=0`$ holds for any numeric matrix because of the $`sl(2)`$ symmetry of the $`R`$-matrix. The transfer matrix is defined as $`t(\lambda ):=tr_{(0)}T(\lambda )`$ where $`tr_{(0)}`$ means trace in the auxiliary space. $`t(\lambda )`$ is a generating functional of integrals of motion since: $`[t(\lambda ),t(\mu )]=0`$. For the present case the transfer matrix can be chosen as a polynomial in $`\eta `$: $`t(\lambda )=_{l=g}^h\eta ^lt_l(\lambda )`$; then the coefficients of $`[t(\lambda ),t(\mu )]=_{l=2g}^{2h}\eta ^lC_l(\lambda ,\mu )`$ vanish at any $`\eta `$-power. The following assumption is crucial for our purposes: The entries of the matrices $`K`$ depend on $`\eta `$

$$K_{Xij}=K_{Xij}^{(0)}+K_{Xij}^{(1)}\eta +\mathrm{},X=\{B,S\}.$$
(5)
The parameter $`\eta `$ is usually called “quantum parameter” since it controls the limit how to recover the classical scattering matrix $`r(\lambda )`$ from the matrix $`R(\lambda )`$. In this sense our matrices $`K`$ in (5) describe “quantum systems” which we couple to the the boson and to the spin at the boundary (twist matrices that are independent on $`\eta `$ might be considered as classical boundaries). In brief, the main idea of our procedure is to play with the boundaries $`K`$ in such a way that $`C_l(\lambda ,\mu )=0[t_l(\lambda ),t_l(\mu )]=0`$ for certain $`l`$. In order to obtain the model we are interested in, the degree and the coefficients of the polynomials are fixed such that: i) $`t_m(\lambda )`$ describes an integrable model, then $`t_l(\lambda )`$ must be $`\text{ }\mathrm{C}`$-numbers for all $`l<m`$; ii) the model results containing the counter-rotating terms; iii) the obtained operator is Hermitian. All these conditions translate in a system of equations for the entries of the matrices $`K`$; we shall see that these parameters will be the coupling constants of the Hamiltonian. It turns out to be sufficient to consider the entries of the $`K`$ matrices to be linear in $`\eta `$. Such entries are restricted to

$`K_{B\mathrm{\hspace{0.17em}11}}^{(0)}=K_{B\mathrm{\hspace{0.17em}22}}^{(0)},K_{S\mathrm{\hspace{0.17em}11}}^{(0)}=K_{S\mathrm{\hspace{0.17em}22}}^{(0)},K_{B\mathrm{\hspace{0.17em}21}}^{(0)}=K_{B\mathrm{\hspace{0.17em}12}}^{(0)},K_{S\mathrm{\hspace{0.17em}12}}^{(0)}=K_{S\mathrm{\hspace{0.17em}21}}^{(0)},K_{S\mathrm{\hspace{0.17em}12}}^{(0)}={\displaystyle \frac{K_{B\mathrm{\hspace{0.17em}21}}^{(0)}K_{S\mathrm{\hspace{0.17em}22}}^{(0)}}{K_{B\mathrm{\hspace{0.17em}11}}^{(0)}}},`$ (6)
$`K_{S\mathrm{\hspace{0.17em}11}}^{(1)}={\displaystyle \frac{K_{B\mathrm{\hspace{0.17em}22}}^{(1)}K_{S\mathrm{\hspace{0.17em}22}}^{(0)}}{K_{B\mathrm{\hspace{0.17em}11}}^{(0)}}}+{\displaystyle \frac{K_{B\mathrm{\hspace{0.17em}11}}^{(0)}}{K_{B\mathrm{\hspace{0.17em}12}}^{(0)}}}(K_{S\mathrm{\hspace{0.17em}21}}^{(1)}K_{S\mathrm{\hspace{0.17em}12}}^{(1)}){\displaystyle \frac{K_{S\mathrm{\hspace{0.17em}11}}^{(0)}}{K_{B\mathrm{\hspace{0.17em}12}}^{(0)}}}(K_{B\mathrm{\hspace{0.17em}12}}^{(1)}K_{B\mathrm{\hspace{0.17em}21}}^{(1)})+{\displaystyle \frac{K_{S\mathrm{\hspace{0.17em}22}}^{(0)}}{K_{B\mathrm{\hspace{0.17em}11}}^{(0)}}}K_{B\mathrm{\hspace{0.17em}11}}^{(1)}K_{S\mathrm{\hspace{0.17em}22}}^{(1)}.`$ (7)
We take $`t_1(\lambda )`$ as the Hamiltonian

$`t_1(\lambda )H=W(\lambda )S^z+\lambda (U+V)a^{}a+2\left[Y+\sqrt{UV}(\mathrm{\Delta }\lambda )\right]S^x+Z\lambda (a+a^{})`$ (8)
$`U(aS^++a^{}S^{})+V(a^{}S^++aS^{})2\sqrt{UV}(a+a^{})S^z,`$ (9)
where the couplings are

$`W(\lambda )={\displaystyle \frac{(U+V)Z}{\sqrt{UV}}}+(VU)(\mathrm{\Delta }\lambda {\displaystyle \frac{Y}{\sqrt{UV}}}),`$ (13)
$`U=K_{B\mathrm{\hspace{0.17em}22}}^{(0)}K_{S\mathrm{\hspace{0.17em}22}}^{(0)},V=K_{B\mathrm{\hspace{0.17em}22}}^{(0)}(K_{S\mathrm{\hspace{0.17em}21}}^{(0)})^2/K_{S\mathrm{\hspace{0.17em}22}}^{(0)}`$
$`Y=K_{B\mathrm{\hspace{0.17em}21}}^{(0)}\left(K_{B\mathrm{\hspace{0.17em}11}}^{(1)}K_{S\mathrm{\hspace{0.17em}21}}^{(0)}/K_{B\mathrm{\hspace{0.17em}22}}^{(0)}K_{S\mathrm{\hspace{0.17em}21}}^{(1)}\right)K_{B\mathrm{\hspace{0.17em}12}}^{(0)}K_{S\mathrm{\hspace{0.17em}22}}^{(1)}K_{B\mathrm{\hspace{0.17em}22}}^{(0)}K_{S\mathrm{\hspace{0.17em}12}}^{(1)},`$
$`Z=Y+K_{B\mathrm{\hspace{0.17em}22}}^{(0)}\left(K_{S\mathrm{\hspace{0.17em}12}}^{(1)}+K_{S\mathrm{\hspace{0.17em}21}}^{(1)}\right)+2K_{B\mathrm{\hspace{0.17em}22}}^{(1)}K_{S\mathrm{\hspace{0.17em}21}}^{(0)}.`$
The coupling constants obey (13) for the model to be integrable; parameters $`\mathrm{\Delta }`$, $`\lambda `$, $`X`$, $`Z`$ can be set freely; the quantity $`UV`$ must be positive. Nevertheless the rotating and counter-rotating terms can be adjusted to have the same sign by acting on the operators: $`S^{}S^{}e^{i\frac{\pi }{2}}`$ and $`aae^{i\frac{\pi }{2}}`$; the third, fourth and last terms in Eq. (9) are transformed accordingly. We observe that the simultaneous presence of rotating and counter-rotating terms preserves the integrability only if a further term $`(a+a^{})S^z`$ appears in the model (the term $`a+a^{}`$ can be transformed out by a translation: $`aa+\xi `$ with $`\xi =Z/(U+V)`$; the coefficients of $`S^x`$ and $`S^z`$ are shifted by $`\xi (UV)`$ and $`2\xi \sqrt{UV}`$ respectively). Restricting the boundary conditions: $`K_B=K_S`$ (non-diagonal) induces the further constraint $`U=V`$. In this case our model reduces to the Rashba Hamiltonian in a constant magnetic field  (the bosonic number labeling the Landau levels; see also Ref.).

The constants of motion of (9) are $`t_1(0)`$ and $`_\lambda t_1(\lambda )`$ (only two of $`H`$, $`t_1(0)`$, $`_\lambda t_1(\lambda )`$ are independent). $`_\lambda t_1(\lambda )`$ can be easily diagonalized : $`𝒰_\lambda t_1(\lambda )𝒰^1|\varphi =((U+V)(n+m)M_{})|\varphi `$. In the $`|\varphi `$ basis the Hilbert space of the Hamiltonian blocks into invariant subspaces labelled by the bosonic number $`n0`$ and with $`S^z|\varphi =m|\varphi `$; $`M_{}Z^2/(U+V)`$. This will be used to classify the excitations in the Bethe equations (24).

Exact eigenvalues. To diagonalize the model (9) we define an auxiliary inhomogeneous spin problem. We use the property that a spin $`j`$-$`su(2)`$ can be contracted to the Weyl-Heisenberg algebra  through the singular limit $`\epsilon \mathrm{}`$ of a Dyson-Maleev transformation

$$\{\eta J^{},\frac{1}{\eta \epsilon ^2}J^+,J^z\}\{a^{},a,a^{}a\frac{\epsilon ^2}{2}\};$$
(14)
such a limit corresponds to $`j=\epsilon ^2/2\mathrm{}`$. The bosonic Lax matrix is thus expressed as limit of a spin-$`j`$ Lax matrix:

$$L_B(\lambda )=\underset{\epsilon \mathrm{}}{lim}(\frac{1}{\eta \epsilon })K(\epsilon )\sigma _y\sigma _zL_J(\lambda \nu )\sigma _y,$$
(15)
where $`K(\epsilon )=diag\{\eta \epsilon ,(\eta \epsilon )^1\}`$ and

$$L_J(\lambda \nu )=\left(\begin{array}{cc}\lambda \nu \eta J^z& \eta J^+\\ \eta J^{}& \lambda \nu +\eta J^z\end{array}\right),$$
(16)
the “inhomogeneity” parameter being set to

$$\nu =\eta \epsilon ^2/2+\eta ^1+\mathrm{\Delta }.$$
(17)
Thus the monodromy matrix Eq (4) can be written as $`T(\lambda )=lim_ϵ\mathrm{}T_a(\lambda )`$ where $`T_a`$ is an auxiliary monodromy matrix defined as

$$T_aK_JL_J(\lambda \nu )\sigma _yK_SL_S(\lambda ),$$
(18)
with $`K_J1/(\eta \epsilon )K_BK(ϵ)\sigma _y\sigma _z`$. $`t_a=tr_0\{T_a\}`$ can be diagonalized adapting the Baxter method  for off-diagonal twisted spin-$`s`$ chain . In the present case the “chain” consists of only two sites; the twist matrices are distinct and containing both diagonal and off-diagonal entries. The details of the calculations will be reported elsewhere. The Baxter equation reads

$`t_a(\lambda )Q(\lambda )=R^{()}(\lambda s\eta )(\lambda \nu j\eta )Q(\lambda +\eta )+R^{(\pm )}(\lambda +s\eta )(\lambda \nu +j\eta )Q(\lambda \eta ),`$ (19)
where $`R^{(\pm )}\frac{1}{2}\left[tr(K_J\sigma _yK_S)\pm \sqrt{(tr(K_J\sigma _yK_S))^2+4det(K_JK_S)}\right]`$. The quantities $`Q(\lambda )`$ are $`(2s+2j+2)\times (2s+2j+2)`$ matrices and constructed in the standard way; they fulfill $`[Q(\lambda ),Q(\mu )]=[Q(\lambda ),t_a(\lambda )]=0`$. Eq (19) fixes the eigenvalue $`\tau _a(\lambda )^{(\pm )}`$ of the transfer matrix

$`\tau _a(\lambda )^{(\pm )}=R^{()}(\lambda s\eta )(\lambda \nu j\eta )`$ $`{\displaystyle \underset{i=1}{\overset{2s+2j}{}}}{\displaystyle \frac{\lambda \lambda _i+\eta }{\lambda \lambda _i}}`$ (21)
$`+R^{(\pm )}(\lambda +s\eta )(\lambda \nu +j\eta ){\displaystyle \underset{i=1}{\overset{2s+2j}{}}}{\displaystyle \frac{\lambda \lambda _i\eta }{\lambda \lambda _i}},`$
where the variables $`\lambda _j`$ are solutions of the equations

$$\frac{R^{()}}{R^{(\pm )}}\frac{(\lambda _ks\eta )(\lambda _k\nu j\eta )}{(\lambda _k+s\eta )(\lambda _k\nu +j\eta )}=\underset{\genfrac{}{}{0pt}{}{i=1}{ik}}{\overset{2s+2j}{}}\frac{\lambda _k\lambda _i\eta }{\lambda _k\lambda _i+\eta },k=1,\mathrm{},2(s+j).$$
(22)
To obtain the solution of the bosonic problem we perform the algebraic contraction at the level of Eqs. (22). This can be done expanding $`R^{(\pm )}=_mR_{2m}^{(\pm )}/(\epsilon \eta )^{2m}`$ and taking into account of Eq. (17); a generalizations of the Bethe equations found in  are obtained. By a further expansion $`R_{2m}^{(\pm )}=_{n=0}^2R_{2m,n}^{(\pm )}\eta ^n`$ in the latter equations, the eigenvalues $``$ of the model (9) can be obtained as their linear terms in $`\eta `$:

$`=_0(\lambda )+{\displaystyle \underset{j=1}{\overset{2s+n}{}}}{\displaystyle \frac{1}{\lambda \lambda _j}}\left[\lambda x_0^{}{\displaystyle \underset{lj}{\overset{2s+n}{}}}{\displaystyle \frac{1}{\lambda _j\lambda _l}}\lambda (\lambda \mathrm{\Delta })y_0^{}\lambda x_1^++sx_0^{}\right],`$ (23)
where $`_0(\lambda )=(\lambda \mathrm{\Delta })(\lambda y_1^++sy_0^{})\lambda x_2^{}+sx_1^+`$; with $`x_n^\alpha =R_{0,n}^{()}\alpha R_{0,n}^{(\pm )}R_{2,n}^{()}`$ and $`y_n^\alpha =R_{0,n}^{(\pm )}+\alpha R_{0,n}^{()}`$ $`\alpha =\{+,\}`$. The quantities $`\lambda _k`$ are fixed by

$$\frac{sx_0^{}}{\lambda _k}\left(\lambda _k\mathrm{\Delta }\right)y_0^{}=x_0^{}\underset{lk}{\overset{2s+n}{}}\frac{1}{\lambda _k\lambda _l}+x_1^+,$$
(24)
where $`k=1,\mathrm{},2s+n`$; the quantum number $`n`$ labels the excitations as discussed before. For generic $`2s+n`$ the equations above can be solved numerically. Alternatively the quantities $`\lambda _k`$ can be obtained as roots of the polynomial $`P(\lambda )=_{s=1}^{2s+n}(\lambda \lambda _s)`$ satisfying

$$\lambda P^{\prime \prime }(\lambda )+\frac{2}{x_0^{}}\left[y_0^{}\lambda ^2+(y_0^{}\mathrm{\Delta }+x_1^+)\lambda sx_0^{}\right]P^{}(\lambda )\left[\zeta (2s+n)\lambda \right]P(\lambda )=0$$
(25)
where $`\zeta `$ is fixed by imposing that $`\lambda =0`$ is a simple root of $`\lambda P(\lambda )`$: $`\zeta =sx_0^{}\frac{P^{}(0)}{P0}`$ .

Conclusions. By the Quantum Inverse scattering method we have constructed integrable T-C models with twisted boundary conditions. The twist matrices are generic in the sense that they contain both diagonal and non-diagonal entries. They are responsible for the presence of rotating and counter-rotating terms in the Hamiltonian. The spectrum is computed through the Baxter method. As far as we know this method is applied to spin-boson systems for the first time; the subtleties related to the bosonic limit, recovered for infinite spin length are dug out. Integrability and exact solution can be obtained provided that a further term $`(a+a^{})S^z`$ is considered. Interestingly enough we found a global rotation of the spin/bosonic degrees of freedom such that the rotating terms (alternatively, the counter-rotating terms) are compensated out. We conjecture that “true” counter-rotating terms in the Tavis-Cummings model could be inserted considering dynamical boundaries: $`K_X=K_X(\lambda )`$; alternatively one should consider $`XYZ`$ symmetry of the scattering matrix. These terms serve to a reliable description of certain systems in quantum optics or to model the spin-orbit interaction in heterostructures where the simultaneous Rashba and Dresselhaus terms are important. Our paper could pave the way to costruct integrable Hamiltonians for such physical situations. As immediate application, we notice that the Hamiltonian (9) describes the quantum circuit of Fig.(1). Two coupled dc-Superconducting Quantum Interference Devices (SQUIDs) are coupled inductively. The primary device $`p`$ is intended built with large Josephson junctions to be described by a classical SQUID Hamiltonian (whose degrees of freedom are, then, bosonic) flowed by the current $`I_p`$ (this circuit plays the role of the $`LC`$ resonant circuit of Ref. ); the secondary SQUID $`s`$, with small junctions, is accommodated inside the primary and pierced by the magnetic flux: $`\mathrm{\Phi }=\varphi _{ext}+L_pI_p`$ ($`L_p`$ is the inductance of the circuit). Thus the secondary is a quantum SQUID controlled by the classical one.

The effective Josephson coupling of the quantum SQUID depends on the flux $`E_J^s(\mathrm{\Phi })E_J^s(\varphi _{ext})+\stackrel{~}{L}_pI_p`$ . This kind of setups are intensively studied as controllable flux-qubits to data-bus transferring in many protocols of quantum computation. The circuit Hamiltonian is

$$_{circuit}=\omega _pa^{}a2E_J^s(\varphi _{ext})S^x2\stackrel{~}{L}_p(a+a^{})S^xiM(aa^{})S^y+2V_C(a+a^{})S^z$$
(26)
where $`\omega _p`$ is the “frequency” of the primary SQUID (or the natural frequency of the resonant circuit), $`V_C`$ is due to the capacitive coupling between the SQUID’s; $`E_J^\alpha ,E_C^\alpha ,\alpha =\{p,s\}`$ are related respectively to the Josephson and the charging energies of the junctions and $`M`$ is the mutual inductance; the gate voltage $`V_g`$ is tuned to the charge degeneracy point. For generic circuit-parameters the dynamics of the qubit is intricated by the presence of the counter-rotating terms, making the device not reliable in the communication protocols. Our calculation suggests how the circuit parameters can be tuned to reproduce our model (9); for it the dynamics is not altered by the presence of the counter-rotating terms. Using this trick the qubit dynamics can be effectively “protected” at any frequency $`\omega _p`$. The relation between the circuit-parameters and coefficients in the Hamiltonian (9) is: $`\{\omega _p,E_J^s,2\stackrel{~}{L}_p,M,V_C\}\{\lambda (U+V),2Y,U+V,UV,2\sqrt{UV}\}`$ ($`Z=0`$ is set for simplicity), implying that $`V_C`$ should be tuned to $`V_C=\sqrt{M^24\stackrel{~}{L_p}^2}/2`$ to make un-effective the counter-rotating terms.

###### Acknowledgements.

Discussions with G. Falci and A. Osterloh are acknowledged. We are grateful to D. Averin, P. Kulish, E. Paladino, R. Fazio and R. Richardson for helpful comments. While completing this paper, equations similar to (19)-(22) appeared in .