text:
# References

Theoretical Investigation of One-Dimensional Cavities in Two-Dimensional Photonic Crystals

S. Foteinopoulou 1 and C.M. Soukoulis 1,2

1Ames Laboratory-USDOE and Department of Physics and Astronomy
Iowa State University, Ames, Iowa 50011

2Research Center of Crete-FORTH and Department of Physics,
University of Crete, Heraklion, Crete 71110, Greece

We study numerically the features of the resonant peak of one-dimensional (1-D) dielectric cavities in a two-dimensional (2-D) hexagonal lattice. We use both the transfer matrix method and the finite difference time-domain (FDTD) method to calculate the transmission coefficient. We compare the two methods and discuss their results for the transmission and quality factor Q of the resonant peak. We also examine the dependence of Q on absorption and losses, the thickness of the sample and the lateral width of the cavity. The Q- factor dependence on the width of the source in the FDTD calculations is also given.

I. Introduction

In semiconductor crystals the presence of the periodic potential affects the properties of the electrons. Likewise in photonic crystals (PCs), that are periodic dielectric arrangements in one, two or three dimensions the properties of the photon can be controlled. For certain frequency regions known as “the photonic band gap”, propagation of electromagnetic (EM) waves is prohibited in the PC . A defect present in an otherwise periodic crystal may introduce one or more propagation modes in the band gap and forms a cavity. Cavities of different geometries can exist in a photonic crystal and the resulting resonant mode can be classified into two general types indicating whether it drops from the “air” (higher) band or raises from the “dielectric” (lower) band . The first type is associated with defects corresponding to the removal of dielectric material while the second corresponds to defects involving the addition of dielectric material. The features of the resonant mode will depend on both the bulk crystal and the cavity characteristics and can so be tuned accordingly .

An important virtue of the PC cavity is that it can control the atomic spontaneous emission, while metallic cavities for the related frequency range are generally lossy . Spontaneous emission is important for a number of semiconductor devices . Resonant cavities have undergone extensive studies , both theoretical and experimental. A variety of applications incorporating a PC cavity have been suggested or reported , such as optical laser components , optical filters , single mode light emiting diode (SMLED)  and optical imaging .

In this paper we focus on the theoretical study of an 1-D dielectric cavity in between a hexagonal patterned region of air cylinders in a dielectric matrix for the case of H polarization (magnetic field along the cylinder axis). We will present the cavity mode frequency versus the cavity length ($`L_c`$) (see Fig. 1), as well as the quality factor (Q) of the modes versus their corresponding frequency. The quality factor is defined as $`\lambda _p`$/$`\delta `$$`\lambda `$ where $`\lambda _p`$ the frequency of the resonant peak and $`\delta \lambda `$ the width of the resonance at transmission half of its maximum value. In our case the quality factor Q of the resonance is relatively small ($``$200), and other ways of determining the Q factor (such as from the energy decay of the resonance ) give similar results. The dependence of the resonant mode on the size of the system (in the lateral direction and in the direction of propagation) will be also shown. In our calculations, two different numerical techinques (the transfer matrix and the FDTD) are used and the results of the two methods are compared and discussed. With both techniques we calculate the transmission through the structure and from that we obtain all the relevant features of the resonant mode.

In Section II we describe briefly the two calculational methods. In Section III we present and discuss our calculational results concerning the quality factor and position of the cavity resonance. We also compare the latter with the values obtained from the experiment previously performed on the structure under study . In Section IV we give a summary of our results.

II. Calculational methods

The first technique we use is the transfer matrix method (TMM) developed by Pendry . In the TMM a grid lattice is used to discretize the space, and the structure is divided into finite blocks along the propagation direction. With the use of the Maxwells equations, that are solved on the grid lattice, the electric and magnetic field can be integrated throughout the blocks and so the respective transmission coefficient can be calculated. The final transmission will result by combining the ones of the individual blocks. In the TMM the modeled structure is finite in the propagation direction (y) but infinite in the x,z directions (where z is the direction of the air cylinders), and the incoming electromagnetic wave is a plane wave. The structure is embedded in a medium with dielectric constant equal to the background dielectric constant to simulate the experiment .

The second method we use is the finite difference time domain (FDTD) technique . The real space is discretized on a grid lattice and Maxwell’s equations are solved in time domain. The electric and magnetic fields are updated in every point of the grid lattice in finite time steps. The structure is infinite in the direction of the cylinders (z direction) but has finite size in both x and y direction and is embedded in a finite sized (in x and y) dielectric slab with dielectric constant equal to that of the matrix medium. Liao absorbing boundary conditions  are applied at the walls of the slab to avoid reflections. The source is a pulse, located close to x=0, y=0 ( y being the direction of propagation), has finite length and generates fields with a trapezoidal (extented with gaussian tails) profile in space. To calculate the transmission, a line detector is placed along the latteral direcion right after the structure. The component that is perpendicular to the detector of the Poynting vector (for the Fourier transformed fields) is taken and averaged over the detector. This is normalized with the respective value for the same detector but positioned close to the source in the absense of the structure (to avoid corrupting the data with reflections) and yields the transmission of the structure.

In both methods the quality factor is calculated from the data of transmission versus frequency data from $`\lambda _p`$/$`\delta `$$`\lambda `$ where $`\lambda _p`$ the frequency of the resonant peak and $`\delta \lambda `$ the width of the resonance at transmission half of its maximum.

The two methods differ not only in the calculational approach but also in the characteristics of the source of the incoming EM waves, and in the lateral size of the sample that is simulated. We also performed calculations using the FDTD method but applying Bloch (periodic) boundary conditions in the lateral (x) direction to have conditions similar to those for the transfer matrix. The source that is used in the latter case is still a pulse but with a plane-wave front.

III. Results

In Fig. 1 the cross section of the cavity structure under study with the x-y plane can be seen. It corresponds to air cylinders ($`ϵ`$a= 1.0) in a GaAs background ($`ϵ`$b= 11.3) (the value of the GaAs dieletric constant for $``$1000 nm is taken in the calculations for simplicity). The whole cavity structure is embedded in GaAs ($`ϵ`$= 11.3). The symmetry directions of the bulk crystal are also shown. With $`L_c`$, we will refer to the cavity width that corresponds to the length of the dielectric defect introduced along the propagation direction. With $`L_w`$, we will refer to the lateral width of the cavity (which is the same as the sample’s lateral width). $`N_c`$ will be the number of rows in each side of the cavity ($`N_c=4`$ for both calculations and experiment). The thickness of the sample then along the propagation direction would be $`L_{th}=2N_cb+L_c`$, where 2b=$`\sqrt{3}`$a ( but is approximated with 1.7a in the TMM for numerical reasons). (Such an approximation introduces an error of at most 3$`\%`$ in the position of the resonance and does not affect the quality factor). The radius R of the air holes is R$``$0.2803a (a being the lattice constant) that corresponds to a filling factor of $``$0.285. (Actually in the TMM because of the approximation mentioned above the simulated structure will have a slightly larger filling factor). No specific value for the lattice constant a is assumed. All the band gap and transmission properties in a photonic crystal scale with the lattice constant a, so all the subsequent results will be presented in dimensionless units of frequency (reduced frequency u=a/$`\lambda `$) and length.

From the transfer matrix results for the various $`L_c`$ we observed that the spectral gap when the cavity is introduced is wider than the one for the periodic crystal. Such widening of the gap was also observed in  when a defect is introduced by removal of 2 rows of sites in a periodic 2-D system of dielectric rods in air background. TMM calculation have been performed for cavity structures with dimensionless $`L_c/a`$ ranging from $``$0.15 to 2.0. At a certain frequency different resonances (peaks) can occur corresponding to different $`L_c`$ values. Each of these peaks is characterized by an order, that basically indicates the order at which a resonant peak appears at this specific frequency while $`L_c`$ is increasing.

In Fig. 2 every cavity width (in units of a) is plotted as a function of the corresponding reduced frequency of the observed peak and three different curves are recovered for the three orders of the resonant peaks. Starting from the lower to the higher curve in Fig. 2 the peaks are the first-, second- and third- order, respectively. It will be seen later that peaks of the same frequency but different order can have different features (e.g., quality factor). It is evident from Fig. 2 that for certain $`L_c`$ values more than one resonant peaks can appear within the gap. The experimental results  for the resonant peaks for various lattice constants (a=200 nm, 210nm, 220nm, 230nm, 240nm) are also included in the figure. The agreement is generally good but there seems to be a small discrepancy that increases with the frequency. This can be attributed to the fact that we have taken the GaAs dielectric constant not to vary with the frequency (for simplicity the value at $``$1000 nm is taken in our calculations) while this is not the case in the actual system. This difference in the dielectric will alter slightly the position of the peaks . This is consistent with the fact that the discrepancy at the higher end of the reduced frequencies seems to be a little larger for the smaller lattice constant (the corresponding wavelength differs more from the value of a 1000 nm). There are also some experimental peaks for a=200 nm that are falling above the theoretical curves, as for example, the first order peaks with $`L_c`$/a from $``$0.60 to 0.85. We have not numerically calculated the Qs for these, because they lie too close to the band edges and the determination of their quality factor can be vaque and so we have not included them in Fig. 2.

