text:
# H I Recycling: Formation of Tidal Dwarf Galaxies

## 1. Recycling the H I gas

The VLA and other synthesis arrays have, for twenty years and more, revealed the spectacular distribution of the atomic hydrogen (H i) in interacting systems. Long tidal tails that are even more prominent than their optical counterparts, bridges, and ring–like structures are among the many weird and wonderful features which show up in high resolution maps (see Schiminovich et al., this volume). Because in general the H i in disk galaxies extends well beyond the optical $`R_{25}`$ radius, this atomic hydrogen reacts very efficiently to any external perturbation as evidenced by the fact that most of the neutral gas is found outside the colliding disks. In systems like Arp 105 and NGC 7252 (see Fig. 1), up to 90% of the gas visible at 21 cm is situated in the intergalactic medium (IGM) and little, if any gas remains within the parent spirals. In addition to gas being actively removed, a large part of the original H i has likely been funnelled inward where it was transformed into another phase (see review by Struck, 1999).

The fate of the tidally expelled gas will largely depend on its density and location with respect to the interacting galaxies. Whereas the clouds closest to the interacting (or merging) pair will fall back at timescales of a few Myr (Hibbard & Mihos, 1995), the most distant ones will become gravitationally unbound. They might slowly diffuse and enrich the IGM with heavy elements or the stars might even form the basis of the faint background light in the intracluster medium (ICM) if the interacting system belongs to a cluster. Under certain conditions, however, tidal debris may be recycled. If self–gravity is sufficiently large, expelled clouds will condense and collapse again to build new star–forming objects. Offspring as massive as magellanic dwarf irregular galaxies has been observed around several interacting systems (see for examples Fig. 1). Generally situated at the tip of 50–100 kpc long tidal tails, they are referred to as Tidal Dwarf Galaxies (TDGs). The total gas fraction that could end up in a TDG is yet largely unknown and it is hoped that numerical simulations might one day get a handle on this. Observationally, H i clouds towards TDGs as massive as $`510^9\mathrm{M}_{}`$ have been measured.

## 2. Forming tidal dwarf galaxies

Not only H i clouds are expelled during tidal interactions, but any material that was originally (rotating) in a disk, in particular the stars. Therefore TDGs are mixed bags, composed of young stars formed in situ from collapsing gas clouds and an older population pulled out from their parent galaxies. The latter component may be unimportant in systems involving gas–rich early type galaxies (i.e. NGC 5291, Duc & Mirabel 1998; see also Fig. 1). Instabilities in the gaseous — H i — component seem to be the driving factor in the formation of TDGs. Star formation occurs at rates which might reach 0.1 M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$. Given that TDGs contain huge H i reservoirs, typically $`10^9\mathrm{M}_{}`$, it is to be expected that the relative importance of an older stellar population will decrease with time as SF proceeds. Appreciable quantities of molecular gas have also recently been detected in TDGs by Braine et al. (2000). They suggest that this gas has been formed in situ out of the collapsed H i.

## 3. Identifying tidal dwarf galaxies

Tidal objects will survive provided that they have a potential well that is deep enough to sustain themselves against internal or external disruption. When identifying TDGs, it is therefore important to check that they are not simply the agglomerated debris of a collision but are self–gravitating entities (Duc et al. 2000). One should hence try to distinguish those tidal features that are kinematically decoupled from their host tails, the kinematics of which is governed by streaming motions. Because of the difficulty of obtaining high sensitivity, high resolution H i data, and problems related to projection effects along the line of sight, this is a difficult task that requires a careful examination of the available datacubes, such as those provided by synthesis radio observations (for the neutral gas component), and Fabry–Perot or any integral field instrument (for the ionised component). Addressing the stellar kinematics would be even more challenging. So far, evidence for such self–gravitating clouds have been found in the interacting systems Arp 105 (Duc et al. 1997; see Fig. 2), NGC 5291 (Duc & Mirabel 1998) and perhaps Arp 245 (Duc et al. 2000).

It is one thing to pick out TDGs which are still linked to the tidal tails out of which they formed. But how can one recognise old tidal dwarf galaxies that would have lost their physical connection with their parent galaxies as the tidal tail linking it vanishes with time? One might make use of three special properties of TDGs. First of all, their metallicity: due to the fact that they are recycled material, it is much higher than that of classical dwarf galaxies of the same luminosity (Duc et al. 2000). The oxygen abundance of recently formed TDGs averages out at one third of solar, i.e., the abundance of spirals at or slightly beyond the optical $`R_{25}`$ radius. Secondly, it is expected that TDGs contain little if any dark matter if DM in their progenitors was distributed in a large halo, as is traditionally assumed (Barnes & Hernquist, 1992). If however, DM is present in disks, for instance in the form of cold molecular gas (Pfenniger et al., this volume) DM dominated tidal tails should form and TDGs would have a substantial DM content. Finally, the stellar population of TDGs should at least be bimodal: a fraction of their stars has originally come from the parent disks, whereas another part has been formed in situ. Reconstructing the star formation history of a galaxy using for instance its color–magnitude diagram might reveal a tidal origin. These methods have already been used to identify some older TDG candidates, both in the field as well as in clusters (Hunter et al. 2000; Duc et al., in prep). But so far, the overall fraction of dwarf galaxies of tidal origin remains unknown.

### Acknowledgments.

We are grateful to our many collaborators involved in this multiwavelength study, and in particular to Felix Mirabel, Volker Springel, Barbara Pichardo, Peter Weilbacher and Jonathan Braine. We thank the organisers for a very successful and pleasant conference.

## References

Barnes, J.E. & Hernquist, L. 1992, Nature, 360, 715

Braine, J., Lisenfeld, U., Duc, P.–A. & Leon, S. 2000, Nature 403, 867

Dickey, J.M. 1997, AJ, 113, 1939

Duc, P.–A., Brinks, E., Wink, J., & Mirabel, I.F. 1997, A&A, 326, 537

Duc, P.–A. & Mirabel, I.F. 1998, A&A, 333, 813

Duc, P.–A., Brinks, E., Springel, V., Pichardo, B., Weilbacher, P., & Mirabel, I.F., 2000, AJ, in press (astro–ph/0006038)

Hibbard, J.E. & Mihos, J.C. 1995, AJ, 110, 140

Hibbard, J.E., Guhathakurta, P., van Gorkom, J.H. & Schweizer, F. 1994, AJ, 107, 67

Hunter, D.A., Hunsberger, S.D. & Roye, E.W. 2000, ApJ, in press (astro–ph/0005257)

Malphrus, B.K., Simpson, C.E., Gottesman, S.T. & Hawarden, T.G. 1997, AJ, 114, 1427

Struck, 1999, Physics Reports, 321, 1


format: text

text:
# Blob ejection from advection-dominated accretion flow: observational consequences

## 1. Introduction

The process of accretion onto a black hole is thought to be an efficient way to release gravitational energy. However for low enough accretion rates the released gravitational energy around a black hole is stored as entropy of accreting gas rather than being radiated. It is then expected that the internal energy of accreting gas will be very high and the temperature may reach the virial value. This is the well-known optically thin ADAF (Narayan & Yi 1994), which is quite different from that of standard accretion disks (Shakura & Sunyaev 1973).

Recently, based on the standard model of accretion disks, Tagger & Pellat (1999) found an accretion-ejection instability. The essential physics for this instability is that the extracted angular momentum is not emitted radially but ends up as an Alfvén wave which can eject material in the direction perpendicular to the disk plane. Such an instability may be stronger in ADAFs, because the instability grows with disk thickness although the detail of the instability in ADAFs remains open (Tagger 2000, private communication). Such a process provides another potential channel to release the gravitational energy in black hole accreting systems. The interesting properties of K$`\alpha `$line emission in radio-loud quasars lead to the suggestion that the central engine of radio-loud quasars may be powered by optically thin ADAFs (Sambruna et al. 2000a, b). If such an instability works in an ADAF, some blobs should be ejected. Then a natural question arises: what is the fate of a blob ejected from the optically thin ADAF? Comparison of the energy density of the ejected blob with that of the inter-cloud medium (ICM) in broad line regions (BLR) will give a clue to this question. The typical values of density and temperature of the ICM in the BLR are thought to be $`5.0\times 10^7`$K and $`5.0\times 10^6`$ cm-3, respectively, in radio-loud quasars (Netzer 1991). If a blob is ejected from an ADAF, what are its observational consequences both in line and continuum emission? These are the main goals of the present paper.

## 2. The Fate of a Blob Ejected from an ADAF

The physical properties of the ejected blob should be obtained from the global solution of an ADAF, which has already been done by Narayan et al. (1997). Since we are primarily interested in the fate of a blob, we can avoid the details of global structure of an ADAF. For simplicity, we deal with the two limiting cases to constrain the flow, i.e., fully and partially advected structure with free-free emission. The actual disk temperature and density functions in practice will be somewhere in between those given by these two limiting cases.

### 2.1. Optically Thin ADAF

For the first limiting case, namely that of full advection, we use the self-similar solution of Narayan & Yi (1994), and get temperature $`T_{_{\mathrm{ADAF}}}2\times 10^{12}r^1`$K, and density $`n_{_{\mathrm{ADAF}}}1.1\times 10^{11}\alpha ^1\dot{m}M_8^1r^{1.5}`$cm-3, where $`\alpha `$ is the viscosity parameter, $`\dot{m}`$ is the dimensionless accretion rate by $`\dot{M}_{\mathrm{Edd}}=1.4\times 10^{26}M_8`$(g/s)(the accretion efficient $`\eta =0.1`$), $`M_8`$ is the mass of the central black hole in units of $`10^8M_{}`$, and $`r=R/R_g`$ is the dimensionless radius with $`R_g=2GM_{}/c^2`$ and $`M_{}`$ is the mass of black hole.

We consider the other limiting case of an accretion disk, assuming a standard accretion disk (Shakura & Sunyaev 1973) holds. In such a case the flow cools via free-free emission. In this work we only consider the consequences of a single-temperature treatment as a first approximation. The height-averaged equations of an accretion disk are used since a more sophisticated approach is not warranted for this level of treatment. Using the $`\alpha `$-viscosity law, the dissipated energy production rate $`Q^+=\alpha p\mathrm{\Omega }__\mathrm{K}\gamma __\mathrm{K}`$, where $`\gamma __\mathrm{K}=d\mathrm{ln}\mathrm{\Omega }__\mathrm{K}/d\mathrm{ln}R`$ and $`\mathrm{\Omega }__\mathrm{K}=(GM_{}/R^3)^{1/2}`$ is the Keplerian angular velocity, is balanced by the free-free cooling rate $`Q^{}=\mathrm{\Lambda }_0\rho ^2T^{1/2}`$ ($`\mathrm{\Lambda }_0=1.24\times 10^{21}`$). Suppose that the vertical hydrostatic balance holds, i.e., $`H^2=c_s^2/\mathrm{\Omega }__\mathrm{K}^2=kT/m_p\mathrm{\Omega }__\mathrm{K}^2`$, where $`H`$ is the half thickness of a disk, $`k`$ is the Boltzmann constant, and $`m_p`$ is the proton mass. The $`\varphi `$-direction equation reads $`\dot{M}l__\mathrm{K}=4\pi \alpha R^2pH`$, where $`l__\mathrm{K}`$ is the Keplerian angular momentum. For simplicity, we have neglected the inner boundary condition that no heat generation occurs at the inner boundary. This is a reasonable approximation for the region not very close to the inner radius. Temperature and density are then given by

