text:
# Associative Recall in Non-Randomly Diluted Neuronal Networks

## 1 Introduction

A number of mathematical concepts and tools have been extensively applied in order to construct models and simulations of neuronal behavior that could shine some light into the intricate workings of the brain. In physics, particular importance and interest have been placed on Hopfield models because of their potential for associative recall and close relationship with statistical mechanics . However, while most such approaches have considered fully connected structures, experimental evidences indicate that biological neuronal networks involve a combination of intense local interconnections (responsible, among other things, for lateral inhibition ) coexisting with long range mappings which often exhibit topographical structure, in the sense that the connections preserve the spatial adjacency of the signal. Studies addressing partially connected Hopfield networks on lattices have appeared sporadically in the literature . More recently , the Barabási-Albert  connecting scheme, where the “rich gets richer”, was used to define Hopfield architectures with symmetric couplings whose potential for associative recall was assessed with respect to the number of connections between the constituent cells. The obtained results indicated that, despite the relatively sparse and heterogeneous connections, for suitable parameters such architectures preserved a good deal of the associative recall when compared to fully connected counterparts, bridging the gap to further investigations.

The present work reports the extension of the models described in  in three aspects that are specially relevant from the biological point of view, namely the asymmetry that neuron $`i`$ may influence neuron $`k`$ even if neuron $`k`$ has no influence on neuron $`i`$, the spatial distribution of cells along the input space (assumed to be a one-dimensional vector), and the spatial distribution of the connections between those cells. More especifically, as the network is constructed, each added cell is placed preferentially at regions of the input space that are denser in cells, in such a way that dense regions tend to become denser. Being diluted, not every bit of the input pattern will be connected to a neuron. Also, the connections between a cell $`i`$ and the rest of the network are established in such a way that cells that are closer to each other have higher probability to interconnect, which is achieved through Monte Carlo simulation assuming normal weighting centered at cell $`i`$. The potential of the networks for associative recall is experimentally assessed with respect to several parameter configurations, with special attention given to the effect of the uniformity of the cell spatial distribution over the obtained performance.

The article starts by revising some key aspects of biological neuronal structures and proceeds by presenting the obtained results and respective discussion, with special attention placed on the biological interpretation and implication of the observed behavior in terms of the parameter configurations.

## 2 Spatial Organization in Biological and Computational Networks

First, in what we call the asymmetric Hopfield-Barabási-Albert model, we put the Hopfield model onto the scale-free Barabási-Albert  network. In contrast to spin glasses and many physics systems, the interactions between biological neurons are not symmetric: Neuron $`k`$ may influence neuron $`i`$ even if neuron $`i`$ has no influence on neuron $`k`$. In other words, the synaptic strengths $`J(i,k)`$ are not symmetric, i.e. $`J(i,k)`$ is not equal to $`J(k,i)`$. We thus construct our model in such a way that every neuron $`i`$ selects itself plus exactly $`m`$ other neurons (repeated connections are allowed) as neighbours. This fixed number $`m+1`$ of influencing neighbours greatly simplifies the algorithm.

Second, in what we call the geometrical model, we try to incorporate the geometry instead of only the topology of neural nets. Another characteristic feature exhibited by the mammalian cortex is the spatial structuring permeating both the input signal and neuronal architectures. For instance, the visual field is mapped onto the retina space, which innervates, through the lateral geniculate nucleus, into the primary visual cortex in such a way that the spatial adjancency of the input patterns is maintained. The importance of such *topographical mappings* is corroborated by the fact that they are found ubiquitously along the cortex, and not only in its primary sensory regions. It should be observed that such mappings are not isometric. One important neuronal pheonomenon closely related to topographical mappings is known as *lateral inhibition*, where the activity at a given cell is inhibited by its spatial neighbours along the input space. In spite of its key role in cortical organization, relatively few neuronal networks have explicitly taken spatial organization into account, Kohonen’s SOM model  being one of the pioneering exceptions. At the same time, much of the efforts devoted to Hopfield models have systematically overlooked the spatial structure of input and topographical mappings. The situation is similar in the biological community, where only more recently (e.g. ) attention was focused on this important organizational principle. The fact that neuronal cells that are closer together tend to have increased chance to interconnect (e.g. ) has also been often overlooked. Indeed, cortical organization seems to combine local, intensive connections with long distance, sparser, mappings between distinct neuronal modules. Such local/global organization may also underly other spatial scales structures in the brain, possibly leading to scale-free characteristics over limited spatial scale intervals. In the current work, we propose a novel variation of the Hopfield model for associative recall that explicitly takes into account the spatial organization of the neuronal cells and the localized and global spatial distribution of their connections.

The input space is assumed to be the one-dimensional vector $`I(i)`$, $`i=1,2,\mathrm{},`$ $`L`$, to which a total of $`N`$ neuronal cells are attached. The number of connections made by each cell is fixed and equal to $`m`$. The network is diluted, in the sense that $`R=N/L<1`$. Figure 1 illustrates a typical realization of the assumed neuronal architecture for $`N=7`$, $`m=3`$ and $`L=15`$, with the respective asymmetric weight matrix shown in Equation ( 1). Each line $`i`$ shows the connections of neuron $`i`$ to other neurons by a ’1’. The construction of the network involves the two following steps: (a) neurons are placed at spatial positions $`p0,1,\mathrm{},L`$ according to a probability density function $`h(p)`$ (cells can only be added at empty sites); and (b) $`m`$ not necessarily distinct synaptical connections are established between the just added cell and its neighbours according to a normal (Gaussian) density function $`g_\sigma (d)`$, where $`d`$ is the distance between the position $`p`$ of the reference cell and each of the other cells in the network and $`\sigma `$ is the standard deviation of the normal distribution. It is assumed throughout this article that $`L=1000`$.

$$J=\left[\begin{array}{ccccccc}1& 1& 1& 0& 1& 0& 0\\ 1& 1& 0& 1& 1& 0& 0\\ 0& 1& 1& 1& 1& 0& 0\\ 1& 0& 1& 1& 0& 1& 0\\ 0& 0& 0& 1& 1& 1& 1\\ 0& 0& 1& 1& 0& 1& 1\\ 0& 0& 0& 1& 1& 1& 1\end{array}\right]$$
(1)
A non-stationary (along time) density function $`h(p)`$ is obtained at each interaction during the network construction by adding the function $`f(i)=1`$ if $`|ip|a`$ and $`0`$ otherwise to the previous instance of $`h(p)`$, where $`p`$ is the position of the most recently added cell, and renormalizing $`h(p)`$. Standard Monte Carlo sampling is then performed over the respective distribution function so as to select the position $`p`$ for a new cell. The probability function $`h(p)`$ is initially set as uniform. Monte Carlo simulation is also employed to select $`m`$ cells for symmetric connections with the current cell. The connections are performed as follows: (i) a vector $`v(i)`$ is built such that $`v(i)=1`$ at the positions $`i`$ where neurons exist and $`0`$ otherwise; (ii) this vector is weighted according to a normal density function centered at $`p`$, the position of the current cell; and (iii) Monte Carlo is used to select $`m`$ sites of $`v`$. It is observed that all cells are always self-connected. Figure 2 illustrates typically obtained connection matrices. The reinforced diagonal structure of such matrices favors conditioning and reflects the locality of the synapses at specific spatial scales defined by $`\sigma `$.

For both the Hopfield-Barabási-Albert and the geometrical model, the neurons $`i,i=1,2,\mathrm{}N,`$ in our models are either firing ($`S_i=1`$) or silent ($`S_i=1`$) and are updated according to

$$S_i\mathrm{sign}(\underset{k}{}J_{ik}S_k)$$

with synaptic strengths $`J_{ik}=_\mu \xi _i^\mu \xi _k^\mu `$ (Hebb rule) if $`i`$ and $`k`$ are connected. Here $`\xi _i^\mu =\pm 1,\mu =1,2,\mathrm{}P,`$ are $`P`$ random bit-strings called patterns, and one of them is supposed to be recalled by this updating rule if it is presented to it in a perturbed fashion: Associative memory. The quality of recall is measured by the overlap $`\mathrm{\Psi }=_iS_i\xi _i^1`$ if, without loss of generality, the first pattern is supposed to be recovered.

## 3 Results

For the asymmetric Hopfield-Barabási-Albert model, $`N`$ binary neurons are added to the initial core of $`m+1`$ neurons, and $`P`$ binary patterns are stored. Initially, the first of these patterns is presented to the network with ten percent of the sites flipped (overlap = 0.8), and after a few iterations the dynamics comes to a fixed point for which the overlap $`\mathrm{\Psi }`$ is determined. Fig. 3 shows as in the symmetric case that for $`PmN`$ the desired pattern is fully restored. However, Fig. 4 shows that the scaling law: $`\mathrm{\Psi }=f(m/P)`$ for infinite $`N`$ is not fulfilled well, in contrast to the symmetric case. The distribution of the number of neurons, influencing $`q`$ sites each, follows roughly a $`1/q^3`$ law for small $`q`$ but not for large $`q`$ (not shown).

For the geometrical model, simulations were performed considering combinations of $`m=25`$ and $`50`$ with $`a=10`$ and $`100`$. Parallel updating is adopted for Hopfield recovery and the number of trained patterns is $`P=0.05N`$. Relative overlaps are used , i.e. $`R=\mathrm{\Psi }/N`$. The relative overlaps are calculated taking into account only the neurons attached to the input pattern.

Figure 5 shows the correlation integral $`C(r)`$ , namely the number of neurons with distance along the input space smaller or equal to $`r`$ divided by $`N^2`$. Figures 5 to 7 show examples of typically obtained density and respective distribution functions obtained for $`a=10`$ (a) and $`100`$ (b) and considering 10 different realizations of the spatial distribution of cells for each parameter settings.

## 4 Discussion

