text:
# Spin-down of Relativistic Stars with Phase Transitions and PSR J0537-6910

## 1 Introduction

The properties of the highest-density region of matter in the interior of relativistic stars have been studied to date by several different theoretical models, which yield very different macroscopic properties of compact stars. Many new interactions have been proposed to occur at the highest densities attained, ranging from the occurrence of hyperons, to quark deconfinement, to the formation of kaon condensates and H-matter (see Weber, 2001; Heiselberg 2002, for recent reviews). Observational constraints on pulsar properties are currently still too weak to allow to distinguish which (if any) of all the proposed theoretical descriptions of the interior of relativistic stars is the correct one. A solution to the problem of determining the correct equation of state in the central region of a compact star must come from new, more accurate observations of their properties and any new observational method that will help in this direction is more than welcome.

One new observational method, that could be interesting in the case that the equation of state (EOS) features a deconfinement phase transition to quark matter, has been proposed by Glendenning, Pei and Weber (1997), (hereafter GPW). Specifically, when a rapidly rotating pulsar spins down, its central density increases with time. At a certain central density, a mixed quark-hadron phase can appear (or can already be present in the pulsar). For some chosen parameters of the equation of state and of the baryonic mass of the pulsar, a pure quark core appears when the pulsar reaches very high densities. In GPW. a particular spin-down sequence was studied, that in the nonrotating limit reached nearly at the maximum mass allowed by the chosen EOS. For this sequence, it was noticed that, when the pure quark core appears, the star undergoes an brief era of spin-up (brief, compared to the pulsar’s lifetime) and its braking index shows an anomalous behaviour (it becomes singular at a certain rotation rate). This change in the braking index has been proposed in GPW to be potentially observable. The anomalous behaviour of the spin evolution is traced back to a “backbending” of the moment of inertia (as a function of angular velocity).

Other authors (Heiselberg & Hjorth-Jensen 1998, Chubarian et al. 2000) have further investigated the observational consequences in the models proposed in GPW and in Glendenning and Weber (2000). In all these studies the relativistic slow-rotation approximation is used.111After this work was completed, we learned of a new study using rapidly rotating relativistic stars, with a numerical code related to ours, (Cheng, Yuan & Zhang, 2002), in which the backbending behaviour of the moment of inertia is shown to be sensitive to the assumed properties of the crust. In the present work, we re-investigate the deconfinement phase-transition of spinning-down pulsars, using fully relativistic, rapidly rotating models and find that, when the quark-core appears (in a sequence of models similar to the one used in GPW), the backbending in the moment of inertia is generically absent for normal pulsars, but very strong for supramassive millisecond pulsars. Furthermore, our numerical results suggest that in the limiting region between normal and suprammasive pulsars, the backbending behaviour is sensitive to truncation errors of the applied numerical scheme. The moment of inertia is a very sensitive function of angular velocity and must be computed with high accuracy. Truncation errors, such as a limited number of points in the EOS table, or the use of the slow-rotation approximation can be large enough to systematically affect the computation of the moment of inertia.

The case of supramassive millisecond pulsars is very interesting, as the spin-up era lasts from the onset of the quark-core appearence to their collapse to black holes, i.e. for their remaining lifetime. For supramassive millisecond pulsars created via accretion-induced spin-up, the spin-up era could last for the most part of their lifetime, if a quark-hadron phase transition occurs. Thus, the first time-derivative of the rotational period can be used as a tool for observing quark-hadron phase transitions and one does not have to rely on the measurement of the braking index. We emphasize that the extended spin-up era in supramassive millisecond pulsars is not restricted for the case of a quark-hadron phase transition, but it can also occur for other types of phase transitions, such as a phase transition to a pion condensate. An example is EOS M (Pandharipande and Smith, 1975) for which an extended spin-up with angular momentum loss (for supramassive, constant baryonic mass sequences) has been found in Cook, Shapiro and Teukolsky (1994). On the other hand, the absence of spin-up in the currently known millisecond pulsar population favours the absence of large phase transitions or implies that the EOS is such that millisecond pulsars created via accretion-induced spin-up are probably not supramassive.

In previous studies of spinning-down pulsars that were limited to the slow-rotation approximation (including up to $`O(\mathrm{\Omega }^2)`$ rotational effects) expressions for the evolution of angular velocity and braking index were derived (see Glendenning, 1997). However, we find that these expressions are incomplete, as changes in the gravitational mass of the star (other than changes in the kinetic energy) were not taken into account in the assumed spin-down law. Here, we adopt a general-relativistic version of the spin-down law, which dictates that the energy lost in the form of e.g. magnetic dipole radiation or gravitational waves is not lost to the expense of the kinetic energy only, but to the expense of the total mass-energy of the star. Hence, making use of a relation due to Bardeen (1970), which expresses the first law of thermodynamics along sequences of constant baryonic mass, we derive expressions for the angular velocity and braking index evolution, which are valid for any rotation rate. In addition, our adopted spin-down law is also valid during a spin-up era, while it’s Newtonian, slow-rotation limiting expression (which is usually adopted in the literature) is not valid during spin-up.

As an application of our numerical investigations, we compute a possible value for the initial period of pulsar J0537-6910, assuming that the high-density EOS exhibits a quark-deconfinement phase transition and that the star spins down only due to magnetic dipole radiation. Without the assumption of a phase transition, previous studies suggested that the initial spin of the pulsar would be extremely small (for a braking index of $`n=3`$) (Marshall et al. 1998). Here we show that the presence of the phase transition can yield a more reasonable initial period for this pulsar, without having to assume a braking index smaller than $`n=3`$. In our example, we obtain a significantly larger initial period of 6ms, which is in better agreement with theoretical expectations. Furthermore, this pulsar has displayed a large number of giant glitches (Gotthelf et al. 2001). We propose that these could be due to the appearance of a pure quark core in the center of the star, which results in an increased rate of contraction of the star during its spin-down. In contrast, rapidly rotating pulsars do not exhibit such giant glitches, when their central density is less than the critical density for the appearance of the quark core. The fact that pulsar J0537-6910 displays such glitches more frequently than other Crab-like222In this paper, “Crab-like” only refers to the young age of a pulsar and to its observational association with a supernova remnant. pulsars may be related to the fact that it is the most rapidly spinning of all Crab-like pulsars and thus a pure quark core could have appeared more recently in its lifetime.

The outline of the paper is as follows: In Sect. 2, we describe the chosen equation of state and its numerical refinement, while in Sec. 3 we obtain the equilibrium properties of a specific sequence of spinning-down compact stars. The behaviour of the moment of inertia and gravitational mass along the sequence is described in Sec. 4 and thermodynamic consistency is used as an independent test in Sec. 5. In Sec. 6 expressions for the evolution of the angular velocity and the braking index are derived, and in Sec. 7 the initial period of PSR J0537-6910 is estimated. Finally, in Sec. 8, we comment on the expected event rate of the considered signal. We conclude with a discussion of our results in Sec. 9.

## 2 Equation of State and Numerical Refinement

The EOS used in GPW features a mixed quark-hadron phase, in between pure hadron and pure quark phases. The mixed phase is obtained by applying the Gibbs criteria for two-component systems, abandoning local (but not global) charge-neutrality (Glendenning 1992). This yields a second-order phase transition between each region. However, it has been suggested (see Heiselberg & Hjorth-Jensen, 1999) that the mixed phase can exist only if the interface tension between quarks and hadrons is not too large. But, on the other hand, the interface tension may not be a free parameter, but it could have the freedom to adjust (via a minimization prinicple) to the nature of the two phases (Glendenning 2002), see Christiansen, Glendenning & Schaffner-Bielich (2000) for a relevant example.

The specific EOS we use is the quark deconfinement EOS given in Table 9.2 of Glendenning (1997) (same as in GPW) in which, however, we replace the pure quark phase by an analytic expression (see discussion below). The pure confined hadronic phase exists below baryon number density $`0.246\mathrm{fm}^3`$ while pure deconfined quark matter exists above $`0.862\mathrm{fm}^3`$ (1fm=$`10^{13}`$cm). In between these two phases, a mixed quark-hadron phase exists (see Fig. 1, which displays the pressure $`P`$ vs. energy density $`ϵ`$ relation for this EOS).

When computing nonrotating models with the original tabulated EOS, using the 2-D code by Stergioulas & Friedman (1995), we noticed that for central energy densities near the central energy density where the pure quark phase appears, the gravitational mass vs. central energy density plot shows an unusual behaviour (see Fig. 2). There is a local maximum followed by a local minimum and the mass increases again, before it reaches the maximum mass allowed by the general-relativistic radial instability (see Shapiro & Teukolsky, 1983). Such a behaviour is not present, for exactly the same tabulated EOS, when using a simple 1-D code that integrates the well-known radial TOV equations of structure for non-rotating models. Thus, the oscillatory behaviour of the results obtained with the 2-D code for nonrotating models is a numerical effect. In the 2-D code, the angular behaviour of all variables is expanded as a sum of Legendre polynomials $`P_l(\mu )`$ of argument $`\mu =\mathrm{cos}\theta `$. At discontinuities, such as the surface of the star (where e.g. the velocity profile is discontinuous), the accuracy of the code is reduced due to the Gibbs phenomenon (see Nozawa et al. 1998). Similarly, when the pure quark phase appears at the center of the star, discontinuities appear in the numerical solution, as the EOS features an abrupt change of slope (see Fig. 1). In addition, as the EOS table only contains a limited number of points, the interpolation methods (such as 4-point Lagrange interpolation or Hermite polynomial interpolation) that are used in our code, do not avoid an oscillatory behaviour in the interpolated quantities, near the discontinuity in the slope of the EOS. It appears that, in the 2-D code, such an interpolation error enhances the Gibbs phenomenon at the center of the star, when the pure quark core appears, resulting in the oscillatory behaviour in Fig. 2.

To solve this numerical problem, we replace the part of the EOS that represents the pure quark phase by an analytic expression. In fact, we find that this part can be represented rather accurately by the following linear relation:

$`P[\mathrm{dynes}/\mathrm{cm}^2]`$ $`=`$ $`1.85240\times 10^{35}`$
$`+2.885209\times 10^{20}(ϵ/c^21.69100\times 10^{15})`$
where $`P`$ is pressure, $`ϵ`$ is energy density and $`c`$ is the speed of light (in CGS units). In this way, the tabulated EOS is used only in the region below the pure quark phase and no interpolation is done through the phase transition, thus eliminating the problem of the oscillatory behaviour. For nonrotating models, this refined tabulated/analytic EOS now reproduces the same mass vs. central energy density relation with the 2-D code, as the original tabulated EOS with the 1-D code (see Fig. 2). Notice that for this comparison we avoid interpolating through the region where the pure quark appears, by choosing points either before or after the phase transition only. This increases the accuracy of computed models (as compared to interpolating through all points in the table) even with the 1-D code.

Thus, we have verified that replacing the pure quark region of the original tabulated EOS by the analytic expression (2) eliminates numerical problems in the 2-D code, while still representing the same EOS (with sufficiently high accuracy).

## 3 Spin-down Sequence