$$T_{\mathrm{ss}}=8.3\times 10^{11}r^{3/4}\alpha ^1\dot{m}^{1/2}\mathrm{K},$$
(1)
$$n_{\mathrm{ss}}=3.47\times 10^{13}r^{15/8}\alpha ^{1/2}\dot{m}^{1/4}M_8^1\mathrm{cm}^3.$$
(2)
Then the actual temperature $`T`$ and number density $`n`$ should be $`T_{\mathrm{ss}}TT_{_{\mathrm{ADAF}}}`$ and $`n_{_{\mathrm{ADAF}}}nn_{\mathrm{ss}}`$. We can see that even the free-free cooling disk is of much higher temperature (within $`10^4R_g`$) than the one producing recombination lines.

### 2.2. Blob Ejected from Optically Thin ADAF

In ADAFs the magnetic field strength is close to the equipartition value (Narayan & Yi 1995). The self-similar solution for such magnetized ADAFs has been found by Soria et al. (1997). It is interesting to note that there is an accretion-ejection instability within several gravitational radii in magnetized geometrically thin disks (Tagger & Pellat 1999). Such an accretion-ejection instability in ADAFs would be stronger than that in thin disks because its growth rate will increase with increasing disk thickness (Tagger 2000, private communication), though the instability in ADAFs needs to be studied in detail in future. Since most of the dissipated energy will be mainly deposited on protons, the internal energy density ($`U=\frac{3}{2}n_{\mathrm{ss}}kT_{\mathrm{ss}}`$) of accreting gas can be obtained by equations (1) and (2):

$$U=5.9\times 10^9r^{\frac{21}{8}}\alpha ^{\frac{1}{2}}\dot{m}^{\frac{3}{4}}M_8^1\mathrm{erg}/\mathrm{cm}^3.$$
(3)
We can see from above equation that the internal energy density is very sensitive to the ejection location of a disk, i.e., $`Ur^{21/8}`$. Thus the ejection at large radii is less significant. It is then assumed that a blob with size $`R_0(R_e)`$ is ejected from an ADAF at radius $`R_e`$ via such an instability, and that the most likely ejection radius is within $`R_e=10R_g`$. Also we assume that changes in temperature and density are little during the ejection.

The ejected blob will expand rapidly, and strongly interact with its surroundings. This behavior is very similar to the fireball of $`\gamma `$-ray bursts (Piran 1999). The interaction between the blob and its surroundings is very complex. There are two reasons for this complexity, because of two kinds of interactions between the ejected blob and its surrounding medium, i.e., expansion and transverse motions in the medium. In this paper we pay our main attention to the interaction between the blob expansion and it surrounding medium. Neglecting the detailed process of expansion (Sedov 1969) due to the high internal energy of the blob, we have $`d(\rho _{_{\mathrm{SS}}}v^2/2)/dt=dU/dt`$, for the expansion velocity, which is based on the energy conservation. Similarly to the fireball of $`\gamma `$-ray bursts, the expansion will be decelerated when the swept mass by expansion is equal to the blob initial mass $`M_{\mathrm{blob}}`$, namely, $`\frac{4}{3}\pi R_2^3n_2m_p=M_{\mathrm{blob}}`$, where $`R_2`$ is its final radius, $`R_2=R_0+v\tau v\tau `$, with $`\tau `$ being the timescale of expansion. These equations give $`v=(3kT_1/m_p)^{1/2}`$ and $`\tau =(R_0/v)(n_1/n_2)^{1/3}`$, where $`T_1`$ and $`n_1`$ are the temperature and the number density of the ejected blob at the initial state, respectively. The expansion proceeds with Mach number

$$^2\left(\frac{v}{c_{s2}}\right)^2\frac{T_1}{T_2},$$
(4)
where $`c_{s2}`$ is the sound speed of its surroundings. Typically $`^22\times 10^3`$, for $`T_2=5.0\times 10^7`$K and $`T_1=10^{11}`$K, and strong shock will be formed. It is expected that the temperature of the post shock wave will go down to $`T_1/^2`$ (Sedov 1969). Thus the ejected blob will be efficiently cooled down, so that the condition for some recombination lines, e.g., FeXXV or FeXXIV ions, will be satisfied (Raymond & Smith 1977). The separation of the ejected blob from the central black hole would be $`d=\beta c\tau `$, where $`\beta c`$ is the ejection velocity of the blob. There are mainly three effects strongly affecting the observed profile of emission lines: 1) spherical expansion causes the broadening of emission line profile, 2) the gravitational redshift (perhaps not strong), and 3) Doppler shift due to the bulk motion of the blob. The profile of an emission line from the relativistic outflow close to a rotating black hole has been calculated in our recent paper (Wang et al. 2000). It is expected that a line from such an ejected blob will have a more complex profile. The evolution of the blob will be mainly determined by free-free cooling after its swept mass of the surrounding medium equals to the initial mass of the blob. The density is then approximated as $`n_e\left[n_1R_0^3+n_2(v\tau )^3\right]/(v\tau )^3=2n_2=1.0\times 10^7`$ (cm-3). The timescale is $`\tau _{\mathrm{ff}}6.6T_8^{1/2}n_{e,7}^1`$ (yr), where $`T_8=T/10^8`$ and $`n_{e,7}=n_e/10^7`$. It is thus expected that the blob with temperature $`10^8`$K will cool down at 2pc from the nucleus, if the ejection occurs at almost light speed.

Next we briefly describe the possibility of nonthermal emission from a blob; blast waves are expected to be generated during expansion and to accelerate a group of thermal electrons to be relativistic. The conversion of the blob kinetic energy to nonthermal particles involves complicated plasma physics and shock acceleration processes. However, we can avoid dealing with the detailed acceleration physics by simply assuming that some fraction of the internal energy is transformed into a group of nonthermal electrons (e.g., Blandford & Eichler 1987). Here we assume that fraction $`\xi `$ of all the internal energy is channelled into the nonthermal electrons. From the energy conservation, the total number of relativistic electrons can be estimated from $`\xi R_0^3(n_1kT_{\mathrm{ss}})=(v\tau )^3n_0m_ec^2\gamma ^{1p}𝑑\gamma `$, where we assumed that shock acceleration forms a population of electrons with power-law index $`p`$ as $`n_0\gamma ^p`$. The total number of relativistic electrons in one blob is

$$N_e=\frac{4}{3}\pi \xi n_1R_0^3\left(\frac{kT_{\mathrm{ss}}}{m_ec^2}\right)\left(\frac{2p}{1p}\right)\left(\frac{\gamma _{\mathrm{max}}^{1p}\gamma _{\mathrm{min}}^{1p}}{\gamma _{\mathrm{max}}^{2p}\gamma _{\mathrm{min}}^{2p}}\right).$$
(5)
The fraction of the relativistic electrons to the total number of the initial blob is about $`\xi kT_{\mathrm{ss}}/\gamma _{\mathrm{min}}m_ec^2=0.012\xi _2`$ for $`p>2`$, where $`\xi _2=\xi /0.01`$ and $`\gamma _{\mathrm{min}}=50`$. This strongly implies that the radiating electrons may originate from the present mechanism, which may be a proposing mechanism responsible for the acceleration of relativistic electrons in some of radio-loud quasars. It is expected the nonthermal emission spectrum and light curves should be similar to some properties predicted by Li & Kusunose (2000) and Kusunose et al. (2000).

It is thus inevitable for the ejected blob that the high-temperature drops due to the strong interaction between the blob and its surroundings, and as a natural consequence there is a group of electrons to be accelerated to relativistic energy responsible for the continuum and some recombination lines from the cooled blob.

## 3. Applications

It has been generally accepted that accretion onto super massive black holes leads to the release of gravitational energy in active galactic nuclei (Rees 1984). However, the status of the accretion disk in different kinds of AGNs remains uncertain. Ion-supported tori may power the central engine in radio-loud quasars (Rees et al. 1982). A famous low luminosity AGN NGC 4258 is generally thought as representative of ADAFs (Gammie et al. 1999). However, recent observations of K$`\alpha `$by ASCA show the origination is not from an accretion disk (Reynolds et al. 2000). Yaqoob et al. (1999) detected a highly Doppler blue-shifted K-emission line in PKS 2149–306. We thus apply our model to the two objects.

NGC 4258: The mass of the central black hole is well determined to be $`M_{}=3.6\times 10^7M_{}`$ by Miyoshi et al. (1995). Radio observation set the transition radius from the standard disk to the ADAF to be 100$`R_g`$ (Herrnstein et al. 1999), and the accretion rate is $`1.0\times 10^3\dot{M}_{\mathrm{Edd}}`$ (Gammie et al. 1999). Taking the typical value of $`\alpha =0.1`$, we have $`T_1=4.7\times 10^{10}`$K and $`n_1=7.2\times 10^{10}`$cm-3, if the ejection takes place at 10$`R_g`$, where the viscous heating rate is highest. The initial dimension of the ejected blob is taken to be $`R_0=10R_g=10^{14}`$ cm. The additional component with luminosity $`L_X^A=2.0\times 10^{40}`$ ergs/s in X-ray shows that the tenuous plasma is of temperature 0.5 keV (Reynolds et al. 2000). This sets constraints on the density of the environment in the vicinity of the nucleus via the Thomson scattering depth $`\tau _{_{\mathrm{Th}}}=\sigma _{_{\mathrm{Th}}}n_2\mathrm{}`$, where $`\sigma _{_{\mathrm{Th}}}`$ is the Thomson scattering cross section and $`\mathrm{}`$ is the dimension of its surroundings. Using the X-ray luminosity of the additional component contributed from the tenuous plasma $`L_X^A=\frac{4}{3}\pi \mathrm{}^3\mathrm{\Lambda }n_2^2T_e^{1/2}`$ ($`\mathrm{\Lambda }=3.4\times 10^{27}`$), we have the number density of the surroundings to be $`n_2=10^8(\tau _{_{\mathrm{Th}}}/0.1)^3`$cm-3. From equation (4) we have the Mach number $`=97`$. The timescale of cooling due to expansion is then given by $`\tau =2.36\times 10^5`$s, which is much shorter than that of free-free cooling. The expanded radius of the blob is $`R_2=9.0\times 10^{14}`$cm. The temperature will drop from $`4.7\times 10^{10}`$K to $`T=T_1/^25.0\times 10^6`$K. During the cooling some recombination lines will be produced (Raymond & Smith 1977), such as FeXXVI$`\lambda 1.78`$Å, FeXXV$`\lambda 1.85`$Å, FeXXV$`\lambda 1.86`$Å, and oxygen line OVIII$`\lambda 18.97`$, and silicon Si XIV$`\lambda 6.18`$Å, etc. If the blob moves with Doppler factor of 10 (Gammie et al. 1999), we expect to observe the highly Doppler blue-shifted line at 64$`𝒟_{10}`$ keV, where $`𝒟=10𝒟_{10}`$ is the Doppler factor. Such lines should be detected by future observations of INTEGRAL.