For the Hopfield-Barabási-Albert model, the asymmetric version gave results similar to the symmetric version , provided self-interaction in taken into account in both cases. For the geometrical model, as indicated by Fig. 5, the considered scheme for positioning the neurons along the input space led to similar power law behavior for $`a=10`$ and $`100`$ for the considered correlation lags. At the same time, as is clear from Fig. 6, the value of $`a`$ strongly influences the characteristic spatial distribution of neurons, in the sense that shorter clusters of cells, indicated by the clumps in Fig. 6, are obtained for smaller values of $`a`$. As expected, higher values of $`m`$ tend to produce more widespread distribution of synaptic connections, leading to the wider dispersions around the main diagonal of the matrices in Fig. 2(b) and (d) and substantially superior potential for associative recall, shown in Fig. 7. This is a consequence of the fact that a larger number of connections has direct impact over the amount of memory in the network. At the same time, the network performance tends to improve as $`a`$ is increased (see Fig. 7). A more uniform spatial distribution of the neurons is obtained for larger values of $`a`$, with the limit case $`a\mathrm{}`$ tending to the uniform distribution. It is clear that, at least for the considered kind of stimuli and noise (uniform distribution of signal changes), the effect of irregular distribution of neurons along the input pattern space tends to deteriorate associative recall. The important biological implication is that the development (ontogeny) of real neuronal systems coping with uniformly distributed stimuli and noise has to incorporate mechanisms for ensuring uniform distribution of neuronal placement. In addition to topographical maps, other possible biological mechanisms that could be involved in such spatial tuning may involve controled apoptosis (i.e. programmed neuronal cell death) and chemotatic migration modulated by density. It is important to bear in mind that different results could be obtained for distinct input stimuli and noise models. For instance, stimuli characterized by specific correlation lengths could be better processed by networks with neurons and connections presenting congruent spatial characteristis.

All in all, the present work has considered extensions of the classical Hopfield network which are characterized by enhanced biological plausibility as far as asymetric connections and spatial distribution of cell bodies and synapses are concerned. Such models were analysed according to the shape-function paradigm, i.e. by assessing how the model parameters influenced the associative recal performance, and a series of perspectives has been opened for further investigation. A particularly interesting issue being currently addressed regards the characterization of the effect of the individual morphology of real neuronal cells over the network behavior, as well as the influence of stimuli and noise spatial features over the respective network performance.

### Acknowledgments

Luciano da F. Costa is grateful to the Human Frontiers Science Program, FAPESP (99/12765-2) and CNPq (468413/00-6 and 301422/92-3) for financial help. D. Stauffer thanks H. Sompolinsky for suggesting asymmetric networks and GIF for travel support.


format: text

text:
# Stream eclipse mapping with ’fire-flies’

## 1. Introduction

Eclipse mapping is an inversion technique that can be used to reconstruct the distribution of bright material in the accretion streams of eclipsing polars. The eclipsing nature of these systems means that for a short duration the stream is the dominant contributor to the observed brightness, the white dwarf and accretion region being hidden from view by the secondary. The technique of eclipse mapping was first applied to accretion disks in the early 1980s, and the method was later applied to the bright eclipsing polar HU Aqr by Hakala (1995). Because it was found that the accretion stream contributed around half the optical emission in HU Aqr, this system was ideal for the application of eclipse mapping to determine the distribution of this bright stream material in the binary system.

Hakala (1995) developed this method assuming a one-dimensional stream confined to the orbital plane. This was represented as an arc connecting the inner Lagrangian point ($`L_1`$) to the white dwarf. Developments to this method included a two part trajectory, and models using this have been successfully applied to the eclipse light curves of the bright eclipsing polars HU Aqr and UZ For (Bridge et al. 2002; Vrielmann & Schwope 2001; Harrop-Allin, Potter & Cropper 2001; Kube, Gänsicke & Beuermann 2000; Harrop-Allin et al. 1999b). This trajectory comprises a ballistic free-fall component from the $`L_1`$-point (Lubow & Shu 1975), followed at some distance from the white dwarf by a magnetically confined component. The confined part of the trajectory originates in the threading region from where the stream material is channeled by the field lines to the white dwarf surface. The field lines are assumed to have a dipolar geometry and be centred on the white dwarf. However, such a clear division between two trajectories is likely to be an inadequate assumption (e.g. Heerlein, Horne & Schwope 1999) and the assumption of a ballistic free-fall component may not be appropriate in all cases (e.g. the low state of HU Aqr; Harrop-Allin et al. 2001). A further development is a move to three-dimensional streams (Vrielmann & Schwope 2001; Kube et al. 2000), thus attempting to include the effects of uneven heating of the stream and features related to the eclipse by the stream of the white dwarf and accretion region.

All these previous model techniques have the common feature that the stream trajectory is fixed prior to the application of the model. Hakala, Cropper & Ramsay (2002) introduced a technique which makes no assumptions about the location of the stream material prior to the modeling process. In principle, emission is allowed from anywhere within the Roche lobe of the white dwarf. They have successfully applied this technique to synthetic data sets, but here we apply it to observed light curves for the first time.

We model the ‘white’ light curves of the eclipsing polars EP Dra and HU Aqr (see Bridge et al. 2002, 2003 for a discussion of the light curves), observations of which were obtained with the superconducting tunnel junction (STJ) camera S-Cam 2 (see e.g. Perryman et al. 2001). The S-Cam 2 camera records the location on the array, time of arrival and energy of each incident photon. This allows the study of the energy dependence of short timescale variations observed in polars.

## 2. Model technique

We discuss briefly here the main points of the model (for a detailed discussion see Hakala et al. 2002). The model consists of a number of ‘fire-flies’, where each fly is an individual bright emission point. A fly is created with an initially random location within the Roche lobe of the white dwarf, but is subsequently free to move during the evolution of the model. The fly has a brightness ($`F_{fly}`$) which varies with the angle $`\alpha `$ between the observer, the fly and the white dwarf. This is defined as $`F_{fly}(\alpha )=F_0+Acos(\alpha )`$ where $`F_0`$ is the minimum brightness of a fly and $`A`$ is the amplitude of angular dependence. The brightness of those flies not eclipsed by the secondary at each phase are summed to form the model light curve. Many swarms of such flies are created and for each a model light curve can be generated.

The goodness of fit of an individual model light curve to the observed light curve is then defined using a fitness function, $`R`$, which consists of a $`\chi ^2`$ term and a regularisation term: $`R=\chi ^2+\lambda S_{reg}`$. The best fit solution is ‘bred’ through a number of generations using a genetic algorithm (GA; see e.g. Charbonneau 1995 for a review) that is regulated by a self-organising map (SOM; Kohonen 1990). The regularisation places a number of sections of a smooth curve through the fly distribution (see Hakala et al. 2002), this helps to constrain the model where generally there are more degrees of freedom than data points. The two end points of the curve are defined as the $`L_1`$-point and the white dwarf surface, thus providing two physically realistic constraints on the stream trajectory. The model then prefers those swarms with minimal distance to the curve, resulting in smooth tube-like fly distributions. Hakala et al. (2002) also demonstrated the use of the technique with no regularisation term. Here there is no computational pressure on the flies to create any form of tube and they are not constrained to lie between the white dwarf and the $`L_1`$-point. In this case the model will define those regions of the Roche lobe volume where emission can and cannot be.

## 3. Results

The light curves were obtained in October 2000 at the William Herschel Telescope, La Palma. Preliminary results for fits to the ‘white’ light curves of EP Dra and HU Aqr are shown in Figure 1, the light curves are binned into 3 s time intervals, with the model fits superimposed as a solid line. We used 500 swarms each with 200 flies, and the models are run for 500 generations. This provides a sufficient number of swarms to produce a unique solution and enough iterations to reach convergence. Figures 2 and 3 show the corresponding fly distributions.

The fitting process is necessarily a trade-off, mediated by the value of $`\lambda `$, between fitting real gradient changes and avoiding statistical noise features. The presence of noise in the data will have the effect of broadening the fly distribution, but we can attempt to minimise this by adjusting the $`\lambda `$ term to place more emphasis on the regularisation. However, while increasing $`\lambda `$ has the effect of smoothing the distribution of flies into a more tube-like structure, we must be careful not to place too much emphasis on $`\lambda `$ as this will have the effect of smoothing out any real brightness variations.

Figure 2 shows that the emission points in HU Aqr are located mostly out of the orbital plane. The extent in volume of the emission points, and location at different distances in the z-direction may indicate that material is threaded along many different field lines, and may be the cause of the lack of evidence for a bright threading region at any one place in the orbital plane. In contrast to HU Aqr, the emission points in EP Dra are located close to the orbital plane, possibly indicating accretion along field lines that meet the white dwarf surface at low latitudes. The emission in EP Dra is also distributed over a larger volume than HU Aqr, particularly the region where we expect coupling to field lines to occur. To some extent this is a product of a noisier light curve, but it could also indicate material located in a large stagnation region, where it collects before coupling to the white dwarf field lines.

## 4. Discussion

We have applied this more objective technique to real data in an attempt to make the least number of assumptions about the location of bright material in these eclipsing binary systems. We find that previous assumptions are largely consistent with our results. For both EP Dra and HU Aqr most of the emission is located in the region where we expect material to be confined to the magnetic field lines of the white dwarf. In particular both models show an enhanced concentration of emission points towards the white dwarf. This may be indicative of where stream material is heated by X-ray emission from the accretion region on the white dwarf. However, the contribution from the accretion region, seen as the large amplitude change at ingress and egress, is fixed for the duration of the model. Therefore, if the input value is too low the model will assign flies to the white dwarf to boost the contribution at the ingress and egress of the accretion region. The number of flies near to the white dwarf may be indicative of this effect, rather than a real stream brightness enhancement.

An important advantage of placing no constraint on the stream trajectory, is the lack of sensitivity to input parameters. The HU Aqr data presented here has been modeled previously with a version of the technique of Harrop-Allin et al. (1999a). This method was found to be acutely sensitive to the exact value chosen for the distance to the threading region, and also the location of the accretion region on the white dwarf surface, which fixes the geometry of the magnetic field lines (Bridge et al. 2002). This meant that the application to the data was restricted, and results had to be interpreted accordingly. As this new technique does not rely on the choice of initial parameters for the stream trajectory, we reduce these uncertainties.

The stream trajectory used previously by Bridge et al. (2002) to model the HU Aqr data is consistent with that found from this new technique. The results of this modeling identified brightness enhancements towards the white dwarf and in the threading region. However with this new technique we do not see much enhancement where we expect material to be threaded by the field lines. The technique of Harrop-Allin et al. (1999a) was also applied to other observations of HU Aqr (Harrop-Allin et al. 1999b). They found that there was generally enhanced brightness towards the white dwarf and in the threading region, and for some cycles the threading region was significantly brighter. By comparison to both these methods, Vrielmann & Schwope (2000) did not find any enhancement towards the white dwarf, but did identify a brighter threading region. However, caution should be exercised as they used line emission only, while here we use the total emission (line and continuum). The results from modeling different light curves of HU Aqr may indicate some evolution of the stream brightness with time. However, the results of applying different models to the same data indicate that the results may be dependent upon the model parameters, as found for the technique of Harrop-Allin et al. (Bridge et al. 2002).