We compare now the two computational methods ( TMM and FDTD) that are used to simulate the cavity structure. There are basically three fundamental differences between the two methods. The first one is that the transfer matrix technique is a time independent method while in the FDTD the fields are solved in time domain. Also, in the TMM, the structure is infinite in x by virtue of the periodic boundary conditions applied along this direction, while in the FDTD the system is finite in the x-y plane bounded by absorbing boundary conditions. Finally, the incoming EM fields in the TMM framework are extended plane waves incident on the whole x-z boundary (corresponding to a source infinite in width) while in the FDTD a pulse emitted from a source with finite width is considered. Also the way the transmission coefficient is obtained is different in the two methods. In the TMM it is calculated from the transfer matrix using the field values while in the FDTD it is calculated from the energy ratio that passes through the crystal along the propagation direction. In order to make the comparison between the two methods in this particular FDTD calculation, the 2-D space is discretized as in transfer matrix which gives 2b=1.7a instead of $`\sqrt{3}`$a . The results for one particular cavity width are shown in Fig. 3. Good agreement is found between the two methods regarding the positions of the resonant peak and the band gap edges but there is a small discrepancy in the quality factor of the peak.

It is interesting therefore to see how the quality factor of the resonant peak is affected by the finiteness of the size of the system and the width of the source, in our FDTD calculations. This can be seen in Fig. 4(a) and 4(b) where the quality factor versus the lateral width of the sample ($`L_w`$) and versus the source width are shown respectively. In Fig. 4(a) the ratio of the source width over the lateral extent of the structure $`L_w`$ (defined in Fig. 1) is kept constant (equal to 0.5). It can be seen that as the lateral width approaches 20a, Q saturates at a value$``$170. In Fig. 4(b) we present results of the Q versus the source width for a system with a lateral width of $`L_w=40a`$ . The width $`L_w`$ is large enough for the quality factor to have reached it’s saturation value (which is related to the lateral size) so that one can see the dependence on the source width only. It is clearly seen that for a small source width the quality factor drops significantly. From the data shown in Fig. 4(a) and 4(b) it can therefore be argued that the finiteness of the source or the system’s latteral size or both give 2-D characteristics in what would essentially be a 1-D resonance for the infinite system and result in a reduced quality factor. By an infinite system we mean infinite latteral size and that the incoming waves have plane wave front. In our system when the latteral size is small the resonance is supported by a smaller number of scatterers and is forced to terminate at the edges. Also, when the system is large but the source is small the resonace occupies a smaller fraction of the cavity area and has magnitude that decreases as one moves away from the center of the resonance. The resonance can decay in time not only along the propagation direction y as the true 1-D resonance does but also along the latteral direction. In both cases (small latteral width or small source) the result would be a smaller quality factor. Now, for the case when both the system and the source is very large the resonance can be extended in the cavity region and be very close to the 1-D resonance. As seen in Fig. 4(b) for source width equal to 30a and $`L_w=40a`$ the quality factor is $`174`$. That value is close to $`180`$ which is the value of the quality factor as obtained from the FDTD that models the infinite system (Bloch boundary conditions) (dotted horizontal line in Fig. 4). The solid horizontal lines in Fig. 4(a) and 4(b) give the quality factor for the infinite system as obtained from the TMM. The respective value is $``$230. We have also looked at the dependence of the quality factor on the latteral width and kept the width of the source constant and equal to 2a. What we found was that the quality factor increases as the latteral width increases, reaches a maximum value , then startes to decrease and eventually saturates. This is the result of a combined effect. As described previously when the latteral width is small this causes the resonant peak to terminate at the edges. As the latteral width increases the resonance occupies the maximum area it can –according to the specific source width. Further increasing the latteral width will not cause the resonance to occupy more area of the cavity along the latteral direction x. It only allows it to decay in time along the x direction as well and therefore the quality factor starts to lessen. Another component that can cause the quality factor to decrease with decreasing source width is the angle span of the source that increases with decreasing source width. This would cause the peak to shift and broaden (towards the higher frequencies). We have numerically checked that the magnitude for the off normal components of the incoming wave is relatively small. Also, we have not observed significant shift in frequency with decreasing source width in the FDTD results. We have observed though that the height of the transmission peak decreases as the source width decreases. It approaches the value of one as the source width approaches the lateral size of the system (for the large system $`L_w=40a`$). It becomes one for the infinite system modeled either through the TMM or the FDTD with Bloch boundary conditions in the lateral direction of the crystal.

Fig. 5(a) shows the transmission calculated with the TMM method versus the reduced frequency u=a/$`\lambda `$ for a 1-D cavity of width $`L_c`$/a=0.985 and $`N_c`$=4, for different imaginary parts in the air dielectric constant. Introducing imaginary parts with values of 0.02-0.05 in the air dielectric constant in the TMM system, has been suggested as a mechanism to model out of plane losses in . This leads as well to a reduced Q factor value. $`N_c`$ is the number of rows of air cylinders in each side of the 1-D dielectric cavity. As expected the transmission peak decreases with Im $`ϵ`$ while the width of the transmission increases with Im $`ϵ`$. We have also examined how the size of the system along the propagation direction affects the characteristics of the resonance. In Fig. 5(b) the dimesionless linewidth $`\mathrm{\Delta }`$/a \- which is defined as $`\delta \lambda `$/a (where $`\delta \lambda `$ is the width in wavelength of the resonance at the half maximum value) - is plotted versus $`N_c`$ (the number of rows of cylinders on each side of the dielectric cavity) in semi-logarithmic scale. As mentioned previously $`N_c`$ is related with the thickness of the system $`L_{th}`$ through the relation $`L_{th}=2N_cb+L_c`$ ($`L_{th}`$ is shown in Fig. 1). A linear relation in the semi-log scale between $`\mathrm{\Delta }`$/a with $`N_c`$ is obtained, meaning that the quality factor Q, which is inversely propotional to $`\mathrm{\Delta }`$/a would increase exponentially with $`N_c`$. That is consistent with results for other defect cavities . In Sigalas et. al. , where an air defect in a dielectric array is studied, it was shown that when absorption is introduced in the dielectric the quality factor saturates and does not increase any further as the size of the system increases. Such saturation with the introduction of an imaginary part in the dielectric of the holes is observed in our case also, as can be seen in Fig. 5(b). We also observe (as in ) that the saturation value of the linewidth $`\mathrm{\Delta }`$ is smaller the larger the imaginary part in the dielectric.

In Fig. 6 we present the results achieved with the TMM of the transmision height versus $`N_c`$ for 1-D cavity of width $`L_c`$/a=0.985 for three different values of Im $`ϵ`$ of the air dielectric. Notice that for Im $`ϵ`$=0 the height is always one, i.e perfect transmision. However for nonzero Im $`ϵ`$, the transmission height drops as $`N_c`$ increases. Two experimental values are shown for almost the same width of $`L_c`$/a as the one in the calculations. Notice that this suggests that the Im$`ϵ`$ that can fit the experimental data  can have a value close to 0.04.

For the ideal case with no losses present the quality factor was calculated by the TMM for various cavity widths $`L_c`$ that led to peaks that span most of the bandgap. In Fig. 7 the quality factor is plotted versus the reduced frequency of the peak. (Peaks too close to the edges of the bandgap were not used because the determination of their Q would be vague as stated earlier in this paper). For every frequency in the plot there are three peaks that are characterized by different orders and correspond to different $`L_c`$ values (see Fig. 2) . So three different curves are recovered when grouping the peaks according to order. Higher order peaks are characterized by a higher value of the quality factor. For every order, the quality factor is maximum for a frequency close to the center of the band gap and reduces as the frequency approaches either edge of the band gap. The calculations above were performed with a grid lattice dividing the lattice constant a into ten intervals. We have investigated the dependence of the quality factor with the accuracy of the TMM given by the number N of intervals the lattice constant is divided to, i.e. N is the number of grid points. It was done for three different peaks. One close to the band gap center and the others close to the lower and higher frequency band edge respectively. It is seen from Fig. 8(a) that the quality factor increases with increasing accuracy and eventually saturates at a value. The same holds for the quality factor calculated from the FDTD method. This is shown on Fig. 8(b) for cavity spacer value Lc/b=1. The FDTD data of Fig. 8(b) are obtained with Bloch boundary conditions applied in the lateral (x) direction. This way one has an infinite source with a plane wave front as in the TMM. The saturation values of the Q-factor for the TMM and FDTD are $``$230 and $``$250 respectively. This difference might be due the different way of calculating the transmission coefficient in the two methods as was mentioned before in Section III.

IV. Conclusions

We have studied the properties of 1-D dielectric cavity structure in a 2-D hexagonal array with the transfer matrix method and the finite difference time domain method. Both methods agree in the position of the defect peak and both yield good agreement (for the peak’s position) with the experiment  as well. The quality factor though shows sensitivity to a lot of parameters such as the size of the system (both lateral and along the propagation direction), the type of the incoming EM fields and losses out of the plane of periodicity. In that context the two numerical methods (FDTD and TMM) were compared.