We study several evolutionary sequences of pulsars, that differ from each other in the assumed constant baryonic mass. For a baryonic mass of $`M_0=1.555M_{}`$ we obtain a sequence that terminates at the maximum mass nonrotating model allowed by the EOS. This sequence is the limiting sequence between normal pulsars and supramassive pulsars and is very similar to the sequence studied in GPW. We assume that a compact star with a baryonic mass of $`M_0`$ is rapidly rotating. The compact star subsequently spins down due to e.g. magnetic braking and/or gravitational radiation emission, and the sequence ends either at the nonrotating limit (for normal pulsars), or the pulsar encounters the axisymmetric instability limit to black hole collapse (see e.g. Cook, Shapiro and Teukolsky, 1984, for examples of such sequences for various realistic EOSs). The sequences are constructed using a highly accurate numerical code (rns, Stergioulas & Friedman, 1995) that computes exact numerical solutions of rapidly rotating stars in general relativity, without any approximation (apart from the truncation errors of the numerical scheme). Solutions are obtained on a 2-D numerical grid, with typical grid-size of 800x400 or 1200x400 (radial $`\times `$angular) points. Each model is defined by its central energy density and polar to equatorial axis ratio. Apart from the requirement of high accuracy for each individual model, essential to the construction of a constant-baryon number sequence, is the accuracy with which a model of a given baryonic mass is located. To this end, we construct various models by keeping the central energy density fixed, varying only the axes ratio and use a robust root-finding method to locate a model of given baryonic mass within a specified accuracy. The root-finding method of our choice is Ridders’ method (Ridders, 1979), since the procedure is guaranteed to stay within the initial brackets (avoiding unphysical parameters) and since it has a relatively high order of convergence ($`\sqrt{2}`$). In practice, about 10-12 individual models are required for the root-finding method to converge to a model of given baryonic mass, with a relative accuracy of $`10^9`$ (the accuracy in computing the baryonic mass for an individual model is of the order of $`10^5`$, but this does not limit the root-finding method to converge to a specified value with higher relative accuracy).

## 4 Moment of Inertia and Gravitational Mass

We compute several physical parameters of the equilibrium models, the moment of inertia $`I`$, gravitational mass $`M`$, angular momentum $`J`$ and angular velocity $`\mathrm{\Omega }`$ being the most relevant for our present study. We find that the behaviour of the moment of inertia depends on the baryonic mass of the sequence. Along a normal sequence with $`M_0=1.551M_{}`$, which is very close to the maximum mass nonrotating limit, but still below it, the change in the moment of inertia, as the star spins down, is shown in Fig. 3, which plots $`I`$ as a function of $`\mathrm{\Omega }`$ (solid line). Using our refined tabulated/analytic EOS, described above, the behaviour of the moment of inertia is seen to be regular, i.e. $`I`$ is a single-valued function of $`\mathrm{\Omega }`$. We notice that an apparent backbending appears with our 2-D code for this sequence, when the EOS in the original tabulated form (and interpolation through all points in the table) is used (dashed line in Fig. 3). This is an indication that the backbending of the moment of inertia is sensitive to the available points in the EOS table and/or to truncation of rotational effects on the structure of the relativistic star.

Another important equilibrium property that we monitor along the equilibrium sequence is the gravitational mass. In Fig. 4, the gravitational mass $`M`$ is plotted as a function of central energy density $`ϵ_c`$. The solid line corresponds to the refined tabulated/analytic EOS and shows a monotonous decrease of the gravitational mass during spin-down. This normal behaviour corresponds to the monotonous behaviour of the moment of inertia in Fig. 3. In contrast, when one uses the original tabulated form of the EOS, the gravitational mass appears to oscillate as a function of angular velocity. This oscillatory behaviour is not correct, as it is the result of numerical truncation errors, due to the interpolation through the phase transition in the original tabulated EOS.

The oscillatory behaviour obtained with the original tabulated EOS deserves some further attention. As we will show, if it were the correct behaviour, it would correspond to a part of a sequence being unstable to quasi-radial perturbations: In this sequence, before the quark core appears, the equilibrium models are stable, while the star loses kinetic energy $`T`$ and angular momentum $`J`$, and the radius decreases only due to the decrease in rotational flattening of the star during spin down. When the pure quark core appears at the center of the star, the radius starts decreasing significantly, as the EOS becomes softer. Looking at models with larger central densities along the equilibrium sequence, the angular velocity $`\mathrm{\Omega }`$, angular momentum $`J`$, kinetic energy $`T`$, absolute value of gravitational binding energy $`|W|`$ and gravitational mass $`M`$, all increase. Thus, a star cannot continue evolving along this uniformly rotating sequence to the expense of the star’s available mass-energy. In Fig. 4 (dashed lined) a turning point appears and the gravitational mass increases with increasing central density. This is a typical case of the onset of the instability to the fundamental axisymmetric, quasi-radial mode of oscillation, in rotating relativistic stars, as shown by Friedman, Ipser & Sorkin (1988) (hereafter FIS). The turning-point criterion for the onset of instability was first proved in FIS for a constant-angular-momentum sequence. Along such a sequence, an extremum in the $`M`$ vs. $`ϵ_c`$ curve separates the regions of stable and unstable stars, i.e.

$$\left(\frac{M}{ϵ_c}\right)_J=0\mathrm{onset}\mathrm{of}\mathrm{instability},$$
(2)
for a sequence of uniformly rotating stars obeying a one-parameter $`P=P(ϵ)`$ equation of state. Such a criterion was then proved by Cook, Shapiro & Teukolsky (1992) to also exist for sequences of constant baryonic mass $`M_0`$:

$$\left(\frac{M}{ϵ_c}\right)_{M_0}=0\mathrm{onset}\mathrm{of}\mathrm{instability}.$$
(3)
The instability to the fundamental quasi-radial mode of oscillation is first a secular one, i.e. it proceeds on a viscous timescale. It has been shown, however, that a dynamical instability sets in near the secular instability point (see Shibata, Baumgarte & Shapiro 2000). Thus, if the star would evolve along the dashed curve in Fig. 4, reaching the secular instability point, it would first evolve away from the uniformly-rotating equilibrium sequence, on a viscous timescale, driven by magnetic braking torques and/or gravitational radiation reaction and could then also become dynamically unstable, in which case it would suddenly contract. Eventually, the star would bounce at higher densities, when reaching again a stable branch in the equilibrium sequence (after the second turning point in Fig. 4).

We find that, when computing sequences of constant baryonic mass, the behaviour of the gravitational mass can be used as a guide to distinguish whether the sequence has been computed accurately: a behaviour such as the apparent instability described above, indicates a strong influence of numerical errors in the computation.

Looking at other constant baryonic mass sequences, with baryonic mass between $`M_0=1.551M_{}`$ and $`M_0=1.555M_{}`$ (see Fig. 5), we see that the moment of inertia gradually approaches a vertical behaviour ($`I/\mathrm{\Omega }\mathrm{}`$), as the limiting sequence that divides normal from supramassive stars is approached. A very small backbending (spin-up) behaviour is first observed for the limiting sequence of $`M_0=1.555M_{}`$ (and clearly not for e.g. $`M_0=1.554M_{}`$, see Fig. 6). For larger baryonic masses (supramassive sequences) the backbending behaviour becomes stronger and for e.g. a baryonic mass of $`M_0=1.6M_{}`$ it becomes the dominant behaviour, i.e. once the quark core appears, the star only spins-up with angular momentum loss, for the rest of its lifetime.

To summarize, for the particular EOS that we have studied in this paper, normal pulsars do not feature a spin-up era during phase transition, while supramassive pulsars only spin up after the phase transition. The limiting sequence that divides normal from supramassive pulsars shows a very small era of spin-up and appears to be the dividing line between stars that can spin up and stars that only spin down.

## 5 Thermodynamic Consistency

In the previous sections we demonstrated that the refined tabulated/analytic EOS produces a constant baryonic mass sequence with equilibrium properties that are in sharp contrast to the equilibrium properties computed on the basis of the original tabulated EOS. In this section we will use an independent check to verify that the equilibrium properties obtained with the refined EOS are indeed the physically acceptable ones, while the equilibrium sequence obtained with the tabulated EOS suffers from numerical errors.

The independent check of the accuracy of the computed equilibrium properties is provided by the use of a relation that was first proved by Ostriker and Gunn (1969) in Newtonian theory and then derived in general relativity by Bardeen (1970).

Along a sequence of uniformly rotating models of constant baryonic mass, changes in $`M`$ and $`J`$ are related by

$$c^2dM=\mathrm{\Omega }dJ,$$
(4)
which can be regarded as an expression of the first law of thermodynamics for such sequences.

We have evaluated the above relation numerically along the obtained sequences of equilibrium models (using sufficiently many individual models to ensure an accurate evaluation of the numerical derivative). We define an error indicator, $`\lambda `$, as

$$\lambda =\left|1\frac{c^2dM}{\mathrm{\Omega }dJ}\right|.$$
(5)
The evaluation of $`\lambda `$ along the evolutionary sequence of $`M_0=1.551M_{}`$ (Fig. 7), shows that the above relation is satisfied with an accuracy of $`1`$% for models with central energy densities smaller than the energy density at which the pure quark-core appears. After the appearance of the pure quark core, numerical errors result in huge values of $`\lambda `$, when using the original tabulated EOS. With the use of the refined tabulated/analytic EOS, these errors are reduced significantly to an acceptable level and diminish as the resolution is increased.

The above check confirms that the refined EOS indeed produces physically acceptable results, while the original tabulated EOS produces results that are dominated by numerical error in the region in question.

## 6 Behaviour of the Observable Braking Index

If a compact star is detected as a pulsar and one can measure the first two time-derivatives $`\dot{\mathrm{\Omega }}`$, $`\ddot{\mathrm{\Omega }}`$ of its angular velocity, then one can define an observational braking index

$$n=n(\mathrm{\Omega })\frac{\mathrm{\Omega }\ddot{\mathrm{\Omega }}}{\dot{\mathrm{\Omega }}^2}.$$
(6)
In the Newtonian, slow-rotation limit, the spin-down of a pulsar is usually modeled in the form of a power law (see Shapiro & Teukolsky. 1983), namely assuming that only kinetic energy is lost and that the rate of loss is proportional to some power of the angular velocity of the star, i.e.

$$\frac{dT}{dt}=\kappa \mathrm{\Omega }^\alpha ,$$
(7)
where $`\kappa <0`$ and $`\alpha >0`$ are usually assumed to be real constants. With these assumptions, the braking index is equal to $`n=\alpha 1`$. For example, it is assumed that, for magnetic braking, $`\dot{T}\mathrm{\Omega }^4`$, which yields an expected braking index $`n=3`$, while, for gravitational wave emission, $`\dot{T}\mathrm{\Omega }^6`$, which yields $`n=5`$. For slowly rotating pulsars, the above considerations are, in most cases, appropriate. However, for rapidly rotating pulsars one has to take into account the rotational flattening of the star. To a first approximation, this can be done considering rotational effects up to order $`O(\mathrm{\Omega }^2)`$. Glendenning (1997) gives the rotationally corrected $`\dot{\mathrm{\Omega }}`$ and braking index $`n(\mathrm{\Omega })`$ for a sequence of uniformly rotating stars, assuming the spin-down law of Equation (7). As we will show here, the expressions given in Glendenning (1997) are incomplete, in the sense that the derivation is not fully consistent to $`O(\mathrm{\Omega }^2)`$, but misses additional contributions of the same order.