In HU Aqr we see no evidence for emission from where we would expect material to be following a ballistic trajectory. In EP Dra however, there is some indication that there is emission from along a ballistic trajectory (in Figures 2 and 3 the dashed curve represents the ballistic trajectory, seen most clearly in the left-hand panels). The location of material close to the orbital plane in EP Dra causes accretion at low latitude and is consistent with the suggested location of the accretion region from the cyclotron models of Schwope & Mengel (1997). Low latitude accretion may be indicative of a more complicated magnetic field geometry, with a quadrupole more likely to cause accretion at lower latitudes (Wu & Wickramasinghe 1993).

## 5. Summary and further work

While the previous techniques all assumed a stream trajectory, here we use a less constrained model that makes few a priori assumptions about the location of material. This shows that the trajectory is largely as expected, with the model providing good fits to the light curves of HU Aqr and EP Dra. The different results obtained from modeling the same data imply some dependence upon the technique, although this new technique is less sensitive to input parameters than that of Harrop-Allin et al. (1999a). A more detailed investigation of the effects of the input parameters and the regularisation curve is needed to fully understand any effects. We see most of the emission from material close to the white dwarf in both systems, and it is unclear why we see no emission from where we expect material to be following a ballistic trajectory in HU Aqr.

The colour information available from S-Cam 2 means that we can use the model to determine the colour dependence of the stream emission. In particular this will help explore the temperature variations along the stream by emphasising which parts of the stream are brightest at bluer wavelengths, and therefore hotter, and where the stream is redder and hence cooler. We can also apply the model to the consecutive cycles of HU Aqr and EP Dra, to examine any evolution of the stream structure and temperature.

## References

Bridge C.M. et al. 2002, MNRAS, 336, 1129

Bridge C.M. et al. 2003, MNRAS, accepted

Charbonneau, P. 1995, ApJS, 101, 309

Hakala P. 1995, A&A, 296, 164

Hakala P., Cropper M., Ramsay G. 2002, MNRAS, 334, 990

Harrop-Allin M.K. 1999, PhD thesis, Univ. London

Harrop-Allin M.K. et al. 1999a, MNRAS, 302, 362

Harrop-Allin M.K. et al. 1999b, MNRAS, 308, 807

Harrop-Allin M.K., Potter S.B., Cropper M. 2001, MNRAS, 326, 788

Heerlein C., Horne K. & Schwope, A.D. 1999, MNRAS, 304, 145

Kohonen T. 1990, Proc. of the IEEE, 78, 1464

Kube J., Gänsicke B.T. & Beuermann K. 2000, A&A, 356, 490

Lubow, S.H. & Shu, F.H. 1975, ApJ, 198, 383

Perryman M.A.C. et al. 2001, MNRAS, 324, 899

Schwope A.D & Mengel S. 1997, AN, 318, 25

Vrielmann S. & Schwope A.D. 2001, MNRAS, 322, 269

Wu, K & Wickramasinghe, D.T. 1993, MNRAS, 260, 141


text:
# The Traveling Salesman Problem for Cubic Graphs

## 1 Introduction

The traveling salesman problem and the closely related Hamiltonian cycle problem are two of the most fundamental of NP-complete graph problems . However, despite much progress on exponential-time solutions to other graph problems such as chromatic number  or maximal independent sets , the only worst-case bound known for finding Hamiltonian cycles or traveling salesman tours is that for a simple dynamic program, using time and space $`O(2^nn^{O(1)})`$, that finds Hamiltonian paths with specified endpoints for each induced subgraph of the input graph (D. S. Johnson, personal communication). Hamiltonian cycle and TSP heuristics without worst case analysis have also been studied extensively . Therefore, it is of interest to find special cases of the problem that, while still NP-complete, may be solved more quickly in the worst case than the general problem.

In this paper, we consider one such case: the traveling salesman problem in graphs with maximum degree three, which arises e.g. in computer graphics in the problem of stripification of triangulated surface models . Bounded-degree maximum independent sets had previously been considered  but we are unaware of similar work for the traveling salesman problem in bounded degree graphs. More generally, we consider the forced traveling salesman problem in which the input is a multigraph $`G`$ and set of forced edges $`F`$; the output is a minimum cost Hamiltonian cycle of $`G`$, containing all edges of $`F`$. A naive branching search that repeatedly adds one edge to a growing path, choosing at each step one of two edges at the path endpoint, and backtracking when the chosen edge leads to a previous vertex, solves this problem in time $`O(2^n)`$ and linear space; this is already an improvement over the general graph dynamic programming algorithm. One could also use algorithms for listing all maximal independent sets in the line graph of the input, to solve the problem in time $`O(3^{n/2})1.732^n`$. We show that more sophisticated backtracking can solve the forced traveling salesman problem (and therefore also the traveling salesman and Hamiltonian cycle problems) for cubic graphs in time $`O(2^{n/3})1.260^n`$ and linear space. We also provide a randomized reduction from degree four graphs to degree three graphs solving the traveling salesman problem in better time than the general case for those graphs. We then consider the problem of listing all Hamiltonian cycles. We show that all such cycles can be found in time $`O(2^{3n/8})1.297^n`$ and linear space. We can also count all Hamiltonian cycles in the same time bound. Our proof implies that every degree three graph has $`O(2^{3n/8})`$ Hamiltonian cycles; we do not know whether this bound is tight, but we exhibit an infinite family of graphs with $`2^{n/3}`$ Hamiltonian cycles per graph.

## 2 The Algorithm and its Correctness

Our algorithm is based on a simple case-based backtracking technique. Recall that $`G`$ is a graph with maximum degree 3, while $`F`$ is a set of edges that must be used in our traveling salesman tour. For simplicity, we describe a version of the algorithm that returns only the cost of the optimal tour, or the special value None if there is no solution. The tour itself can be reconstructed by keeping track of which branch of the backtracking process led to the returned cost; we omit the details The steps of the algorithm are listed in Table 1. Roughly, our algorithm proceeds in the following stages. Step 1 of the algorithm reduces the size of the input without branching, after which the graph can be assumed to be cubic and triangle-free, with forced edges forming a matching. Step 2 tests for a case in which all unforced edges form disjoint 4-cycles; we can then solve the problem immediately via a minimum spanning tree algorithm. Finally (steps 36), we choose an edge to branch on, and divide the solution space into two subspaces, one in which the edge is forced to be in the solution and one in which it is excluded. These two subproblems are solved recursively, and it is our goal to minimize the number of times this recursive branching occurs.

All steps of the algorithm either return or reduce the input graph to one or more smaller graphs that also have maximum degree three, so the algorithm must eventually terminate. To show correctness, each step must preserve the existence and weight of the optimal traveling salesman tour. This is easy to verify for most cases of steps 1 and 36. Case 1i performs a so-called $`\mathrm{\Delta }`$-$`Y`$ transformation on the graph; case analysis (Figure 1, left) shows that each edge of the contracted triangle participates in a Hamiltonian cycle exactly when the opposite non-triangle edge also participates.

###### Lemma 1

Let $`a`$, $`b`$, and $`c`$ be three edges forming a triangle in a graph $`G`$ in which every vertex has degree at most three, and let $`d`$ be a non-triangle edge incident to $`b`$ and $`c`$ (and opposite $`a`$). Then every Hamiltonian cycle of $`G`$ that contains $`a`$ also contains $`d`$, and vice versa.

###### Proof

Let $`T`$ be the set of endpoints of the triangle. Then a Hamiltonian cycle must enter and exit $`T`$ at least once. Each entry-exit pair has two edges with one endpoint in the cycle, but there are at most three such edges overall in $`G`$, so the cycle can only visit $`T`$ once: it must enter $`T`$ at one of its vertices, use triangle edges to cover a second vertex, and exit $`T`$ from the third vertex. The three ways the cycle might do this are shown in , and it is easy to see from the figure that the property specified in the lemma holds for each case.

This lemma justifies the correctness of adding non-triangle edges to $`F`$ when they are opposite forced triangle edges, and of adding the triangle edges’ weights to the weights of the opposite non-triangle edges. To justify the correctness of contracting the triangle, we use the same case analysis to prove another lemma:

###### Lemma 2

Let $`a`$, $`b`$, and $`c`$ be three edges forming a triangle in a graph $`G`$ in which every vertex has degree at most three, and let graph $`G^{}`$ be formed by removing those three edges from $`G`$ and replacing their endpoints by a single supervertex. Then the Hamiltonian cycles of $`G^{}`$ are in one-to-one correspondence with the Hamiltonian cycles of $`G`$. If $`C^{}`$ is the set of edges in a cycle in $`G^{}`$ corresponding to a cycle $`C`$ in $`G`$, then $`C^{}=C\{a,b,c\}`$.

We now turn to step 1j. This step concerns a 4-cycle in $`G`$, with edges in $`F`$ forcing the Hamiltonian cycle to enter or exit on two opposite vertices (Figure 1, right). If a Hamiltonian cycle enters and exits a cycle in $`G`$ only once, it does so on two adjacent vertices of the cycle, so the 4-cycle of this case is entered and exited twice by every Hamiltonian cycle, and the step’s addition of edges to $`F`$ does not change the set of solutions of the problem.

It remains to prove correctness of step 2 of the algorithm.

###### Lemma 3

Suppose that $`G`$, $`F`$ can not be reduced by step 1 of the algorithm described in Table 1, and that $`GF`$ forms a collection of disjoint 4-cycles. Then step 2 of the algorithm correctly solves the forced traveling salesman problem in polynomial time for $`G`$ and $`F`$.

###### Proof

Let $`C_i`$, $`H_i`$, $`H`$, and $`G^{}`$ be as defined in step 2 of the algorithm. Figure 2(left) depicts $`F`$ as the thick edges, $`C_i`$ as the thin edges, and $`H_i`$ and $`H`$ as the thin solid edges; Figure 2(middle) depicts the corresponding graph $`G^{}`$.

We first show that the weight of the optimal tour $`T`$ is at least as large as what the algorithm computes. The symmetric difference $`T(FH)`$ contains edges only from the 4-cycles $`C_i`$. Analysis similar to that for step 1j shows that, within each 4-cycle $`C_i`$, $`T`$ must contain either the two edges in $`H_i`$ or the two edges in $`C_iH_i`$. Therefore, $`T(FH)`$ forms a collection of 4-cycles which is a subset of the 4-cycles in $`GF`$ and which corresponds to some subgraph $`S`$ of $`G^{}`$. Further, due to the way we defined the edge weights in $`G^{}`$, the difference between the weights of $`T`$ and of $`FH`$ is equal to the weight of $`S`$. $`S`$ must be a connected spanning subgraph of $`G^{}`$, for otherwise the vertices in some two components of $`FH`$ would not be connected to each other in $`T`$. Since all edge weights in $`G^{}`$ are non-negative, the weight of spanning subgraph $`S`$ is at least equal to that of the minimum spanning tree of $`G^{}`$.