Acknowledgments

The authors would like to thank Mario Agio for providing the FDTD code and helpful discussions. We would also like to acknowledge H. Benisty, M. Rattier, and M. Sigalas for helpful discussions. Ames Laboratory is operated by the U.S Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work is supported by the IST project PCIC of the European Union and the NSF grant INT-0001236.


format: text

text:
# Redistribution of Spectral Weight in Spin-1/2-Doped Haldane Chains

\[

## Abstract

We study the evolution of the dynamical spin structure factor in a spin-1 antiferromagnetic Heisenberg chain which is randomly doped with spin-1/2 moments. Using stochastic series expansion Quantum Monte Carlo simulations combined with the Maximum Entropy method, we monitor the crossover from the spin-1 chain, dominated by a single resonance at the Haldane energy, to the spin-1/2 chain with a continuum of spinon states which diverges at zero-frequency. Upon increasing the doping level, spectral weight is rapidly transferred from the Haldane peak to lower energies. If the exchange couplings between the spin-1/2 substituents and the spin-1 sites of the host are sufficiently small, finite-frequency bound states are observed below the spin gap for small doping concentrations.

\]

It is well established that antiferromagnetic (AF) spin-1 Heisenberg chains have a low-energy spin gap $`\mathrm{\Delta }0.41J`$, leading to properties significantly different from quasi-long-range-ordered spin-1/2 AF Heisenberg chains.  The effects of non-magnetic impurities on such low-dimensional quantum spin liquids have been the focus of considerable recent attention. In particular, it was shown that open ends in segmented Haldane chains induce low-energy states below the spin gap, and that vacancies in AF Heisenberg ladders liberate free spin degrees of freedom in their immediate vicinity. At sufficiently low energies, RKKY-like interactions among these “pruned” spins lead to the formation of networks described by effective models of random spins interacting via random exchange couplings. The effects of substitutions with non-vanishing magnetic moments tend to be more complex, depending on the spins of the dopants as well as on their exchange couplings with the host material.  We address this issue by studying the evolution of the dynamical spin structure factor of spin-1 Heisenberg chains when doped with static spin-1/2 moments, and by monitoring the crossover from a completely gapped Haldane spectrum to that of a spin-1/2 chain upon complete substitution.

There have been several experimental and theoretical studies of the effects of spin-1/2 impurities in the spin-1 chain. For example, the substitution of $`\mathrm{Cu}^{2+}`$ for $`\mathrm{Ni}^{2+}`$ in the compound $`\mathrm{Ni}(\mathrm{C}_2\mathrm{H}_8\mathrm{N}_2)\mathrm{NO}_2(\mathrm{ClO}_4)`$ (NENP) introduces spin-1/2 impurities on the original Ni sites, which have a weak ferromagnetic coupling to the $`S=1`$ chain-end spins. On the other hand, substituting Ca2+ for Y3+ in $`\mathrm{Y}_2\mathrm{BaNiO}_5`$ introduces carriers in the chains which become localized on the oxygen sites between the $`\mathrm{Ni}^{2+}`$ ions. This breaks the superexchange path and replaces it by a direct exchange coupling to the spin-1/2 impurities. In this latter compound, inelastic neutron scattering experiments found a substantial increase of the spectral function below the Haldane gap, indicating the creation of states below the energy of the spin gap.

Numerical diagonalization and density matrix renormalization group (DMRG) techniques were used to study the emergence of such impurity-induced low-energy states in the Heisenberg model. The main conclusion of these studies is that localized excited states exist below the Haldane gap as long as the exchange couplings $`J^{}`$ between the spin-1/2 and spin-1 sites are small compared to the couplings $`J`$ between the original spin-1 sites $`J`$, namely if $`J^{}/J0.45`$. Other works have considered mobile holes.  In addition to enhanced low-energy spectral weight in the doped compound, these studies also observe low-temperature spin-glass behavior, due to impurity-induced frustrating ferromagnetic inter-chain interactions.

Here, we concentrate on the static case discussed in Refs.  and apply stochastic series expansion Monte Carlo (SSE QMC) simulations in combination with the maximum entropy method to obtain the spectral properties of Haldane chains randomly doped with spin-1/2 sites. This approach allows us to study the evolution of the dynamical spin structure factor upon increasing the spin-1/2 concentration far beyond the limit of low impurity concentration. The doped spin-1 AF Heisenberg chain pictured in Fig. 1 is modeled by a mixed-spin Heisenberg model,

$$H=\underset{i=1}{\overset{N}{}}J_i𝐒_i𝐒_{i+1},$$
(1)
where $`𝐒_i`$ denotes a spin-1 or spin-1/2 moment. In general, three inequivalent types of bonds exist in a mixed-spin chain: (i) bonds between two spin-1 sites, (ii) bonds between two spin-1/2 sites, and (iii) bonds between spin-1 and spin-1/2 sites. In the mixed spin model of Eq. 1, we use $`J`$ as the original coupling of the Haldane chain with $`J_i=J`$ if $`S_i`$ and $`S_{i+1}`$ are both spin-1, and otherwise $`J_i=J^{}`$, such that the couplings between spin-1 and spin-1/2 sites are $`J^{}`$. Since we are interested in the generic properties of this problem, in our simulations we concentrate on representative parameter sets.

The significance of $`J^{}`$ has been emphasized in Ref.  for the case of isolated spin-1/2 impurities. Consider first the two extreme limits $`J^{}/J=0`$ and $`J^{}/J=\mathrm{}`$. If the spin-1/2 dopants are only weakly coupled to the spin-1 moments of the host, an effective Hamiltonian for the three-site cluster formed by the spin-1 chain-end spins and the impurity spin can be derived, correctly describing the emergence of bound states below the Haldane gap of the spin-1 chain. In the strong coupling limit, the coupling to the remainder of the chain destabilizes the local excitations. A critical coupling ratio $`(J^{}/J)_{crit}=0.45`$ separates these two regimes.  Therefore, in the following we discuss separately the cases of weak coupling, which we illustrate with the parameter choice $`J^{}=0.1J`$, and strong coupling, for which we take $`J^{}=J`$.

In order to observe the evolution of the dynamical spectrum with doping concentration, we have performed SSE QMC simulations and measured the imaginary-time spin-spin correlation function

$$S(𝐤,\tau )=\frac{1}{N}\underset{i,j}{}\mathrm{exp}[i𝐤(𝐫_i𝐫_j)]S_i^z(\tau )S_j^z(0),$$
(2)
considering the dominant AF wave vector $`𝐤=\pi `$ by using the Monte Carlo estimator derived in Ref. . We studied chains of up to 500 sites for up to 500 realizations of the random distribution of spin-1/2 substituents down to temperatures $`T0.01J`$. The maximum entropy analytical continuation procedure was used to obtain the spectral function $`S(\pi ,\omega )`$ from the impurity-averaged data for the imaginary-time correlation functions.

Let us first consider the case $`J^{}/J=0.1`$, where the emergence of localized states below the spin gap is expected from DMRG calculations in the single impurity problem. In Fig. 2 we show the evolution of the spectral function $`S(\pi ,\omega )`$ upon increasing the concentration $`\rho `$ of spin-1/2 dopants. For the pure spin-1 chain, $`\rho =0`$, a dominant resonance is observed at the Haldane gap energy $`\omega =\mathrm{\Delta }=0.41J`$, consistent with earlier results using similar methods. Upon introducing spin-1/2 sites, the spectral weight of the Haldane peak is strongly suppressed and its width broadens. In addition, a significant amount of spectral weight appears at low but finite energies, indicating the emergence of impurity-induced states below the gap edge. In fact, the shape of the spectral density for $`\rho =0.01`$ is consistent with Ref. , i.e. bound states with energies of approximately $`\alpha J^{}/2`$ and $`\alpha 3J^{}/2`$ appear in the low-energy spectrum, where $`\alpha =1.064`$ is a renormalization factor for the effective spin-1/2 chain-end spins $`𝐒_{L,R}`$ next to the impurity site.  These localized states are the lowest excited states of the effective Hamiltonian $`H_{eff}=\alpha J^{}(𝐒_L𝐒^{}+𝐒^{}𝐒_R)`$, where $`𝐒^{}`$ denotes the impurity spin. The maximum entropy method does not allow us to resolve the detailed structure of the spectrum, but clearly indicates the increase of spectral weight at this energy scale with increasing impurity concentration. Upon a further increase of the spin-1/2 concentration beyond $`\rho 0.05`$, a zero-frequency peak emerges out of these low-energy contributions. Simultaneously, the Haldane peak appears to move towards larger frequencies . The spectrum at larger impurity concentrations ($`\rho >0.2`$) exhibits features similar to those of the pure spin-1/2 chain, i.e., a divergence at $`\omega =0`$ and no remaining peak structure at the energy of the Haldane gap.