The energy lost in the form of electromagnetic or gravitation radiation is not only to the expense of the star’s kinetic energy (which would be the case only in the $`O(\mathrm{\Omega })`$ slow-rotation approximation) but to the expense of the star’s total mass-energy (gravitational mass). In the Newtonian limit, this principle has been applied to the spin-down of rapidly rotating neutron stars by Finn & Shapiro (1990) and by Spyrou & Stergioulas (2001). In relativistic stars, the total mass-energy is

$$Mc^2=M_0c^2+U+T+W,$$
(8)
(see Friedman, Ipser & Parker 1986) where, additionally, $`U`$ is the internal energy. We thus adopt the following natural generalization of Equation (7) in general relativity:

$$c^2\frac{dM}{dt}=\kappa \mathrm{\Omega }^\alpha .$$
(9)
The kinetic energy in relativity is defined as $`T=(1/2)J\mathrm{\Omega }`$, so that

$$\frac{dT}{dt}=\frac{1}{2}\left(c^2\frac{dM}{dt}+\frac{d\mathrm{\Omega }}{dt}J\right).$$
(10)
Expressing this in terms of the moment of inertia, which is defined as in Friedman et al. (1986), $`I=J/\mathrm{\Omega }`$, and with the help of Equation (4), we find

$$\frac{dT}{dt}=\left(c^2\frac{dM}{dt}\frac{1}{2}\mathrm{\Omega }^2\frac{dI}{dt}J\right).$$
(11)
Equation (11) shows that, in a rapidly rotating star, under the above assumptions, the time derivatives of the kinetic energy and gravitational mass differ by an $`O(\mathrm{\Omega }^2)`$ term that is proportional to the time-derivative of the moment of inertia. Notice that this expression is exactly valid for any rotation rate (even for stars rotating at the mass-shedding limit) and not only in an $`O(\mathrm{\Omega }^2)`$ approximation. The relations derived in Glendenning (1997) for the braking index and time-derivative of the angular velocity in the $`O(\mathrm{\Omega }^2)`$ approximation, are missing the contribution from the $`O(\mathrm{\Omega }^2)`$ term on the r.h.s. of Equation (11).

With the above assumptions, it is easy to show that Equations (5.18) and (5.20) in Glendenning (1997) are replaced by

$$\dot{\mathrm{\Omega }}=\frac{\kappa }{I}\mathrm{\Omega }^{\alpha 1}\left(1+\frac{\mathrm{\Omega }I^{}}{I}\right)^1,$$
(12)
where $`I^{}=dI/d\mathrm{\Omega }`$ and

$$n(\mathrm{\Omega })=\alpha 1\frac{2I^{}\mathrm{\Omega }+I^{\prime \prime }\mathrm{\Omega }^2}{I+I^{}\mathrm{\Omega }}.$$
(13)
We have verified that Equation (13) is the correct one, by numerically computing the braking index as obtained above and comparing it to the direct computation of the braking index, using the definition Equation (6). The two numerical results agree well even close to the mass-shedding limit. In contrast, the expression for $`n(\mathrm{\Omega })`$ in Glendenning (1997) shows a large disagreement with the result obtained from the definition of $`n(\mathrm{\Omega })`$ for rapidly rotating stars. This difference between the two expressions is significant also in the case of supramassive sequences.

In Fig. 8, the observational braking index $`n(\mathrm{\Omega })`$ is shown as a function of angular velocity (computed along the normal spin-down sequence with $`M0=1.551M_{}`$), assuming that the spin-down is only due to magnetic braking. The computed $`n(\mathrm{\Omega })`$ is roughly $`n3`$ before the pure quark core appearance. During the phase transition, the braking index rises to roughly $`n10`$ and then decreases to a value of roughly $`n2`$, before returning to the canonical value of $`n3`$. This behaviour is qualitatively and quantitatively different from the singular behaviour in the braking index obtained with our code when one uses the original tabulated EOS (see inset in Fig. 8).

## 7 Initial Period of PSR J0537-6910

PSR J0537-6910 was discovered in SNR N157B by Marshall et al. (1998). It has a spin period of only 16ms, the fastest of all Crab-like pulsars, and a period derivative of $`\dot{P}=5.13\times 10^{14}\mathrm{ss}^1`$. Both the characteristic age $`\tau =P/(2\dot{P})`$ and the age estimate for SNR N157B suggest that the pulsar is 5,000 years old. An upper limit on the age of the pulsar based on observations of $`H\alpha `$-emitting gas in N157B (Chu et al. 1992), is 20,000 years. Marshall et al. show that, if the braking index is $`n=3`$ and the pulsar age is 5,000 years, then the pulsar must have been born with an extremely small (even sub-millisecond) period. Such small initial periods in pulsars associated with supernova remnants have not been detected to date. There are theoretical arguments, based e.g. on the possible occurrence of gravitational-wave induced instabilities in rapidly rotating relativistic stars (see e.g. Stergioulas 1998, Andersson & Kokkotas 2001 for reviews) that suggest a pulsar initial period (shortly after the compact star is born) of at least several ms. How can this theoretical expectation be compatible with the observed properties of PSR J0537-6910 ? One possible solution would be that the actual braking index for this pulsar is closer to $`n=2`$ (for the upper limit on age of 20,000 years, a very low braking index of $`n1`$ is required).

Here we show that an initial period of several ms (6 ms in our example) can be obtained for a braking index of $`n=3`$, if the appearance of a pure quark core is considered. We use the same EOS and the $`M0=1.551M_{}`$ constant baryonic mass sequence mentioned in previous sections and numerically integrate Equation (12). We find that the age of the pulsar can be expressed in the following integral form

$$t_{\mathrm{age}}=\frac{1}{\kappa }_{\mathrm{\Omega }_i}^\mathrm{\Omega }\frac{I+\mathrm{\Omega }dI/d\mathrm{\Omega }}{\mathrm{\Omega }^{\alpha 1}}𝑑\mathrm{\Omega },$$
(14)
where $`\mathrm{\Omega }_i`$ is the initial angular velocity, or, equivalently,

$$t_{\mathrm{age}}=\frac{1}{\kappa }_{J_i}^J\mathrm{\Omega }^{1\alpha }𝑑J,$$
(15)
where $`J_i`$ is the initial angular momentum. For numerical integration, it is more convenient to introduce a dimensionless central energy density $`\stackrel{~}{ϵ}_c=(ϵ_c/c^2)/(10^{15}\mathrm{g}/\mathrm{cm}^3)`$ as the integration parameter, and evaluate

$$t_{\mathrm{age}}=\frac{1}{\kappa }_{\stackrel{~}{ϵ}_c(i)}^{\stackrel{~}{ϵ}_c}\mathrm{\Omega }^{1\alpha }\frac{dJ}{d\stackrel{~}{ϵ}_c}𝑑\stackrel{~}{ϵ}_c.$$
(16)
The result of the numerical integration is shown in Fig. 9, which displays the age of the pulsar as a function of the central energy density. For a current spin of 16ms and a pulsar age of 5,000 years, the initial spin is obtained to be 6ms. Thus, even an $`n=3`$ braking index can yield an initial spin that is in agreement with theoretical expectations, provided the occurrence of a pure quark core is taken into account. In the above example, a particular choice of EOS and baryonic mass was made. Different choices would lead to different initial spin period estimates. However, the qualitative effect of the presence of a pure quark core is to increase (with respect to an EOS without a phase transition) the computed initial spin period.

## 8 Event Rate

An important parameter for the successful detection of a phase-transition signal in the braking index of spinning-down pulsars is the event rate, which is determined by several important factors. Let’s assume that pulsars are born with an initial period of only a few milliseconds. Since normal pulsars have a magnetic field strength of $`10^{12}G`$, they would quickly spin down to much larger periods due to magnetic dipole radiation. Thus, if the baryonic mass of a newly-born pulsar is such that the phase transition appears at a short rotational period, then the phase transition signal will not last for an extended period of time to have a realistic chance of being detected. Notice that no young, Crab-like pulsar with period shorter than $`16`$ms has been observed, which is indicative of the very small event rate that should be expected for such an observation.

On the other hand, if pulsars are born with rotational periods larger than, say, 6ms, then a phase transition can occur during their lifetime only if their baryonic mass falls within an extremely small range of values. To illustrate this point, we plot in Figure 10 the constant baryonic mass sequence considered in Sec. 3 and the corresponding sequence of nonrotating models. The vertical dashed line specifies the central density at which the pure quark core appears. It crosses the constant-baryonic-mass sequence at a model with rotational period of 6ms and gravitational mass of $`1.419M_{}`$. Notice that the gravitational mass of a nonrotating star with central energy density equal to that of the quark core appearance ($`1.691\times 10^{15}\mathrm{g}/\mathrm{cm}^3`$) is $`1.414M_{}`$. Thus, a pure quark core can appear in a spinning down pulsar with rotational period less than 6ms only if its mass falls within $`\mathrm{\Delta }M0.0025M_{}`$ of the value of $`1.4165M_{}`$ (the the sequence considered in this example). The probability for this to happen is obviously negligible.

The above arguments show that the only observed population remaining, in which one could hope for a significant event rate for the considered signal, are old millisecond pulsars, spun-up by accretion, with low magnetic field and spin-down rate. The actual event rate will then depend on the mass-distribution for this population and on the precise value of the central energy density at which the quark core appears. None of these parameters are currently known with sufficient accuracy to compute a reliable value for the event rate. In fact, the observational determination of the braking index is itself a difficult task and no braking index has been measured for any of the known millisecond pulsars, to date. However, should a millisecond pulsar become supramassive, due to mass-accretion, then there is a realistic chance for detecting the presence of a quark core in observations of the first time-derivative, $`\dot{P}`$, of the rotational period. Once the quark core appears and if the pulsar is sufficiently massive, then $`\dot{P}<0`$ for the rest of the pulsar’s lifetime.

## 9 Discussion

Observational constraints for the properties of compact stars are currently still too weak to determine the correct equation of state of high-density matter. The presence of a quark/hadron deconfinement phase transition has been suggested by Glendenning, Pei and Weber (1997) (GPW) to be observable in the braking index of pulsars. In this method, the observational signal comes from a spin-up era, originating from a backbending behaviour of the moment of inertia during the pure quark core appearance.

We have investigated the deconfinement phase-transition of spinning-down pulsars, using fully relativistic, rapidly rotating models. We find that, for normal pulsars, the star does not go through a spin-up era and the pure quark core appears without a backbending behaviour in the moment of inertia. The backbending behaviour first occurs for the limiting sequence of constant baryon mass that terminates at the maximum mass nonrotating model allowed by the EOS. In contrast, for sequences with larger mass (supramassive sequences) the backbending behaviour becomes dominant and pulsars only spin up, after a quark core appears.

From our numerical results, it is evident that the moment of inertia along a spin-down sequence is very sensitive to the appearance of phase transitions and must be computed with great care. In particular, we emphasize that low resolution and the restriction to a limited number of points in the EOS table (with interpolation through all points) can easily spoil the correct computation of the behaviour of the moment of inertia, at least in the present 2-dimensional code used in this paper.