In the other direction, one can show by induction that, if $`T^{}`$ is any spanning tree of $`G^{}`$, such as the one shown by the thick edges in Figure 2(middle), and $`S^{}`$ is the set of 4-cycles in $`G`$ corresponding to the edges of $`T^{}`$, then $`S^{}(FH)`$ is a Hamiltonian cycle of $`G`$ with weight equal to that of $`FH`$ plus the weight of $`T^{}`$ (Figure 2(right)). Therefore, the weight of the optimal tour $`T`$ is at most equal to that of $`FH`$ plus the weight of the minimum spanning tree of $`G^{}`$.

We have bounded the weight of the traveling salesman tour both above and below by the quantity computed by the algorithm, so the algorithm correctly solves the traveling salesman problem for this class of graphs.

We summarize our results below.

###### Theorem 2.1

The algorithm described in Table 1 always terminates, and returns the weight of the optimal traveling salesman tour of the input graph $`G`$.

## 3 Implementation Details

Define a step of the algorithm of Table 1 to be a single execution of one of the numbered or lettered items in the algorithm description. As described, each step involves searching for some kind of configuration in the graph, and could therefore take as much as linear time. Although a linear factor is insignificant compared to the exponential time bound of our overall algorithm, it is nevertheless important (and will simplify our bounds) to reduce such factors to the extent possible. As we now show, we can maintain some simple data structures that let us avoid repeatedly searching for configurations in the graph.

###### Lemma 4

The algorithm of Table 1 can be implemented in such a way that step 3, and each substep of step 1, take constant time per step.

###### Proof

The key observation is that most of these steps require finding a connected pattern of $`O(1)`$ edges in the graph. Since the graph has bounded degree, there can be at most $`O(n)`$ matches to any such pattern. We can maintain the set of matches by removing a match from a set whenever one of the graph transformations changes one of its edges, and after each transformation searching within a constant radius of the changed portion of the graph for new matches to add to the set. In this way, finding a matching pattern is a constant time operation (simply pick the first one from the set of known matches), and updating the set of matches is also constant time per operation.

The only two steps for which this technique does not work are step 1c and step 1d, which each involve finding a cycle of possibly unbounded size in $`G`$. However, if a long cycle of forced edges exists, step 1e or step 1f must be applicable to the graph; repeated application of these steps will eventually either discover that the graph is non-Hamiltonian or reduce the cycle to a single self-loop. So we can safely replace step 1c and step 1d by steps that search for a one-vertex cycle in $`F`$, detect the applicability of the modified steps by a finite pattern matching procedure, and use the same technique for maintaining sets of matches described above to solve this pattern matching problem in constant time per step.

To aid in our analysis, we restrict our algorithm so that, when it can choose among several applicable steps, it gives first priority to steps which immediately return (that is, step 1a and steps 1c1e, with the modifications to step 1c and step 1d described in the lemma above), and second priority to step 1f. The prioritization among the remaining steps is unimportant to our analysis.

## 4 Analysis

By the results of the previous section, in order to compute an overall time bound for the algorithm outlined in Table 1, we need only estimate the number of steps it performs. Neglecting recursive calls that immediately return, we must count the number of iterations of steps 1b, 1f1j, and 36.

###### Lemma 5

If we prioritize the steps of the algorithm as described in the previous section, the number of iterations of step 1f is at most $`O(n)`$ plus a number proportional to the number of iterations of the other steps of the algorithm.

###### Proof

The algorithm may perform at most $`O(n)`$ iterations of step 1f prior to executing any other step. After that point, each additional forced edge can cause at most two iterations of step 1f, merging that edge with previously existing forced edges on either side of it, and each step other than step 1f creates at most a constant number of new forced edges.

The key idea of the analysis for the remaining steps is to bound the number of iterations by a recurrence involving a nonstandard measure of the size of a graph $`G`$: let $`s(G,F)=|V(G)||F||C|`$, where $`C`$ denotes the set of 4-cycles of $`G`$ that form connected components of $`GF`$. Clearly, $`sn`$, so a bound on the time complexity of our algorithm in terms of $`s`$ will lead to a similar bound in terms of $`n`$. Equivalently, we can view our analysis as involving a three-parameter recurrence in $`n`$, $`|F|`$, and $`|C|`$; in recent work  we showed that the asymptotic behavior of this type of multivariate recurrence can be analyzed by using weighted combinations of variables to reduce it to a univariate recurrence, similarly to our definition here of $`s`$ as a combination of $`n`$, $`|F|`$, and $`|C|`$. Note that step 1f leaves $`s`$ unchanged and the other steps do not increase it.

###### Lemma 6

Let a graph $`G`$ and nonempty forced edge set $`F`$ be given in which neither an immediate return nor step 1f can be performed, and let $`s(G,F)`$ be as defined above. Then the algorithm of Table 1, within a constant number of steps, reduces the problem to one of the following situations:

* a single subproblem $`G^{},F^{}`$, with $`s(G^{},F^{})s(G,F)1`$, or
* subproblems $`G_1,F_1`$ and $`G_2,F_2`$, with $`s(G_1,F_1),s(G_2,F_2)s(G,F)3`$, or
* subproblems $`G_1,F_1`$ and $`G_2,F_2`$, with $`s(G_1,F_1)s(G,F)2`$ and $`s(G_2,F_2)s(G,F)5`$.

###### Proof

If step 1b, step 1g, step 1h, or step 1j applies, the problem is immediately reduced to a single subproblem with more forced edges, and if step 1i applies, the number of vertices is reduced. Step 2 provides an immediate return from the algorithm. So, we can restrict our attention to problems in which the algorithm is immediately forced to apply steps 36. In such problems, the input must be a simple cubic triangle-free graph, and $`F`$ must form a matching in this graph, for otherwise one of the earlier steps would apply.

We now analyze cases according to the neighborhood of the edge $`yz`$ chosen in step 3. To help explain the cases, we let $`yw`$ denote the third edge of $`G`$ incident to the same vertex as $`xy`$ and $`yz`$. We also assume that no immediate return is performed within $`O(1)`$ steps of the initial problem, for otherwise we would again have reduced the problem to a single smaller subproblem.

* In the first case, corresponding to step 3a of the algorithm, $`yz`$ is adjacent to a 4-cycle in $`GF`$ which already is adjacent to two other edges of $`F`$. Adding $`yz`$ to $`F`$ in the recursive call in step 4 leads to a situation in which step 1j applies, adding the fourth adjacent edge of the cycle to $`F`$ and forming a 4-cycle component of $`GF`$. Thus $`|F|`$ increases by two and $`|C|`$ increases by one. In step 5, $`yz`$ is removed from $`F`$, following which step 1b adds two edges of the 4-cycle to $`F`$, step 1f contracts these two edges to a single edge, shrinking the 4-cycle to a triangle, and step 1i contracts the triangle to a single vertex, so the number of vertices in the graph is decreased by three.
* In the next case, $`yz`$ is chosen by step 3b to be adjacent to forced edge $`xy`$, and neither $`yz`$ nor $`yw`$ is incident to a second edge in $`F`$. If we add $`yz`$ to $`F`$, an application of step 1f removes $`yw`$, and another application of step 1b adds the two edges adjoining $`yw`$ to $`F`$, so the number of forced edges is increased by three. The subproblem in which we remove $`yz`$ from $`F`$ is symmetric. This case and its two subproblems are shown in Figure 3.
* If step 3b chooses edge $`yz`$, and $`z`$ or $`w`$ is incident to a forced edge, then with $`y`$ it forms part of a chain of two or more vertices, each incident to exactly two unforced edges that connect vertices in the chain. This chain may terminate at vertices with three adjacent unforced edges (Figure 4, left). If it does, a similar analysis to the previous case shows that adding $`yz`$ to $`F`$ or removing it from $`G`$ causes alternating members of the chain to be added to $`F`$ or removed from $`G`$, so that no chain edge is left unforced. In addition, when an edge at the end of the chain is removed from $`G`$, two adjacent unforced edges are added to $`F`$, so these chains generally lead to a greater reduction in size than the previous case. The smallest reduction happens when the chain consists of exactly two vertices adjacent to forced edges. In this case, one of the two subproblems is formed by adding two new forced edges at the ends of the chain, and removing one edge interior to the chain; it has $`s(G_1,F_1)=s(G,F)2`$. The other subproblem is formed by removing the two edges at the ends of the chain, and adding to $`F`$ the edge in the middle of the chain and the other unforced edges adjacent to the ends of the chain. None of these other edges can coincide with each other without creating a 4-cycle that would have been treated in the first case of our analysis, so in this case there are five new forced edges and $`s(G_2,F_2)=s(G,F)5`$.
* In the remaining case, step 3b chooses an edge belonging to a cycle of unforced edges, each vertex of which is also incident to a forced edge (Figure 4, right). In this case, adding or removing one of the cycle edges causes a chain reaction which alternately adds and removes all cycle edges. This case only arises when the cycle length is five or more, and if it is exactly five then an inconsistency quickly arises causing both recursive calls to return within a constant number of steps. When the cycle length is six or more, both resulting subproblems end up with at least three more forced edges.

Note that the analysis need not consider choices made by step 3c of the algorithm, as $`F`$ is assumed nonempty; step 3c can occur only once and does not contribute to the asymptotic complexity of the algorithm. In all cases, the graph is reduced to subproblems that have sizes bounded as stated in the lemma.

###### Theorem 4.1

The algorithm of Table 1 solves the forced traveling salesman problem on graphs of degree three in time $`O(2^{n/3})`$.

###### Proof

The algorithm’s correctness has already been discussed. By Lemmas 3, 4, 5, and 6, the time for the algorithm can be bounded within a constant factor by the solution to the recurrence

$$T(s)1+\mathrm{max}\{s^{O(1)},T(s1),2T(s3),T(s2)+T(s5)\}.$$

Standard techniques for linear recurrences give the solution as $`T(s)=O(2^{s/3})`$. In any $`n`$-vertex cubic graph, $`s`$ is at most $`n`$, so expressed in terms of $`n`$ this gives a bound of $`O(2^{n/3})`$ on the running time of the algorithm.