Turning to the isotropic case, $`J^{}=J`$, DMRG calculations for single spin-1/2 impurities do not exhibit any low-energy states below the Haldane gap. Instead, an extended state with an energy of the order of $`\mathrm{\Delta }`$ is reported. We performed SSE QMC simulations for the isotropic case and obtained the dynamical structure factors shown in Fig. 3. Concentrating first on the low-energy properties of this spectral function, we can confirm that in the limit of small impurity concentrations $`\rho <0.06`$ no indications of low-lying impurity states are observed, consistent with the DMRG results. Upon further increasing the impurity concentration, a peak appears at zero energy, which develops into the zero-frequency singularity characteristic of the gapless spinon continuum of states in the pure spin-1/2 Heisenberg chain, i.e. $`S(\pi ,\omega )\sqrt{\mathrm{ln}\omega }/\omega `$. Analogous to the discussion of the weak-coupling limit, spectral weight is shifted from the Haldane peak to lower energies. However, in the isotropic case no bound states are observed at finite frequencies below the spin gap. Furthermore, the depletion of the Haldane resonance is less rapid than in the weak-coupling limit, indicating that the presence of weaker bonds has a more dramatic effect on the redistribution of spectral weight than the presence of spin-1/2 substituents.

In order to quantify this statement, the dependence of the spectral weight $`I`$ at the Haldane energy is plotted in Fig. 4 as a function of the doping concentration. Here, a local Gaussian fit of the spectral function centered around the peak position was used to estimate the spectral weight from the height and the width of the approximant. It is found that in the isotropic case remnants of the Haldane peak can be resolved at much higher spin-1/2 concentrations than in the limit of weak coupling. For example, at $`\rho =0.5`$, the Haldane resonance is still visible in Fig. 3 ($`J^{}=J`$), whereas in Fig. 2 ($`J^{}=0.1J`$) this feature cannot be resolved beyond $`\rho =0.2`$. This critical value for $`J^{}=0.1J`$ is obtained by a linear extrapolation of $`I`$ (dashed line shown in Fig. 4). In the case of isotropic couplings the spectral weight of the Haldane peak initially shows a rapid decrease, similar to the weak-coupling limit. However, beyond a doping concentration $`\rho 0.06`$ at which zero-frequency singularities emerge in the dynamical structure factor (Fig. 3), $`I`$ decays at a significantly slower rate. It is obviously more difficult to deplete the extended states at the energy of the Haldane gap for isotropic couplings than dissolving the sub-gap bound states in the weak-coupling limit.  A linear extrapolation of $`I`$ (dot-dashed line in Fig. 4) suggests a critical doping concentration $`\rho 0.8`$ for isotropic couplings. Beyond this concentration the dynamical structure factor resembles that of the pure spin-1/2 AF Heisenberg chain.

In conclusion, we have used a combination of stochastic series expansion Quantum Monte Carlo simulations and Maximum Entropy analytical continuation techniques to study the effects of random spin-1/2 substitutions on the dynamical response function of spin-1 AF Heisenberg chains. In the limit of small doping concentrations, the dynamical spin structure factor is found to agree very well with results of earlier DMRG and exact diagonalization studies. If the AF couplings between the spin-1/2 dopants and the spin-1 sites of the host are weak, bound states are formed below the Haldane gap $`\mathrm{\Delta }0.41J`$. On the other hand, in the case of isotropic coupling constants there are extended states at the energy of the spin gap. Upon increasing the doping concentration, a proliferation of low-energy states is observed at the expense of the spectral weight of the dominant resonance at $`\omega =\mathrm{\Delta }`$ of the pure spin-1 AF Heisenberg chain. In the limit of weak coupling, the spectral weight of the Haldane peak disappears rapidly. However, for isotropic couplings the transfer of spectral weight to lower energies is more gradual. Here, the gapless spinon continuum spectrum of the pure spin-1/2 AF Heisenberg chain is recovered beyond a critical doping concentration $`\rho 0.8`$.

We thank B. Normand for useful discussions, and acknowledge financial support by the Petroleum Research Fund, Grant No. ACS-PRF 35972-G6 and by the National Science Foundation, Grant No. DMR 0089882.


text:
# Acknowledgments

## Acknowledgments

These notes are based both on recent independent work and on recent work in collaboration with Tsvi Piran and with Luisa Doplicher and Soonkeon Nam. Conversations with Luisa, Soonkeon and Tsvi, and with Roberto Aloisio, Aurelio Grillo and Jurek Kowalski-Glikman, are gratefully acknowledged.


text:
# The Blazar Main Sequence

## 1 Introduction

The Blazars are commonly perceived as a class of Active Galactic Nuclei marked out from the rest of their kin by a number of features: strong, high frequency radio emission from compact cores, often with superluminal expanding components (see Jorstad et al. 2001); powerful $`\gamma `$ rays extending to energies $`h\nu 10^2`$ MeV and beyond (see Mukherjee 1999); rapid multi-wavelength variability (see Böttcher 2000); high and variable optical polarization (see Yuan et al. 2001). Comprehensive reviews and extensive references concerning these sources may be found in Urry & Padovani (1995), Padovani & Urry (2001, henceforth PU01).

The above features defining the Blazar class are widely explained in terms of the view originally proposed by Begelman, Blandford & Rees (1984). This holds the radiation from all these sources to be produced in a collimated, relativistic jet of particles with bulk Lorentz factor $`\mathrm{\Gamma }520`$; when observed at small angles of order $`\mathrm{\Gamma }^1`$ the jet produces the “blazing” effects that provide the class denomination. The intrinsic luminosity emitted by the jet over the entire solid angle, and total jet kinetic power are reviewed by Ghisellini (1999a), whose definitions and beaming factors we will adopt in the following.

Within the Blazar class two main subclasses are usually identified, namely, the BL Lac Objects and the Flat Spectrum Radioloud Quasars (FSRQs). The former differs from the latter subclass on several accounts, that we shall discuss in this paper: a) lack of emission lines with EW $`>\text{ }5`$ A, and lack of blue-UV bumps; b) total outputs $`L<\text{ }10^{46}`$ erg s-1, compared with FSRQs often exceeding this limit and in some cases approaching $`10^{48}`$ erg s-1 (see specifically Maraschi 2001); c) spectral energy distributions peaking at frequencies from optical to X-ray and in the range of $`10`$ GeV, compared with the FSRQ peaks at far IR to optical frequencies and around $`10^1`$ GeV; d) no signs of cosmological evolution, compared with the strong evolution clearly displayed by the FSRQs in common with the other quasars.

These distinctive features stand out of the many selections entangling this field, and even though intermediate objects are currently found (Sambruna et al. 2000, Perlman et al. 2001), the FSRQ/BL Lac strong bimodality calls for explanation; this is the scope of the present work, that expands some preliminary results reported by D’Elia & Cavaliere (2001).

We will base on the accreting black hole (BH) paradigm (see Begelman, Blandford & Rees 1984), in the framework of a 4-parameter classification similar to that discussed by Blandford (1990, 1993); this comprises the BH mass $`M`$, its angular momentum $`J`$, the accretion rate $`\dot{M}`$, and the viewing angle $`\theta `$ from the axis of the jet (if any). In the scheme we contemplate, radio quiet AGNs correspond to slowly rotating BHs with low (Seyfert galaxies) or full (quasars) $`\dot{M}`$ compared to the Eddington rates; the quasars require large BH masses $`M`$ of order $`10^9M_{}`$ to attain proportionally high Eddington luminosities $`L_E10^{47}M_9`$ erg s-1.

Radio loudness constitutes a debated issue. Radio loud AGNs are often observed to reside in very bright galaxies. On the one hand, these sources have been inferred to correlate with large BH masses; for instance, Laor (2000) stresses the case for radio loudness being a consequence of large $`M`$. On the other hand, the majority of the optically selected quasars including very bright ones are radio quiet (see for recent data Goldschmidt at al. 1999, Hamilton et al. 2001); in addition, Scarpa & Urry 2001 find that the hosts of quasars strong in the radio are no different from giant ellipticals in general, including hosts of radio quiet quasars. Finally, Ho (2001) from a comprehensive database finds no simple relation between radio loudness and $`M`$. Similar data from more limited observations spurred the view (see Blandford 1993) that the other BH parameter, namely the spin, should be involved as a necessary, intrinsic condition for directing and launching the jets associated with radio emissions. In the present paper we focus on this view.

In fact, we will consider the Blazars in the context of supermassive holes with $`M10^9M_{}`$ close to the maximal observed values for MDOs (see Gebhardt et al. 2000, Ferrarese & Merritt 2000), with high angular momentum $`JJ_{max}`$ close to the maximal $`J_{max}=GM^2/c`$, observed nearly pole on within angles $`\theta <\text{ }\mathrm{\Gamma }^1`$ (see Sambruna, Maraschi & Urry 1996). We shall concentrate on discussing how all Blazar properties can be arranged into a Blazar main sequence (BMS) in terms of levels of the accretion rate $`\dot{m}=\dot{M}c^2/L_E`$. Levels $`\dot{m}\eta ^1110`$ (the radiative efficiency being up to about 0.4 around a Kerr hole) mark the FSRQs, while values $`\dot{m}<\text{ }10^2`$ mark the BL Lacs.