The expected event rate of a phase-transition signal in the pulsar braking index is negligibly small for the population of normal pulsars. For the population of old millisecond pulsars, a reliable estimate for the event rate cannot be obtained, as it involves several uncertain factors and, in addition, the observational determination of braking indices in millisecond pulsars has not become technically possible, yet. However, if millisecond pulsars become supramassive, then the presence of a quark core could be detected by observing a negative first time-derivative of the rotational period. On the other hand, the absence of any such observation in the current population of known millisecond pulsars indicates that either large phase transitions do not occur, or that the EOS is sufficiently stiff to prevent accreting relativistic stars from becoming supramassive. A detailed study of phase transitions in supramassive models of relativistic stars, constructed with various EOSs, could provide interesting constraints on the equation of state at very high densities.

In order to describe the spin down of a pulsar due to magnetic dipole or gravitational radiation, we adopt a general-relativistic version of the usual power-law assumed in the slow-rotation limit. The energy loss is assumed to be to the expense of the total mass-energy of the star and not simply to the expense of the star’s kinetic energy. In this way, we derive new expressions for the evolution of the angular velocity and braking index during spin down. As an application, assuming magnetic dipole braking only, the initial spin period of PSR J0537-6910 is obtained to be 6ms (for the EOS and spin-down sequence considered in this paper). The presence of the deconfinement phase transition can thus give a more reasonable initial period for this pulsar, compared to previously obtained estimates.

PSR J0537-6910 exhibits the highest rate of large glitch events (6 in 3 years) (see Gotthelf et al. 2001) of all known Crab-like pulsars. At the same time, it is the most rapidly rotating Crab-like pulsar. If a pure quark core appears at small spin periods, then the quark-core fraction in such pulsars continuously changes during spin down, most probably contributing to the occurrence of large glitches. The fact that PSR J0537-6910 exhibits the highest glitch rate could be related to the fact that the phase transition may have happened more recently in its lifetime (according to the spin-down sequence used here) and its spin-down behaviour is still influenced by the gradual change of the quark core fraction during spin down. This suggestion has to be examined in more detail, as it could be a generic feature of phase transitions in the high-density EOS.

## Acknowledgements

We thank Nils Andersson, Pawel Haensel and John L. Friedman for useful discussions and Emanuele Berti for a careful reading of the manuscript. We are grateful to Leszek Zdunik and Pawel Haensel for comparing numerical results of several evolutionary sequences, before publication. Finally, we are indebted to the referee, Prof. N. Glendenning, for critical comments that substantially improved the final version of this paper. This work has been supported by the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137), KBN-5P03D01721 and the Greek GSRT Grant EPAN-M.43/2013555.


format: text

text:
# Dilution Effects on Two-Dimensional Heisenberg Antiferromagnets with Non-Magnetic Spin-Gapped Ground State

## 1 Introduction

Dilution effects on low-dimensional antiferromagnets with a spin-gapped ground state exhibit various interesting phenomena inherent to low-dimensional quantum spin systems. In the first inorganic spin-Peierls compound CuGeO3, the dimerized state is realized at low temperature and the system has the spin-gapped ground state.  When a small amount of the non-magnetic impurity Zn or Mg is substituted for Cu, the spin-gapped ground state of the pure system changes to the antiferromagnetic (AF) long-range ordered state,  thereby the lattice dimerization is preserved. The impurity-induced AF long-range order (LRO) has been observed also for other spin-gapped materials such as the two-leg ladder SrCu2O3  and the Haldane compound PbNi2V2O8 The non-magnetic-impurity-induced AF-LRO has been numerically investigated by simulating the site-diluted AF Heisenberg model with the spin-gapped ground state.  In the site-diluted system, the AF-LRO is induced even at 1 % of site dilution for the AF Heisenberg model of the coupled alternating chains,  suggesting the critical concentration of non-magnetic impurities for AF-LRO to come up is actually 0 (see also Fig. 1 below). When spins are removed from sites, the effective spins are induced, whose magnitudes are peaked at the neighboring site connected to the removed site by a strong bond. They have the same extent as the correlation length of the pure system. Furthermore, there exists an effective coupling $`\stackrel{~}{J}_{mn}(1)^{|r_mr_n+1|}\mathrm{exp}[l/\sqrt{\xi _\mathrm{p}^x\xi _\mathrm{p}^y}]`$ between the effective spins,  which preserves the staggerdness with respect to the original lattice and decays exponentially as the distance between the effective spins, $`l=|r_mr_n|`$, increases. Here $`\xi _\mathrm{p}^{x,y}`$ are the anisotropic correlation lengths of the pure system in the two directions on a square lattice. Thus the AF-LRO is induced at an infinitesimal concentration of dilution for the two-dimensional (2D) site-diluted system.

In the present work we have investigated ground-state phase transitions in the bond-diluted 2D spin-1/2 AF Heisenberg model of the coupled alternating chains on a square lattice by means of the quantum Monte Carlo (QMC) simulation with the continuous-imaginary-time loop algorithm. In the bond-diluted case, in contrast to the site-diluted one, a finite critical concentration of dilution, $`x_\mathrm{c}`$, exists, above which the system has the AF-LRO. In the extremely-low-concentration region, the system has a spin gap independent of the concentration of the diluted bond, $`x`$. The spin gap $`\mathrm{\Delta }`$ is proportional to the squared interchain coupling $`J^2`$ significantly smaller than the original spin gap in the pure system, $`\mathrm{\Delta }_\mathrm{p}`$. This spin gap $`\mathrm{\Delta }`$ is caused by singlet pairs of two edge spins on the diluted bond by the AF coupling $`\stackrel{~}{J}_{\mathrm{af}}`$ through the 2D shortest paths. Furthermore, when $`x`$ increases, $`\mathrm{\Delta }`$ starts to decrease, and tends to vanish at $`x`$ which is smaller than $`x_\mathrm{c}`$. This may suggest the existence of the quantum Griffiths phase in this system.

## 2 Model and method