## 5 Degree Four

It is natural to ask to what extent our algorithm can be generalized to higher vertex degrees. We provide a first step in this direction, by describing a randomized (Monte Carlo) algorithm: that is, an algorithm that may produce incorrect results with bounded probability. To describe the algorithm, let $`f`$ denote the number of degree four vertices in the given graph. The algorithm consists of $`(3/2)^f`$ repetitions of the following: for each degree four vertex, choose randomly among the three possible partitions of its incoming edges into two sets of two edges; split the vertex into two vertices, with the edges assigned to one or the other vertex according to the partition, and connect the two vertices by a new forced edge (Figure 5). Once all vertices are split, the graph has maximum degree 3 and we can apply our previous forced TSP algorithm.

It is not hard to see that each such split preserves the traveling salesman tour only when the two tour edges do not belong to the same set of the partition, which happens with probability $`2/3`$; therefore, each repetition of the algorithm has probability $`(2/3)^f`$ of finding the correct TSP solution. Since there are $`(3/2)^f`$ repetitions, there is a bounded probability that the overall algorithm finds the correct solution. Each split leaves unchanged the parameter $`s`$ used in our analysis of the algorithm for cubic graphs, so the time for the algorithm is $`O((3/2)^f2^{n/3})=O((27/4)^{n/3})`$. By increasing the number of repetitions the failure probability can be made exponentially small with only a polynomial increase in runtime.

A hitting set technique due to Beigel (personal communication, 1995) allows this algorithm to be derandomized. We group the degree four vertices arbitrarily into groups of $`k`$ vertices per group. Within each group, a single choice among the three possible expansions of each degree-four vertex can be described by a $`k`$-digit ternary word $`w\{0,1,2\}^k`$. If $`x=x_0x_1\mathrm{}x_{k1}`$ is a ternary word, let $`D(x)=\{w_0w_1\mathrm{}w_{k1}\{0,1,2\}^k\text{for all }i\text{,}w_ix_i\}`$. As two of the three expansions of each vertex preserve the correct traveling salesman solution, the set of expansions of the group that preserve the TSP solution is a set of the form $`D(x)`$ for some unknown $`x`$. We can find one such expansion if we test all words in a hitting set: that is, a set $`S`$ of ternary words such that $`SD(x)\mathrm{}`$ for all $`x`$. In order to find a fast deterministic algorithm, we seek a hitting set that is as small as possible.

As each set $`D(x)`$ has the same cardinality $`2^k`$, the hitting set problem for the sets $`D(x)`$ has a simple fractional solution: assign weight $`2^k`$ to each of the $`3^k`$ possible ternary words. The total weight of this solution is therefore $`(3/2)^k`$. By standard results relating integer to fractional solutions of the hitting set (or equivalently set cover) problem, this implies that there is a hitting set of cardinality at most $`(3/2)^k\mathrm{ln}𝒟=O(k(3/2)^k)`$, where $`𝒟=3^k`$ is the cardinality of the family of sets $`D(x)`$.

We use this hitting set $`H`$ as part of a deterministic search algorithm that tests each choice of a member of $`H`$ for each group of $`k`$ vertices. For each of these $`|H|^{n/k}`$ choices, we expand the degree four vertices in each group as specified by the choice for that group, and then apply our degree-three TSP algorithm on the expanded graph. At least one choice hits the set $`D(x)`$ in each group of expansions preserving the TSP, so the best TSP solution among the $`|H|^{n/k}`$ expansions must equal the TSP of the original graph. By choosing a suitably large constant $`k`$, we can achieve time $`O((27/4+ϵ)^{n/3})`$ for any constant $`ϵ>0`$.

We omit further details as this result seems unlikely to be optimal.

## 6 Listing All Hamiltonian Cycles

Suppose we want not just a single best Hamiltonian cycle (the Traveling Salesman Problem) but rather a list of all such cycles. As we show in Table 2, most of the steps of our traveling salesman algorithm can be generalized in a straightforward way to this cycle listing problem. However, we do not know of a cycle listing analogue to the minimum spanning tree algorithm described in step 2 of Table 1, and proven correct in Lemma 3 for graphs in which the unforced edges form disjoint 4-cycles. It is tempting to try listing all spanning trees instead of computing minimum spanning trees, however not every Hamiltonian cycle of the input graph $`G`$ corresponds to a spanning tree of the derived graph $`G^{}`$ used in that step. Omitting the steps related to these 4-cycles gives the simplified algorithm shown in Table 2. We analyze this algorithm in a similar way to the previous one; however in this case we use as the parameter of our analysis the number of unforced edges $`U(G)`$ in the graph $`G`$. Like $`s(G)`$, $`U`$ does not increase at any step of the algorithm; we now show that it decreases by sufficiently large amounts at certain key steps.

###### Lemma 7

If $`G`$ initially has no parallel edges or self-loops, none are introduced by the steps of this algorithm.

###### Proof

The only step that adds a new edge to the graph is step 1e, which replaces a two-edge path with a single edge. If the graph has no parallel edges before this step, the step cannot create a self-loop. It can create a multiple adjacency if the two path edges belong to a triangle, but as part of the step we immediately detect and eliminate this adjacency.

###### Lemma 8

Let a graph $`G`$ be given in which neither a backtrack nor step 1e can be performed, let $`F`$ be nonempty, and let $`U(G)`$ denote the number of unforced edges in $`G`$. Then the algorithm of Table 2, within a constant number of steps, reduces the problem to one of the following situations:

* a single subproblem $`G^{}`$, with $`U(G^{})U(G)1`$, or
* two subproblems $`G_1`$ and $`G_2`$, with $`U(G_1),U(G_2)U(G)4`$, or
* two subproblems $`G_1`$ and $`G_2`$, with $`U(G_1)U(G)3`$ and $`U(G_2)U(G)6`$.

###### Proof

If step 1c or step 1d applies, the problem is immediately reduced to a single subproblem with fewer unforced edges. So, we can restrict our attention to problems in which the algorithm is immediately forced to apply step 2. In such problems, the input must be a simple cubic triangle-free graph, and $`F`$ must form a matching in this graph, for otherwise one of the earlier steps would apply.

We now perform a case analysis according to the local neighborhood of the edge $`yz`$ chosen in step 1d. As in Lemma 6, let $`yw`$ denote the third edge of $`G`$ incident to the same vertex as $`xy`$ and $`yz`$.

* In the first case, neither $`yz`$ nor $`yw`$ is incident to a second edge in $`F`$. In the subproblem in which we add $`yz`$ to $`F`$, an application of step 1e removes $`yw`$, and another application of step 1c adds the two edges adjoining $`yw`$ to $`F`$. Thus, in this case, the number of unforced edges is reduced by four. The subproblem in which we remove $`yz`$ from $`F`$ is symmetric. This case and the two subproblems it produces are shown in Figure 3.
* If $`z`$ or $`w`$ is incident to a forced edge, then with $`y`$ it forms part of a chain of two or more vertices, each of which is incident to exactly two unforced edges that connect vertices in the chain. This chain may terminate at vertices with three adjacent unforced edges (Figure 4, left). If it does, a similar analysis to the previous case shows that adding $`yz`$ to $`F`$ or removing it from $`G`$ causes alternating members of the chain to be added to $`F`$ or removed from $`G`$ as well, so that no chain edge is left unforced. In addition, when an edge at the end of the chain is removed from $`G`$, the two adjacent unforced edges are added to $`F`$. The smallest reduction in unforced edges happens when the chain consists of exactly two vertices adjacent to forced edges. In this case, one of the two subproblems is formed by adding two new forced edges at the ends of the chain, and removing one edge interior to the chain; it thus has $`U(G_1)=U(G)3`$. The other subproblem is formed by removing the two edges at the ends of the chain, and adding to $`F`$ the edge in the middle of the chain and the other unforced edges adjacent to the ends of the chain. Thus it would seem that this subproblem has $`U(G_2)=U(G)7`$, however it is possible for one unforced edge to be adjacent to both ends of the chain, in which case we only get $`U(G_2)=U(G)6`$.
* In the remaining case, we have a cycle of four or more unforced edges, each vertex of which is also incident to a forced edge (Figure 4, right). In this case, adding or removing one of the cycle edges causes a chain reaction which alternately adds and removes all cycle edges, so both resulting subproblems end up with at least four fewer unforced edges.

Thus, in all cases, the graph is reduced to two subproblems that have numbers of unforced edges bounded as in the statement of the lemma.

###### Theorem 6.1

For any simple graph $`G`$ with maximum degree 3, set $`F`$ of forced edges in $`G`$, and assignment of weights to the edges of $`G`$ from a commutative semiring, we can list all Hamiltonian cycles in $`G`$ in $`O(2^{3n/8})`$ and linear space.

###### Proof

By the previous lemma, the number of calls in the algorithm can be bounded within a constant factor by the solution to the recurrence

$$T(u)1+\mathrm{max}\{T(u1),2T(u4),T(u3)+T(u6)\}.$$

Standard techniques for linear recurrences give the solution as $`T(u)=O(2^{u/4})`$. In any $`n`$-vertex cubic graph, $`u`$ is at most $`3n/2`$, so expressed in terms of $`n`$ this gives a bound of $`O(2^{3n/8})`$ on the number of operations. As in the previous algorithm, by the appropriate use of simple data structures we can implement each step of the algorithm in constant time per step.

In order to achieve this time bound, we must output each Hamiltonian cycle using an implicit representation that changes by a constant amount in each step. Explicitly listing the vertices in order in each Hamiltonian cycle would take $`O(2^{3n/8}+Kn)`$ where $`K`$ is the number of Hamiltonian cycles produced by the algorithm. Appendix A presents an implementation of this algorithm in the Python language, using such an implicit representation. This implementation was able to find a Hamiltonian cycle of a 400-vertex 3-regular graph, dual to a triangulated torus model (Figure 6) in approximately two seconds on an 800 MHz PowerPC G4 computer. However, on a slightly larger 480-vertex graph dual to a triangulated sphere model, we aborted the computation after 11 hours with no result.

###### Corollary 1

We can count the number of Hamiltonian cycles in any cubic graph in time $`O(2^{3n/8})`$ and linear space.

###### Proof

We simply maintain a counter of the number of cycles seen so far, and increment it each time we output another cycle. The average number of bits changed in the counter per step is $`O(1)`$, so the total additional time to maintain the counter is $`O(2^{3n/8})`$.