Our plan is as follows. In §2 we will stress how the optical properties of the Blazars can be understood in terms of different levels of $`\dot{m}`$. In §3 we show that high rates $`\dot{m}1`$ can produce the huge thermal and non-thermal outputs featured by many FSRQs, while conditions where $`\dot{m}1`$ lead to the moderate, but mainly non-thermal outputs of the BL Lacs. In §4 we will argue that for lower values of $`\dot{m}`$ along the BMS the typical particle and photon energies are expected to increase; a secondary physical parameter relevant at low $`\dot{m}`$ causes considerable variance within the BL Lac subclass, adding to the effects from the orientation. In §5 we show how in the same framework we understand the widely different evolutionary properties observed in FSRQs and in BL Lacs. Finally, in §6 we summarize and discuss our conclusions.

## 2 Optical Properties

In the standard view (see Frank, King & Raine 1992; Peterson 1997) all thermal outputs are produced in, or around the accretion disk, with only weak anisotropies and with power given by $`L_{th}\eta \dot{m}L_E`$.

Here we take up this view, and just stress how from the broad optical emission lines and the blue-UV bump in particular we derive a first indication toward a BMS in terms of levels of $`\dot{m}`$. The bump is produced by power reprocessed in a superposition of black body emissions from the inner rings of the accretion disk, or in a hot corona (see Sun & Malkan 1989, Siemoginowska et al. 1995); the lines are produced by the high UV-soft X-ray continuum reprocessed by distant gas “clouds” (see Netzer 1990). In both cases the density of the reprocessing material is expected to scale up with increasing $`\dot{m}`$.

The FSRQs share with other quasars conspicuous bumps and prominent broad emission lines; the strength of these features are consistent with high accretion rates $`\dot{m}1`$.

On the other hand, the weakness or absence of such features in the BL Lac spectra can be understood in terms of $`\dot{m}<\text{ }10^2`$. In fact, at such levels of $`\dot{m}`$ not only a weaker ionizing/exciting continuuum is emitted by the disk, but also lower gas densities are expected in the disk and its environment out to pc and larger scales; these conditions concur to account for the weak or absent emission lines.

The other effect that in the BL Lacs acts to swamp any residual bump and to further reduce the EW of intrinsically weak emission lines is constituted by the blazing continuum from the jet when observed at a small angle, see Vagnetti, Giallongo & Cavaliere (1991); Georganopoulos (2000).

We recall that low values of $`\dot{m}`$ related to radiosources had been proposed by Begelman et al. (1984). Specifically, the different power levels and optical spectra of the FRII and FRI radiosources have been previously related to high and low values of $`\dot{m}`$, respectively, by Baum, Zirbel, & O’Dea 1995; this point will be further discussed in §6.

## 3 The Blazar Power

As to the jets, we assume them to be powered by variants of the mechanism originally proposed by Blandford & Znajek (1977, BZ) for direct extraction of rotational energy from a Kerr hole via the Poynting-like flux associated with the surrounding magnetosphere. The BZ power scales as

$$L_KB_h^2r_h^2,$$
$`(1)`$
in terms of the magnetic field $`B_h`$ threading the hole horizon at $`r_h`$.

Variants are necessary in view of the limitations to $`L_K`$ recently stressed by Modersky & Sikora (1996); Ghosh & Abramowicz (1997); Livio, Ogilvie & Pringle (1999). Such variants involve either high strengths of $`B_h`$, or the similar electrodynamic output of the disk from a larger radius $`r_d`$.

In the first variant, very strong $`B_h`$ as advocated by Meier (1999) are required to account for the huge outputs of some FSRQs. Fields up to $`B_h^2/8\pi \rho c^2`$ in the plunging orbit region have been argued by Krolik (1999); Armitage, Reynolds & Chiang (2000) and Paczynski (2000) discuss why such field values are unlikely in a thin disk. In thick disks the status of such enhanced fields is still unsettled.

The second variant takes up from recent discussions (Ghosh & Abramowicz 1997; Livio, Ogilvie & Pringle 1999) that have stressed the continuity of $`B_h`$ with the field $`B_d`$ at the inner rim of a standard $`\alpha `$-disk; in turn, $`B_d`$ is bounded by the maximum pressure, following $`B_d^2/4\pi <P_{max}`$. Depending on $`\dot{m}`$, the inner disk is dominated by gas pressure (GPD) or by radiation pressure (RPD); the pressure scales as $`P_{max}(\alpha M_9)^{9/10}\dot{m}_4^{4/5}`$ in the first, and as as $`P_{max}(\alpha M_9)^1`$ in the second regime. Correspondingly, in the first regime the maximum power extractable from the hole grows with $`\dot{m}`$, and in the jet frame the emission reads (see Modersky & Sikora 1996; Ghosh & Abramowicz 1997)

$$L_K=10^{44}M_9^{11/10}\dot{m}_4^{4/5}(J/J_{max})^2ergs^1(GPD)$$
$`(2)`$
in terms of the Kerr hole mass $`M=M_9\mathrm{\hspace{0.17em}10}^9M_{}`$; but for $`\dot{m}>\text{ }10^3`$ the power $`L_K`$ saturates to

$$L_K=\mathrm{2\hspace{0.17em}10}^{45}M_9(J/J_{max})^2ergs^1(RPD).$$
$`(3)`$
This is because the radius $`r_c`$ defining the RPD region grows with $`\dot{m}^{16/21}`$ (Novikov & Thorne 1973), and it turns out to exceed the radius $`r_{ms}`$ of the last stable orbit when $`\dot{m}>\text{ }10^3`$ holds; in the latter conditions an inner RPD region actually exists.

On the other hand, the electrodynamic power $`L_d`$ contributed by the disk and itself strongly anysotropic scales as

$$L_dB_d^2r_d^2;$$
$`(4)`$
here the effective disk radius $`r_d`$ is bounded at several times $`r_h`$ either by the value of $`r_c`$ or by the radial decline in the disk of the efficiency $`\eta GM\dot{M}/2r`$. We recall that the ratio $`L_d/L_K`$ had originally been estimated by Blandford & Znajek 1977 at values $`>\text{ }5`$; recent reappraisals (e.g., by Livio, Ogilvie & Pringle 1999) have evaluated it at $`L_d/L_K(B_d/B_h)^2(r_d/r_h)^{3/2}`$.

We draw the following implications concerning the BMS. Many BL Lac jets can live on accretion rates $`\dot{m}10^2`$, since the power levels $`L_K`$ given by eq. (3) are adequate even including the kinetic energy in the jets (Celotti, Padovani & Ghisellini 1997, Maraschi 2001). A larger but comparable contribution comes from the disk, so that $`L_KL_d>\text{ }L_{th}`$ holds. The disk is coupled to the rotating hole both dynamically and magnetically (see Livio, Ogilvie & Pringle 1999) through the magnetic field lines threading the hole and tethered in the disk.

On the other hand, many FSRQs feature total outputs $`L>10^{46}`$ erg s-1, with specific sources exceeding $`10^{47}`$erg s-1, see Tavecchio et al. (2000), and Maraschi (2001), who include a high proton contribution to the jet energy and adopt average values of $`\mathrm{\Gamma }10`$ (see Sikora 2001). Such luminosities require not only a saturated $`L_K`$ (see eq. 3), but also a dominant component up to $`L_d10^2L_K`$ from the disk corresponding to $`B_d`$ a few times larger than $`B_h`$ and $`r_d5r_h`$; this implies an extended region dominated by radiation pressure, hence conditions where $`\dot{m}1`$, consistent with the requirements from the optical properties discussed in §2.

We note that in FSRQ conditions both $`L_d`$ and the thermal emission $`L_{th}`$ are fed by the gravitational power supply of order $`GM\dot{M}/2r_d`$. The share between $`L_d`$ and $`L_{th}`$ depends on the power spectrum of the magnetic field inhomogeneities, that is, on the power share between two scales: the large scale, coherent vs. the small scale, turbulent component. A dominant share of the former, such as to yield a large $`L_d`$ and a large associated transfer of angular momentum outward, may require modifications of the standard $`\alpha `$-disk model as discussed by Salvati (1997).

The hole output $`L_K`$, even when it constitutes a minor component to the total, is likely to be important to provide on the jet axis a “high-velocity spine” instrumental for the outward jet propagation, see Livio (2000), Chiaberge et al. (2000).

## 4 The Blazar SEDs

Another feature of the Blazars is constituted by their extended spectral energy distribution (SED). This shows a first peak commonly interpreted as synchrotron radiation (see Urry & Padovani 1995). In the FSRQ spectra the peak frequencies $`\nu _{peak}`$ lie mainly in the range from the far IR to the optical band; in the BL Lacs the peaks are generally shifted to higher frequencies, and in some objects up to the hard X rays (Costamante et al. 2001).

An additional high energy component to the radiation (see PU01) is observed into the 10 GeV range for the FSRQs, and out to several TeVs for the BL Lacs 111Note that the $`\gamma `$-ray spectra at the sources may extend to higher energies than observed at Earth after the intervening absorption due to the IR background. If so, the following argument concerning coherent acceleration would be reinforced. This is likely to be inverse Compton emission produced by substantial numbers of GeV or of $`10^210^3`$ GeV electrons, respectively, scattering off external or synchrotron seed photons as discussed by Ghisellini (1999a).