The bond-diluted AF Heisenberg model of the coupled alternating chains on a square lattice investigated in the present work is described by the Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{i,j}{}}ϵ_{(2i,j)(2i+1,j)}𝐒_{2i,j}𝐒_{2i+1,j}+\alpha {\displaystyle \underset{i,j}{}}ϵ_{(2i+1,j)(2i+2,j)}𝐒_{2i+1,j}𝐒_{2i+2,j}`$ (2.1)
$`+`$ $`J^{}{\displaystyle \underset{i,j}{}}ϵ_{(i,j)(i,j+1)}𝐒_{i,j}𝐒_{i,j+1},`$
where 1 and $`\alpha `$ ($`0\alpha 1`$) are the AF intrachain alternating coupling constants, $`J^{}`$ ($`0`$) the AF interchain coupling constant, and $`𝐒_{i,j}`$ the quantum spin operator with $`S=1/2`$ at site ($`i,j`$). Randomly quenched bond occupation factors {$`ϵ_{(i,j)(k,l)}`$} independently take either 1 or 0 with probability $`1x`$ and $`x`$, respectively. The pure system described by (2.1) with $`ϵ_{(i,j)(k,l)}=1`$ for all bonds is in either a spin-gapped phase or an AF long-range ordered phase at zero temperature depending on the strengths of $`\alpha `$ and $`J^{}`$ In the $`\alpha =0.5`$ system which we examine in the present work, there exists the spin-gapped phase for $`J^{}<J_\mathrm{c}^{}0.55`$.

The QMC simulations with the continuous-imaginary-time loop algorithm  are carried out on $`L\times L`$ square lattices with the periodic boundary conditions. For each sample with a bond-diluted configuration, $`10^310^4`$ Monte Carlo steps (MCS) are spent for measurement after $`50010^3`$ MCS for thermalization. Sample average for dilution is taken over $`1010^3`$ samples depending on $`L`$, $`x`$, and the temperature $`T`$.

## 3 Results and discussions

The staggered magnetization $`M_\mathrm{s}(x)`$ at $`T=0`$ is evaluated by

$$M_\mathrm{s}^2(x)=\underset{L\mathrm{}}{lim}\underset{T0}{lim}\frac{3}{L^4}\underset{i,j}{}(1)^{|ij|}S_i^zS_j^z,$$
(3.1)
where the bracket $`\mathrm{}`$ denotes both the thermal and random averages. The concentration dependences of $`M_\mathrm{s}`$ in the site-  and bond-diluted systems with $`\alpha =J^{}=0.5`$ are shown in Fig. 1. In the site-diluted case, the value of $`M_\mathrm{s}`$ becomes finite even at 1 % of site dilution. This is caused by the strong correlation between the effective spins due to the effective coupling $`\stackrel{~}{J}_{mn}`$ mentioned in $`\mathrm{\S }`$ 1. In the bond-diluted case, on the other hand, there is the finite critical concentration of $`x_\mathrm{c}0.05`$ even for $`\alpha =J^{}=0.5`$, which is the system near the critical point $`J_\mathrm{c}^{}0.55`$ for $`\alpha =0.5`$ and $`x=0`$. Further dilution would yield the disappearance of $`M_\mathrm{s}`$ near the percolation thresholds.

Next we examine the spin gap $`\mathrm{\Delta }`$ between the ground state and the first excited state to reveal the nature of the disordered phase below $`x_\mathrm{c}`$ in the bond-diluted system. We adopt the second-moment method,  by which the value of $`\mathrm{\Delta }^{(2)}`$ defined as

$$(\mathrm{\Delta }^{(2)})^1=\underset{L\mathrm{}}{lim}\underset{\beta \mathrm{}}{lim}\frac{\beta }{2\pi }\sqrt{\frac{S(\pi ,\pi ,0)}{S(\pi ,\pi ,2\pi /\beta )}1},$$
(3.2)
is generally considered to be a good estimate of $`\mathrm{\Delta }`$. In Eq. (3.2) $`\beta =1/T`$ and $`S(\pi ,\pi ,\omega )`$ is the dynamic structure factor at momentum ($`\pi `$, $`\pi `$) described by $`S(\pi ,\pi ,\omega )=_0^\beta C(\tau )\mathrm{e}^{i\omega \tau }d\tau `$, where $`C(\tau )`$ is the correlation function of the staggered magnetization in the imaginary-time direction:

$$C(\tau )=\frac{1}{L^2\beta }\underset{i,j}{}_0^\beta d\tau ^{}(1)^{|ij|}S_i^z(\tau ^{})S_j^z(\tau ^{}+\tau )$$
(3.3)
with $`S_i^z(\tau )`$ being the $`z`$ component of the spin operator at the imaginary time $`\tau `$. Its $`x`$-dependence for the system with $`\alpha =0.5`$ and $`J^{}=0.3`$ is shown in Fig. 3 below. For a small $`x(<0.01)`$ $`\mathrm{\Delta }^{(2)}`$ rapidly decreases with increasing $`x`$. However as we will argue below this is an artifact of the second-moment method except for exactly at $`x=0`$.

Let us assume that, in the low-concentration limit ($`x0`$), there exists a mid-gap $`\mathrm{\Delta }`$ induced by dilution inside the original spin gap, $`\mathrm{\Delta }_\mathrm{p}`$, of the pure system. The former excitation mode contributes to $`C(\tau )`$ in proportion to $`x`$ so that it may less contribute to $`\mathrm{\Delta }^{(2)}`$ as compared to the back ground mode of $`\mathrm{\Delta }_\mathrm{p}`$. More explicitly, $`C(\tau )`$ in this situation and in the limit $`L\mathrm{}`$ may be written as

$$C(\tau )=(1ax)(\mathrm{sinh}\frac{\beta \mathrm{\Delta }_\mathrm{p}}{2})^1\mathrm{cosh}[(\tau \frac{\beta }{2})\mathrm{\Delta }_\mathrm{p}]+ax(\mathrm{sinh}\frac{\beta \mathrm{\Delta }}{2})^1\mathrm{cosh}[(\tau \frac{\beta }{2})\mathrm{\Delta }],$$
(3.4)
where $`a`$ is a non-trivial constant depending on $`\alpha `$ and $`J^{}`$. Substituting Eq. (3.4) for Eq. (3.2), we obtain the following relation between $`\mathrm{\Delta }^{(2)}`$, $`\mathrm{\Delta }_\mathrm{p}`$ and $`\mathrm{\Delta }`$:

$$\mathrm{\Delta }^{(2)}(x)=\mathrm{\Delta }\sqrt{\frac{1+\frac{1ax}{ax}\frac{\mathrm{\Delta }}{\mathrm{\Delta }_\mathrm{p}}}{1+\frac{1ax}{ax}(\frac{\mathrm{\Delta }}{\mathrm{\Delta }_\mathrm{p}})^3}}.$$
(3.5)
To confirm that the above assumption is in fact the case, we simulate the systems in which only one strong bond is removed, i.e., $`x=1/2L^2`$. In Fig. 2 the $`x`$ dependences of the QMC data of $`\mathrm{\Delta }^{(2)}`$ are shown for $`J^{}=0.3`$, 0.25, 0.2, 0.15 and $`\alpha =0.5`$. They are well fitted to Eq. (3.5), the lines in the figure. This implies that the spin gap $`\mathrm{\Delta }`$, or the mid-gap, is independent of $`x`$ at least in the $`x`$-range shown in Fig. 2. By this fitting the value of $`\mathrm{\Delta }`$ is obtained as shown in the inset of Fig. 2, in particular $`\mathrm{\Delta }=0.138(5)`$ for $`J^{}=0.3`$ as indicated in Fig. 3. The existence of the spin gap proportional to $`J^2`$ is interpreted that two edge spins induced at both ends of a diluted bond reform a singlet pair due to the effective coupling $`\stackrel{~}{J}_{\mathrm{af}}`$ ($`J^2`$) through the 2D shortest paths.  Since two edge spins on the diluted bond are on the different sublattices, the interaction is always AF.

As $`x`$ increases, the distinction between a mid-gap and the gap $`\mathrm{\Delta }_\mathrm{p}`$ becomes obscure. Instead, the value of $`\mathrm{\Delta }^{(2)}`$ gives us an upper bound of the spin gap $`\mathrm{\Delta }`$ of our interest, whose $`x`$-dependence is schematically shown by the dotted line in Fig. 3. Presumably $`\mathrm{\Delta }^{(2)}`$ is a good estimate for $`\mathrm{\Delta }`$ in the range $`x>0.02`$.

In Fig. 3 we also show the $`x`$ dependence of the spatial correlation length $`\xi ^x(x)=lim_L\mathrm{}\xi ^x(L;x)`$ at $`T0`$. The critical concentration $`x_\mathrm{c}`$ is about 0.11, which is estimated by the finite-size-scaling analysis using the QMC data of $`\xi ^x(L;x)`$. Its details will be reported elsewhere. In marked contrast to $`1/\xi ^x`$, $`\mathrm{\Delta }^{(2)}(\mathrm{\Delta }`$ in the range $`x>0.02`$) seems to vanish at a certain value $`x_\mathrm{c}^{}`$ smaller than $`x_\mathrm{c}`$, which strongly suggests the existence of the quantum Griffiths phase at $`x_\mathrm{c}^{}<x<x_\mathrm{c}`$. Thus the bond-diluted system has the more peculiar properties than the site-diluted system.

The $`x`$ dependence of the Néel temperature in a bond-random system CuGe1-xSixO3 has been investigated experimentally,  and has turned out to be similar to that of the site-diluted system Cu1-xMgxGeO3. Besides the problem whether the substitution of Si for Ge simply corresponds to the bond dilution or not, we have to investigate the systems at finite but low temperatures, explicitly taking into account the interlayer coupling. Experimentally, analyses on systems with further smaller $`x`$ and at lower temperatures are desirable.

## 4 Summary

We investigated dilution effects on the spin-1/2 quantum AF Heisenberg model with the spin-gapped ground state. When the strong bonds are removed from the lattice for the spin-gapped ground state, two effective interactions with the couplings $`\stackrel{~}{J}_{mn}`$ and $`\stackrel{~}{J}_{\mathrm{af}}`$ are induced between the effective spins. The competition between $`\stackrel{~}{J}_{mn}`$ and $`\stackrel{~}{J}_{\mathrm{af}}`$ leads to peculiar quantum phase transitions. In the low-concentration region, $`\stackrel{~}{J}_{\mathrm{af}}`$ between edge spins mainly contributes to magnetic properties and, therefore, the new spin gap due to singlet pairs reformed by $`\stackrel{~}{J}_{\mathrm{af}}`$ is induced within the original spin gap. The value of the spin gap decreases as $`x`$ increases and tends to vanish at $`x_\mathrm{c}^{}`$ which is less than $`x_\mathrm{c}`$ or just at $`x_\mathrm{c}`$. If the former is the case, the quantum Griffiths phase would exist. Thus there are quite different magnetic properties at $`T=0`$ between the site- and bond-diluted systems.

## Acknowledgements

Most of numerical calculations in the present work have been performed on the SGI 2800 at Institute for Solid State Physics, University of Tokyo. The program is based on ’Looper version 2’ developed by S.T. and K. Kato and ’PARAPACK version 2’ by S.T. The present work is supported by the “Research for the Future Program” (JSPS-RFTF97P01103) of Japan Society for the Promotion of Science. S.T. acknowledges support of the Swiss National Science Foundation.


text:
# Self-induced transparency and giant nonlinearity in doped photonic crystals

## I Introduction

Photonic crystals (PCs) can exhibit an interplay between Bragg reflections, which block the propagation of light in photonic band gaps (PBGs) pbg-bib ; yablon ; john ; joannopoulos ; defect\_a ; defect\_b , and the dynamical modifications of these reflections by nonlinear light-matter interactions kofman ; john\_wang ; pl ; Ze95 ; Scal94 ; Scal94a . A very interesting situation arises when foreign atoms or ions—dopants—with transition frequencies within the PBG are implanted in the PC kofman ; john\_wang ; pl . Then light near one of these frequencies resonantly interacts with the dopants and is concurrently affected by the PBG dispersion. Consequently, highly nonlinear processes with a rich variety of unusual PC-related features are anticipated.

Our aim in recent years has been to identify those regimes of nonlinear optical propagation in doped PCs that allow transmission of extremely weak pulses, while filtering out undesirable noise, and are therefore highly advantageous for optical communications, data storage and processing, near or at the quantum limit. These requirements are satisfied by novel regimes surveyed in this paper that have been theoretically discovered and investigated by us: a) Self-induced transparency (SIT) solitons propagating inside or near a PBG at a frequency that is near resonant with the transition frequency of the dopant Kuri2001 ; Kozh95 ; Kozh98a ; Opat99 . This peculiar form of gap solitons (GSs) is immune to resonant absorption even for a small number of photons, and may also possess two- or three-dimensional (2D or 3D) localization in the form of light bullets (LBs)Blaa00a  (Sec. II). b) Cross-coupling of SIT and electromagnetically-induced transparency (EIT) pulses in PCs. We put forward a new regime in Sec. III: a strong modulation of the phase of a weak pulse subject to EIT by a control pulse in the form of an SIT GS moving at the same slow velocity. Thereby giant cross-phase modulation can be formed between the GS and the EIT pulses. In Sec. IV we summarize and discuss our findings in this paper.

## II Self-induced transparency (SIT) gap solitons

### II.1 Background

A GS is usually understood as a self-localized moving or standing (quiescent) bright region, where light is confined by Bragg reflections against a dark background. The soliton spectrum is tuned away from the Bragg resonance by the nonlinearity at sufficiently high field intensities. The first type of GS had been predicted Chri89 ; Acev89 ; Feng93 ; Ster94 , and later observed Eggl96 , in a Bragg grating possessing Kerr-nonlinearity. A principally different mechanism of GS formation has been theoretically discovered by our group in a periodic array of thin layers of resonant two-level atoms (TLA) separated by half-wavelength nonabsorbing dielectric layers, i.e., a resonantly absorbing Bragg reflector (RABR) Kuri2001 ; Kozh95 ; Kozh98a ; Opat99 . As opposed to the $`2\pi `$-solitons arising in SIT, i.e., resonant field–TLA interaction in uniform media McCa67 ; McCa69 , their GS counterparts in a RABR may have an arbitrary pulse area Kozh95 ; Kozh98a . It must be stressed that stable, moving or standing, GS solutions have been consistently obtained only in a RABR with thin active TLA layers. By contrast, the case of a periodic structure uniformly doped with active TLA calls either for a solution of the wave (Maxwell) equation without the spatial slowly-varying envelope approximation (SVEA), or for a solution of an infinite set of coupled Bloch equations for all spatial harmonics of the atomic polarization (Fourier components) Kuri2001 . Therefore, an attempt Akoz98  to obtain a self-consistent solution for a uniformly-doped periodic structure by imposing the SVEA, or by arbitrarily truncating the infinite hierarchy of equations for the harmonics of atomic population inversion and polarization to its first two orders, is generally unjustified. In fact, it can be shown numerically to fail for many parameter values.

In the simplest case of a uniform (bulk) medium, when the driving field is resonant with the atomic transition, the TLA Bloch equations can be easily integrated and the Maxwell equation then reduces to the sine-Gordon equation

$$\frac{^2\theta }{\zeta \stackrel{~}{\tau }}=\mathrm{sin}\theta $$
(1)
for the pulse area $`\theta =_{\mathrm{}}^t\mathrm{\Omega }𝑑t^{}`$, i.e., the time integral of the Rabi frequency $`\mathrm{\Omega }`$. Equation (1) is written in terms of the dimensionless variables $`\stackrel{~}{\tau }=(tn_0z/c)/\tau _0`$ and $`\zeta =n_0z/c\tau _0`$, where $`\tau _0=\frac{n_0}{\mu }\sqrt{\frac{\mathrm{}}{2\pi \omega _c\varrho _0}}`$, is the cooperative resonant absorption time, $`\varrho _0`$ being the TLA density (averaged over $`z`$), $`\mu `$ the dipole moment of the TLA transition at frequency $`\omega _0`$ and $`n_0`$ is the refraction index of the host medium. This sine-Gordon equation is known to have solitary-wave solutions, which propagate without attenuation or distortion with a conserved pulse area of $`2\pi `$ McCa67 ; McCa69 . These SIT solitons have the form:

$$\mathrm{\Omega }(\zeta ,\stackrel{~}{\tau })=(\tau _0)^1A_0\text{sech}[\beta (\zeta v\stackrel{~}{\tau })],$$
(2)
where the pulse width $`\beta `$ is an arbitrary real parameter uniquely defining the amplitude $`A_0=2/\beta `$ and group velocity $`v=1/\beta ^2`$ of the soliton. In what follows, Eq. (2) will be compared with an SIT GS in a RABR.

### II.2 SIT in RABR: The Model

Let us assume Kuri2001 ; Kozh95 ; Kozh98a ; Opat99  a one-dimensional (1D) periodic modulation of the linear refractive index $`n^2(z)=n_0^2[1+a_1\mathrm{cos}(2k_cz)]`$. The periodic grating gives rise to a PBG with a central frequency $`\omega _c=k_cc/n_0`$ and gap edges at $`\omega _{1,2}=\omega _c\left(1\pm a_1/4\right)`$. The electric field $`E`$ of a pulse propagating along $`z`$ can be expressed by means of the dimensionless quantities $`\mathrm{\Sigma }_\pm 2\tau _o\mu \mathrm{}^1(_F\pm _B)`$, where $`_F`$ and $`_B`$ denote the slowly varying amplitudes of the forward and backward propagating fields, respectively, as

$$E(z,t)=\mathrm{}(\mu \tau _0)^1\left\{\mathrm{Re}[\mathrm{\Sigma }_+(z,t)e^{i\omega _ct}]\mathrm{cos}k_cz\mathrm{Im}[\mathrm{\Sigma }_{}(z,t)e^{i\omega _ct}]\mathrm{sin}k_cz\right\}.$$
(3)
We further assume that very thin TLA layers (much thinner than $`1/k_c`$), whose resonance frequency $`\omega _0`$ is close to the gap center $`\omega _c`$, are placed at the maxima of the modulated refraction index (Fig. 1). They are located at positions $`z_j`$ such that

$$e^{ik_cz_{2j}}=1,e^{ik_cz_{2j+1}}=1,$$
(4)
i.e., the TLA density is described by $`\varrho =\varrho _0\lambda /2_j\delta (zz_j)`$, where $`\lambda =2\pi /k_c`$ is the wavelength.

The Bloch equations for the slowly varying polarization envelope $`P`$ and inversion $`w`$ in the even numbered layers can be obtained (in the slowly varying envelope approximation) by substituting for the Rabi frequency $`\mathrm{\Omega }=\tau _0^1(\mathrm{\Sigma }_+\mathrm{cos}k_cz+i\mathrm{\Sigma }_{}\mathrm{sin}k_cz)`$ and applying Eq. (4) at the positions of these layers:

$`{\displaystyle \frac{P}{\tau }}`$ $`=`$ $`i\delta P+\mathrm{\Sigma }_+w,`$ (5)
$`{\displaystyle \frac{w}{\tau }}`$ $`=`$ $`\mathrm{Re}\left(\mathrm{\Sigma }_+P^{}\right).`$ (6)
Combining Eqs. (5) and (6), one can eliminate the TLA population inversion: $`w=\sqrt{1|P|^2}`$. The remaining equation, together with the Maxwell equations for $`\mathrm{\Sigma }_\pm `$ (driven by $`P`$), form a closed system,

$`{\displaystyle \frac{^2\mathrm{\Sigma }_+}{\tau ^2}}{\displaystyle \frac{^2\mathrm{\Sigma }_+}{\zeta ^2}}`$ $`=`$ $`\eta ^2\mathrm{\Sigma }_++2i(\eta \delta )P2\sqrt{1|P|^2}\mathrm{\Sigma }_+,`$ (7)
$`{\displaystyle \frac{^2\mathrm{\Sigma }_{}}{\tau ^2}}{\displaystyle \frac{^2\mathrm{\Sigma }_{}}{\zeta ^2}}`$ $`=`$ $`\eta ^2\mathrm{\Sigma }_{}2{\displaystyle \frac{P}{\zeta }},`$ (8)
$`{\displaystyle \frac{P}{\tau }}`$ $`=`$ $`i\delta P\sqrt{1|P|^2}\mathrm{\Sigma }_+,`$ (9)
where $`\tau t/\tau _0`$, $`\zeta \left(n_0/c\tau _0\right)z`$ and $`\delta (\omega _0\omega _c)\tau _0`$ are the dimensionless time, coordinate, and detuning, respectively, and $`\eta =l_{\mathrm{abs}}/l_{\mathrm{refl}}=a_1\omega _c\tau _0/4`$ is the dimensionless modulation strength, which can be expressed as the ratio of the TLA absorption distance $`l_{\mathrm{abs}}=\tau _0c/n_0`$ to the Bragg reflection distance $`l_{\mathrm{refl}}=4c/(a_1\omega _cn_0)`$. We emphasize that the above equations are obtained using the SVEA, which is valid under the assumption that the Bragg reflection does not appreciably change the pulse envelope over a distance of a wavelength, $`l_{\mathrm{refl}}\lambda `$, whence $`a_12/\pi `$.

To reach general understanding of the dynamics of the model, one should first consider the spectrum produced by the linearized version of Eqs. (7)–(9), which describes weak fields in the limit of infinitely thin TLA layers. Setting $`\mathrm{\Sigma }_+=Ae^{i(\kappa \zeta \chi \tau )}`$, $`\mathrm{\Sigma }_{}=Be^{i(\kappa \zeta \chi \tau )}`$, $`w=1`$, and $`P=Ce^{i(\kappa \zeta \chi \tau )}`$, we obtain from the linearized equation (9) that $`C=i(\delta \chi )^1A`$. Substituting this into Eqs. (7) and (8), we arrive at the dispersion relation for the wavenumber $`\kappa `$ and frequency $`\chi `$ in the form

$$(\chi ^2\kappa ^2\eta ^2)(\chi \delta )\times \left\{(\chi \delta )\left[\chi ^2\kappa ^2(2+\eta ^2)\right]+2(\eta \delta )\right\}=0.$$
(10)
Different branches of the dispersion relation generated by Eq. (10) are shown in Fig. 2. The roots $`\chi =\pm \sqrt{\kappa ^2+\eta ^2}`$ (corresponding to the solid lines in Fig. 2) originate from the driven equation (8) and represent the dispersion relation of a Bragg reflector with the gap $`|\chi |<\eta `$, that does not feel the interaction with the active layers. Important roots of Eq. (10) are those of the expression in the curly brackets, shown by the dashed and dash-dotted lines in Fig. 1. These roots correspond to nontrivial spectral features: bright or dark solitons in the indicated (shaded) bands.

The frequencies corresponding to $`k=0`$ are $`\chi _0=\eta `$ and $`\chi _{0,\pm }=(\eta \delta )/2\pm \sqrt{2+(\eta +\delta )^2/2}`$, while at $`k^2\mathrm{}`$ the asymptotic expressions for different branches of the dispersion relation are $`\chi =\pm k`$ and $`\chi =\delta +2\left(\eta \delta \right)k^2`$. Thus, the linearized spectrum always splits into two gaps, separated by an allowed band, except for the special case, $`\eta =\eta _0\delta /2+\sqrt{1+\delta ^2/4}`$, when the upper gap closes down. The upper and lower band edges are those of the periodic structure, shifted by the induced TLA polarization in the limit of a strong reflection. They approach the SIT spectral gap for forward- and backward-propagating waves Mant95  in the limit of weak reflection. The allowed middle band corresponds to a polaritonic excitation (collective atomic polarization) in the periodic structure.

### II.3 Standing (quiescent) self-localized pulses

We seek the stationary solutions of Eqs. (7) and (9) corresponding to bright solitons in the form

$$\mathrm{\Sigma }_+=e^{i\chi \tau }𝒮(\zeta ),P=ie^{i\chi \tau }𝒫(\zeta )$$
(11)
with real $`𝒫`$ and $`𝒮`$. Substituting this into (9), we eliminate $`𝒫`$ in favor of $`𝒮`$ and obtain an equation for $`𝒮(\zeta )`$,

$$\frac{d^2𝒮}{d\zeta ^2}=(\eta ^2\chi ^2)𝒮2𝒮\frac{(\eta \chi )\mathrm{sign}(\chi \delta )}{\sqrt{(\chi \delta )^2+𝒮^2}}.$$
(12)
It then follows Kozh98a  that bright solitons can appear in two frequency bands $`\chi `$, the lower band being $`\chi _1<\chi <\mathrm{min}\{\chi _2,\eta ,\delta \}`$, and the upper band being $`\mathrm{max}\{\chi _1,\eta ,\delta \}<\chi <\chi _2`$, where $`\chi _{1,2}1/2[\delta \eta \sqrt{(\eta +\delta )^2+8}]`$ are the boundary frequencies. The lower band exists for all values $`\eta >0`$ and $`\delta `$, while the upper one only exists for the weak-reflectivity case $`\delta >\eta 1/\eta `$. On comparing these expressions with the spectrum shown in Fig. 2, we conclude that part of the lower gap is always empty from solitons, while the upper gap is completely filled with stationary solitons in the weak-reflectivity case, and completely empty in the opposite limit.

In an implicit form, the solution of Eq. (12) reads

$$𝒮(\zeta )=2|\chi \delta |(\zeta )\left(1^2(\zeta )\right)^1,$$
(13)
with

$$|\zeta |=\sqrt{2\left|\frac{\chi \delta }{\chi \eta }\right|}\left[(1_0^2)^{1/2}\mathrm{tan}^1\sqrt{\frac{_0^2^2}{1_0^2}}+(2_0)^1\mathrm{ln}\left(\frac{_0+\sqrt{_0^2^2}}{}\right)\right],$$
(14)
and $`_0^2=1|(\chi +\eta )(\chi \delta )|/2`$. This zero-velocity (ZV) gap soliton is always single-humped and its amplitude, found from Eq. (14), is given by

$$𝒮_{\mathrm{max}}=4_0/\sqrt{\left|\chi +\eta \right|}.$$
(15)
To calculate the electric field in the antisymmetric $`\mathrm{\Sigma }_{}`$ mode, we substitute $`\mathrm{\Sigma }_{}=ie^{i\chi \theta }𝒜(\zeta )`$ into Eq. (8) and obtain

$$𝒜^{\prime \prime }+\left(\chi ^2\eta ^2\right)𝒜=2𝒫^{},$$
(16)
which can be easily solved by the Fourier transform, once $`𝒫(\zeta )`$ is known. We note that, depending on the parameters $`\eta `$, $`\delta `$ and $`\chi `$, the main part of the soliton energy can be carried either by the $`\mathrm{\Sigma }_+`$ or the $`\mathrm{\Sigma }_{}`$ mode.

The most drastic difference of these new solitons from the well-known SIT solitons in Eq. (2) is that the area of the ZV soliton (integrated over $`\zeta `$) is not restricted to $`2\pi `$, but, instead, may take an arbitrary value. This basic new feature shows that the Bragg reflector can enhance (by multiple reflections) the field coupling to the TLA, so as to make the pulse area effectively equivalent to $`2\pi `$. In the limit of the small-amplitude and small-area solitons, $`_0^21`$, Eq. (14) can be easily inverted, the ZV soliton becoming a broad sech-like pulse:

$$𝒮2|\chi \delta |_0\mathrm{sech}\left(\sqrt{2\left|\frac{\chi \eta }{\chi \delta }\right|}_0\zeta \right).$$
(17)
In the opposite limit, $`1_0^20`$, i.e., for vanishingly small $`|\chi +\eta |`$, the the soliton is characterized by a broad central part with a width $`\left(1_0^2\right)^{1/2}`$ and its amplitude (15) becomes very large. Thus, although the ZV soliton has a single hump, its shape is, in general, strongly different from that of the traditional nonlinear-Schrödinger (NLS) sech pulse.

### II.4 Moving solitons

One could expect a translational invariance of the ZV solitons (13) on physical grounds. Hence, a full family of soliton solutions should have velocity as one of its parameters. This can be explicitly demonstrated in the limit of the small-amplitude large-width solitons \[cf. Eq. (17)\]. We search for the corresponding solutions in the form $`\mathrm{\Sigma }_+(\zeta ,\tau )=𝒮(\zeta ,\tau )\mathrm{exp}(i\chi _0\tau )`$, $`P(\zeta ,\tau )=i𝒫(\zeta ,\tau )\mathrm{exp}(i\chi _0\tau )`$ \[cf. Eqs. (11)\], where $`\chi _0`$ is the frequency corresponding to $`k=0`$ on any of the three branches of the dispersion relation (10) (see Fig. 2), and the functions $`𝒮(\zeta ,\tau )`$ and $`𝒫(\zeta ,\tau )`$ are assumed to be slowly varying in comparison with $`\mathrm{exp}(i\chi _0\tau )`$. Under these assumptions, we arrive at the following asymptotic equation for $`𝒮(\zeta ,\tau )`$:

$$\left[2i\frac{\chi _0(\chi _0\delta )^2\eta +\delta }{(\chi _0\delta )^2}\frac{}{\tau }+\frac{^2}{\zeta ^2}+\frac{\chi _0\eta }{(\chi _0\delta )^3}|𝒮|^2\right]𝒮=\left(\eta ^2\chi _0^2+2\frac{\chi _0\eta }{\chi _0\delta }\right)𝒮.$$
(18)
Since this equation is of the NLS form, it has the full two-parameter family of soliton solutions, including the moving ones Newe92 .

In order to check the existence and stability of the moving solitons numerically, the following procedure has been used Kozh98a : Eqs. (7) and (9) were simulated for an initial configuration in the form of the ZV soliton multiplied by $`\mathrm{exp}(ip\zeta )`$ with some wavenumber $`p`$, in order to ‘push’ the soliton. The results demonstrate that, at sufficiently small $`p`$, the ‘push’ indeed produces a moving stable soliton \[Fig. 3(a)\]. However, if $`p`$ is large enough, the multiplication by $`\mathrm{exp}(ip\zeta )`$ turns out to be a more violent perturbation, splitting the initial pulse into two solitons, one quiescent and one moving \[Fig. 3(b)\].

### II.5 Light bullets (spatiotemporal solitons) in PCs

The advantageous properties of SIT GS can be supplemented by immunity to transverse diffraction, i.e., simultaneous transverse and longitudinal self-localization of light in a PC: multi-dimensional spatio-temporal solitons or “light bullets” (LBs) Silb90  have been analytically and numerically predicted by our group to exist and be stable, not only in uniform 2D and 3D SIT media Blaa00 , but also in 2D or 3D periodic structures, wherein SIT solutions combining LB and GS properties are demonstrated Blaa00a . Our objective is to consider the propagation of an electromagnetic wave with a frequency close to $`\omega _c`$ through a 2D PC doped by thin TLA layers. The forward- and backward-propagating components satisfy equations that are a straightforward generalization of the 1D equations (7) and (8)

$`i{\displaystyle \frac{^3\mathrm{\Sigma }_+}{\tau x^2}}+i{\displaystyle \frac{^3\mathrm{\Sigma }_{}}{\zeta x^2}}+{\displaystyle \frac{^2\mathrm{\Sigma }_+}{\tau ^2}}{\displaystyle \frac{^2\mathrm{\Sigma }_+}{\zeta ^2}}+\eta {\displaystyle \frac{^2\mathrm{\Sigma }_+}{x^2}}+\eta ^2\mathrm{\Sigma }_+2{\displaystyle \frac{P}{\tau }}2i\eta P=0,`$ (19)
$`i{\displaystyle \frac{^3\mathrm{\Sigma }_{}}{\tau x^2}}+i{\displaystyle \frac{^3\mathrm{\Sigma }_+}{\zeta x^2}}+{\displaystyle \frac{^2\mathrm{\Sigma }_{}}{\tau ^2}}{\displaystyle \frac{^2\mathrm{\Sigma }_{}}{\zeta ^2}}\eta {\displaystyle \frac{^2\mathrm{\Sigma }_{}}{x^2}}+\eta ^2\mathrm{\Sigma }_{}+2{\displaystyle \frac{P}{\zeta }}=0,`$ (20)
where the Fresnel number $`F>0`$, which governs the transverse diffraction in the 2D and 3D propagation, was incorporated into $`x`$ denoting the transverse coordinate. The equations for the polarization $`P`$ and inversion $`w`$ are the same as Eqs. (5) and (6).

We search for analytical LB solutions of Eqs. (19), (20), (5) and (6), by the following ansatz that reduces in 1D to the exact moving GS Kuri2001 ; Kozh95

$`\mathrm{\Sigma }_+`$ $`=`$ $`A_0\sqrt{\text{sech}\mathrm{\Theta }_1\text{sech}\mathrm{\Theta }_2}e^{i(\kappa \zeta \chi \tau )+i\pi /4},`$ (21)
$`\mathrm{\Sigma }_{}`$ $`=`$ $`\mathrm{\Sigma }_+/v,`$ (22)
$`P`$ $`=`$ $`\sqrt{\text{sech}\mathrm{\Theta }_1\text{sech}\mathrm{\Theta }_2}\{(\mathrm{tanh}\mathrm{\Theta }_1+\mathrm{tanh}\mathrm{\Theta }_2)^2+`$ (23)
$`{\displaystyle \frac{\delta \eta }{4\eta }}C^4[(\mathrm{tanh}\mathrm{\Theta }_1\mathrm{tanh}\mathrm{\Theta }_2)^22(\text{sech}^2\mathrm{\Theta }_1+\text{sech}^2\mathrm{\Theta }_2)]^2\}^{1/2}e^{i(\kappa \zeta \chi \tau )+i\nu },`$
$`w`$ $`=`$ $`\left[1|P|^2\right]^{1/2},`$ (24)
with $`\mathrm{\Theta }_1(\tau ,\zeta )\beta (\zeta v\tau )+\mathrm{\Theta }_0+Cx`$, $`\mathrm{\Theta }_2(\tau ,\zeta )\beta (\zeta v\tau )+\mathrm{\Theta }_0Cx`$, the phase $`\nu `$ and coefficients $`\mathrm{\Theta }_0`$ and $`C`$ being real constants, while the other parameters are defined as $`A_0=2\sqrt{\delta /\eta 1}`$, $`\beta =\sqrt{\delta /\eta +1}`$, $`v=\sqrt{(\delta \eta )/(\delta +\eta )}`$, $`\kappa =\sqrt{\delta ^2\eta ^2}`$, and $`\chi =\delta `$.

The ansatz (21)- (24) satisfies Eqs. (19) and (20) exactly, while Eqs. (5) and (6) are satisfied to order $`\sqrt{\delta /\eta 1}C^2`$, which requires that $`\sqrt{\delta /\eta 1}C^21`$. The ansatz applies for arbitrary $`\eta `$, admitting both weak ($`\eta 1`$) and strong ($`\eta >1`$) reflectivities of the Bragg grating, provided that the detuning remains small with respect to the gap frequency. Comparison with numerical simulations of Eqs. (19), (20), (5) and (6), using Eqs. (21)-(24) as an initial configuration, tests this analytical approximation and shows that it is indeed fairly close to a numerically exact solution; in particular, the shape of the bullet remains within 98% of its originally presumed shape after having propagated a large distance, as is shown in Fig. 4.

Three-dimensional (3D) LB solutions with axial symmetry have also been constructed in an approximate analytical form and successfully tested in direct simulations, following a similar approach Blaa00a . Generally, they are not drastically different from their 2D counterparts described above.

### II.6 Information transmission by SIT GSs and LBs

The efficiency of information transmission is characterized either by channel (information) capacity $`C=W\mathrm{ln}(I_s/I_n)`$, where $`W`$ is the bandwidth and $`I_s/I_n`$ is the ratio of the signal-to-noise intensities, or by the data transmission density $`D=NM`$, where $`N`$ is the number of bits per channel and $`M`$ is the number of accessible channels. Both $`C`$ and $`D`$ can be very high in the case of a SIT GS or LB for the following reasons: (a) The bandwidth $`W`$ is large, being limited by the PBG width, which can be very large in the optical domain. At the same time, noise is very effectively suppressed by the Bragg reflection and by the absence of diffraction losses in the case of a LB. (b) The maximal transmission density $`D`$ can be estimated Mossberg  as the ratio of the accessible bandwidth, in our case the PBG width (in excess of $`10^{13}`$ s<sup>-1</sup> in the optical domain), to the spontaneous linewidth ($`10^6`$ s<sup>-1</sup> for rare-earth ions). Hence, these modes of transmission can be very effective for optical communications.

## III Cross coupling between electromagnetically-induced and self-induced transparency pulses

### III.1 EIT in bulk media: Background

Electromagnetically induced transparency (EIT) is based on the phenomenon of coherent population trapping eit\_a ; eit\_b , in which the application of two laser fields to a three-level atomic system creates the so-called “dark state”, which is stable against absorption of both fields. Consider a four-level atomic system whose level configuration is depicted in Fig. 5(a). The phase shift and absorption of an optical field $`E_i`$ are given by the real and imaginary parts of its complex polarizability $`\alpha _i`$. In the absence of the “control” field $`E_b`$, the usual EIT spectrum \[Fig. 5(b)\] for the weak probe field $`E_a`$ exhibits vanishing phase shift and absorption \[$`\text{Re}(\alpha _a)=\text{Im}(\alpha _a)=0`$\] at the two-photon Raman resonance $`\omega _a=\omega _{21}+\omega _d`$, where $`\omega _a`$ and $`\omega _d`$ are the frequencies of the probe and driving fields, respectively, and $`\omega _{ij}`$ is the frequency of the atomic transition $`|i|j`$. An off-resonant control field $`E_b`$ with the frequency $`\omega _b`$ such that $`|\mathrm{\Delta }_b|=|\omega _b\omega _{43}|\gamma _4`$, where $`\gamma _i`$ is the decay rate of the corresponding atomic level, induces an ac Stark shift of level $`|3`$ and thereby shifts the EIT spectrum \[Fig. 5(b)\].

Due to the steepness of the dispersion curve in the vicinity of the Raman resonance, $`|_{\omega _a}\text{Re}(\alpha _a)||_{\omega _a}\text{Im}(\alpha _b)|`$, this Stark shift leads to a large phase shift along with small absorption of the probe field $`E_a`$:

$$\varphi _a=\text{Re}(\alpha _a)z\frac{\alpha _0\gamma _2|\mathrm{\Omega }_b|^2}{2\mathrm{\Delta }_b|\mathrm{\Omega }_d|^2}z,\text{Im}(\alpha _a)=\frac{\gamma _4\text{Re}(\alpha _a)}{2\mathrm{\Delta }_b}\text{Re}(\alpha _a),$$
(25)
where $`\alpha _0`$ is the resonant absorption coefficient of the medium at the frequency $`\omega _{21}`$ and $`\mathrm{\Omega }_i=\mu _{ij}E_i/\mathrm{}`$ is the Rabi frequency of the corresponding field ($`\mu _{ij}`$ the dipole matrix element on the respective transition). This is the essence of the so-called giant Kerr cross-phase modulation of a probe field by a control field, introduced first by Schmidt and Imamoğlu imam . Later Harris and Yamamoto harris\_a  have predicted that a resonant control field $`E_b`$ with $`|\mathrm{\Delta }_b|<\gamma _4`$ can destroy the coherence between the two ground levels $`|1`$ and $`|3`$, which leads to a two-photon absorption $`\text{Im}(\alpha _{a,b})=\frac{\alpha _0\gamma _2|\mathrm{\Omega }_{b,a}|^2}{\gamma _4|\mathrm{\Omega }_d|^2}`$, that is, the medium absorbs two fields simultaneously, but does not absorb one field alone. This is the essence of a probe-photon switch, conditional on the presence of control photons.

The main limitation of the above schemes imam ; harris\_a ; harris\_b  stems for the fact that the effective interaction length is limited by the mismatch between the group velocity of the slowly propagating $`E_a`$ field, $`v_g^{(a)}\frac{2|\mathrm{\Omega }_d|^2}{\alpha _0\gamma _3}c/n_0`$ and that of the nearly-free propagating $`E_b`$ field, $`v_g^{(b)}c/n_0`$. For weak (few-photon) pulses, this mismatch ultimately limits the maximal phase shift or absorption of the probe in the presence of the control field.

### III.2 Simultaneous EIT and SIT in RABR

In this section we propose a new implementation of the cross-phase modulation, in which the group velocities of both fields can be matched, allowing one to obtain any desired phase shift of the probe field with a weak control field. To this end, we consider the same configuration as in Sec. II, leading to SIT GS and LB solutions: a PC periodically doped by thin layers of atoms at the maxima of its refractive index. However, the multi-level structure of the atoms is now playing a role: it is shown in Fig. 6, along with the polarizations and propagation directions of the fields involved. The states $`|1`$, $`|3`$ and $`|5`$ are the degenerate Zeeman components with $`M_F=1,0,+1`$, respectively, of the atomic ground level having total angular momentum $`F=1`$. Similarly, the states $`|4`$ and $`|6`$ are the degenerate Zeeman components with $`M_F=1,0`$, respectively, of the excited level having angular momentum $`F=1`$. Finally, the state $`|2`$ corresponds to the single Zeeman component with $`M_F=0`$ of another excited level having $`F=0`$. Such a level scheme is found, e.g., in alkali atoms, where the ground level is $`S_{1/2},F=1`$ and the two excited levels are $`P_{1/2},F=1`$ and $`P_{3/2},F=0`$. Due to the dipole selection rules, the $`\pi `$-polarized driving field couples the states with $`\mathrm{\Delta }M=0`$, the $`\sigma _+`$-polarized $`E_a`$ field couples the states with $`\mathrm{\Delta }M=1`$ and the $`\sigma _{}`$-polarized $`E_b`$ field couples the states with $`\mathrm{\Delta }M=1`$.

We assume that initially all the atoms are optically pumped into the states $`|1`$ and $`|5`$, which then acquire equal populations 1/2. Hence, the sequence of transitions $`|1|2|3|4`$ repeats that of Fig. 5(a), realizing the cross-phase modulation scheme of Sec. IIIIII.1. The frequency of the $`E_a`$ field is far from the band gap frequencies of the PC, while the frequency of the $`E_b`$ field is within the band gap.

As was shown in Sec. II, PC structures doped with the near-resonant TLAs can support standing and slowly moving SIT GSs, whose pulse area (integrated over $`z`$) can take an arbitrarily small value. In the present setup, the transition $`|5|6`$ realizes that near-resonant TLA, allowing for slow propagation of the $`E_b`$ field through the PC.

Let us write the propagation equation for the slowly moving SIT soliton in the form

$$\mathrm{\Omega }_{\text{SIT}}(z,t)=\mathrm{\Omega }_b\text{sech}\left(\frac{tv_g^{(b)}z}{2\beta }\right),$$
(26)
where $`\mathrm{\Omega }_b`$ is the peak Rabi frequency, $`v_g^{(b)}`$ is the group velocity and $`\beta =(2\alpha _b)^1`$, with $`\alpha _b`$ being the absorption coefficient of the active medium at the carrier frequency $`\omega _b`$ of the soliton. The temporal width of the pulse is given by $`\tau _b=2\beta /v_g^{(b)}`$. The area of the $`E_b`$ pulse

$$\theta _b=_{\mathrm{}}^{\mathrm{}}\mathrm{\Omega }_{\text{SIT}}𝑑t=\frac{\mathrm{\Omega }_b}{v_g^{(b)}\alpha _b}\pi ,$$
(27)
is then inversely proportional to the group velocity of the pulse. Hence, the SIT condition (Sec. IIIIII.1) $`\theta _b=2\pi `$ imposes a unique relation between the Rabi frequency of the SIT soliton and its group velocity:

$$v_g^{(b)}=\frac{\mathrm{\Omega }_b}{2\alpha _b}.$$
(28)
The absorption-free propagation of the SIT soliton is limited to $`z<v_g^{(b)}/\gamma _6`$, where $`\gamma _6`$ is the decay rate of the upper atomic state $`|6`$.

Our aim is to match the group velocities of the $`E_a`$ field subject to EIT and the $`E_b`$ field having the form of a slow SIT gap soliton: $`v_g^{(a)}=v_g^{(b)}`$. This requires that $`|\mathrm{\Omega }_d|^2=\mathrm{\Omega }_b\alpha _a\gamma _2/4\alpha _b`$, i.e., an appropriate choice of the driving field Rabi frequency $`\mathrm{\Omega }_d`$, for a given Rabi frequency $`\mathrm{\Omega }_b`$ of the control field. Such velocity matching of the two copropagating weak fields would maximize their interaction.

One possibility to launch the required slow SIT soliton is to irradiate the PC by a laser beam at a small angle $`\psi `$ relative to the periodicity direction $`z`$, $`\psi Dn_0v_g^{(b)}/(Lc)1`$, where $`D`$ and $`L`$ are, respectively, the transverse and longitudinal dimensions of the structure. This choice of $`\psi `$ ensures that the $`z`$-component of the beam, which forms the SIT soliton and propagates in the PC with the group velocity $`v_g^{(b)}`$ over the distance $`L`$, will traverse the structure during the same time as the transverse component of that beam, which covers the distance $`D`$ with the velocity $`(c/n_0)\mathrm{sin}\psi `$.

We have checked that for the parameter values corresponding to dopant atoms (or ions) with the mean density $`N=10^{13}`$ cm<sup>-3</sup> (surface density of $`4\times 10^8`$ cm<sup>-2</sup> in the thin layers), $`\mathrm{\Delta }_b=30\gamma _4`$, $`|\mathrm{\Omega }_b|10^6`$ rad/s and $`|\mathrm{\Omega }_d|4\times 10^6`$ rad/s, we obtain $`\pi `$ phase shift of the $`E_a`$ field over a distance $`z4`$ cm, while the absorption probability remains less that 10%.

One possible difficulty of our scheme is that, with the parameters listed above, the temporal width of the $`E_b`$ field is $`\tau _b10^6`$ s, and the interaction time is of the order of $`10^3`$ s, while the lifetime of the SIT soliton is of the order of the decay time of the excited atomic state $`\gamma _6^110^7`$ s. One can cope with this problem by employing the atomic level scheme shown in Fig. 7, which allows one to launch Raman solitons. We irradiate the system with an additional strong cw field $`E_R`$, which couples the $`F=1`$ excited state with the $`F=2`$ metastable ground states. The fields $`E_R`$ and $`E_b`$ are largely detuned from the fast decaying excited states $`|4`$ and $`|6`$ by an amount $`\mathrm{\Delta }\gamma _{4,6}`$. Then, upon adiabatically eliminating the states $`|4`$ and $`|6`$, we obtain that the Rabi frequency $`\mathrm{\Omega }_b`$ of the control field in Eq. (26) is simply replaced by $`\mathrm{\Omega }_b\mathrm{\Omega }_R/\mathrm{\Delta }`$. The lifetime of the SIT soliton is given now by the lifetime of the $`F=2`$ ground states, which can be very large, reaching in some instances a fraction of a second! In addition, such a setup allows one to launch slow Raman GSs perlin , and thus circumvent the difficulty of launching a standing (ZV) or slowly moving GSs, which must otherwise overcome the high reflectivity of the PC boundaries.

## IV Conclusions

In this paper we have focused on properties of solitons in a doped PC or RABR, combining a periodic refractive-index superlattice (Bragg reflector in 1D or 2D) and a periodic set of thin active layers (consisting of TLAs resonantly interacting with the field). We have demonstrated that the system supports a vast family of bright GSs, whose properties differ substantially from their counterparts in periodic structures with either cubic or quadratic off-resonant nonlinearities. Depending on the initial conditions, these can be either standing (ZV) or slowly moving stable solitons that exhibit SIT irrespective of their photon number (pulse energy) for an appropriate group velocity. A multidimensional version of this model corresponds to a periodic set of thin active layers placed at the maxima of a 2D- or 3D-periodic refractive index. It has been found to support stable propagation of spatiotemporal solitons in the form of 2D- and 3D-localized LBs.

The best prospect of realizing a PC which is adequate for observing the GSs and LBs is to use thin layers of rare-earth ions Grei99  embedded in a spatially-periodic semiconductor structure Khit99 . The TLAs in the layers should be rare-earth-ions with the density of $`10^{15}10^{16}`$ cm<sup>-3</sup>, and large transition dipole moments. The parameter $`\eta `$ can vary from 0 to 100 and the detuning is $`10^{12}10^{13}`$ s<sup>-1</sup>. Cryogenic conditions in such structures can strongly extend the dephasing time $`T_2`$ and thus the soliton’s or LB’s lifetime, well into the $`\mu `$sec range Grei99 , which would greatly facilitate the experiment.

In a 2D PC, LBs can be envisaged to be localized on the time and transverse-length scales, respectively, $`10^{12}`$ s and $`1\mu `$m. The incident pulse has uniform transverse intensity and the transverse diffraction is strong enough. One needs $`d^2/l_{\mathrm{abs}}\lambda _0<1`$, where $`l_{\mathrm{abs}}`$, $`\lambda _0`$ and $`d`$ are the resonant-absorption length, carrier wavelength, and the pulse diameter, respectively. For $`l_{\mathrm{abs}}10^3`$ m and $`\lambda _010^4`$ m, one thus requires $`d<10^4`$ m, which implies that the transverse size of the PC must be a few $`\mu `$m.

We have considered here (Sec. IIIIII.2) the cross-coupling of optical beams in a PC or RABR. We have pointed out, for the first time, the advantageous features of the cross coupling between EIT and SIT pulses, which is capable of producing extremely strong correlations between the two pulses. With doping parameters as above, and driving and control fields with Rabi frequencies of the order of $`10^6`$ rad/s, we can obtain a phase shift of $`\pi `$ for the weak probe pulse over a distance of few cm. This is much larger than any corresponding phase shift (for similar control fields) in other media.

We strongly believe that the highly promising payoff expected from the construction of suitable structures justifies the experimental challenge they pose. If and when the schemes proposed above are experimentally realized, they may prove to be useful for producing ultrasensitive nonlinear phase shifters or logical photon switches for both classical and quantum information processing or communication, owing to the unique advantages of the doped PCs over conventional EIT schemes imam ; harris\_a ; harris\_b ; lukin  or high-Q cavities cavity\_a ; cavity\_b :

## Acknowledgments

This work was supported by the EU (ATESIT) Network, the US-Israel BSF and the Feinberg Fellowship (D.P.).