PKS 2149–306: The highly Doppler shifted iron K-emission line has been detected by ASCA in this high redshift ($`z`$=2.345) radio-loud quasar (Yaqoob et al. 1999), although the current data does not allow to determine its unambiguous profile. We attempt to apply our present model to this source. The central mass can be obtained from the full-width-at-half-maximum (FWHM) and the luminosity of emission line CIV (Peterson 1998). With the measurements of $`v_{_{\mathrm{FWHM}}}=6400`$km/s (Wilkes 1986), and $`L(\mathrm{CIV})=1.2\times 10^{46}`$erg/s (Wilkes et al. 1983) (which we measured from the spectrum), we have the mass of the central black hole, $`M_{}=3.4\times 10^9M_{}`$, which corresponds to Eddington luminosity $`L_{_{\mathrm{Edd}}}=4.3\times 10^{47}`$erg/s. If we follow the mean spectrum of radio-loud quasars (Elvis et al. 1994), we get the bolometric luminosity $`L_{\mathrm{bol}}5.0\times 10^{47}`$erg/s based on the continuum (Siebert et al. 1996). Then the accretion rate is roughly $`\dot{m}=L_{\mathrm{bol}}/(𝒟^4L_{_{\mathrm{Edd}}})2.0\times 10^2`$, where Doppler factor $`𝒟=2.65`$; here we directly used the Doppler factor of the iron K emission line as jet’s. It suggests that this object may be powered by an ADAF. The transition radius from the standard disk to the ADAF can be obtained by the approximate formula $`r_{\mathrm{tr}}2.6\times 10^3\alpha ^4\dot{m}^2`$ for the case that half the released energy is advected \[their eq. (4.1) of Narayan & Yi 1995\]. We get $`r_{\mathrm{tr}}400`$ for $`\alpha =0.1`$. If the ejection takes place at $`10R_g`$ with initial radius $`R_0=10R_g10^{16}`$cm, then the density is $`n_1=1.7\times 10^8`$cm-3 and the temperature is $`T_1=2.3\times 10^{11}`$K. Taking the typical values of the parameters, $`n_2=5.0\times 10^6`$cm-3 and $`T_2=5.0\times 10^7`$K, for ICM in the BLR (Netzer 1991), we have the typical expansion timescale $`\tau 2.8\times 10^6`$s, and the temperature will drop to $`10^{78}`$K, where Mach number $`^2=4.7\times 10^3`$ from equation (4). The blob’s distance from the center is about $`\beta c\tau 8.4\times 10^{16}`$cm$`80R_g`$ ($`\beta 1`$). It is thus expected that some high energy recombination lines will appear in $`2.8\times 10^6`$s since the ejection of a blob.

There should be two steps to eject a blob. The first is the blob formation with the timescale of $`\tau _\mathrm{f}`$, which is approximate to the viscous timescale (timescale to accumulate matter), namely $`\tau _\mathrm{f}\alpha ^1\tau __\mathrm{K}`$ \[$`\tau __\mathrm{K}=1/\mathrm{\Omega }__\mathrm{K}=(R_g/c)r^{3/2}=10^6`$s at $`r=10`$\]. In the more accurate global solution the radial velocity is higher than the above one. The second is the relativistic ejection with timescale $`\tau _\mathrm{e}`$. It seems reasonable to assume that the timescale of ejection is much shorter than the timescale of blob formation. The duration for a blob to produce the iron recombination lines is approximately $`\tau _{\mathrm{ff}}=6.6`$yr for $`T_e=10^8`$K and $`n_e=1.0\times 10^7`$cm-3. Therefore the number of the blobs producing the recombination line of iron K-emission would be $`N_b\tau _{\mathrm{ff}}/(\tau _\mathrm{e}+\tau _\mathrm{f})\alpha \tau _{\mathrm{ff}}/\tau __\mathrm{K}20`$ if the ejection takes place at 10$`R_g`$ (It should be noted that this number is not the total of the ejected blobs). We assume the cosmic abundance of iron in the blob. Considering the boosting effects of flux due to the relativistic motion, the intrinsic flux is given by $`F_{\mathrm{in}}=F_{\mathrm{obs}}𝒟_b^4`$, where $`F_{\mathrm{obs}}`$ is the observed flux (erg s-1 cm-2) and $`𝒟_b`$ is the Doppler factor of outflow. The energy of K-emission line with emissivity $`j__\mathrm{K}`$ within the time interval $`\mathrm{\Delta }t`$ is $`\mathrm{\Delta }E\frac{4}{3}\pi N_bj__\mathrm{K}n_e^2R_2^3\mathrm{\Delta }t`$, which can be compared with the observable $`\mathrm{\Delta }E=4\pi d_L^2F_{\mathrm{obs}}\mathrm{\Delta }t/𝒟_b^4`$, where $`n_e`$ is the number density of electrons in the blob. With the help of energy conservation, we can predict the observed flux: $`F_{\mathrm{obs}}=N_\mathrm{b}𝒟_b^4j__\mathrm{K}n_e^2R_2^3/3d_L^20.8\times 10^{13}`$erg s-1 cm-2, where we take the emissivity of iron recombination lines, $`j__\mathrm{K}`$=1.4$`\times 10^{24}`$erg s-1 cm-3 (Raymond & Smith 1977), and $`𝒟_b=2.65`$ (Yaqoob et al 1999). We find this estimated value is close to the observed flux, $`F_{\mathrm{obs}}3.6_{2.5}^{+2.5}\times 10^{13}`$erg s-1 cm-2 (Yaqoob et al. 1999). Here we do not rule out the contribution of iron K$`\alpha `$fluorescence emission to the observed flux. The total number of relativistic electrons in one blob is $`N_e1.5\times 10^{54}`$, if we take the minimum energy of electrons $`\gamma _{\mathrm{min}}=50`$. These electrons will be responsible for the nonthermal emission. The detailed radiation processes and spectra with the evolution of electron’s energy distribution will be carried out in a separate paper.

## 4. Conclusions and Discussions

If a blob ejection takes place in such an ADAF the fate of the blob is an interesting topic. In the ADAF regime the released gravitational energy from the viscous dissipation is restored as the entropy of the accreting gas. We have explored the observational consequences of the blob ejected from an optically thin ADAF. We found it inevitably drastically expands into its surroundings with Mach number of $`10^2`$. The hot blob cools down rapidly via the interaction between the blob and its surroundings via the blob’s expansion. Consequently a group of thermal electrons will be accelerated to be relativistic via strong shocks, leading to nonthermal radiation from the blob. Simultaneously there will appear some of recombination lines from the blob. If the ejected blob moves, these lines will be highly Doppler shifted. We think the observed iron K emission line in PKS 2149–306 results from the ejected blob from an optically thin ADAF. We also made predictions for NGC 4258 for future observations of INTEGRAL.

Recent studies on the K$`\alpha `$line properties in radio-loud quasars show that the line is weak and narrow in radio-loud quasars, but strong and broad in radio-quiet quasars, which are quite similar to that in Seyfert galaxies (Reynolds 1997, Sambruna et al. 1999a, Wozniak et al. 1998, Reeves & Turner 2000, Fabian et al. 2000). The ADAF thus may power the central engines of radio-loud quasars (Sambruna et al. 1999b). The proposed model is consistent with the existing data of PKS 2149-306 and supported by the recent observations of K$`\alpha `$line in other radio-loud quasars. The thermal X-ray emission lines from an ADAF have been studied by Narayan & Raymond (1999), who suggest these lines may be detected by the current X-ray instruments (Narayan & Raymond 1999, Perna et al. 2000). If the future observations (high space resolution of radio telescopes and high energy resolution of X-ray telescopes) confirm our predictions, the highly Doppler blue-shifted iron K emission line may shed light on jet formation and be very helpful to the construction of central engine models of radio-loud quasars.

It is likely for the ejected blob to encounter BLR clouds when it passes through the BLR. One may expect to see changes in the optical/UV emission line spectrum in the simultaneous multiwavelength observations. For a radio-loud quasar it would be difficult to test the X-ray continuum from the ADAF itself since the nonthermal emission dominates its continuum by the Doppler boosting effects. The spatial distribution of ejected blobs is expected to be perpendicular to the accretion disk plane. Their distributions may lead to interesting observables, for example, the line profile will be broadened, if they occupy a considerable solid angle.

The proposed model is based on two important assumptions. First, we assumed a blob is ejected from an ADAF, but we did not specify the working mechanism responsible for the ejection. This leads us not to accurately predict the number of blobs. The accretion-ejection instability found by Tagger & Pellat (1999) could be a promising mechanism and we will study it in future. The second involves the interaction between the blob with its surroundings. If the density of blob’s surroundings is too low, the fate of the ejected blob would be different from the proposed model.

The authors are grateful to an anonymous referee for the helpful comments and suggestions on the manuscript clarifying several points. The useful discussions with T.P. Li, F.J. Lu and J.L. Qu are acknowledged. J.M.W. is supported by ’Hundred Talents Program of CAS’. This work is financed by the Special Funds for Major State Basic Research Projects and by NSFC.


text:
# A transfer matrix approach to the enumeration of plane meanders

## 1 Introduction

Meanders  form a set of combinatorial problems concerned with the enumeration of self-avoiding curves crossing a line through a given number of points. Meanders are considered distinct up to any smooth deformation. This problem seems to date back at least to work of Poincaré on differential geometry . Since then it has from time to time been studied by mathematicians in various contexts such as the folding of a strip of stamps  or folding of maps . More recently it has been related to enumerations of ovals in planar algebraic curves  and the classification of 3-manifolds . During the last decade or so has it has received considerable attention in other areas of science. In computer science meanders are related to the sorting of Jordan sequences  and have been used for lower bound arguments . In physics meanders are relevant to the study of compact foldings of polymers , properties of the Temperley-Lieb algebra , matrix models , and defects in liquid crystals and $`2+1`$ dimensional gravity .

The difficulty in the enumeration of most interesting combinatorial problems is that, computationally, they are of exponential complexity. That is to say, the time it takes to calculate the first $`n`$ terms in the generating function grows asymptotically as $`\lambda ^n`$, where $`\lambda >1`$ is the growth rate. Initial efforts at computer enumeration of meanders were based on direct counting. Independently, Koehler  and Lunnon , studied the number of ways of folding a strip of stamps (or a map) of length $`n`$, and published result up to, respectively, $`n=16`$ and 24. Lando and Zvonkin  studied closed meanders, open meanders and multi-component systems of closed meanders, and calculated the number of open meanders up to $`n=26`$ and the number of closed meanders up to $`n=14`$. The calculation of the number of closed meanders was subsequently extended up to $`n=16`$ by Pratt . Di Francesco et al. studied semi-meanders, a problem equivalent to the stamp folding problem , and extended the calculation to $`n=29`$ , in addition they studied numerous other problems including that of multi-component systems of semi-meanders.

Few exact and mathematically rigorous results have been obtained for any of the many meander problems. However, in a recent paper it was conjectured that some of the meander problems can be related to a gravitational version of a certain loop model . From the conformal field theory of the model, conjectures were proposed for the exact critical exponent of closed and open meanders, $`\alpha =(29+\sqrt{145})/12=3.4201328\mathrm{}`$, as well as the exponent for semi-meanders, $`\overline{\alpha }=1+\sqrt{11}(\sqrt{29}+\sqrt{5})/24=2.0531987\mathrm{}`$. This work has recently been extended to multi-component systems of closed and semi-meanders  and to various other geometries. Conjectures were given for the critical exponents as functions of the loop-fugacity $`q`$. These were checked numerically  and found to be correct within numerical error. In a recent paper  we analysed extended series for the meander generating functions. Using the numerical technique of differential approximants  we obtained accurate estimates for the exponents and found that the conjecture for the semi-meander exponent is unlikely to be correct, while the conjecture for closed meanders is just consistent with the results from the analysis.