The apparent inverse scaling of $`\nu _{peak}`$ with the radiative $`L`$ in moving from FSRQs to (quiescent) BL Lacs has been pointed out by Fossati et al. (1998), and interpreted as a cooling sequence due to the radiative energy losses being faster in the more powerful sources, up to controlling the maximal electron energies. Concerning the acceleration, one view considers weak electric fields over relatively large distances; these occur in the internal shock scenario where, however, the typical energies produced are of order $`2\mathrm{\Gamma }m_pc^2<20`$ GeV and fall short of the top energies required in some BL Lacs as discussed by Ghisellini (1999b).

So we are lead to consider as an addition or as an alternative particularly relevant to BL Lacs, the electron acceleration due to higher electric fields acting over shorter distances. We shall see that this leads to an acceleration sequence also linked with the BMS; this is because in the BL Lacs lower levels of $`\dot{m}`$ allow for higher electron energies and so produce a SEDs peaked at higher frequencies. Meanwhile, a 2nd intrinsic parameter produces the specific scatter of $`\nu _{peak}`$ within this subclass.

Such high $`E`$ fields naturally arise around the BZ magnetosphere in association with considerable energy transport via “Pointing flux” along the jet, with the magnetic field decreasing outwards rather slowly, like $`B=B_d(r_{ms}/r)^{(1+p)}`$ with $`p<1`$; typically $`p=1/4`$, in a similar manner to Blandford & Payne (1982). The force-free condition $`EB=0`$ governing the BZ magnetosphere has to break down at the flow boundaries, and to give way to electric fields $`E<\text{ }B`$; this is likely to occur time-dependently and inhomogeneously within sheets or filaments located at average distances $`r10^{16}÷10^{18}`$ cm. But such $`E`$ fields are electrodynamically screened out beyond distances exceeding a few times $`d=c/\omega _p(\gamma /n)^{1/2}`$, where $`\omega _p=(4\pi e^2n/m\gamma )^{1/2}`$ is the plasma frequency of the screening particles with Lorentz factor $`\gamma `$ and density $`n`$.

Lower bounds to the densities may be estimated from the powers emitted in the form of synchrotron or inverse Compton radiation, that is, $`L\gamma ^2UR^3n`$ where $`U`$ is the total energy density in quiescent conditions. We obtain $`n1`$ cm-3 for BL Lacs with $`L10^{44}`$ erg s-1, and $`n10^3`$ cm-3 for FSRQs with $`L10^{47}`$ erg s-1. Correspondingly, we expect many electrons to attain high energies, up to

$$\gamma _{max}<\text{ }eBd/m_ec^210^8B_4M_9^{1.25}d_{10}r_{17}^{1.25},$$
$`(5)`$
where we have used $`d=10^{10}d_{10}`$ cm, and $`r=10^{17}r_{17}`$ cm; we have also expressed as $`B_d=10^4B_4`$ G the field values at $`r_{ms}\mathrm{1.5\hspace{0.17em}10}^{14}M_9`$ cm, which decrease into the emitting region like $`B(r)=B_d(r_{ms}/r)^{1.25}`$ as recalled above.

In eq. (5) there is room for variations of $`B`$ and $`n`$; the former may be down by a factor 10 or the latter up by factors $`10^2`$ (with $`M`$ again set at $`10^9M_{}`$), and eq. (5) still yields $`\gamma 10^7`$. This is enough to produce inverse Compton photons up to $`10`$ TeV, the highest energy photons observed to now from the BL Lacs. In FSRQs densities $`10^3`$ larger give electron energies some 30 times smaller, consistently with the parameters derived by Böttcher & Dermer 2001 from advanced spectral modeling. Given $`\gamma _{max}`$, we expect less energetic but more numerous electrons to arise naturally when the particles escape before full acceleration or $`E<B`$ holds, and when the inhomogeneities of $`B`$ are considered; the latter varies outwards following $`Br^{1.25}`$ as said, and this condition by itself produces an energy range by factors close to $`10^3`$ within a distance range $`r10^{17\pm 1}`$ cm.

The electron energies given by eq. (5) are proportional to $`dn^{1/2}`$, and so scale like

$$\gamma _{max}(\gamma ^2UR^3)^{1/2}L^{1/2}.$$
$`(6)`$
This results in $`\gamma _{max}L^{1/2}`$ if we adopt the empirical relation $`\gamma ^2U`$ const reported by Ghisellini (1999a) for the FSRQs. On the other hand, the same author finds that the $`\gamma `$ values for the BL Lacs exceed this relation, but may be included on modifying it to read $`\gamma ^2UU^{0.2}`$; from this we obtain $`\gamma _{max}L^{3/5}`$. Thus we expect the peaks of the non-thermal emissions, and more clearly the syncrotron $`\nu _{peak}\gamma _{max}^2`$, to anticorrelate strongly with the disk luminosity following $`\nu _{peak}L^{6/5}`$ or steeper, when comparing the FSRQ with the BL Lac subclass.

In closer detail, the emitting particles are likely to constitute the high energy subset of the total particle population; but similar results obtain if all electrons are described as in Ghisellini (1999a) by a power-law energy distribution $`n(\gamma )\gamma ^a`$ from $`\gamma _{min}10^3`$; with $`a2.5`$ we obtain $`\nu _{peak}\gamma _{max}^2L^{4/3}`$, while a steeper relation $`\nu _{peak}L^2`$ holds with $`a2`$. Concerning the full SED we refer to the complete modeling of the emitting regions by Böttcher & Dermer 2001, who treat quantitatively the spectral shapes radiated by synchrotron and inverse Compton emissions in the jets, considering also the cooling limit to the electron energies provided by the high photon densities in FSRQs; they evaluate specific parameters agreeing with the orders of magnitude we estimate.

The other outcome from our view is that within the BL Lac subclass we also expect the second parameter $`L_K/L_d`$ to affect the scaling laws provided by eq. (6). This is because the total radiative luminosity of the BL Lacs $`L=L_d+L_K`$ comprises a disk and a Kerr hole component with comparable powers. Of these, only $`L_dGM\dot{M}/2d_{ms}`$ is directly related to $`\dot{m}`$ and to $`n`$ while $`L_K`$, when considerable, tends to be independent of them as shown by eq. (3). Thus at a given $`L_d`$ (hence at given $`\dot{M}`$ and $`n`$) higher total luminosities obtain in BL Lacs where $`L_K`$ is important. In other words, $`\nu _{peak}`$ scales inversely with the component $`L_d`$ only, and it is bound to show scatter specific to the BL Lacs when one tries to correlate it with the total luminosity. As a consequence, we do not expect the scaling to extend unblurred from the so called LBLs to the HBLs (low and high energy peaked BL Lacs, respectively, see Padovani & Giommi 1995). In fact, considerable scatter of $`\nu _{peak}`$ at given $`L`$ is increasingly found within the BL Lac subclass, see Giommi et al. (2001).

On the other hand, for the FSRQs $`LL_d`$ holds, that is, $`L_K`$ is irrelevant as to the total power, and so is its contribution to the scatter.

The specific BL Lac scatter discussed here adds to other components common to all Blazars due to different values of $`B,M`$, and viewing angles (see Sambruna, Maraschi & Urry 1996, Georganopoulos 2000), and to variations of $`\nu _{peak}`$ and $`L`$ arising during flares.

To this point we have shown how all spectral properties of FSRQs and BL Lacs are expected to consistently vary along the BMS. In the next section we discuss why we expect the Blazar evolutions to be also related to different levels of $`\dot{m}`$.

## 5 The Blazar Evolutions

Strong cosmological evolution is closely shared by the FSRQs with the rest of the quasars (see Wall & Jackson 1997, Goldschmidt et al. 1999), and shows up, e.g., in their steep number counts; the BL Lacs, instead, show no signs of a similar behavior (Giommi, Menna & Padovani 1999; Padovani 2001).

Here we argue that such an evolutionary peculiarity is related to the following two independent circumstances: i) At levels $`\dot{m}<\text{ }10^2`$ all BL Lacs shine for several Gyrs with only a mild decrease in luminosity (a very slow “luminosity evolution”, LE), and this implies counts close to Euclidean. ii) A genetic link is conceivable from FSRQs to BL Lacs, that is, some of the former may switch from a short-lived regime of high $`\dot{m}`$ to the long-lived regime of low $`\dot{m}`$; so the BL Lacs undergo also some negative “density evolution”(negative DE), and their counts are further flattened. Since point i) has been discussed by Cavaliere & Malquori 1999, we just recall their results below, and concentrate on point ii).