A preliminary version of this paper used a similar recursion for counting Hamiltonian cycles, but returned the counts from each recursive subproblem, incurring an additional polynomial factor overhead for the arithmetic involved.

## 7 Graphs with Many Hamiltonian Cycles

The following result follows immediately from Theorem 6.1:

###### Corollary 2

Every simple $`n`$-vertex graph with maximum degree three has at most $`2^{3n/8}`$ Hamiltonian cycles.

For cubic multigraphs, one can show a $`2^{n/2}`$ bound on the number of Hamiltonian cycles by choosing one distinguished cycle, and corresponding the other cycles to subsets of the remaining $`n/2`$ edges. This bound is tight as can be seen by the example of an $`n`$-gon with alternating single and double bonds. We do not know whether our $`2^{3n/8}`$ bound is tight, but we exhibit in Figure 7 the construction for an infinite family of graphs with $`2^{n/3}`$ Hamiltonian cycles per graph. Each graph in the family is formed by connecting $`n/6`$ $`6`$-vertex gadgets into a cycle, where each gadget is formed by removing an edge from the graph $`K_{3,3}`$. A Hamiltonian cycle of a graph formed in this way must pass once through each gadget, and there are four possible paths through each gadget, so the total number of Hamiltonian cycles is $`2^{n/3}`$. We note that the $`n/6`$ edges connecting pairs of gadgets in this graph are all forced to occur in any Hamiltonian cycle, so in terms of the number $`u`$ of unforced edges the number of Hamiltonian cycles is $`4^{u/8}=2^{u/4}`$, matching the worst-case bound in terms of $`u`$ for our cycle-listing algorithm.

We used our cycle listing algorithm to search, unsuccessfully, for a better gadget among all possible cubic graphs on 20 or fewer vertices, using tables of these graphs provided online by Gordon Royle . Our experiments support the conjecture that, in any $`n`$-vertex cubic graph, each edge participates in at most $`2^{n/3}`$ Hamiltonian cycles. This bound is achieved by the graphs of Figure 7 as well as by similar graphs that include one or two four-vertex gadgets formed by removing an edge from $`K_4`$, so if true this bound would be tight. The truth of this conjecture would imply that any cubic graph has $`O(2^{n/3})`$ Hamiltonian cycles.

## Appendix A Implementation of Cycle Listing Algorithm

We present here an implementation of the $`O(2^{3n/8})`$ algorithm for listing all Hamiltonian cycles in a degree three graph, in the Python programming language . The yield keyword triggers Python’s simple generator protocol , which creates an iterator object suitable for use in for-loops and similar contexts and returns it from each call to HamiltonianCycles. A more elaborate version of the implementation, which backtracks when it discovers that the current graph is not biconnected or its unforced edges have no perfect matching, and includes code for testing the algorithm on several simple families of 3-regular graphs, is available for download at http://www.ics.uci.edu/~eppstein/PADS/CubicHam.py.


text:
# The Early Light Curve of the Optical Afterglow of GRB 021211

## 1 INTRODUCTION

The majority of gamma-ray burst (GRB) optical afterglow (OA) observations have been obtained several hours after the bursts, owing mainly to the delay in getting X-ray images and/or the difficulties in localizing the GRBs in real time. Out of the $`40`$ detected GRB OAs, only GRB 990123 (Akerlof et al. 1999), GRB 021004 (Fox 2002), and GRB 021211 (Fox & Price 2002) were successfully detected within 10 min after the initial burst of gamma rays.

In this Letter, we report on early-time optical observations of GRB 021211 with the Katzman Automatic Imaging Telescope (KAIT). KAIT is a 0.76-m robotic telescope at Lick Observatory that is dedicated to searching for and observing supernovae and monitoring other celestial events. It is equipped with an Apogee AP7 back-illuminated CCD camera having a resolution of $`0\stackrel{}{\mathrm{.}}8`$ pixel-1 and a total field of view of $``$ $`6\stackrel{}{\mathrm{.}}8\times 6\stackrel{}{\mathrm{.}}8`$. More information on KAIT can be found in Li et al. (2000) and Filippenko et al. (2001), and the KAIT GRB alert system is described in detail elsewhere (Li et al. 2003).

## 2 OBSERVATIONS

On 2002 December 11, 11:18:34.03 UT (defined as $`t=0`$ in this paper), a bright GRB triggered the HETE-2 satellite (Crew et al. 2002). The gamma-ray light curve shows a single peak lasting $``$ 2.5 s. The flight localization of the GRB was reported in a GRB Coordinates Network (GCN; Barthelmy et al. 1994) notice at $`t=22`$ s.

The KAIT GRB alert program received the GCN notice with the position of GRB 021211 at $`t=32`$ s. At the time, KAIT was conducting followup observations of one of its supernova discoveries, SN 2002he, and was 478 s into a 600-s $`B`$-band exposure (a rare, long exposure; most are only 20 s). KAIT performed some checks to ensure that the GRB was observable, stopped the autoguider, rotated the filter wheel so that no filter was in the way of the light path, and began to slew across much of the sky to the GRB position at $`t=49`$ s. Currently, KAIT cannot stop an exposure in progress, so the shutter remained open during the long move to the opposite side of the sky, fortunately imaging only the dark, slowly rotating dome. The dome slit was centered on the GRB position at $`t=108`$ s, with 48 s left in the ongoing exposure. The effective start time for the GRB exposure is $`t=105`$ s.111We have conducted tests to repeat the same pointing procedure, and concluded that the dome cleared the telescope 3 s earlier than centering on the GRB position. This exposure, with the SN 2002he field and the GRB 021211 field superimposed, is shown in Figure 1a.

KAIT then began its pre-programmed sequence of GRB images. It started with a 10-s exposure at $`t=171`$ s (Fig. 1b), continued with a 40-s exposure, five $`3\times 20`$ s grid images (for a total of 15 individual exposures of the GRB), two 60-s exposures, one 120-s guided exposure, one 300-s guided exposure, and four 120-s guided exposures that together make a mosaic covering $`12\stackrel{}{\mathrm{.}}7\times 12\stackrel{}{\mathrm{.}}7`$ around the GRB position. One of the $`3\times 20`$ s grid images is shown in Figure 1c. Such grid images, in which the telescope is moved by $`36^{\prime \prime }`$ in right ascension between consecutive exposures, save two readout times (10 s each).

KAIT received three more GCN notices on GRB 021211 in the next 3 hr, and repeated the same sequence each time. The notices arrived at 11:53:30 UT, 12:25:24 UT, and 13:29:33 UT, which correspond to 2096 s, 4010 s, and 7859 s after the burst, respectively. KAIT was conducting the regular supernova search (20-s exposures) when these notices arrived, and was able to start the respective GRB exposure 35 s, 23 s, and 45 s after receiving the GCN notice. A total of 61 images (including the one superimposed on the field of SN 2002he) were obtained, with 101 individual exposures of GRB 021211. Unfortunately, there were passing clouds during the night, and the short exposures in the last three sequences do not go sufficiently deep to detect the GRB.

Since our observations were obtained without filters, it is important to transform the measured magnitudes to a standard photometric system, so that comparisons with observations obtained elsewhere can be made. Details of transforming the unfiltered magnitudes to the standard Cousins $`R`$ are described by Li et al. (2003); see also Riess et al. (1999). These papers demonstrate that when the unfiltered observations are treated as though they were obtained with a broadband filter, and color terms are applied to correct for the difference in the throughput curves, the observed unfiltered magnitude can be transferred to a standard passband with a precision of $``$ 5%. The color terms are determined by observing Landolt standard stars in $`BVRI`$ and without filters during photometric nights; a value of 0.27$`\pm `$0.05 mag is found for transforming KAIT unfiltered magnitudes to Cousins $`R`$. Knowledge of the ($`VR`$) color of GRB 021211 is also needed for the transformation, and we assume ($`VR`$) = 0.2 mag, similar to that of GRB 990123 (Fruchter et al. 1999). A different choice of ($`VR`$) for GRB 021211 will make our measurements systematically brighter or fainter. For example, if ($`VR`$) = 0.5 mag were used, our measurements would need to be made 0.08 mag brighter.

We used the PSF-fitting technique to reduce our data. For the grid images, the PSF for an individual exposure is constructed by using stars in that exposure only. Absolute photometric calibration of the GRB 021211 field is provided by Henden (2002). Our final photometry in the $`R`$ band is listed in Table 1. In total, 28 useful data points were measured from the observations, 18 of which were obtained within 10 min after the burst.

## 3 RESULTS

Figure 2 displays the early-time light curve of GRB 021211 from the KAIT observations. We assume no host galaxy reddening to GRB 021211, and adopt a Galactic extinction of $`A_R=0.065`$ mag (Schlegel, Finkbeiner, & Davis 1998) for the figure. The data show a steep initial decline, well fit by a power-law decay, with indexes measured as $`1.56\pm 0.02`$ for the first 4 points, and $`1.58\pm 0.02`$ for the first 12 points. An apparent break is seen at $`t0.1`$–0.2 hr, when the GRB switches to a slower power-law decay, with indexes measured as $`1.01\pm 0.05`$ from the last 7 points and $`0.94\pm 0.06`$ from the last 6 points. The power-law indexes reported here are consistent with, but supersede, those measured by Chornock et al. (2002) from the preliminary reduction of the KAIT data (Li et al. 2002).

To better describe the temporal evolution of the OA of GRB 021211, we fit the data with the sum of two power-law decays:

$$F_\nu =0.5F_{\nu ,b}[(t/t_b)^{\alpha _1}+(t/t_b)^{\alpha _2}],$$
(1)
where $`t_b`$ is the time of the break, and $`F_{\nu ,b}`$ is the flux at $`t_b`$. We obtained the following values for the parameters using a $`\chi ^2`$-minimizing procedure: $`\alpha _1`$ = 1.82$`\pm `$0.02, $`\alpha _2`$ = 0.82$`\pm `$0.11, and $`t_b`$ = 0.18$`\pm `$0.03 hr (= 10.2 min). These parameters are consistent with those found by Fox et al. (2003) with independent later-epoch data. The solid line in the upper panel of Figure 2 shows this fit, while the lower panel shows the residuals of the fit. The reduced $`\chi ^2`$ of the fit is 3.34 for 25 degrees of freedom, indicating that eq. (1) is not a good model for the data at the 8$`\sigma `$ confidence level,222We have additionally tried to fit the data with smoothly broken power-law models of varying sharpness, but these did not change the results significantly. and suggest that the short timescale variations in the residuals seen in the lower panel of Figure 2 (and especially before $`t_b`$) are real features.