The purpose of this paper is to give a detailed description of the new and improved algorithm used to derive the series studied in . While the algorithm still has exponential complexity, the growth rate is much smaller than that experienced with direct counting, and consequently the calculation can be carried much further. The algorithm is easily modified to enumerate various multi-component systems of closed meanders, semi-meanders or open meanders. In particular we have extended the calculation for closed meanders up to $`n=24`$, for open meanders up to $`n=43`$, and for semi-meanders up to $`n=45`$.

In section 2 we shall briefly describe some meander problems and define the meandric numbers. Section 3 contains a detailed description of the transfer matrix algorithm for the enumeration of closed meanders and an outline of generalisations to other problems. Finally we give our conclusions in section 4.

## 2 Definitions of meanders

A closed meander of order $`n`$ is a closed self-avoiding curve crossing an infinite line $`2n`$ times (see figure 1). The meandric number $`M_n`$ is simply the number of such meanders distinct up to smooth transformations. Note that each meander forms a single connected curve. The number of closed meanders is expected to grow exponentially, with a sub-dominant term given by a critical exponent, $`M_nCR^{2n}/n^\alpha `$. The exponential growth constant $`R`$ is often called the connective constant. The generating function is expected to behave as

$$(x)=\underset{n=1}{\overset{\mathrm{}}{}}M_nx^nA(x)(1R^2x)^{\alpha 1},$$
(1)
and hence have a singularity at $`x_c=1/R^2`$ with exponent $`\alpha 1`$.

We can extend the definition to multi-component systems of closed meanders, where we allow configurations with several disconnected closed meanders. The meandric numbers $`M_n^{(k)}`$ are the number of meanders with $`2n`$ crossings and $`k`$ components, and we thus obtain the more general generating function:

$$(x,q)=\underset{n=1}{\overset{\mathrm{}}{}}\underset{k=1}{\overset{n}{}}M_n^{(k)}x^nq^k.$$
(2)
Obviously, $`M_n=M_n^{(1)}`$, and $`(x)=(x,q)/q|_{q=0}`$. An example of a multi-component meander system with several pieces is shown in figure 2. In addition one could also look at further restrictions on these systems. Two obvious examples would be to exclude configurations where meanders are nested within one another or where meanders can be separated from one another, i.e., all smaller meanders are completely enclosed within a larger meander.

An open meander of order $`n`$ is a self-avoiding curve running from west to east while crossing an infinite line $`n`$ times (see figure 3). The number of such curves is $`m_n`$ and we can define a generating function for this problem in analogy with (1). It should be noted  that $`M_n=m_{2n1}`$, and hence the critical exponent is identical to that of closed meanders and the connective constant is $`R`$.

Finally, instead of looking at crossings of an infinite line one could consider a semi-infinite line and allow the curve to wind around the end-point of the line . A semi-meander of order $`n`$ is a closed self-avoiding curve crossing the semi-infinite line $`n`$ times. The number of semi-meanders of order $`n`$ is denoted by $`\overline{M}_nC^{}\overline{R}/n^{\overline{\alpha }}`$ and we define a generating function as in (1). In this case a further interesting generalisation is to study the number of semi-meanders $`\overline{M}_n(w)`$ which wind around the end-point of the line exactly $`w`$ times. Again we could also study systems of multi-component semi-meanders according to the number of independent meanders. Two semi-meanders are shown in figure 4.

## 3 Enumeration of meanders

The method used to enumerate meanders is based on the transfer matrix algorithm devised by Enting  in his pioneering work on the enumeration of self-avoiding polygons. Derrida  used a similar algorithm to study self-avoiding walks. The transfer matrix technique involves drawing a boundary line perpendicular to the infinite line. The intersection between the boundary and a given meander results in a set of loop-ends. Each loop-end is connected (to the left of the boundary) to at most one other loop-end. In the case of closed meanders the matching is perfect and each loop-end is connected to exactly one other loop-end. In the case of open meanders there is in addition one and only one loop-end which is free and therefore not connected to any other loop-ends, as illustrated in figure 5. In addition to the information describing the configuration of loop-ends, and how they are connected, we need to know where the infinite line is situated within the loop-ends. This can be done simply by specifying how many loop-ends lie beneath the infinite line. For each such configuration we keep count of all the possible (partially completed) meanders which gives rise to that particular configuration of loop-ends. Meanders can then be enumerated by successive moves of the boundary line, so that exactly one crossing is added with each move. An extra crossing is added either by putting in a new loop across the infinite line or by taking an existing loop-end immediately above/below the line and dragging it to the other side.

These general remarks hold for all the meander enumeration problems. In the following we give a detailed description of the algorithm used in the enumeration of closed connected meanders. Afterwards we describe how to generalise the method to other meander problems.

### 3.1 Enumeration of closed connected meanders

In the enumeration of closed connected meanders the two major constraints which must be observed are self-avoidance and the constraint that all meanders constructed during the calculation must consist of a single connected component. As we move the boundary line, the partially completed meanders will consist of a number of disjoint loop segments, which must be connected to each other if a valid closed meander is to be produced. A pair of loops can be placed relative to one another in two distinct ways, namely, side by side or nested, as shown in the left panel of figure 6. In each case it is possible to connect the loop-ends so as to form a single loop (middle panel) or so as to form graphs with two separate components (right panel). It is connections equivalent to these latter cases which we must avoid. So the constraint, which must be observed in order to avoid separate components, is that a loop can be closed on itself only if the boundary intersects no other loops.

To avoid loops closing on themselves we have to label the loop-ends so we can keep track of how they are connected to one another. One scheme would be to assign a unique label to each loop as indicated in figure 5. However, due to the two-dimensional nature of the problem and since the curve making up a meander is self-avoiding, there is a scheme better suited to a computer algorithm. Each configuration of loop-ends can be represented by an ordered set of states $`\{\sigma _i\}`$, where