The evolution of the bright optically selected (and mainly radio quiet) quasars is strong, see Boyle et al. (2000), and so is in X-rays, see Della Ceca et al. (1994); Miyaji, Hasinger & Schmidt (2000). For the bright sources, this behavior is widely traced back to the exhaustion in the host galaxy of the circumnuclear gas stockpile usable for further accretion onto the central black hole. Such an exhaustion is caused by previous accretion episodes and by ongoing star formation (Cattaneo, Haehnelt & Rees 1999; Kauffmann & Haehnelt 2000; Cavaliere & Vittorini 2000, CV00); the ensuing decrease of the average accretion rate $`\dot{m}`$ is fast, and causes the average bolometric luminosities to scale down rapidly for $`z<2.5`$ (a strong LE). The radio loud quasars share this evolutionary behavior, constituting a fraction (not yet quantitatively understood, see §1, and the discussions by Moderski, Sikora & Lasota 1998 and by Livio 2000) of order 10% of the total; the FSRQs in particular are widely held to be the same sources seen pole on.

But as the rates $`\dot{m}`$ decrease we also expect a considerable fraction of the underlying hole-disk systems to eventually switch from being overwhelmingly fueled by accretion to being partly fed by the hole rotational supply $`E_K`$; the latter had been stockpiled during the phases of accretion of coherent angular momentum $`J`$ along with the mass (Bardeen 1970; Moderski, Sikora & Lasota 1998). In other words, we expect some radio loud sources to switch from the FSRQ to the BL Lac mode. The moderate outputs typical of this mode can be sustained for several Gyrs by the dynamically and magnetically coupled system constituted by Kerr hole and inner disk; so the decrease of BL Lac luminosities is expected to be slow (Cavaliere & Malquori 1999), with time scales in the range estimated to a first approximation from $`\tau _LE_K/L_K5÷10`$ Gyr, rather uniform within the population.

In closer detail, CV00 (see also refs. therein) trace the accretion episodes feeding the quasars back to interactions of the host galaxies with their companions in a group. These events destabilize the host gas, and trigger accretion over times of order $`\tau `$ a few $`10^1`$ Gyr that comprise the duration of a close interaction and a few galactic dynamical times for the inflow to develop and subside; some $`35`$ repetitions are expected per host galaxy since $`z2.5`$. The frequency of such interaction episodes diminishes in time, and gives rise to a positive, weak DE of the quasar population over scales $`\tau _D`$ several Gyrs. Meanwhile, the efficiency of such episodes drops as anticipated above, due to the exhaustion of the host gas; the latter is halved on the scale $`\tau _L3`$ Gyr. The result is that the average luminosities drop over the same scale, in other words, a strong LE is produced in the quasar population.

On this basis CV00 compute the luminosity function (LF) for the optically selected objects at $`z<2.5`$ to read

$$N_F(L,z)\frac{\tau N_G(z)}{\tau _r(z)L_b(z)}\frac{1}{l^{1.2}+l^{3.2}}.$$
$`(7)`$
Here $`N_G(z)`$ is the space density of groups of galaxies, $`\tau _r^1(z)(1+z)^{3/2}`$ is the average interaction frequency for a host galaxy in a group, and $`lL/L_b(z)`$ is normalized using the break luminosity that scales as $`L_b(z)(1+z)^3`$ in the critical universe; we use this simple scaling for numerical evaluations, since CV00 show that the observables do not vary much in the $`\mathrm{\Omega }_o+\mathrm{\Omega }_\lambda =1`$ cosmology. For the FSRQs we will use here eq. (7) with the proportionality factor of order $`\mathrm{5\hspace{0.33em}10}^2`$.

Our main point here is that the episodes of powerful FSRQ activity feeding on high $`\dot{m}`$ will last some $`10^1`$ Gyr and die out over similar scales; the sources may be no longer reactivated after a “last interaction”. This leaves behind a rapidly spinning hole (and a relatively quiescent though bright galaxy, see Scarpa 2001) in about 1/2 of the objects. Thus the DE timescale $`\tau _D`$ for the deaths of the bright FSRQs also constitutes the timescale for BL Lac births; we evaluate it from eq. (7) integrated over $`LL_b`$, to obtain

$$\tau _D^1=\frac{d}{dt}ln[_1^{\mathrm{}}N(l,t)𝑑l]=\frac{d}{dt}ln[N_G(t)/\tau _r(t)],$$
$`(8)`$
with the second equality holding because only the time-dependent amplitude in eq. (7) is relevant here. We obtain values of $`\tau _D`$ between 7 and 5 Gyr, depending on the form for small $`z`$ of $`N_G(z)`$ which may scale as $`(1+z)`$ on assuming a Press & Schechter (1974) mass distribution, or as $`(1+z)^2`$ following Lacey & Cole (1993).

One sign of evolution is provided by the integrated source counts, which may be evaluated from the relation (see Cavaliere & Maccacaro 1990)

$$N(>S)S^{3/2}[1C(S_0/S)^{1/2}+0(S^1)].$$
$`(9)`$
This simple relation is limited to high-intermediate fluxes. On the other hand, it has the advantage of exposing how the evolutionary properties of the two Blazar subclasses depend on the two population time scales $`\tau _D`$ and $`\tau _L`$ (the latter being long for the BL Lacs, and short for the FSRQs); these appear in the coefficient

$$C=3D_0l^2[2(1+\alpha )1/H_0\tau _D(\beta 1)/H_0\tau _L]/4R_Hl^{3/2}.$$
$`(10)`$
For the BL Lacs counted in the radio band a relevant value is $`\tau _L10`$ Gyr (marking the slow LE of the BL Lacs). Other quantities involved are: the flat spectral index $`\alpha =0.3`$ in the GHz range; the slope $`\beta 2.5`$ of the radio LF; the normalized moments $`l^n`$ of the LF; the typical distance to high flux BL Lacs $`D_0=(L_0/4\pi S_0)^{1/2}0.03R_H`$ in Hubble units.

On using only the slow LE on the scale $`\tau _L8`$ Gyr, eq. (10) with the Hubble constant $`H_0=65`$ km/s Mpc yields $`C=\mathrm{\hspace{0.17em}0.1}`$; then from eq. (9) the result is $`N(>F)F^{1.6}`$ as in Cavaliere & Malquori (1999).

But BL Lac births imply also negative DE to occur, described by the scale $`\tau _D\mathrm{\hspace{0.17em}6}`$ Gyr evaluated below eq. (8). When also this scale is introduced in eq. (10), the values of $`C`$ increases to $`C0.1`$. Now the result is $`N(>S)S^{1.5}`$ at high-intermediate fluxes, consistent with the data by Giommi et al. (1999), see fig. 1; at 5 GHz similarly flat counts are both expected (see Jackson & Wall 2001) and observed (see the preliminary data by Padovani 2001), while the FSRQ counts observed are quite steep (see Wall & Jackson 1997, Jackson & Wall 2001). We add that no evidence of evolution has been observed by Caccianiga et al. (2001) from the values of $`V_e/V_a`$ in a large sample of BL Lacs from the REX survey.

The count slopes we evaluate for the BL Lacs are in sharp contrast with the result from a similar evaluation for the FSRQs. The parameters appropriate for them are recalled above and read: $`\tau _D=+\mathrm{\hspace{0.17em}6}`$ Gyr, $`\tau _L=3`$ Gyr, while $`D_00.3R_H`$ holds and again $`\beta 2.5`$. Eq. (10) then gives $`C8`$, and in turn eq. (9) yields radio counts $`N(>S)`$ steep well above $`S^{1.5}`$ and consistent with the data, see fig. 2.

We thus see that also the evolutionary properties vary systematically along the BMS. Actually, the faint end of the LF is expected to be flattened by the varying beaming effects associated with a distribution of $`\mathrm{\Gamma }`$ values (Urry & Padovani 1995); this holds both for the FSRQs and the BL Lacs, and will contribute to the flattening of the counts at their respective faint ends, but it does not affect the intermediate range where sharply different slopes stand out.

## 6 Discussion and Conclusions

We have argued that Blazars constitute basically similar sources comprising a supermassive, nearly maximally spinning Kerr hole and its accretion disk; but large, systematic and correlated changes are bound to arise in all observational properties when the accretion rate levels vary from $`\dot{m}1`$ that mark the FSRQs to $`\dot{m}1`$ that mark the BL Lacs. Such systematic changes basically arise because lower values of $`\dot{m}`$ imply depleted densities in and around the disk; by the same token, these cause both higher energies of the emitting particles with less reprocessing of the emitted photons, and lower total luminosities over longer time scales leading to less evolution.

Our view (in short, BMS for Blazar main sequence) outlined here lead us to specifically predict the following correlated features, also summarized in Table 1.

a) Optical properties: in §2 we recall that high values of $`\dot{m}`$ yield strong thermal features, namely, emission lines and big blue bump. On the other hand, we stress that low values of $`\dot{m}`$ imply intrinsic weakness or drowning of all these features, especially in the presence of a strong beamed continuum.

b) Luminosities: in §3 we argue – in view of the limitations to the power directly extractable from a rotating hole – that disks with high values of $`\dot{m}`$ are required to produce outputs as strong as those observed in FSRQs. In BL Lacs, on the other hand, low values of $`\dot{m}`$ sustain intermediate or low disk luminosities, then comparable to those directly extractable from a Kerr hole after eq. (3). This yields a limiting power of about $`10^{46}`$ erg s-1; it is intrinsic to the present view that the top BL Lac outputs should not considerably exceed such values, as in fact observed (see Maraschi 2001).