If we exclude the first data point, the model fit to eq. (1) becomes acceptable (reduced $`\chi ^2`$ = 1.03). However, we have no a priori reason to suspect the first data point is anomalous. For the first point to be consistent with the simple double-power-law model, (1) the flux-weighted mean exposure time would have to be 135.0 s, which is 4.5 s after the midpoint of the exposure, (2) the photometry error bar would have to be increased to 0.08 mag, or (3) the GRB would have to have been 0.11 mag brighter than our measurement. The first option requires the GRB to be brightening during the exposure, and is ruled out as the earlier RAPTOR image (Wozniak et al. 2002) shows the GRB already declining. The second option is also not viable, because the OA is bright and well-detected. The third option is possible if the color of the GRB underwent a dramatic change at early times — e.g., the $`VR`$ color was 0.6 mag during the first exposure, and became bluer to 0.2 mag (adopted in the photometric reductions) starting from the second exposure.

## 4 DISCUSSION

The OA of GRB 021211 was first identified by Fox & Price (2002) in images taken with the Palomar 1.2-m Oschin Schmidt telescope, beginning about 20 min after the burst. It was caught at early times by two robotic telescopes besides KAIT. The RAPTOR team obtained one image of GRB 021211 at $`t`$ = 65 s (Wozniak et al. 2002), the second-earliest detection ever of an afterglow (only GRB 990123 was detected earlier at $`t`$ = 23 s; Akerlof et al. 1999), and Super-LOTIS observations began at $`t`$ = 143 s (Park, Williams, & Barthelmy 2002). Attempts to detect the afterglow of GRB 021211 in other energy bands all failed: these include the prompt GeV/TeV emission by Milagro (McEnery 2002), the sub-mm search by Hoge et al. (2002) on Dec. 12, and the radio searches by Berger & Frail (2002) on Dec. 12 and by Rol & Strom (2002) on Dec. 12 and 17. Late-time spectroscopy at the location of GRB 021211 shows underlying emission consistent with a host galaxy at probable redshift 1.006 (Vreeswijk et al. 2003).

GRB 021211 is only the second GRB with prompt optical data, so it is of interest to compare its light curve with that of GRB 990123 (Akerlof et al. 1999). In Figure 3, the data for GRB 990123 were collected from various sources as listed in Fruchter et al. (1999); only the last 5 ROTSE points are included. We have also adopted the time zero point for GRB 990123 from Fruchter et al. (1999). The data for GRB 021211 reported in the GCN Circulars are also plotted. We fit the same power-law model \[eq. (1)\] to the GRB 990123 data within 2 d after the burst and obtained the following parameters: $`\alpha _1`$ = 1.80$`\pm `$0.11, $`\alpha _2`$ = 1.06$`\pm `$0.05, and $`t_b`$ = 0.15$`\pm `$0.04 hr (= 9.0 min). The reduced $`\chi ^2`$ of the fit is 1.10 for 20 degrees of freedom, which indicates eq. (1) is an adequate model to fit the data of GRB 990123. The fit parameters of GRB 990123 are very similar to those of GRB 021211, except for perhaps $`\alpha _2`$, as the late-time decline of GRB 990123 is a bit steeper (as shown also by the residuals in the lower panel of Figure 3). Although the temporal evolution of both GRBs is consistent with a break at $`t`$ 10 min, the well-sampled KAIT data of GRB 021211 are the first to directly demonstrate the break and show a smooth transition between the two power-law decays.

Despite the overall similarity between the temporal evolution of the OA of the two bursts, GRB 990123 is brighter than GRB 021211 by about 3 mag at similar epochs, even though the former is probably at a larger redshift ($`z`$ = 1.600; Hjorth et al. 1999). Other important differences between the two bursts are that (1) GRB 990123 is a much longer burst (lasting $``$ 63.3 s; Galama et al. 1999) than GRB 021211 (lasting only $``$ 2.5 s); (2) GRB 990123 is exceptionally bright in gamma rays, and has a fluence that is in the top 0.3% of BATSE bursts, while the gamma-ray fluence of GRB 021211 is in the middle of HETE-2 bursts reported by Barraud et al. (2003); and (3) the afterglow of GRB 990123 was successfully observed in optical, infrared, sub-mm, mm, and radio passbands, while the afterglow of GRB 021211 is detected only in the optical and infrared. In particular, a radio flare was observed for GRB 990123 (Kulkarni et al. 1999; Galama et al. 1999) and interpreted by Kulkarni et al. (1999) as evidence of reverse-shock emission. The apparent faintness of the OA of GRB 021211 compared with the OA of GRB 990123 suggests that the afterglow of GRB 021211 in other energy bands may be fainter as well, contributing to the reasons why it was not detected.

The current standard model for GRBs and their afterglows is the fireball shock model (for a review see Mészáros 2002). Based on this model, we expect a prompt multi-wavelength flash, contemporaneous with the gamma-ray emission and reaching $`V`$ magnitudes ranging between 9 and 18 for typical GRB distances (Mészáros & Rees 1997; Sari & Piran 1999a; Sari & Piran 1999b). The prompt optical emission of GRB 990123 is not proportional to the gamma-ray flux, nor can it be understood as the low-energy extrapolation of the GRB spectrum (Galama et al. 1999). This indicates that the prompt optical emission and the GRB originate from different regions in the fireball. It has been popular practice to ascribe the origin of the GRB to the internal shocks (Rees & Mészáros 1994), and the long-term afterglow to the forward component of the external shocks. The prompt optical emission is explained as a consequence of the reverse component of the external shocks (Mészáros & Rees 1997; Sari & Piran 1999a). The decay rate of the optical flux from the reverse shocks is much faster than that of the forward shocks, so the emission of the latter dominates after tens of minutes. If this interpretation is correct, GRB 990123 and GRB 021211 would be the only two bursts in which all three emitting regions have been convincingly seen, and demonstrate that afterglow observations may give us insight into the geometry of the central engine of GRBs, as well as the structure of the relativistic ejecta.

The interpretation that the prompt and afterglow emissions come from different physical origins also implies that the two components are largely independent, so if both are declining as power laws, eq. (1) should be an adequate model to fit the data. The rather poor fit as discussed in §3 is thus surprising. Besides the possible explanation of a dramatic color change at early times (§3), another alternative is that there are small-scale variations in the reverse-shock emission superimposed on a power-law decline. It is also possible that the afterglow emission does not begin its power-law decay until after undergoing a brightening phase to a peak (e.g., Fox et al. 2003), in which case eq. (1) would not be an accurate model for the data. A full treatment of these possibilities is beyond the scope of the current Letter. We also emphasize that early-time photometry with high temporal resolution, preferably in multiple passbands, is necessary to determine which explanation is correct.

The small-scale variations in the reverse-shock emission, if real, are characterized by a timescale of $`\mathrm{\Delta }tt`$. This suggests a process that involves the whole observable fireball surface. It will be difficult to achieve these variations in the reverse shock with external density inhomogeneities, and a refreshed shock (an extra energy input from the inner engine) seems to be a more promising explanation (e.g., Panaitescu, Mészáros, & Rees 1998). This would indicate a difference between the variability in GRB 021211 and GRB 021004, whose variability is believed to be caused by external density inhomogeneities (e.g., Lazzati et al. 2002).

The OA of GRB 021211 declined 4 mag in half an hour, and was fainter than 20 mag just over 2 hr after the burst. These observations may shed some light on the mystery of the so-called “dark bursts,” in which no OA is seen (Fynbo et al. 2001; Berger et al. 2002). Explanations for the non-detections of OAs when deep, prompt searches exist are that they are (1) obscured by dust in their host galaxies, (2) at high redshifts, (3) intrinsically faint, and (4) fast decliners. Out of these, extinction by dust is often favored, since GRBs are suggested to be associated with massive star formation, and there are GRBs whose radio and infrared afterglows are detected but not the OA (e.g., GRB 980329; Yost et al. 2002). The relative faintness of the OA of GRB 021211 compared to the OA of GRB 990123, and its rapid optical decay, suggest that some dark bursts may be intrinsically dim, and not obscured by dust. The faintness of such OAs, coupled with the fast decline typical of OAs, will make them difficult to observe at high redshifts. Thus, the diversity in the underlying properties of OAs must be observationally determined before substantive inferences can be drawn from the statistics of dark bursts.

The KAIT observations of GRB 021211 demonstrate the importance of rapid turnaround in localizing GRB positions and intensive followup observations. Future prospects for optical observations of GRBs are bright: with continuing contributions from HETE-2, and the planned mission of Swift, fast and accurate GRB positions will be delivered automatically to the ground for rapid followup by KAIT and other robotic systems.

We thank S. Barthelmy for his help setting up the KAIT GRB alert system, and the referee for useful comments. The work of A.V.F.’s group at U. C. Berkeley is supported by NSF grant AST-9987438. KAIT was made possible by generous donations from Sun Microsystems, Inc., the Hewlett-Packard Company, AutoScope Corporation, Lick Observatory, the NSF, the University of California, and the Sylvia and Jim Katzman Foundation. S.J. thanks the Miller Institute for Basic Research in Science at U. C. Berkeley for support through a research fellowship.


text:
# Supermassive black holes in radio-loud AGN

## 1 Introduction

Over the last decade the combination of HST and ground-based data has led to the measurement of the black-hole masses of some 40 inactive galaxies in the nearby Universe via gas and stellar dynamics. These black-hole mass measurements have led to the discovery of correlations between black-hole mass and the properties of the surrounding host galaxy, namely bulge luminosity and stellar velocity dispersion. In parallel to the discovery of these correlations in inactive galaxies, the last few years has also seen a new industry emerge in estimating the black-hole masses of active galaxies. In the first section of this review I will discuss the recent advances in our understanding of the connection between black-hole mass and host-galaxy properties in both active and inactive galaxies at low redshift. In the second section I will proceed to discuss how black-hole mass estimation in active galaxies has been used to investigate the connection between black-hole mass and radio power.

## 2 The black-hole masses of nearby inactive galaxies

Within the context of the study of AGN black-hole masses, the most important aspect of the recent increase in the numbers of nearby quiescent galaxies with black-hole mass measurements has been the discovery of two correlations between black-hole mass and the properties of the surrounding host galaxy. Given the recent demonstration that at low redshift the host galaxies of all powerful AGN are indistinguishable from nearby inactive ellipticals (Dunlop et al. 2003), irrespective of radio luminosity, these correlations can now be applied to estimate the central black-hole masses of active galaxies.