$$\sigma _i=\{\begin{array}{cc}\hfill 0& \text{lower end of a loop},\hfill \\ \hfill 1& \text{upper end of a loop}.\hfill \end{array}$$
(3)
It is easy to see that this encoding uniquely describes which loop-ends are connected. In order to find the upper loop-end, matching a given lower end, we start at the lower end and work upwards in the configuration counting the number of ‘0’s and ‘1’s we pass (the ‘0’ of the initial lower end is not included in the count). We stop when the number of ‘1’s exceeds the number of ‘0’s. This ‘1’ marks the matching upper end of the loop. It is worth noting that there are some restrictions on the possible configurations. Firstly, every lower loop-end must have a corresponding upper end, and it is therefore clear that the total number of ‘0’s is equal to the total number of ‘1’s. Secondly, as we look through the configuration starting from the bottom the number of ‘0’s is never smaller than the number of ‘1’s. Those familiar with algebraic languages will immediately recognise that each configuration of labelled loop-ends forms a Dyck word (see ).

#### 3.1.1 The transfer matrix algorithm

The total configuration of loop-ends and their placement relative to the infinite line can thus be described by a pair of integers $`(h,S)`$, where $`h`$ is the number of loop-ends below the infinite line and $`S`$ is the integer whose binary representation corresponds to the configuration of loop-ends. We shall call such a $`(h,S)`$-pair a signature, and in practise we represent it by a 64-bit integer with the first 6 bits coding $`h`$ and the remaining bits coding $`S`$. In the following we shall often explicitly write out the binary representation, $`\{b_0b_1\mathrm{}b_n\}`$ of $`S`$, and use the notation $`\{S_1S_2\}`$ to mean a configuration of loop-ends obtained by concatenating the binary strings $`S_1`$ and $`S_2`$.

The algorithm for enumerating closed connected meanders is:

1. Set the maximum order $`n`$ of the meandric numbers we wish to calculate. Start with the signature $`(1,\{01\})`$ with a count of 1, that is one loop crossing the infinite line. Set the number of crossings $`n_c=1`$. Set the meandric number $`M_1=1`$.
2. Next move the boundary one step ahead and add one more crossing, $`n_c=n_c+1`$, to all existing signatures. A new crossing is added either by putting in an additional loop or by dragging an existing loop-end (immediately above or below) across the infinite line.

A new loop is inserted simply by taking an existing source signature $`(h,\{S_1S_2\})`$, where $`S_1`$ ($`S_2`$) is the configuration of loop-ends below (above) the infinite line, and creating a new target signature $`(h+1,\{S_101S_2\})`$. The count of the source is added to the count of the target. This first type of move is illustrated in figure 7 where in moving to position 2 we generate the target signature $`(2,\{0011\})`$ from the source $`(1,\{01\})`$. Additional loops are also inserted while moving to positions 4 and 7.

As we cross the infinite line with an existing loop-end we either leave it at that or we may be allowed to connect the loop-end to the loop-end on the other side of the infinite line.

From the source $`(h,S)`$ we produce the the two targets $`(h1,S)`$ and $`(h+1,S)`$, by dragging a loop-end from below to above and vice versa. Both these updates are allowed provided $`h10`$ and $`h+12m_1`$, where $`m_1`$ is the number of ‘1’s in $`S`$. For $`h=0`$ only the target $`(1,S)`$ is allowed and for $`h=2m_1`$ only the target $`(h1,S)`$ is allowed. Examples of these moves are given in figure 7 when moving to positions 3, 5, 8 and 11.

There are four distinct cases depending on whether the loop-ends below and above are of type ‘0’ or ‘1’.

In this case we connect a lower loop-end from below the line to a lower loop-end above the line. From a source $`(h,\{S_100S_2\})`$ we generate the target $`(h1,\{S_1\widehat{S_2}\})`$, where the symbol $`\widehat{S_2}`$ indicates that this string is changed via further processing. This is so because by connecting the two lower loop-ends an upper loop-end elsewhere in the old configuration $`S_2`$ becomes a lower loop-end in the new configuration $`\widehat{S_2}`$. An example of this type of relabelling is shown in figure 7 where in the move to position 9 we see that the signature $`(2,\{000111\})`$ before the step becomes the configuration $`(1,\{0011\})`$ after the step. That is, the upper end of the third loop before the step becomes the lower end of the second loop after the step. In general the nesting of loops could be more complicated and the general rule for the relabelling of the configuration is as follows: When connecting two ‘0’s we work upwards in the configuration, counting the number of ‘0’s and ‘1’s we pass until the number of ‘1’s exceeds the number of ‘0’s. This ‘1’ is the matching end of the inner loop and it should now be changed to a ‘0’, thus becoming the lower end of the outer loop (drawing a few further figures should make this relabelling clearer).

In this case we connect an upper loop-end from below to a lower loop-end above. So from a source $`(h,\{S_110S_2\})`$ we generate the target $`(h1,\{S_1S_2\})`$.

This is never allowed since it would result in a closed loop and thus generate graphs with separate components. The only exception is when there are no other loop-ends in $`S`$, but this case is dealt with in 3.

In this case we connect an upper loop-end from below the line to an upper loop-end above the line. From a source $`(h,\{S_111S_2\})`$ we generate the target $`(h1,\{\widehat{S_1}S_2\})`$. The rule for the relabelling of $`S_1`$ $`\widehat{S_1}`$ is similar to the case ‘00’, but we work downwards in the string $`S_1`$ until we find the unmatched lower loop-end, which is then changed to an upper loop-end.

Note that all of the above moves may be allowed. So from a given source we can generate up to four targets, by adding a new loop, by dragging a loop-end from below to above the infinite line, doing the reverse, or by connecting two loop-ends across the infinite line. As we move along and generate new target signatures their counts are calculated by adding up the count for the various source signatures which could generate that target. For example the target $`(2,\{0011\})`$ is generated from the sources $`(1,\{01\})`$, $`(1,\{0011\})`$, $`(3,\{0011\})`$, and $`(3,\{001011\})`$, by, respectively, putting in an additional loop, moving a loop-end below the line, moving a loop-end above the line and connecting two loop-ends across the line.
3. If $`n_c`$ is odd then set $`j=(n_c+1)/2`$ and extract the meandric number $`M_j`$ as the count of the signature $`(1,\{01\})`$. This is the only case in which we are allowed to close a loop. Doing so obviously adds one more crossing.
4. If $`n_c<2N1`$ go to 2.

Not all possible signatures that can be generated in a calculation to order $`n`$ are actually required. The main restriction is that no meanders should have more than $`2n`$ crossings. Since each move adds one more crossing and reduces the number of loop-ends above/below the infinite line by at most one, it is clear that for a given signature we have to add at least $`n_a=\mathrm{max}(h,2m_1h)`$ additional crossings in order to produce a closed meander. Thus if, for a given signature, $`n_c+n_a>2n`$, we can discard the signature since it would contribute only to a meandric number exceeding the order to which we wish to carry out the calculation. Further savings of a factor of almost 2 is obtained by using the symmetry with respect to reflection in the infinite line. A further factor of approximately 2 is obtained as follows. Note that $`n_a`$ is the minimum number of additional crossings and that for some signatures further crossings are needed. The most obvious case is when the loop-ends above and below the infinite line are connected to one another. In this case we cannot connect the two ends and first we have to move one of the loop-ends across the line. So when $`h=m_1`$, this results in at least two extra crossings. In the general case one can readily write an algorithm to count the actual number of additional crossings required, and as stated above this results in a saving of close to 2 in the number of signatures one need retain.

### 3.2 Generalisations to other meander problems

#### 3.2.1 Multi-component systems of closed meanders

As we noted above connecting a ‘0’ below the line to a ‘1’ above the line results in a closed loop. Failure to observe the restriction on this closure would result in graphs with disconnected components, either one closed meander over another or one closed meander within another. Obviously these are just the types of graphs required in order to enumerate multi-component systems of closed meanders. So by noting that each such closure adds one more component it is straightforward to generalise the algorithm to enumerate systems of closed meanders. The only major change is that, rather than just storing the number of partially completed meanders, for each signature we have to store a generating function, that is a polynomial giving the number of partially completed meanders with $`k`$ components, where $`1kn`$.

#### 3.2.2 Open meanders

Open meanders are a little more complicated. The first part of the necessary generalisation consists in adding an extra piece of information to our signature. We have to keep track of a single free end by specifying its position within the configuration of connected loop-ends. One simple way of doing this is, in analogy with the infinite line, to specify the number, $`h_f1`$, of (connected) loop-ends below the free end, so that $`h_f`$ is the position of the free end as counted from the bottom. So a configuration is now described by a signature $`(h,h_f,S)`$. Naturally, we also have to generalise the algorithm described above. We now start with the signature $`(1,1,0)`$, that is a single free end below the infinite line, and no crossings. The updating rules for adding a new crossing are very similar to the ones described above for closed meanders. One difference is that when a new loop is added (two loop-ends joined) below the free end, $`h_f`$ is increased (decreased) by 2. We also need to consider what happens when joining the free end to a connected loop-end. In this case we have to change the matching end of the connected loop to the new free end in the target signature and we have to change $`h_f`$ accordingly. An example illustrating this is shown in figure 8. The updating rule when the free end does not join with the loop-end on the other side is obviously just to increase (decrease) $`h`$ by 1 as the free end is moved below (above) the line.

#### 3.2.3 Semi-meanders

Finally, in order to enumerate semi-meanders all we need do is change the starting configuration. We now start in a position just before the first crossing of the semi-infinite line with $`w`$ loops nested with one another. By running the algorithm for each $`w`$ from 0 to $`n`$ we can count all semi-meanders with up to $`n`$ crossings. The generalisation to multi-component systems of semi-meanders is the same as for closed meanders.

### 3.3 Computational complexity

Using the new algorithm we have calculated $`M_n`$ up to $`n=24`$ as compared to the previous best of $`n=16`$ obtained by V. R. Pratt . To fully appreciate the advance it should be noted that the computational complexity grows exponentially, that is the time required to obtain $`n`$ term grows asymptotically as $`\lambda ^n`$. For direct enumerations time is simply proportional to $`M_n`$ and thus $`\lambda =lim_n\mathrm{}M_{n+1}/M_n12.26`$. Thus extending the count of the meandric number from 16 to 24 by direct counting would have required approximately $`12.28^85\times 10^8`$ as much CPU time as the calculation of the first 16 terms. The transfer matrix method employed in this paper is far more efficient. In figure 9 we have plotted the maximum number of signatures required in order to calculate the number of closed meanders up to order $`n`$. As can be seen, the number of signatures grows exponentially with $`n`$, and the numerical evidence suggests that the computational complexity is such that $`\lambda 2.5`$, which obviously is a very significant improvement on direct counting. The drawback of the transfer matrix method is that, since we need to store all the different signatures, the memory requirement of the algorithm also grows exponentially with growth rate $`\lambda `$, whereas direct counting algorithms typically have memory requirements which are linear in $`n`$. In fact it is exactly the memory requirement which is the major limitation of the transfer matrix method. The calculations reported in this paper used up to 2Gb of memory and typically took a few days of CPU time.

Di Francesco et al.  carried out a more detailed analysis of the complexity of the algorithm as it applies to multi-component systems of closed meanders and obtained the estimate $`N_{\mathrm{sig}}a^n=(2.582603\mathrm{})^n`$, for the number of signatures, very close to the purely empirical estimate obtained here.

Another way of gauging the improved efficiency is to note that the calculations for semi-meanders carried out in  were “done on the parallel Cray-T3D (128 processors) of the CEA-Grenoble, with approximately 7500 hours $`\times `$ processors.” Or in total about 100 years of CPU time. The equivalent calculations can be done with the transfer matrix algorithm in about 15 minutes on a single processor workstation!

### 3.4 Further particulars

Finally a few remarks of a more technical nature. The number of contributing configurations becomes very sparse in the total set of possible states along the boundary line and as is standard in such cases one uses a hash-addressing scheme . Since the integer coefficients occurring in the expansion become very large, the calculation was performed using modular arithmetic . This involves performing the calculation modulo various prime numbers $`p_i`$ and then reconstructing the full integer coefficients at the end. In calculations involving only single component meanders most of the memory is used storing the signatures so we used primes of the form $`p_i=2^{30}r_i`$. The Chinese remainder theorem ensures that any integer has a unique representation in terms of residues. If the largest integer occurring in the final expansion is $`M`$, then we have to use a number of primes $`m`$ such that $`p_1p_2\mathrm{}p_m>M`$. Up to 3 primes were needed to represent the coefficients correctly. In calculations involving systems of multi-component meanders the main storage requirement comes from the generating functions. In order to save memory we used primes of the form $`p_i=2^{15}r_i`$ so that the residues of the coefficients in the polynomials could be stored using 16 bit integers. In this case we used up to 6 primes.

## 4 Conclusion

We have presented an improved algorithm for the enumeration of plane meanders. The computational complexity of the algorithm for the problem of closed meanders is estimated to be $`2.5^n`$, much better than direct counting algorithms which have complexity $`12.26^n`$. Implementing this algorithm enabled us to obtain closed meanders up to order $`n=24`$, compared to the previous best of $`n=16`$. We also counted the number of open meanders up to order $`n=43`$ and semi-meanders up to order $`n=45`$. From our extended series we obtained precise estimates for the connective constants and critical exponents . This showed that a recent conjecture for the exact value of the semi-meander critical exponent is unlikely to be correct, while the conjectured exponent value for closed and open meanders is just consistent with the results from the analysis.

## E-mail or WWW retrieval of series

The series for the various generating functions so far generated with this algorithm and studied in  can be obtained via e-mail by sending a request to I.Jensen@ms.unimelb.edu.au or via the world wide web on the URL http://www.ms.unimelb.edu.au/~iwan/ by following the instructions.

## Acknowledgements

Financial support from the Australian Research Council is gratefully acknowledged.


text:
# Response of circumnuclear water masers to luminosity changes in an active galactic nucleus

## 1 Introduction

Luminous maser emissions in the $`6_{16}5_{23}`$ 22 GHz line of water have now been detected from more than 20 active galaxies, a list of which is given in Table 2 of the recent review by Moran, Greenhill & Herrnstein (2000). In almost every case where the size of the emission region has been determined interferometrically, the water maser emission has been found to originate primarily within a parsec- or subparsec-sized region that is coincident with the active nucleus of the galaxy, suggesting that active galactic nuclei are the energy sources that ultimately power the observed water maser emission. In several sources, the geometry is strongly suggestive of a circumnuclear accretion disk (e.g. Miyoshi et al. 1995).

Theoretical models for the collisional excitation of water masers (e.g. de Jong 1973, Neufeld & Melnick 1991) place requirements on the physical conditions in the emitting region. Luminous water maser emission is favored by molecular hydrogen densities, $`n(\mathrm{H}_2)`$, in the range $`10^710^{11}\mathrm{cm}^3`$, water abundances, $`x(\mathrm{H}_2\mathrm{O})=n(\mathrm{H}_2\mathrm{O})/n(\mathrm{H}_2)\stackrel{>}{}10^4`$, and gas temperatures $`\stackrel{>}{}400`$ K. A minimum temperature requirement is imposed both by the excitation requirement for the maser transition and by the fact that water abundances are significantly enhanced at temperatures $`\stackrel{>}{}400`$ K by the reaction $`\mathrm{O}+\mathrm{H}_2\mathrm{OH}+\mathrm{H}`$ followed by $`\mathrm{OH}+\mathrm{H}_2\mathrm{H}_2\mathrm{O}+\mathrm{H}`$ (e.g. Elitzur & de Jong 1978; Neufeld, Lepp & Melnick 1995). An additional requirement is that the gas temperature must significantly exceed the temperature of the radiation field at mid- to far-infrared wavelengths: thus some energy source that directly heats the circumnuclear gas in active galaxies is needed to explain the observed water maser emission. X-rays from the active nucleus (Neufeld, Maloney & Conger 1994, hereafter NMC94), shock waves within a circumnuclear accretion disk (Maoz & McKee 1998), and viscous heating (Desch, Wallin & Watson 1998, hereafter DWW98) have all been proposed as the energy source responsible for heating the masing gas.

Detailed models for water masers in X-ray heated gas have been presented by NMC94, who considered the water maser emission as a function of depth into an X-irradiated slab. Those models showed that the physical and chemical conditions within X-ray heated gas are conducive to the excitation of water masers, and that – given reasonable assumptions about the nature of the circumnuclear emission region – X-ray heating could account for the typical luminosities of water masers observed in AGN. Collison and Watson (1995, hereafter CW95) subsequently made the important point that for certain parameters the maser emission can be significantly enhanced by the absorption of non-masing far-infrared water lines by dust, a process not included by NMC94.

Since AGN are known to exhibit significant and rapid variability (e.g. Ulrich, Maraschi, & Urry 1997) at optical, ultraviolet, and particularly X-ray wavelengths (e.g. Mushotzki, Done & Pounds 1993), maser models that invoke X-ray heating naturally raise the question of how the water maser emission would respond to changes in the incident X-ray and bolometric fluxes. This question, to be addressed this Letter, is also strongly motivated by observations of cooordinated maser variability in NGC 1068 that have been reported recently by Gallimore et al. (2000). In §2 below, I present the results of a general parameter study for the equilibrium maser emission from a warm dusty medium. In §3, I consider how the maser emission responds to changes in the illuminating X-ray flux or in the dust temperature within the masing medium. In §4, the results are discussed with specific reference to the water masers in NGC 1068.

## 2 Water maser emission from a warm dusty medium

The collisional excitation of water masers within a warm dusty medium has been modeled by Wallin & Watson (1997) in the limit where the medium is optically thick to far-infrared continuum radiation. As previously pointed out by CW95, this limit often applies to circumnuclear gas in AGN, given the continuum dust opacities typically present (Draine & Laor 1993). The effect of continuum opacity is to place a lower limit on the escape probability for non-masing far-infrared transitions of water, allowing maser emission to occur even at great depths within a circumnuclear disk or gas cloud and enhancing the cooling rate due to far-infrared water emissions (DWW98).

In the limit where the medium is optically thick to far-infrared continuum radiation, the water level populations are determined solely by the local conditions and are dependent on the gas temperature, $`T_{\mathrm{gas}}`$, the dust temperature, $`T_{\mathrm{dust}}`$, the H2 density, $`n(\mathrm{H}_2)`$, and the water abundance $`x(\mathrm{H}_2\mathrm{O})=n(\mathrm{H}_2\mathrm{O})/n(\mathrm{H}_2)`$. Since water abundances $`10^4`$ are predicted for regions of luminous water maser emission, I have constructed a grid of models for fixed $`x(\mathrm{H}_2\mathrm{O})10^4`$ that span a broad region in the three dimensional parameter space defined by $`T_{\mathrm{gas}}`$, $`T_{\mathrm{dust}}`$, and $`n(\mathrm{H}_2)`$.

The water level populations are given by the usual equations of statistical equilibrium (e.g. CW95), with the rates for all radiative processes diminished by an escape probability $`\beta _{ij}^c=[2\alpha _{ij}/(1+2\alpha _{ij})](\mathrm{ln}[e+\pi ^{1/2}\alpha _{ij}^1])^{1/2}`$, where $`\alpha _{ij}`$ is the ratio of dust continuum opacity to line center opacity for the transition from $`i`$ to $`j`$ (Hollenbach & McKee 1979, CW95). For consistency with the previous work of CW95, we assume a continuum absorption coefficient of $`10^{23}n_\mathrm{H}\mathrm{cm}^2`$ for wavelengths $`\lambda 50\mu \mathrm{m}`$ and $`10^{23}n_\mathrm{H}(50\mu \mathrm{m}/\lambda )^2\mathrm{cm}^2`$ for $`\lambda 50\mu \mathrm{m}`$, where $`n_\mathrm{H}`$ is the density of H nuclei. For comparison, the probability of escape from a semi-infinite plane parallel slab is $`\beta _{ij}=(4\pi ^{1/2}\tau _{ij}[\mathrm{ln}\tau _{ij}]^{1/2})^1`$ in the limit where the line center optical depth $`\tau _{ij}1`$ (Hollenbach & McKee 1979). Comparing these expressions, we find that dust absorption dominates the effective escape of far-infrared water emissions for dust optical depths $`(8\pi ^{1/2})^1=0.07`$, corresponding to column densities of H nuclei, $`N_\mathrm{H}\stackrel{>}{}10^{22}\mathrm{cm}^2`$. Following Neufeld & Kaufman (1993, hereafter NK93), I solved the equations of statistical equilibrium for 179 rotational states of ortho-water, adopting the collisional rate coefficients of Green, Maluendes & McLean (1993) for transitions among the lowest 45 states and a simple extrapolation for transitions involving states of higher energy. The assumed line width, defined as velocity shift needed to reduce the opacity from its peak value by a factor e, was $`1\mathrm{km}\mathrm{s}^1`$.

Once the statistical equilibrium equations have been solved, the cooling rate due to water emissions and the rate of maser emission under conditions of saturation may be obtained straightforwardly (e.g. Neufeld & Melnick 1991; NK93). The results I obtained for the maser emissivity are in good agreement with those presented previously by Wallin & Watson (1997). If the gas is in thermal equilbrium, then the required heating rate is equal to the cooling rate, which – for the conditions of relevance here – is dominated by water cooling and gas-grain collisional cooling (Neufeld et al. 1995). Using the gas-grain cooling rate adopted by DWW98 and the water cooling rate obtained above, I computed the equilibrium heating rate, $`H`$, and used it in place of $`T_{\mathrm{gas}}`$ as one of the three dependent variables.

Figure 1a shows a contour plot of the gas temperature, $`T_{\mathrm{gas}}`$, and the maser efficiency, $`ϵ_{\mathrm{sat}}^0=h\nu \mathrm{\Phi }_{\mathrm{sat}}^0/(n_\mathrm{H}H)`$, as function of $`T_{\mathrm{dust}}`$ and $`H/n_\mathrm{H}`$. Here $`\mathrm{\Phi }_{\mathrm{sat}}^0`$ is the rate at which maser photons are emitted per unit volume, the superscript 0 indicating that the results apply to the case where the dust optical depth is large (Watson & Wallin 1997) and the subscript sat indicating that the maser emission is assumed to be saturated. The heating rate, $`H`$, is defined following the notation of Maloney, Hollenbach & Tielens (1996, hereafter MHT96) as the heating rate per hydrogen nucleus (with units erg s-1). The maser efficiency $`ϵ_{\mathrm{sat}}^0`$ is thus the fraction of the input power that emerges as 22 GHz maser radiation. The results shown in Figure 1a apply to a fixed molecular hydrogen density $`n(\mathrm{H}_2)=\frac{1}{2}n_\mathrm{H}=10^9\mathrm{cm}^3`$.

As required by the second law of thermodynamics, Figure 1 shows that $`T_{\mathrm{gas}}T_{\mathrm{dust}}`$ in the limit $`H0`$. There is a minimum value of the heating rate (defined by the curve $`ϵ_{\mathrm{sat}}^0=0`$) below which $`T_{\mathrm{gas}}T_{\mathrm{dust}}`$ is insufficient to support a population inversion . The vertical line at $`H/n_H1.2\times 10^{28}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1`$ represents the maximum X-ray heating rate that can be sustained without the gas being dissociated (NMC94). Thus the region for luminous maser action within an X-ray heated medium lies to the right of the curve $`ϵ_{\mathrm{sat}}^0=0`$, to the left of the line $`H/n_H=1.2\times 10^{28}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1`$, and above the curve $`T_{\mathrm{gas}}=400`$ K (the latter requirement being needed for consistency with the large water abundance of 10-4 that has been assumed). Figure 1b presents analogous results to those in Figure 1a, but now $`ϵ_{\mathrm{sat}}^0`$ and $`T_{\mathrm{gas}}`$ are shown for a fixed dust temperature (500 K) as a function of $`n(\mathrm{H}_2)`$ and $`H/\mathrm{n}_H`$.

## 3 Response of masers to a change in their environment

Results shown in Figure 1 apply to conditions of thermal equilibrium. The effects of changing the heating rate and dust temperature can now be considered. The timescale for the gas to re-establish thermal equilibrium after a change in the heating rate or dust temperature is given by $`\tau _{\mathrm{therm}}=\frac{5}{2}kT/2H=`$ $`0.027(T/1000\mathrm{K})(n(\mathrm{H}_2)/10^9\mathrm{cm}^3)^1([H/n_\mathrm{H}]/10^{28}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1)^1\mathrm{yr}.`$ Let us first consider the effects of changing the heating rate. On timescales greater than $`\tau _{\mathrm{therm}}`$, an increase in the heating rate increases the difference between the gas and dust temperatures and increases the rate of maser emission. The response is then conveniently characterized by the logarithmic derivative $`\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}H=1+(\mathrm{ln}ϵ_{\mathrm{sat}}^0/\mathrm{ln}H)`$, which is plotted in Figure 2 (dashed curve) for $`n(\mathrm{H}_2)=10^9\mathrm{cm}^3`$ and $`T_{\mathrm{dust}}=500`$ K. Except for small values of $`H`$ close to the threshold for maser action, the derivative $`\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}H`$ lies within a factor of 2 of unity, implying that changes in the heating rate are accompanied by comparable fractional changes in the maser emission rate.

The results plotted in Figure 2 are premised on assumptions of constant density and water abundance. The density is only expected to change on a dynamical timescale, which is typically much longer than $`\tau _{\mathrm{therm}}`$. The molecular fraction, $`2n(\mathrm{H}_2)/n_\mathrm{H}`$, only changes on the X-ray ionization timescale, which is given (MHT96) by $`4(n(\mathrm{H}_2)/10^9\mathrm{cm}^3)^1([H/n_\mathrm{H}]/10^{28}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1)^1\mathrm{yr}`$ and is a factor $`200`$ greater than the cooling timescale. Even if the molecular fraction is constant, however, the water abundance can change in response to changes in the temperature. This effect, the inclusion of which lies beyond the scope of the present study, is unlikely to be important except when the gas temperature lies close to the threshold for rapid water production ($`T_{\mathrm{gas}}400`$ K) or when the X-ray ionization rate lies close to the threshold for molecular dissociation.

If the heating rate varies, the gas temperature and maser emissivity $`\mathrm{\Phi }_{\mathrm{sat}}^0`$ cannot respond on timescales any shorter than $`\tau _{\mathrm{therm}}`$. However, $`\mathrm{\Phi }_{\mathrm{sat}}^0`$ can vary on timescales shorter than $`\tau _{\mathrm{therm}}`$ if the dust temperature changes rapidly. The thermal inertia of the dust grains is negligible, so the timescale on which the dust temperature can change is controlled by the diffusion time for the far-infrared continuum radiation111This is the shortest timescale on which the local infrared radiation field can vary; note that the infrared flux as viewed from afar may have a characteristic variability timescale that is considerably longer than that of the local infrared radiation as a result of the different light travel times from the near and far sides of the irradiated disk or torus., given by $`\tau _{\mathrm{rad}}\tau _{\mathrm{FIR}}\mathrm{}/c0.005[N_H/10^{24}\mathrm{cm}^2]^2[\mathrm{n}(\mathrm{H}_2)/10^9\mathrm{cm}^3]^1\mathrm{yr}`$, where $`\mathrm{}`$ is the distance below the irradiated surface of the cloud and $`\tau _{FIR}10(N_H/10^{24}\mathrm{cm}^3)`$ (assumed $`1`$) is the corresponding optical depth at far-infrared wavelengths.

Under circumstances where the dust temperature varies on timescales shorter than $`\tau _{\mathrm{therm}}`$ – as it might in response to variations in the incident bolometric flux from the AGN, for example – the gas temperature remains constant but the maser emissivity varies as a result of changes in the difference between the gas and dust temperatures, $`T_{\mathrm{diff}}=T_{\mathrm{gas}}T_{\mathrm{dust}}`$. The response is now characterized by the derivative $`(\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}T_{\mathrm{dust}})_{T_{\mathrm{gas}}}`$, the subscript $`T_{\mathrm{gas}}`$ indicating that the gas temperature is held constant. The dotted curve in Figure 2 shows the value of this derivative, again for the case where $`n(\mathrm{H}_2)=10^9\mathrm{cm}^3`$ and $`T_{\mathrm{dust}}=500`$ K. Two features of this curve are particularly noteworthy. First, the value of $`(\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}T_{\mathrm{dust}})_{T_{\mathrm{gas}}}`$ is always negative because $`T_{\mathrm{diff}}`$ decreases as $`T_{\mathrm{dust}}`$ increases. Second, $`(\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}T_{\mathrm{dust}})_{T_{\mathrm{gas}}}`$ can become much larger than unity for small values of $`H`$. This behavior occurs because the maser emissivity is primarily controlled by $`T_{\mathrm{diff}}`$; a small fractional change in $`T_{\mathrm{dust}}`$ causes a large fractional change in $`T_{\mathrm{diff}}`$ when $`T_{\mathrm{diff}}T_{\mathrm{dust}}`$.

## 4 Application to NGC 1068

The results presented in Figures 1 and 2 can now be discussed with reference to recent observations of the Seyfert galaxy NGC 1068. Gallimore et al. (2000) have recently reported that two maser features in NGC 1068, symmetrically placed in Doppler velocity about the systemic velocity of the host galaxy, flared simultaneously with a rise time of at most 22 days. Given the emission region structure elucidated by earlier VLBA observations (Greenhill & Gwinn 1997), the velocities of the maser features involved suggest that they are symmetrically located at distance $`d0.7`$ pc on either side of the nucleus. Unless the simultaneity of the flares is entirely happenstance, they were presumably coordinated by a signal that propagated from the nucleus.

The results presented in Figure 2 suggest two possible mechanisms for how the maser flares might be coordinated: (1) variations in the central X-ray luminosity, $`L_X`$, might simultaneously change the heating rates in the masing regions; or (2) variations in the intrinsic bolometric luminosity, $`L_{\mathrm{bol}}`$, might simultaneously change the dust temperature222For the case of saturated masers, an third mechanism is possible: in principle, the observed flares might be coordinated by a varying 22 GHz continuum flux from a compact central source. Since the maser radiation from the flaring regions travels tangentially toward us, the flux we receive from those regions might increase if the radial flux of seed radiation from the central source diminished. This explanation, however, would require that the maser emission is controlled by seed radiation from a compact continuum source; in NGC 1068 the observed distribution of maser emission close to the systemic velocity argues against that being the case.. Unfortunately, the average X-ray and bolometric luminosities are rather uncertain for NGC 1068, since the central source is obscured even for hard X-rays; I adopted the estimates of Pier et al. (1994), who compiled observations obtained over a wide range of wavelengths and assumed that 1% of the total intrinsic luminosity is reflected towards us. These estimates of the luminosities imply a 1 – 100 keV X-ray flux $`7\times 10^5\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ at the distance of the flaring masers and a characteristic dust temperature $`(L_{\mathrm{bol}}/4\pi d^2\sigma )^{1/4}500`$ K. Expressions given by MHT96 then imply an X-ray heating rate $`H5.6\times 10^{20}(N_\mathrm{H}/10^{24}\mathrm{cm}^2)^{0.9}\mathrm{erg}\mathrm{s}^1`$ at column density $`N_\mathrm{H}`$ below the irradiated surface of the maser emission region. For column densities greater than $`N_\mathrm{H}10^{24}\mathrm{cm}^2`$, Thomson scattering becomes important and increases the effective column density appearing in the above equation for $`H`$ by a factor $`\tau _T=6.65\times 10^{25}N_\mathrm{H}\mathrm{cm}^2`$.

In Figure 3, relevant constraints are plotted333The results plotted here apply to a gas cloud that is illuminated at normal incidence such that the column densities for X-ray shielding and for the escape of far-infrared radiation are identical. For oblique illumination at angle $`\theta `$ to the normal, the column density $`N_\mathrm{H}`$ refers to the shielding column density, the timescale $`\tau _{rad}`$ is reduced by a factor sec$`{}_{}{}^{2}\theta `$, and the dust temperature should be reduced by a factor $`(\mathrm{sec}\theta )^{1/4}`$. An additional assumption is that dust absorption dominates the removal of far-infrared line photons, requiring $`N_\mathrm{H}\mathrm{cos}\theta \stackrel{>}{}10^{22}\mathrm{cm}^2`$. as a function of $`n(\mathrm{H}_2)`$ and $`N_\mathrm{H}`$. Molecules are present above and to the right of the line $`H/n_\mathrm{H}=1.2\times 10^{28}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1`$. Maser action can occur below and to the left of the curve $`ϵ_{\mathrm{sat}}^0=0`$. The observed variability can be driven by changes in the X-ray heating rate for parameters to the left of the line $`\tau _{\mathrm{therm}}=0.05`$ yr; or alternatively by changes in the dust temperature for the region to the left of the line $`\tau _{\mathrm{rad}}=0.05`$ yr and to the right of the line $`\tau _{\mathrm{therm}}=0.05`$ yr. In either case, the minimum required H2 density for a variability timescale of 0.05 yr is $`10^8\mathrm{cm}^2`$. In Figure 4, the derivatives $`\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}L_X`$ and $`(\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}L_{\mathrm{bol}})_{T_{\mathrm{gas}}}`$ are shown for several values of the density. Since the dust temperature increases as $`L_{\mathrm{bol}}^{1/4}`$, the latter derivative is given by $`\frac{1}{4}(\mathrm{ln}\mathrm{\Phi }_{\mathrm{sat}}^0/\mathrm{ln}T_{\mathrm{dust}})_{T_{\mathrm{gas}}}`$.

Figures 3 and 4 indicate that the simultaneous flaring of the two maser features in NGC 1068 could plausibly be coordinated either through variations in the heating rate or through variations in the dust temperature. The former mechanism would be unique to X-ray heating whereas the latter might apply in principle even if shock heating or viscous heating were responsible for the maser emission. Observationally, the mechanisms might be discriminated by searching for correlations between the maser emission and X-ray and bolometric luminosities, although in the case of NGC 1068 the fact that we observe only reflected radiation from the nucleus will introduce uncertain delays that will be difficult to remove. Maser variability induced by a varying X-ray heating rate will be characterized by a positive correlation with $`L_X`$ (with a time delay determined by the different light travel times, of course), while variability induced by a varying dust temperature will be characterized by a negative correlation with $`L_{\mathrm{bol}}`$.

It is a pleasure to acknowledge helpful discussions with E. Agol and J. Gallimore. I gratefully acknowledge the support of a National Young Investigator award from the National Science Foundation.


text:
# Many-Body Coulomb Gauge Exotic and Charmed Hybrids

## Abstract

Utilizing a QCD Coulomb gauge Hamiltonian with linear confinement specified by lattice, we report a relativistic many-body calculation for the light exotic and charmed hybrid mesons. The Hamiltonian successfully describes both quark and gluon sectors, with vacuum and quasiparticle properties generated by a BCS transformation and more elaborate TDA and RPA diagonalizations for the meson ($`q\overline{q}`$) and glueball ($`gg`$) masses. Hybrids entail a computationally intense relativistic three quasiparticle ($`q\overline{q}g`$) calculation with the 9 dimensional Hamiltonian matrix elements evaluated variationally by Monte Carlo techniques. Our new TDA (RPA) spectrum for the nonexotic $`1^{}`$ charmed ($`c\overline{c}`$ and $`c\overline{c}g`$) system provides an explanation for the overpopulation of the observed $`J/\psi `$ states. For the important $`1^+`$ light exotic channel we obtain hybrid masses above 2 $`GeV`$, in broad agreement with lattice and flux tube models, indicating that the recently observed resonances at 1.4 and 1.6 $`GeV`$ are of different, perhaps four quark, structure.

Exotic hybrids, hadrons with quantum numbers not possible in simple $`q\overline{q}`$ or $`qqq`$ quark models, have been an elusive, yet signature prediction of Quantum Chromodynamics (QCD). It was therefore quite natural that the recent observation by the E852 collaboration  of two exotic $`J^{PC}=1^+`$ states with masses 1.4 and 1.6 $`GeV`$ would attract widespread interest. Since these states have isospin $`I=1`$ they can not be glueballs (oddballs) and would initially appear to be viable hybrid meson candidates, especially since vintage bag model calculations  predict exotic excitations with explicit gluonic degrees of freedom in this mass range. However, the detailed structure of these states remains uncertain since most contemporary theoretical studies, such as lattice gauge , flux tube , QCD sum rule , perturbative non-relativistic QCD  and constituent models  have focused on heavy quark hybrids. The only two modern light quark hybrid calculations, a lattice gauge  and flux tube , have further compounded this uncertainty by predicting the lightest exotic hybrid mass to be about $`2.0`$ $`GeV`$– significantly above the observed $`1^+`$ states. Because the bag model results are rather dated and lattice calculations are less accurate for light quarks due to extrapolation, it is important to have an additional, alternative hybrid prediction. The purpose of this Letter is to determine if these exotic states can indeed be interpreted as hybrids within a relativisitc many-body constituent approach that has successfully described both conventional meson  and glueball  systems.

Our starting point is the QCD Coulomb gauge Hamiltonian (see for example ref. ) which we simplify to a form amenable for many body calculations

$`H={\displaystyle 𝑑𝐱\mathrm{\Psi }^{}(𝐱)(i\alpha +\beta m)\mathrm{\Psi }(𝐱)}+Tr{\displaystyle 𝑑𝐱(𝚷^a𝚷^a+𝐁_A^a𝐁_A^a)}`$ (1)
$`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐱𝑑𝐲\rho ^a(𝐱)V(|𝐱𝐲|)\rho ^a(𝐲)}.`$ (2)
Here $`\mathrm{\Psi }`$ and $`𝐀`$, are the respective quark and gluon fields, $`𝐁_A^a=\times 𝐀^a`$, and $`\rho ^a=\mathrm{\Psi }^{}T^a\mathrm{\Psi }+f^{abc}𝐀^b𝚷^c`$ is the quark plus gluon color density. The current quark mass, $`m`$, is assigned the values, $`m_u=m_d=5MeV`$ and $`m_c=1200MeV`$ for the $`u,d`$ and $`c`$ flavors, respectively. Confinement and leading canonical interactions are represented by the instantaneous potential, $`V=\frac{\alpha _s}{r}+\sigma r`$, with $`\alpha _s=.2`$, and, $`\sigma =0.135`$ $`GeV^2`$, as determined by the string tension from lattice and Regge fits. We also use a cut-off parameter $`\mathrm{\Lambda }=45GeV`$ to regularize the logarithmic divergent term in the mass gap equation. The model parameters $`\sigma `$, $`\alpha _s`$ and $`\mathrm{\Lambda }`$ are commensurate with our previous pure quark  and gluon  applications which produced reasonable hadronic descriptions including the Regge trajectory slopes for the mesons (e.g. $`\rho `$ tower) and glueballs (pomeron) .

Next we proceed to the many-body diagonalizations but first perform a canonical transformation (BCS rotation) to a new quasiparticle basis

$`\alpha _i^a(𝐤)=\mathrm{cosh}\mathrm{\Theta }_ka_i^a(𝐤)+\mathrm{sinh}\mathrm{\Theta }_ka_i^a(𝐤)`$ (3)
$`B_{c\lambda }(𝐤)=\mathrm{cos}{\displaystyle \frac{\theta _k}{2}}b_{c\lambda }(𝐤)\lambda \mathrm{sin}{\displaystyle \frac{\theta _k}{2}}d_{c\lambda }^{}(𝐤)`$ (4)
$`D_{c\lambda }(𝐤)=\mathrm{cos}{\displaystyle \frac{\theta _k}{2}}d_{c\lambda }(𝐤)+\lambda \mathrm{sin}{\displaystyle \frac{\theta _k}{2}}b_{c\lambda }^{}(𝐤),`$ (5)
where $`\mathrm{\Theta }_k`$, $`\theta _k/2`$ are the BCS angles, further specified below, and $`a(\alpha ),b(B)`$ and $`d(D)`$ are bare (dressed) gluon, quark and antiquark Fock operators, respectively. The indices $`a=1,2\mathrm{}8`$ and $`c=1,2,3`$ denote color while $`\lambda `$ represents spin projection. The new field expansions are

$`A_i^a(𝐱)={\displaystyle \frac{d𝐤}{(2\pi )^3}\frac{1}{\sqrt{2\omega _k}}[\alpha _i^a(𝐤)+\alpha _i^a(𝐤)]e^{i𝐤𝐱}}`$ (6)
$`\mathrm{\Pi }_i^a(𝐱)=i{\displaystyle \frac{d𝐤}{(2\pi )^3}\sqrt{\frac{\omega _k}{2}}[\alpha _i^a(𝐤)\alpha _i^a(𝐤)]e^{i𝐤𝐱}},`$ (7)
for the gluon fields and

$$\mathrm{\Psi }(𝐱)=\underset{c\lambda }{}\frac{d𝐤}{(2\pi )^3}\left[U_{c\lambda }(𝐤)B_{c\lambda }(𝐤)+V_{c\lambda }(𝐤)D_{c\lambda }^{}(𝐤)\right]e^{i𝐤𝐱},$$
(8)
for the fermion field. The rotated Dirac spinors are given in terms of the Pauli spinors, $`\chi `$,

$$U_{c\lambda }(𝐤)=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\sqrt{1+\mathrm{sin}\varphi _k}\chi _{c\lambda }\\ \sqrt{1\mathrm{sin}\varphi _k}\sigma \widehat{𝐤}\chi _{c\lambda }\end{array}\right],V_{c\lambda }(𝐤)=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\sqrt{1\mathrm{sin}\varphi _k}\sigma \widehat{𝐤}\chi _{c\lambda }\\ \sqrt{1+\mathrm{sin}\varphi _k}\chi _{c\lambda }\end{array}\right].$$
(9)
We then find an improved, nontrivial vacuum, $`|\mathrm{\Omega }`$, by minimizing the ground state expectation value of the Hamiltonian variationally with respect to the BCS angles. Actually, the specific variational parameters are the quark gap angle, $`\varphi _k`$, related to the BCS angle by $`tan(\varphi _k\theta _k)=m/k`$, and the gluon self-energy, $`\omega _k`$, satisfying $`\omega _k=ke^{2\mathrm{\Theta }_k}`$. This generates a quark and gluon gap equation (equivalent to the Schwinger-Dyson equation) yielding mass gaps of about 100 $`MeV`$ for the $`u/d`$ quarks and 800 $`MeV`$ for the gluon. The BCS vacuum contains quark and gluon condensates (Cooper pairs) in reasonable agreement with QCD sum rules. For more complete details consult refs. .

We find the two body sector is adequately described by the Tamm-Dancoff approximation (TDA) in which a glueball is represented by Fock states $`g^{}g^{}|\mathrm{\Omega }`$ and mesons by $`q^{}\overline{q}^{}|\mathrm{\Omega }`$. The notable exception is the light pseudoscalar sector, where a more sophisticated, collective approximation (Random Phase approximation or RPA) is needed to correctly reproduce the Goldstone boson nature of the pion due to spontaneous chiral symmetry breaking by our BCS vacuum. This is also further documented in refs. .

Finally we formulate the hybrid meson as $`|hybrid=q^{}\overline{q}^{}g^{}|\mathrm{\Omega }[[q^{}\overline{q}^{}]_8g^{}]_0|\mathrm{\Omega }`$, now involving quark color octet states. The resulting TDA equation for the hybrid mass $`M`$ is

$$hybrid|[H,q^{}\overline{q}^{}g^{}]|\mathrm{\Omega }=Mhybrid|q^{}\overline{q}^{}g^{}|\mathrm{\Omega }.$$
(10)
This projects the hybrid meson wave equation, which is pictorially represented in Fig. 1, onto the three body Fock basis. Unlike our pion application, the $`q\overline{q}`$ pair is now in a color octet and the TDA is sufficient since the Hamiltonian, as well as exact QCD, does not conserve the chiral color octet current (we specifically calculated no difference in the more elaborate, chiral symmetry preserving hybrid RPA calculation, see below). The relevant angular momenta ($`am`$) are the $`q,\overline{q}`$ spins coupled to an intermediate $`S`$, and the gluon spin with its orbital $`am`$ $`L_+`$ (with respect to the $`q\overline{q}`$ cm) coupled to intermediate $`l`$. Coupling $`l`$ with $`L_{}`$ (the orbital $`q\overline{q}am`$) yields $`L`$ which combines with $`S`$ giving the total $`amJ`$. The complete wavefunction in the hybrid cm has form

$`F_{\lambda _g\lambda _q\lambda _{\overline{q}}}^{JPC}(𝐪_+,𝐪_{})={\displaystyle \underset{lL_{}L_+LSm_+m_{}}{}}F_{lL_{}L_+LS}^{JPC}(|𝐪_+|,|𝐪_{}|)Y_{L_+}^{m_+}(\widehat{𝐪}_+)Y_L_{}^m_{}(\widehat{𝐪}_{})`$ (11)
$`(1)^{\lambda _g}L_+m_+1\lambda _g|lm_lL_{}m_{}lm_l|Lm_L{\displaystyle \frac{1}{2}}\lambda _q{\displaystyle \frac{1}{2}}\lambda _{\overline{q}}|Sm_S(1)^{\frac{1}{2}\lambda _{\overline{q}}}Lm_LSm_S|Jm_J,`$ (12)
where $`𝐪_{}`$ and $`𝐪_+`$ are the respective relative momentum of the $`q\overline{q}`$ pair and gluon (with respect to the pair cm). We then impose the transversality condition, $`\widehat{𝐤}\alpha (𝐤)=0`$, from the Coulomb gauge constraint which eliminates states with $`L_+=1`$ and $`l=0`$. For pure $`S`$ waves the lightest hybrid states will then have $`J^{PC}=`$ $`1^+`$, $`0^{++}`$, $`1^{++}`$ and $`2^{++}`$. These are nonexotic states which will mix with conventional mesons and hinder hybrid identification. For exotic states one $`P`$ wave is necessary and we calculate the lightest corresponds to $`L_+=1`$ since the $`L_{}=1`$ excitation is energetically more expensive due to quark repulsion in the octet channel. This generates the exotic states $`1^+`$, $`3^+`$ and $`0^{}`$.

Instead of solving the formidable TDA nonlocal equations (effectively a 12-dimensional problem in momentum space), we evaluate the hybrid mass variationally using an exponential radial wavefunction for each of the two independent momentum variables. In the center of momentum frame the matrix elements reduce to 9-dimensional integrals that we evaluate numerically using the Monte Carlo code VEGAS. We then perform searches for minima on the energy surface in the different angular momentum channels. Our final results and key findings of this Letter are displayed in Figs. 2 and 3.

Note from Fig. 2 the clear agreement between our predictions and the lattice and flux tube results for both light and charmed hybrid states. This agreement sharply contrast with the BNL measurements which strongly suggests that the observed exotic states are not hybrids. To confirm our result is not an artifact of the variational method, we have reproduced our conventional meson and glueball exact TDA spectra to within a few percent. We also varied the least constrained model parameter $`\alpha _s`$, from $`0.2`$ to $`0.4`$. Related, we even performed a more extensive RPA variational calculation and took the chiral limit ($`m_u=m_d0`$), finding only miniscule change in the hybrid mass, consistent with the nonconservation of the chiral color octet charge discussed above. The culmination of our model sensitivity study produced at most a 10 $`\%`$ hybrid mass variation indicated by the box in Fig. 2.

Since four quark states $`q\overline{q}q\overline{q}`$ can also have exotic quantum numbers, one can make simple estimates yielding exotic masses between 1 and 2 $`GeV`$ for quarks in color singlet configurations. This is consistent with a recent unitary quark model calculation  which also concluded that the observed $`1^+`$ states are indeed predominately meson-meson resonances.

We also calculated the lightest non-exotic hybrid (ground state) to have mass slightly above 2 $`GeV`$. It is interesting to note that this state has $`J^{PC}=1^+`$ in contrast to certain heavy hybrid, quenched lattice ground state results  which find near degeneracies among several negative parity states. This was first noted by ref.  who used a similar constituent model in the heavy (static) quark limit. However, the heavy quark lattice simulations do not include quark spin and until this is included the ground state hybrid quantum numbers, as well as related level ordering, remains unclear . Further, another lattice calculation  studying orbital hybrid excitations concluded that the ground state quantum numbers were likely to be $`1^+`$ in agreement with our work. This study noted that the degeneracies in the static framework would be broken by quark spin-orbit effects which shift the exotic levels and also mix nonexotic and conventional $`q\overline{q}`$ states. They were also able to infer the level splitting, yielding the above conclusion, using propogating quarks on the lattice. There is less, but still some, uncertainty in level order within our model as well since our spin dependent interaction is incomplete. Although further level ordering study is necessary, the thrust of this work is the exotic hybrid and improved spin effects will not alter our conclusion that the lightest is above 2 $`GeV`$.

In Fig. 3 we compare our full model spectrum to data for the well studied, believed to be gluon rich, $`1^{}`$, $`J/\psi `$ system to provide an explanation for the anomalous overpopulation of observed states  with respect to quark model predictions. Whereas previous constituent calculations, using only $`S`$ waves, could only account for 3 of the known 6 charmonium levels, we now predict 7 $`c\overline{c}`$ states in addition to 4 $`c\overline{c}g`$ hybrids. Further, ref.  lists an additional charmonium level $`\psi (3836)`$ assigned $`J^{PC}=2^{}`$ which also agrees well with our D wave prediction (not shown). Notice that by simply including $`D`$ waves we have resolved the ”overpopulation” problem. In general all of the $`1^+`$ states, both $`c\overline{c}`$ and $`c\overline{c}g`$, will mix and a more elaborate calculation is in progress. However, our current result is already sufficient to conclude that simple level counting (density of states) arguments will probably not be effective in identifying charmed hybrid states.

Finally, we mention a novel color octet effect leading to an isospin splitting since it only affects the $`I=0`$ states. This is the annihilation process depicted in Fig. 1 corresponding to $`q\overline{q}gq\overline{q}`$ for the $`L_{}=0`$, spin aligned color octet quark pair. Octet quarkonium is the QCD analogue to ortho positronium and the annihilation interaction raises all $`I=0`$ light hybrid states by roughly $`300MeV`$ when the $`q\overline{q}`$ spins are aligned ($`S=1`$).

Summarizing, our large-scale diagonalizations of an effective Coulomb gauge Hamiltonian provide a reasonable, comprehensive description of the meson, lattice glueball and lattice hybrid meson spectra. Further, our composite $`J/\psi `$ spectrum is now also in much better agreement with data, especially in terms of density of states. It is important to note that our quark/gluon unified approach essentially entails only one pre-determined dynamical parameter. Finally, and perhaps most significant, our reaffirmation of lattice and flux tube $`1^+`$ masses indicates that the recently observed exotic states are not hybrids. Based upon preliminary estimates and other independent studies it is more likely that these resonances are four quark states and more rigorous, higher quark Fock state calculations are in progress.

We thank NERSC for providing Cray J-90 CPU time. F. L. E. acknowledges SURA-Jefferson Lab for a graduate fellowship. This work was partially supported by grants DOE DE-FG02-97ER41048 and NSF INT-9807009.