c) SEDs: also the main features of the Blazar spectra can be predicted in this picture as shown in §4. Comparing BL Lacs with FSRQs, in the former sources we expect lower particle densities in the acceleration region, and less effectively screened electric fields. These produce higher electron energies up to some $`10^3`$ GeV, allowed on the other hand by the lower photon densities; the ensuing synchrotron emission will peak at higher frequencies. From b) and c) we derive the inverse correlation $`\nu _{peak}L_d^1`$ or somewhat steeper, see §4. But in BL Lacs specific and considerable scatter is to arise when one tries to correlate $`\nu _{peak}`$ with the total luminosity $`L`$, because this is made up by two comparable components, $`L_K`$ and $`L_d`$.

d) Evolution: comparing FSRQs with BL Lacs, much weaker evolution is expected for the latter, as argued in §5. This is because their output is lower and is partly extracted from the Kerr hole rotational energy, in conditions of low and long lasting accretion rates; so these sources are are long lived, with a slow population evolution of the LE type over scales $`\tau _L`$ of several Gyrs. Using only this scale to predict the source counts from the relations given in §5, the result is $`N(>F)F^{1.6}`$ at fluxes $`F>10^2`$ Jy. In addition, we expect a genetic link to also occur; some FSRQs may make a transition to the BL Lac mode after a “last interaction” of their host with companion galaxies. Such BL Lac births are equivalent to a negative, though weak DE on a scale $`\tau _D6`$ Gyr; this competes with the slow LE to yield counts flatter yet, namely $`N(>F)F^{1.5}`$ at intermediate fluxes, as shown fig. 1; equivalently, we expect values of $`V_e/V_a0.5`$, see Giommi, Menna & Padovani 1999. This scenario implies transitional objects to occur from FSRQs to BL Lacs, with $`\dot{m}<\text{ }10^1`$ and shining with intermediate luminosities for a few $`10^1`$ Gyr; in the optical band these sources will still show thermal features (emission lines and bumps), while their non-thermal SEDs will peak at higher $`\nu _{peak}`$ compared with canonical FSRQs. On the other hand, similar relations computed with the shorter time scales appropriate for the FSRQs yield much steeper counts, as shown in fig. 2.

We note that if the radiosources FRI and II constitute the parent populations of the BL Lacs and the FSRQs, respectively, then a genetic link between BL Lac births and FSRQ deaths implies a similar link between the related radiosurces. Such a link constitutes matter of current debate. We recall from §2 that power levels and optical spectra of the FRII and FRI radiosources have been related to high and low values of $`\dot{m}`$, respectively, by Baum, Zirbel & O’Dea 1995. We stress that overlapping properties and truly transitional objects are being found (see Bicknell 1995, and PU01 with references therein); in addition, similar evolutionary behaviors are observed for the radiosources and for the Blazars related to similar values of $`\dot{m}`$, that is to say strong evolution for the FRII but weak or even negative evolution for the FRI radiosources (Jackson & Wall 2001). So there is scope for discussing “grand unification” schemes of radiosources and Blazars, as pointed out by PU01.

We also note that our BMS agrees well with the parallel work by Böttcher & Dermer 2001; this is based on detailed modeling of the Blazar emitting region, and concerns the full SEDs evolving from FSRQs to BL Lacs as the accretion rates decrease and the black hole environments are depleted in gas, dust and reprocessed photons. In fact, the two papers turn out to be complementary, and concur toward the notion of a sequence or even a genetic link between the two subclasses.

Finally, if we carry the above BMS to its extreme, that is, to residual $`\dot{m}<\text{ }10^4`$, we are led to expect objects with very faint, if any, e.m. emission (further suppressed in the ADAF regime, see Di Matteo et al. 2000), along with nearly unscreened electric fields; these objects would be akin to the cosmic ray accelerators discussed by Boldt & Ghosh (1999).

Under the widespread if arguable assumption that in such conditions a magnetic field $`B10^210^3`$ G can still be held around the hole with the slow upward decrease recalled in §4, particles including protons would be accelerated over long distances, at radii that must exceed the force-free region’s (acceleration at small radii has been shown to be untenable by Levinson 2000). These particles – in fact, cosmic rays with ultra high energies – are predicted by the analogous of eq. (5) but in conditions of unscreened fields, that is, large $`d`$ due to very low $`n`$ (and to consistently large $`\gamma `$). Thus we recover the long recognized energy limit given by

$$_{max}eB_dr_{ms}(r_{ms}/r)^p/p10^{20}B_3M_9^{5/4}r_{16}^{1/4}eV.$$
$`(11)`$
Tens of these accelerators could lie within 50 Mpc. As noted by Boldt & Ghosh (1999), they would evade in the simplest way the GZK cutoff at UHE; as to accounting for the observed particle flux, they need to produce only $`L10^{42}`$ erg s-1 mostly in energetic particles. Intergalactic magnetic fields of nG strength would blur to near isotropy the geometrical memory of the actual sources for most such UHECRs.

Acknowledgements: We are indebted to P. Giommi, L. Maraschi, P. Padovani and M. Salvati for several helpful discussions. We thank our referee for useful and constructive comments. Partial grants from ASI and MIUR are acknowledged.

Table 1 - The Blazar Main Sequence summarized

## Figure Captions

Fig. 1. The BL Lac counts evaluated from eqs. (8) and (9) using $`\tau _D=\mathrm{\hspace{0.17em}7}`$ Gyr (dashed) and $`\tau _D=\mathrm{\hspace{0.17em}5}`$ Gyr (solid). The dotted line represents the Euclidean counts. The data points concern HBLs and are taken from Giommi et al. 1999, see their fig. 5.

Fig. 2. The FSRQ counts evaluated from eqs. (8) and (9) using $`\tau _D=+\mathrm{\hspace{0.17em}6}`$ Gyr and $`\tau _L=\mathrm{\hspace{0.17em}3}`$ Gyr (solid line). The open squares represent FSRQ counts from Wall & Jackson 1997.


text:
# Two symmetry problems in potential theory

## Remarks

###### Proposition 2

Suppose there exists a solution $`uC^2(\overline{\mathrm{\Omega }})`$ to following overdetermined boundary-value problem

$$\begin{array}{c}\mathrm{\Delta }u=r^\alpha \text{in }\mathrm{\Omega }\\ u=0\text{on }\mathrm{\Omega }\\ x_iu_i+c_0r^{2+\alpha }+c_1=0\text{on }\mathrm{\Omega }\text{;}\end{array}$$
$`10`$
where $`\alpha `$, $`c_0`$ and $`c_1`$ are constants. If $`\beta :=(\alpha +2)(c_0(\alpha +N)1)`$ is not a negative integer, then $`\mathrm{\Omega }`$ is an $`N`$-dimensional ball.

###### Demonstration Proof

Introduce the functional $`V=x_iu_i+c_0r^{\alpha +2}+c_1+\beta u`$. Then by direct computation we obtain that V is harmonic. As can easily be seen, the boundary conditions make $`V0`$ on $`\mathrm{\Omega }`$. Thus,

$$\begin{array}{c}\mathrm{\Delta }V=0\text{in }\mathrm{\Omega }\\ V=0\text{on }\mathrm{\Omega }\text{ .}\end{array}$$

Classical Maximum Principles show that $`V0`$ inside $`\mathrm{\Omega }`$, too. Hence, it follows that

$$x_iu_i+c_0r^{\alpha +2}+c_1+\beta u=0\text{in }\overline{\mathrm{\Omega }}\text{.}$$
$`11`$
Rewriting equation (11) one obtains

$$\frac{}{r}\left(r^\beta (uu(0))+\frac{c_0r^{\beta +\alpha +2}}{\beta +\alpha +2}\right)=r^{\beta 1}\left(\beta (uu(0))+r\frac{u}{r}+c_0r^{\alpha +2}\right)=0.$$

This in turn implies that

$$r^\beta (uu(0))+\frac{c_0r^{\beta +\alpha +2}}{\beta +\alpha +2}G(\mathrm{\Theta })$$

for angular variables $`\mathrm{\Theta }`$. Hence,

$$uu(0)=c_0\frac{r^{\alpha +2}}{\beta +\alpha +2}+r^\beta G(\mathrm{\Theta })$$
$`12`$
Now, if $`\beta `$ 0 then $`u`$ cannot be regular at the origin unless $`G(\mathrm{\Theta })0.`$

If $`\beta <0`$ is not a negative integer, then once more we must have $`G(\mathrm{\Theta })0`$ for $`u`$ must satisfy $`\mathrm{\Delta }u=r^\alpha `$ and $`r^\beta G(\mathrm{\Theta })`$ is not harmonic with this value of $`\beta `$. Therefore, the solution $`u`$ is radial and

$$u=u(0)c_0\frac{r^{\alpha +2}}{\beta +\alpha +2}.$$

After using (10) and the vanishing of $`u`$ on $`\mathrm{\Omega }`$, this proves that $`\mathrm{\Omega }`$ is a ball thereby completing the proof of the proposition.