The first of these correlations is between black-hole mass and the bulge luminosity of the host galaxy (Kormendy & Richstone 1992, Magorrian et al. 1998). This $`M_{bh}M_{bulge}`$ correlation has the advantage of being observationally inexpensive, particularly with respect to radio galaxies, but is traditionally thought to suffer from a large associated scatter $`0.5`$ dex (although see below).

The second, more recently discovered, correlation is between black-hole mass and the velocity dispersion of the host galaxy stellar population (Gebhardt et al. 2000, Ferrarese & Merritt 2000). This correlation has been hailed as being the more “fundamental” due to its very low associated scatter ($`0.3`$ dex). However, the $`M_{bh}\sigma `$ correlation, although extremely tight, has several disadvantages with respect to studying the black-hole masses of active galaxies. Firstly, even at low redshift it is not possible to obtain central stellar velocity dispersions for powerful quasars due to the presence of the unresolved nuclear point source. Furthermore, although it is possible to obtain velocity dispersion measurements for radio galaxies at low redshift, this will not be practicable at $`z>1`$ even with 8m-class telescopes. Consequently, before proceeding to discuss a more direct method of estimating black-hole mass in AGN, I will digress slightly to discuss the possibility that the intrinsic scatter in the $`M_{bh}M_{bulge}`$ relation is significantly lower than has previously been appreciated.

### 2.1 The intrinsic scatter in the $`M_{bh}M_{bulge}`$ relation

It has occurred to a number of authors that the high level of scatter associated with the $`M_{bh}M_{bulge}`$ relation could well be an artifact of the difficulties in determining accurate bulge luminosities for late-type galaxies. This issue is addressed in Fig 1 (taken from McLure & Dunlop 2002) which shows a comparison of the $`M_{bh}M_{bulge}`$ and $`M_{bh}\sigma `$ relations for a sample of 18 quiescent galaxies drawn from the list of objects with reliable black-hole mass measurements published by Kormendy & Gebhardt (2001). These 18 objects comprise those objects from the Kormendy & Gebhardt compilation which have pure E-type morphology, and therefore allow an investigation of the influence of incorrectly determined bulge luminosity upon the scatter in the $`M_{bh}M_{bulge}`$ relation. The best-fitting relations are shown as solid lines in both panels of Fig 1 and display an identical level of associated scatter; $`0.3`$ dex. This comparison demonstrates that when objects for which it is problematic to obtain accurate bulge luminosities are excluded, the scatter in the $`M_{bh}M_{bulge}`$ relation is comparable to that of the $`M_{bh}\sigma `$ relation. This fact has obvious consequences for the study of black-hole mass in radio-loud AGN, particularly with respect to radio galaxies for which it is straightforward to obtain relatively accurate bulge luminosities (eg. Bettoni et al. 2003). Furthermore, Graham et al. (2001) and Erwin et al. (2002) have recently demonstrated that, at least for quiescent galaxies, if sufficient care is taken over the modelling of the galaxy surface brightness profiles, then the scatter in both the $`M_{bh}M_{bulge}`$ relation and a correlation between black-hole mass and galaxy light concentration is only $`0.3`$ dex, irrespective of morphological type.

## 3 Estimating black-hole masses of low redshift AGN

Although the black-hole mass - host galaxy correlations discussed above are useful for estimating the black-hole masses of AGN, there is a more direct method of black-hole mass estimation which can be applied to broad-line AGN.

This method is the so-called virial black-hole mass estimator, and utilizes the width of broad emission lines in quasar spectra to directly trace the gravitational potential of the central black-hole on the scale of light-days. The underlying assumption behind the method is that the motion of the line-emitting material is virialized. Under this assumption the width of the broad lines can be used to trace the Keplerian velocity of the broad-line gas, and thereby allow an estimate of the central black-hole mass via the formula : $`M_{bh}=G^1R_{BLR}V_{BLR}^2`$, where $`R_{BLR}`$ is the broad-line region (BLR) radius and $`V_{BLR}`$ is the Keplerian velocity of the BLR gas. Recent evidence to support the Keplerian interpretation of AGN broad-line widths has been presented by Peterson & Wandel (2000) and Onken & Peterson (2002). These authors demonstrate, at least for several well studied Seyfert galaxies, that the FWHM of various emission lines follow the $`V_{BLR}R_{BLR}^{0.5}`$ relation expected for Keplerian motion. It is also worth noting that Peterson (2002) has recently demonstrated that the virial black-hole masses of reverberation mapped Seyfert galaxies lie on the same $`M_{bh}\sigma `$ relation as nearby inactive galaxies.

### 3.1 The $`M_{bh}M_{bulge}`$ relation in low redshift AGN

By utilizing the virial black-hole mass estimator it is possible to investigate whether the host galaxies of AGN do follow the same $`M_{bh}M_{bulge}`$ relation as inactive galaxies in the nearby Universe (Laor 2001, McLure & Dunlop 2002). In Fig 2 the results of a study to investigate the $`M_{bh}M_{bulge}`$ relation of 72 $`z<0.5`$ AGN by McLure & Dunlop (2002) are shown. Also shown is the sample of nearby inactive galaxies with reliable black-hole mass measurements previously discussed in Fig 1. The AGN sample ($`50`$ % radio-loud) consists of objects for which the host-galaxy bulge luminosity was determined via modelling of either HST or high-resolution ground-based imaging. It can clearly be seen that the AGN host galaxies follow a tight (0.4 dex scatter) correlation which is identical to that followed by nearby quiescent galaxies. Indeed, the best-fitting linear relation (dotted line) is equivalent to $`M_{bh}=0.0012M_{bulge}`$, identical to that determined by Merritt & Ferrarese (2001) for nearby quiescent galaxies.

## 4 The radio power - black hole mass connection

In the second section of this review I will proceed to examine recent evidence regarding the connection between black-hole mass and radio power. Figure 3 is taken from a recent paper by Dunlop et al. (2003) and shows 5 GHz radio luminosity plotted against black-hole mass for two samples. The first sample is comprised of 33 AGN, radio galaxies, radio-loud quasars and radio-quiet quasars, who’s bulge/black-hole masses have been determined from the analysis of HST imaging. The second sample features the nearby quiescent galaxies investigated by Franceschini et al. (1998), who found a tight correlation between black-hole mass and 5 GHz radio power of the form: $`P_{rad}M_{bh}^{2.5}`$. In Fig 3 there are three parallel relations plotted, each of the form $`P_{rad}M_{bh}^{2.5}`$, separated from each other by 2.5 decades. The lowest of these three relations appears to represent a lower limit to the radio output of a black-hole of a given mass.

There are two points concerning Fig 3 which are worth highlighting. Firstly, it can be seen that many of the radio-quiet quasars from the Dunlop et al. sample lie up against the apparent radio power lower limit. This is important because the radio-quiet and radio-loud quasars were selected by Dunlop et al. to have identical optical luminosities. Consequently, it can immediately be seen that the “radio-quietness” of these quasars is not due to their central engines being starved of fuel, because they are still producing large optical luminosities. The second point refers to the hypothesized upper envelope plotted in Fig 3. Clearly, if the Franceschini-type relation does represent a lower limit to the radio output of a given black-hole mass, then it is interesting to inquire about the form of the corresponding upper envelope. Although the form of the upper envelope plotted in Fig 3 is speculative, at least based on the data from Dunlop et al., further support is provided by the data from the study of Lacy et al. (2001).

Figure 4 is an adapted version of Fig 2 from Lacy et al. (2001), which again shows the $`P_{rad}M_{bh}`$ plane, but this time populated by objects from the First Bright Quasar Survey (FBQS), the PG quasar survey and several notable low redshift objects. Figure 4 clearly demonstrates that radio power and black-hole mass are related, although a large dynamic range is required in both parameters for the relation to become apparent.

The dotted line in Fig 4 is the best-fitting relation determined by Lacy et al. (not including accretion rate as an additional parameter) and has the form $`P_{rad}M_{bh}^{1.4}`$. However, in order to offer an alternative interpretation of the Lacy et al. data I have annotated Fig 4 with the addition of the upper and lower limits on radio-power suggested by Dunlop et al. (2003). The lower radio power limit describes the Lacy et al. data extremely well, and again shows radio-quiet quasars lying tight up against the relation defined by quiescent galaxies, despite being well supplied with fuel. Furthermore, the Lacy et al. data also obey the upper limit on radio-power over 4 orders of magnitude in black-hole mass and 11 orders of magnitude in radio luminosity. Given that objects appear to exist between radio power boundaries spanning 5 orders of magnitude, it is therefore understandable that studies which adopt the radio-loudness parameter $``$ tend to find little or no correlation between radio-loudness and black-hole mass (eg. Woo & Urry 2002).

One final interesting point is illustrated by Fig 4. Although there is clearly a large scatter associated with the $`P_{rad}M_{bh}`$ relation, it is still true that the most radio-loud quasars are those with the most massive black-holes. For example, Fig 4 shows that quasars with $`P_{rad}>10^{25}`$WHzsr<sup>-1</sup> are virtually exclusively associated with black-holes of mass $`>10^9\text{ }\mathrm{M}_{}`$. One consequence of this fact is that selecting samples based on extreme radio luminosity offers a method for selecting objects with the largest black-holes, and presumably bulge masses, at any given epoch. As a result, it should therefore be possible to cleanly trace the evolution of the $`M_{bh}M_{bulge}`$ relation in a sub-set of the most massive galaxies in the Universe.

## 5 Conclusions

The conclusions regarding the connection between black-hole mass and radio power can be summarizes as follows:

* There does not exist a threshold in black-hole mass above which an AGN must be radio-loud. There is a large overlap in the black-hole mass distributions of radio-loud and radio-quiet quasars, although at a given optical luminosity radio-loud quasars do appear to be biased to higher black-hole masses than their radio-quiet counterparts.
* Black-hole mass and radio luminosity are connected, although an extremely large dynamic range in both parameters is required for this connection to become apparent.
* The position of an object on the $`P_{rad}M_{bh}`$ plane appears to be a function of an additional physical parameter. Although accretion rate may be playing a role, the lack of a tight correlation between radio and optical luminosity suggests that some other parameter, such as black-hole spin, may be involved.
* Selecting the most radio-loud objects, as defined by their radio luminosity alone, is an effective method of isolating the objects with the most massive black-holes, and presumably host galaxies, at any epoch.
* Selecting the most radio-loud objects as a function of redshift provides an ideal method for cleanly studying the evolution of the $`M_{bh}M_{bulge}`$ relation within the sub-set of the most massive galaxies.