text:
# Young Star Clusters: Metallicities, Ages, and Masses

## 1. Why do we want to know GC and YSC masses?

GCs are conventionally believed to be (among) the oldest objects in the Universe, dating back to the times of galaxy formation. They are used to constrain the age of the Universe. LFs of GC systems are believed to be universal enough for determinations of distances to $`>20`$ Mpc and of the Hubble constant. On the other hand, “present-day GCs are the hardiest survivors of a larger original population” (Harris 1991). Hence, their observed LF is not only their LF at formation shifted by stellar evolutionary fading, but might be additionally modified by cluster destruction processes. Dynamical modelling of cluster systems in the Galactic potential shows that destruction processes and timescales strongly depend on cluster masses. While destruction by dynamical friction is more efficient for high mass clusters, tidal shocking and evaporation preferentially destroy low mass clusters.

Bright YSCs are observed in large numbers in interacting galaxies and merger remnants and a burning question with far-reaching implications is if these YSCs are young GCs. Masses derived for these YSCs are much higher than those of open clusters in the Milky Way. Effective radii of YSCs typically are a few pc similar to GC radii. The mass function (MF) of YSC systems in comparison with the MFs of molecular clouds or molecular cloud cores tells us about star and cluster formation processes. The MF of YSCs, as first pointed out by Meurer (1995), may differ in shape from their LF since M/L varies rapidly at young ages and the age spread within a YSC system is not much smaller than its age. The LF of open clusters in the Milky Way, the MFs of molecular clouds and molecular cloud cores, the observed LFs of YSCs, e.g. in NGC 4038/39, NGC 7252, NGC 3256, and of giant HII regions all are power laws with slopes in the range $`\alpha 1.5...1.8`$ (cf. Solomon et al. 1987, Lada et al. 1991, Kennicutt 1989, and reviews by Harris & Pudritz 1994, Elmegreen & Efremov 1997). Yet the LF of old GCs is Gaussian with typical parameters $`\mathrm{M}_\mathrm{V}7.3\mathrm{mag},\sigma (\mathrm{M}_\mathrm{V})1.3\mathrm{mag}`$, their MF is log-normal with typically $`\mathrm{Log}(\mathrm{M}/\mathrm{M}_{})5.5,\sigma 0.5`$ (e.g. Ashman et al. 1995).

Hence the question as to the MF of YSCs has profound implications. If the MF of YSCs were a power law like their LF and if YSC systems are to evolve into something similar to old GC systems, dynamical destruction processes would have to transform the power law MF into a log-normal MF over a Hubble time. If, on the other hand, the MF of YSCs were log-normal (and their LF distorted to a power law by age spread effects), then the star/cluster formation process would have to transform the power law MF of the molecular clouds into the log-normal MF of YSCs. Or, else, might already the MF of molecular clouds/cloud cores in violently star forming mergers (where to my knowledge it has not yet been observed) be different from what it is in quietly star forming ‘normal galaxies’ (where it is observed)?

## 2. Evolutionary synthesis of SSPs

Star clusters are simple stellar populations (SSPs), formed in one short burst of star formation with one metallicity Z. With evolutionary synthesis models, the time evolution of SSPs of various metallicities $`\mathrm{Z}_\mathrm{i}`$ is studied in terms of luminosities $`\mathrm{L}_\lambda `$, colours, spectrum, absorption features, stellar mass loss, and hence M/L by many groups (e.g. Bruzual & Charlot 1993, Worthey 1994, Bressan et al. 1994, F.-v.A. & Burkert 1995, Leitherer et al. 1999, Kurth et al. 1999). Basic parameters of this approach are the IMF and the set of stellar evolutionary tracks used (e.g. from the Padova or Geneva groups). All models agree that the changes in luminosities and colours are rapid in early evolutionary stages and become slower with increasing age. The colour evolution depends significantly on metallicity, already in very early stages, and the fading also depends on metallicity, in particular during the first Gyr (cf. F.-v.A. & Burkert 1995, Kurth et al. 1999).

## 3. YSC observations vs. SSP models

### 3.1. Procedure

If the metallicity of YSCs is known, individual ages can be obtained on the basis of their observed UBVI colours. Unfortunately, metallicity information from spectroscopy is only available yet for a handful of YSCs in NGC 7252 (Schweizer & Seitzer 1993, 1998) and NGC 1275 (Brodie et al. 1998). In all other cases we have to go back to some educated guess. In gas rich galaxy mergers, YSCs form out of the gas from the progenitor (spiral) galaxies, the abundance of which gives a lower limit to YSC abundances. E.g., for the YSCs in NGC 7252, which is a merger of two luminous Sc type spirals, we had predicted $`\mathrm{Z}_{\mathrm{YSC}}>\mathrm{Z}_{\mathrm{ISM}}^{\mathrm{Sc}}{}_{}{}^{1}/_{2}^{}\mathrm{Z}_{}`$ from the ISM abundance evolution in our 1-zone spiral models (F.-v.A. & Gerhard 1994). Spectroscopy of the brightest YSCs yielded $`\mathrm{Z}_{\mathrm{YSC}}\mathrm{Z}_{}`$ with some tentative indication of self-enrichment during the burst from the comparison of Mgb and Fe lines (cf. F.-v.A. & Burkert 1995). Distances are known for all YSC systems. Hence absolute luminosities $`\mathrm{L}_\mathrm{V}`$ can be combined with $`\mathrm{M}/\mathrm{L}_\mathrm{V}`$ values from SSP models of appropriate age and metallicity to derive the masses of YSCs.

### 3.2. Sources of uncertainty

Sources of uncertainty for these mass estimates come both from observations and models. Observational errors on luminosities and colours might be inhomogeneous and probably not independent of each other, dust extinction is not known for individual YSCs, the dust distribution in some systems is inhomogeneous, the metallicity is, at most, known for a small subsample of YSCs and might show an intrinsic scatter, and, finally, the completeness limit need not be homogeneous over the region of YSC observations. Intrinsic model uncertainties are estimated to be $`<0.1`$ mag in (optical) luminosities and colours. Differences between models from various authors are mainly due to differences in the stellar input physics. While serious colour discrepancies at very young ages $`10`$ Myr are seen, e.g. comparing models from Bruzual & Charlot with those of Leitherer et al. , probably due to the inclusion/non-inclusion of emission lines, fairly good agreement is reached among models from various groups for the same metallicity and IMF at all ages $`>60100`$ Myr. E.g., $`\mathrm{\Delta }(\mathrm{V}\mathrm{I})_{12\mathrm{Gyr}}<0.1`$ mag and $`\mathrm{\Delta }(\mathrm{M}/\mathrm{L}_\mathrm{V})_{12\mathrm{Gyr}}<10\%`$.

$`\mathrm{M}/\mathrm{L}_\lambda `$, however, depends on wavelength $`\lambda `$, metallicity Z, and IMF, i.e. on its slope and lower mass limit. In Tab.1, I briefly sketch out these dependencies as obtained from our models using Padova stellar evolutionary tracks for stars in the mass range $`0.160\mathrm{M}_{}`$ (cf. Kurth et al. 1999) and two different IMFs (Salpeter vs. Scalo 1986).

In general and as reported by others before, model $`\mathrm{M}/\mathrm{L}_\mathrm{V}`$ are about twice as large as are $`\mathrm{M}/\mathrm{L}_\mathrm{V}`$ values derived from observations of old GCs, even when a Scalo IMF is assumed. For a Salpeter IMF the discrepancy is higher by another factor of two. Observational $`\mathrm{M}/\mathrm{L}`$ values are obtained by measuring the central velocity dispersion $`\sigma _0`$, central surface brightness $`\mathrm{I}_0`$, and half light radius $`\mathrm{r}_{1/2}`$, and assuming isotropic orbits, no radial gradients in $`\mathrm{M}/\mathrm{L}`$, and no DM halos around GCs. Then

$`\mathrm{M}/\mathrm{L}\frac{\sigma _0^2}{\mathrm{I}_0\mathrm{r}_{1/2}}`$

and typical values quoted for Galactic GCs are $`\mathrm{M}/\mathrm{L}_\mathrm{V}2`$. There are two effects that might invalidate the assumptions going into this derivation. First of all, mass segregation in GCs (cf. Meylan this conf.) will result in an $`\mathrm{M}/\mathrm{L}`$ increasing with radius, as e.g. observed by Côté et al. 1995 for NGC 3201. Second, low mass stars with their high $`\mathrm{M}/\mathrm{L}`$ values are preferentially lost by evaporation (e.g. Gerhard this conf.). While both processes are undoubtedly at work in GCs, attempts to examine quantitatively the validity of the assumptions involved seem to still give ambiguous results. Leonard et al. 1992 use proper motion data and radial velocity measurements for stars in the GC M13 and find that $`50\%`$ of the mass of M13 is in low mass stars and brown dwarfs. Including this unseen mass leads to an increased $`\mathrm{M}/\mathrm{L}4`$, a value well compatible with the models. On the other hand, observations of tidal tails on some GCs in the Milky Way potential are interpreted to indicate that there is not much DM around those GCs, constraining their mass-to-light ratios to $`\mathrm{M}/\mathrm{L}<2.5`$ (Moore 1996).

## 4. YSCs in the Antennae: a 1st example

The Antennae galaxies (= NGC 4038/39) is an interacting pair of two gas rich spirals, probably Sc, of comparable mass. In the ongoing starburst triggered by the interaction a population of bright YSCs is formed, numerous enough to allow for the first time to reasonably define a LF. The LF from WFPC1 observations is a power law with slope $`\alpha 1.8`$ (Whitmore & Schweizer 1995 (WS95)). WFPC2 reobservations are going to be presented by Miller (this conf.).

Assuming a homogeneous metallicity $`\mathrm{Z}1/2\mathrm{Z}_{}`$ lack of individual cluster spectroscopy and in analogy to the YSC system in NGC 7252, I analysed the WFPC1 data of WS95 with evolutionary synthesis models for SSPs. In a first step, the average dereddened $`(\mathrm{V}\mathrm{I})`$ colour of the 550 YSCs is used to derive a mean age of the YSC population of $`(2\pm 2)10^8`$ yr, consistent with Barnes’ (1988) dynamical model for the interaction between NGC 4038 and 4039 and consistent with Kurth’s (1996) global starburst age. SSP models describe both the fading and the reddening of the YSC population as it ages. Assuming a mean age for the YSC population, both the LF and the colour distribution of the YSCs are simply shifted towards fainter magnitudes and redder colours, respectively, without changing shape.

### 4.1. Evolution of the YSC LF

In a second step, individual ages are derived for all the YSCs on the basis of their individual $`(\mathrm{V}\mathrm{I})`$ and $`(\mathrm{U}\mathrm{V})`$ colours. Interestingly, the resulting age distribution not only shows a peak at very young ages $`0410^8`$ yr for the YSCs, but also a contribution from $`12`$ Gyr old GCs from the parent galaxies (Fig. 1a). A small number of apparent interlopers are probably due to inhomogenities in the dust distribution. It is improbable that many of the old GCs we identify are highly reddened YSCs, since they would have to be exeedingly bright intrinsically. The fact that age estimates from $`(\mathrm{U}\mathrm{V})`$ colours, as far as available, agree with age estimates from $`(\mathrm{V}\mathrm{I})`$ supports our metallicity assumption and makes us hope that the average $`\mathrm{E}_{\mathrm{B}\mathrm{V}}`$ is correct for most of the clusters. It is clear, however, that the individual YSC extinctions are a major source of uncertainty in this analysis. With individual YSC ages and observed luminosities, SSP models can be used to predict the time evolution of cluster luminosities. Neglecting any kind of dynamical cluster destruction effects, and only using the young star clusters identified, we obtain the surprising result that by an age of 12 Gyr, when age differences among individual YSCs will be negligible, their LF will have evolved from the presently observed power law into a fairly normal Gaussian GC LF with parameters $`\mathrm{M}_\mathrm{V}=6.9`$ mag and $`\sigma (\mathrm{M}_\mathrm{V})=1.3`$ mag (Fig. 1b). The turn-over then is $`>1`$ mag brighter than the completeness limit which evolves to $`\mathrm{M}_\mathrm{V}=5.7`$ mag, and fainter by $`0.4`$ mag than for typical GC systems. This is a consequence of the enhanced metallicity of the YSC population with respect to that of old GC systems, and in agreement both with observations of old GC systems with a range of metallicities (Ashman et al. 1995) and with our SSP model predictions. With this surprising result we confirm and quantify Meurer’s conjecture that over a Hubble time, age spread effects can transform an observed power law LF of YSCs into the Gaussian LF of old GCs (cf. F.-v.A. 1998 for details).

The bright end of the LF defined by the old GCs from the parent spirals is well described by a Gaussian LF with parameters $`\mathrm{M}_\mathrm{V}=7.3`$ mag and $`\sigma (\mathrm{M}_\mathrm{V})=1.2\mathrm{mag},`$ normalised to the total number of GCs in the Milky Way and Andromeda galaxies together (Fig. 2a). Hence, if the two interacting spirals NGC 4038 and 4039 had a similar number of GCs as those galaxies, and if the bulk of the YSC population really are young GCs, Zepf & Ashman’s (1993) requirement, that the number of secondary GCs formed in mergers should be comparable to the number of primary GCs in the two progenitor spirals, would be fulfilled. This requirement is necessary if the higher specific GC frequency in ellipticals – as compared to spirals – is to be compatible with a spiral-spiral merger origin for those ellipticals.

### 4.2. MF of YSCs

To investigate whether the LFs and MFs of old GC systems are determined by the cluster formation process or whether they are the result of secular dynamical destruction processes, we derive the MF of the very young star cluster system in the Antennae. We restrict ourselves to YSCs brighter than the completeness limit, use individual YSC ages and luminosities, and combine them with model M/L for the respective YSC ages to determine their individual masses. Our tentative conclusion is that for all the 393 YSCs brighter than the completeness limit and with $`(\mathrm{V}\mathrm{I})`$ colours available, the MF is compatible with a Gaussian MF with parameters $`\mathrm{Log}(\mathrm{M}_{\mathrm{YSC}}/\mathrm{M}_{})=5.6,\sigma 0.46`$ $`\mathrm{Mass}(\mathrm{YSC})410^5\mathrm{M}_{}`$ (Fig. 2b). If we include YSCs fainter than the completeness limit, neither the shape nor the parameters of the MF are changed (cf. F.-v.A. 1999 for details).

The uncertainty in the YSC metallicity leads to age uncertainties, which, in turn, lead to uncertainties in model M/L values of the order of 10 % in these early stages. Inhomogenities in the dust distribution do not seem to be very important for the brightest YSCs for which ages from $`(\mathrm{U}\mathrm{V})`$ agree with ages from $`(\mathrm{V}\mathrm{I})`$ colours. We do not know, however, how important they are for those YSCs that are not detected in U. Hence, on the basis of $`(\mathrm{V}\mathrm{I})`$ alone, we cannot quantify in how far an inhomogeneous internal reddening might affect our results. For this reason and since our analysis is based on WFPC1 data, we caution that our conclusions concerning the evolution of the YSC LF and their MF can only be preliminary. As soon as the reduced WFPC2 data presented by Miller will become available to us, we shall repeat our analysis and supplement it with simulations to estimate the effects of differential observational uncertainties.

### 4.3. Implications

If, however, our preliminary result became confirmed, this would imply that the log-normal MF of old GCs is produced by the cluster formation process rather than by secular dynamical evolution of the cluster system. In this context. it seems very interesting to obtain observational information about the molecular cloud mass spectrum in massive interacting galaxies. Jog & Solomon (1992) conjecture that the strongly enhanced ambient pressure in massive gas rich mergers might affect the molecular cloud structure. In Ultraluminous Infrared Galaxies – which all are mergers with strong starbursts – the fraction of gas at very high densities of $`\mathrm{n}10^4\mathrm{and}10^5\mathrm{cm}^3`$, as traced by HCN and CS lines, with respect to gas at $`\mathrm{n}500\mathrm{cm}^3`$, as traced by CO, is indeed observed to be higher by 1 – 2 orders of magnitude (e.g. Solomon et al. 1992).

### 4.4. Discussion

Even with deeper WFPC2 data, however, the MF of the YSCs in the Antennae does not yet seem to be unambiguously settled. Zhang & Fall (1999) use reddening-free colour indices Q1, Q2, solar metallicity Bruzual & Charlot models, consider incompleteness and stellar contamination, and find power law MFs with slopes $`\alpha 2`$ for YSCs in the two age intervals 2.5 – 6.3 Myr and 25 – 160 Myr where Q1 and Q2 give unambiguous results. In an independent analysis of the same data, Miller (this conf.) finds a significant flattening in the observed LF at $`\mathrm{M}_\mathrm{V}10.7`$ mag. Intriguingly, this value is exactly the turnover luminosity implied by our MF at the mean YSC age in the absence of an age spread. Miller’s method of reconstruction, however, leads to a power law YSC MF. Clearly, more work both from the modelling and observational sides is needed before we really understand the YSC MF.

### 4.5. Other systems

For the YSCs in NGC 7252 and NGC 3921, Miller & Fall (1997) and Miller (this conf.) prefer power law MFs. Distances to those systems, however, are larger by factors 3 and 4, respectively, pushing the completeness limit to higher luminosities, even in WFPC2 data. Moreover, the YSCs in these dynamically old merger remnants have higher mean ages (650 – 750 Myr for NGC 7252, 250 – 750 Myr for NGC 3921, as compared to 200 Myr for the Antennae), and hence, are intrinsically fainter already by 0.8 – 1 mag, on average. If we assume that the MFs for the YSCs in NGC 7252 and NGC 3921 were identical to the one we derived for the YSCs in the Antennae, we obtain turn-over luminosities $`\mathrm{M}_\mathrm{V}=10\mathrm{}9.5`$ mag for NGC 7252 and $`\mathrm{M}_\mathrm{V}=9.5\mathrm{}9`$ mag for NGC 3921 at their respective mean YSC ages. In both cases, these turn-over luminosities are close to the completeness limits. The difficulty to disentangle the old GC and the YSC populations also increases with increasing YSC age. We therefore doubt that it is possible to track the MF beyond a possible turnover in NGC 7252 and 3921. In any case, the YSCs in NGC 7252 and 3921 have survived $`10`$ crossing times and from this fact alone are young GCs rather that open ones, as also indicated by their combination of small effective radii and high luminosities (cf. Miller et al. 1997, Schweizer et al. 1996).

### 4.6. Dynamical Evolution

Stellar mass loss is included in our models for the evolution of M/L. Over 12 Gyr, the masses of YSCs will decrease by $`10\mathrm{and}15\%`$ for a Salpeter and Scalo-IMF, respectively. About half of the entire stellar mass loss occurs during the $`1^{\mathrm{st}}`$ Gyr. External dynamical effects from an interaction of the clusters with the potential of the interacting galaxy system seem extremely difficult to model. Comparison of YSC LFs and MFs in an age sequence of interacting galaxies, mergers, merger remnants, and dynamically young ellipticals will allow to “see these processes at work”.

In the Milky Way potential, Vesperini (1998) has shown that an assumed initially log-normal GC MF is conserved in shape and parameters during self-similar evolution over a Hubble time despite the destruction of $`50\%`$ of the cluster population. If, on the other hand, he starts with a power law MF, severe fine-tuning is required for his model parameters to secularly transform it into the observed log-normal MF of old GC systems.

It is clearly important to analyse more YSC systems in order to see if (and in how far) their MFs are universal or might depend on environment. Old GC systems have their turn-over around $`\mathrm{M}_\mathrm{V}7.2`$ mag. With 10m telescopes they are accessible to more than Virgo cluster distances. MOS – e.g. with FORS on the VLT – in combination with HST imaging will allow to determine cluster abundances and hence to more precisely age-date them, and may even provide kinematic information for independent mass estimates.

An open question seems to me if the (globular) cluster formation process in the high metallicity environment of interacting spirals today is the same or not as it was in the Early Universe when the radiation field was stronger and the metallicity lower.

#### Acknowledgments.

I wish to thank the organisers and in particular Ariane Lançon for a stimulating workshop in a very pleasant and warm atmosphere in cold and rainy Strasbourg.

## References

Ashman, K. M., Conti, A., Zepf, S. E., 1995, AJ 110, 1164

Barnes, J. E., 1988, ApJ 331, 699

Bressan, A., Chiosi, C., Fagotto, F., 1994, ApJS 94, 63

Brodie, J. P., Schrd̈er, L. L., Huchra, J. P., et al. 1998, AJ 116, 691

Bruzual, G. A., Charlot, S., 1993, ApJ 405, 538

Côté, P., Welch, D. L., Fischer, P., Gebhardt, K., 1995, ApJ 454, 788

Elmegreen, B. G., Efremov, Y. N., 1997, ApJ 480, 235

Fritze – v. Alvensleben, U., 1998, A&A 336, 83

Fritze – v. Alvensleben, U., 1999, A&A 342, L25

Fritze – v. Alvensleben, U., Burkert, A., 1995, A&A 300, 58

Fritze – v. Alvensleben, U., Gerhard, O. E., 1994, A&A 285, 751 + 775

Harris, W. E., 1991, ARA&A 29, 543

Harris, W.E., Pudritz, R. E., 1994, ApJ 429, 177

Jog, C. J., Solomon, P. M., 1992, ApJ 387, 152

Kennicutt, R. C., 1989, ApJ 344, 685

Kurth, O., 1996, Diploma Thesis, Univ. Göttingen

Kurth, O., Fritze - v. Alvensleben, U., Fricke, K. J., 1999, A&AS 138, 19

Lada, E., Bally, J., Stark, A. A., 1991, ApJ 368,432

Leitherer, C., et al. 1999, ApJS 123, 3

Leonard, P. J. T., Richer, H. B., Fahlman, G. G., 1992, AJ 104, 2104

Meurer, G. R., 1995, Nat 375, 742

Miller, B. W., Fall, S. M., 1997, AAS 119, #115.04

Moore, B., 1996, ApJ 461, L13

Scalo, J. M., 1986, Fundam. Cosm Phys. 11, 1

Schweizer, F., Seitzer, P., 1993, ApJ 417, L29

Schweizer, F., Seitzer, P., 1998, AJ 116, 2206

Schweizer, F., Miller, B. W., Whitmore, B. C., Fall, S. M., 1996, AJ 112, 1839

Solomon, P. M., Rivolo, A. R., Barrett, J., Yahil, A., 1987, ApJ 319, 730

Solomon, P. M., Downes, D., Radford, S. J. E., 1992, ApJ 387, L55

Vesperini, E., 1998, MN 299, 1019

Whitmore, B.C., Schweizer, F., 1995, AJ 109, 960 (WS95)

Worthey, G., 1994, ApJS 95, 107

Zepf, S. E., Ashman, K. M., 1993, MN 264, 611

Zhang, Q., Fall, S. M., 1999, ApJ 527, L81


format: text

text:
# Punctured polygons and polyominoes on the square lattice

## 1 Introduction

A self-avoiding polygon (SAP) can be defined as a walk on a lattice which returns to the origin and has no other self-intersections. Alternatively we can define a SAP as a connected sub-graph (of a lattice) whose vertices are of degree 2. The history and significance of this problem is nicely discussed in . Generally SAP are considered distinct up to a translation, so if there are $`p_n`$ SAP of length $`n`$ there are $`2np_n`$ returns to the origin (the factor of two arising from the two possible directions in which the loop can be traveled). Staircase polygons are a special case of SAP. A staircase polygon can be viewed as the intersection of two directed walks starting at the origin, taking steps either to the right or down and terminating once the walks meet. We define a punctured SAP as a SAP with one or more holes, with the perimeter of each hole itself being a SAP. In other words a punctured SAP is a SAP enclosing one or more SAP. To avoid any possible confusion in our definition of punctured polygons, note that the punctures are disjoint — there are no degree four vertices in the punctured polygons we are considering. A similar definition can be given for punctured staircase polygons. We have shown an example of each case in figure 1. The two principal questions one can ask are “how many polygons, distinct up to a translation, are there of perimeter $`2n`$ (including the perimeter of the holes) with $`k`$ punctures?” and “how many polygons, distinct up to a translation, are there of area $`n`$ with $`k`$ punctures?”

A polyomino is defined as the union of connected cells, where a cell is a unit square with 4 edges (and its interior). Two cells are said to be joined if they share a common edge, and are said to be connected if there exists a sequence of cells joining the two cells such that successive pairs of cells are joined. A $`k`$-punctured polyomino is a polyomino with $`k`$ holes. Unlike the situation for punctured polygons, degree 4 vertices are permitted, but two holes meeting only at a vertex do count as two holes, rather than one, as they are not joined. Punctured polygons are a proper subset of punctured polyominoes. The difference between them is that punctured polygons do not include those polyominoes where there are 2 cells touching each other only at a vertex. An example is shown in figure 2.

Reverting to $`k`$-punctured polygons, for $`k=0`$ we have the simpler — but still unsolved — problem of the number of polygons of perimeter $`2n`$ or area $`n.`$ Both these problems have been extensively studied for several decades. The most recent result for polygon perimeters appears to be  where polygons of perimeter up to 90 steps are given. In that paper our analysis of the polygon perimeter generating function led us to conclude that

$$P^{(0)}(x)=\underset{n}{}p_{2n}^{(0)}x^nA^{(0)}(x)+B^{(0)}(x)(1\mu ^2x)^{2\alpha },$$

where $`p_{2n}^{(0)}`$ is the number of SAP of perimeter $`2n`$, distinct up to translations. The analysis in  yielded a very accurate estimate for the connective constant $`\mu =2.63815853034(10)`$ and confirmed the theoretical prediction $`\alpha =1/2`$ . Furthermore, we obtained estimates for the critical amplitudes $`A^{(0)}(x_c)0.036`$ and $`B^{(0)}(x_c)0.234913,`$ where $`x_c=1/\mu ^2`$. We also concluded that there was no evidence for a non-analytic correction-to-scaling exponent, so that we expect the asymptotic form of the coefficients to behave as:

$$p_{2n}^{(0)}\mu ^{2n}n^{\frac{5}{2}}[a_1+a_2/n+a_3/n^2+a_4/n^3+\mathrm{}].$$

The connective constant $`\mu `$ is of course the same as that for self-avoiding walks on the same lattice .

For polygon areas the most recent published work appears to be  in which the first 20 terms of the area generating function were given and analysed. In that work it was found that

$$A^{(0)}(y)=\underset{n}{}a_n^{(0)}y^nC^{(0)}(y)+D^{(0)}(y)\mathrm{log}(1\kappa y),$$

where $`a_n^{(0)}`$ is the number of SAP of area $`n`$, $`\kappa 3.97087,`$ and the amplitudes $`C^{(0)}`$ and $`D^{(0)}`$ were not estimated. Recently this series has been extended to 26 terms , but in the present work we have devised a new and exponentially faster algorithm, as a result of which we have extended the series to 42 terms, and we present an analysis of this longer series. The connective constant $`\kappa `$ is found to be slightly smaller than that for the related problem of polyominoes .

In the following, we refer to the boundary of a polygon and its interior as a disc and so we will be discussing punctured discs. An unpunctured disc is a SAP.

For punctured discs, the basic problem is, analogously, the calculation of the generating functions

$$P^{(k)}(x)=\underset{n}{}p_{2n}^{(k)}x^nB^{(k)}(x)+C^{(k)}(x)(1(\mu ^{(k)})^2x)^{2\alpha _k},$$
(1)
and

$$A^{(k)}(y)=\underset{n}{}a_n^{(k)}y^nD^{(k)}(y)+E^{(k)}(y)(1\kappa ^{(k)}y)^{\beta _k},$$
(2)
where the superscript $`k`$ refers to the number of holes, or punctures. From the generating function, one wishes to deduce the asymptotic behaviour, believed to be as shown on the r.h.s. of the above equations. The major problem to be investigated is how the behaviour of $`P^{(k)}(x)`$ and $`A^{(k)}(y)`$ changes as $`k`$ is increased. Previous work  has been confined to the study of punctured SAP by area. There it was proved that $`\kappa ^{(k)}=\kappa ^{(0)}=\kappa ,`$ and that if the exponent exists, $`\beta _k=\beta _0+k.`$ These results apply more generally to punctured surfaces, but in this work we are confining ourselves to two dimensions. As far as we are aware, there has been no previous work on the problem of the perimeter generating function of punctured discs.

The problem is interesting for several reasons. The effect of a change in geometry is a much studied topic in lattice statistics, and our study of the change in perimeter exponent with punctures seems to be entirely new. It has only been possible by the algorithms we have designed and implemented, which are exponentially faster than pre-existing algorithms. A number of related models have been studied previously, such as $`c`$-animals  and the behaviour of prime knots in polygons . $`c`$-animals are lattice animals with exactly $`c`$ cycles. In  it was proved that if the number of such animals $`a_n(c)\lambda _c^nn^{\theta _c}`$ as $`n\mathrm{}`$, $`c`$ fixed, then $`\theta _c=\theta _0+c`$ provided that $`\theta _0`$ exists. It had been previously proved  that $`\lambda _c=\lambda _0.`$ The change in connective constant of $`c`$-animals as the number of cycles per vertex changes from zero to non-zero is discussed — among other results — in . Similarly, in a numerical study of knotted polymers , it was conjectured that the exponent $`\alpha `$ depends on the number of prime knots $`n_p`$ that arise in the knot decomposition of a given SAP via the relation $`\alpha (n_p)=\alpha (0)+n_p.`$

Furthermore, there is considerable pedagogical connection between some of these previous studies and our work here, and also between our work and the study of branched polymers. To sketch this connection, we first remark that the number of polyominoes, also called lattice animals, is just the number of (strongly embedded) site animals. This connection is readily seen by placing a site at the centre of every cell of the polyomino, and joining those sites corresponding to joined cells by a bond. In this way, every lattice animal is mapped to a distinct site animal and vice versa. All other models we are considering map similarly to a different subset of strongly embedded site animals. For example, SAP map to site animals whose only cycles are 4-cycles (which may be isolated or joined). Punctured polygons map to site animals with larger cycles, as well as 4-cycles, whereas punctured polyominoes correspond to site animals with more complex restrictions. Thus this study also complements the earlier studies of $`c`$-animals. In those studies, the variation of exponent with the number of cycles is considered, whereas in this study we are varying the types of allowed cycles.

In order to study these and related systems, when an exact solution can’t be found one has to resort to numerical methods. For many problems the method of series expansions is by far the most powerful method of approximation. For other problems Monte Carlo methods are superior. For the analysis of $`P^{(k)}(x)`$ and $`A^{(k)}(y)`$, series analysis is undoubtedly the most appropriate choice. This method consists of calculating the first few coefficients in the expansion of the generating function. Given such a series, using the numerical technique known as differential approximants , highly accurate estimates can frequently be obtained for the critical point and exponents, as well as the location and critical exponents of possible non-physical singularities. Other numerical methods are discussed in , and those used in the present study are described more fully below.

In the next section we will describe the finite lattice method for enumerating punctured polygons. In Section 3 we prove the invariance of $`\mu ^{(k)}`$ as $`k`$ changes, and give an heuristic argument for the exponent shift with $`k.`$ The results of the analysis of the series are presented in Section 4. Results analogous to those known for punctured polygons by area are proved for punctured polyominoes in the Appendix.

## 2 Enumeration of punctured polygons

### 2.1 Enumeration of punctured self-avoiding polygons

The method used to enumerate punctured self-avoiding polygons on the square lattice is a generalisation of the method devised by Enting  for the enumeration of ordinary SAP. In the following we first describe the original method in some detail and show how simple it is to generalize the method to the enumeration of punctured polygons. The first terms in the series for the generating function can be calculated using transfer matrix techniques to count the number of polygons in rectangles $`W+1`$ edges wide and $`L+1`$ edges long. The transfer matrix technique involves drawing a line through the rectangle intersecting a set of $`W+2`$ edges. For each configuration of occupied or empty edges along the intersection we maintain a (perimeter) generating function for loops to the left of the line cutting the intersection in that particular pattern. Polygons in a given rectangle are enumerated by moving the intersection so as to add one vertex at a time, as shown in figure 3. Since the loops are non-intersecting, each configuration can be represented by an ordered set of edge states $`\{n_i\}`$, where

$$n_i=\{\begin{array}{cc}\hfill 0& \text{empty edge},\hfill \\ \hfill 1& \text{lower part of loop closed to the left},\hfill \\ \hfill 2& \text{upper part of loop closed to the left}.\hfill \end{array}$$

Configurations are read from the bottom to the top. So the configuration along the intersection of the polygon in figure 3 is $`\{0112122\}`$. In passing it is worth noting that there are some major restrictions on the possible configurations. Firstly, since all loop-ends are connected to the left of the intersection, every lower loop end must have a corresponding upper end, and it is therefore clear that the total number of ‘1’s is equal to the total number of ‘2’s. Secondly, as we look through the configuration starting from the bottom the number of ‘1’s is never smaller than the number of ‘2’s.

In Table 1 we have listed the possible local ‘input’ states and the ‘output’ states which arise as the kink in the intersection is propagated by one step. Some of these update rules are illustrated further in Fig 4. The first panel represents the input states ‘10’ and ‘01’ and the possible output states are also ‘01’ and ‘10’. The second panel represents the input state ‘11’ as part of the configuration $`\{01122\}`$. In this case we connect the two loop ends, but in doing so we see that the upper part of the second loop before the move becomes the lower part of the one remaining loop after the move. That is the configuration $`\{01122\}`$ becomes $`\{00012\}`$. This relabelling of the other loop-end when connecting two ‘1’s (or two ‘2’s) is denoted by over-lining in Table 1. In general there could be more loops nested in between the two ‘1’s and the corresponding ‘2’ at the other end of the loop. Say for instance we had the configuration $`\{11121222\}`$ and connected the first two ‘1’s then the new configuration of unconnected loop ends would be $`\{00121212\}`$ (drawing a little figure makes this quite clear). The general rule for the relabelling is as follows: When connecting two ‘1’s (‘2’s) we work upward (downward) in the configuration, counting the number of ‘1’s and ‘2’s we pass until the number of ‘2’s (‘1’s) exceeds the number of ‘1’s (‘2’s). This ‘2’ (‘1’) is the other end of the inner loop and it should now be changed to a ‘1’ (‘2’), thus becoming the lower (upper) end of the outer loop (again drawing a few pictures should make this clearer). The weights corresponding to these configuration transformations are simply calculated by counting the number of steps which have been added to the polygon. Note that the input state ‘12’ is special because connecting the two ends results in a closed loop, so this is only allowed if there are no other loops cut by the intersection and the result is a valid polygon, which is then accumulated in the total count for that particular length. Failure to observe this restriction would result in graphs with disconnected components, either one polygon over another or a polygon within another (this latter case is of course of interest when we wish to enumerate punctured polygons). This is illustrated in figure 5 where we show the possible ways a pair of loops can be placed relative to one another (first panels), how the loops can be connected to produce a valid configuration (middle panels) and the ways of connecting loops that lead to invalid graphs containing disconnected components (last panels). In this figure the invalid SAP in the last panel on the bottom could result in a valid punctured SAP. Note also that from the input state ‘00’ we can produce the output state ‘12’ only if there are other loops crossing the intersection (otherwise we would produce disconnected polygons sitting side by side). We refer the interested reader to  for further details regarding the encoding and relabelling of configurations.

Due to the obvious symmetry of the lattice one need only consider rectangles with $`LW`$. In the original approach  valid polygons were required to span the enclosing rectangle in the lengthwise direction. So it is clear that polygons with projection on the $`y`$-axis $`<W`$, that is polygons which are narrower than the width of the rectangle, are counted many times. It is however easy to obtain the polygons of width exactly $`W`$ and length exactly $`L`$ from this enumeration . Any polygon spanning such a rectangle has a perimeter of length at least $`2(W+L)`$. By adding the contributions from all rectangles of width $`WW_{\mathrm{max}}`$ (where the choice of $`W_{\mathrm{max}}`$ depends on available computational resources) and length $`WL2W_{\mathrm{max}}W+1`$, with contributions from rectangles with $`L>W`$ counted twice, the number of polygons per vertex of an infinite lattice is obtained correctly up to perimeter $`4W_{\mathrm{max}}+2`$.

With the original algorithm the number of configurations required as $`W_{\mathrm{max}}`$ increased grew asymptotically as $`3^{W_{\mathrm{max}}}`$ . In a recent improvement of the algorithm  valid polygons were required to span the rectangle in both directions. In other words we directly enumerate polygons of width exactly $`W`$ and length $`L`$. For each configuration of partially completed polygons we keep track of the current minimum number of steps $`N_{\mathrm{cur}}`$ that have been inserted to the left of the intersection and we calculate the minimum number of additional steps $`N_{\mathrm{add}}`$ required to produce a valid polygon that spans a rectangle of size at least $`W\times W`$. If the sum $`N_{\mathrm{cur}}+N_{\mathrm{add}}>4W_{\mathrm{max}}+2`$ the partial generating function for that configuration was discarded because it would make no contribution to the polygon count up to the perimeter lengths we were trying to obtain. Numerical evidence indicated that the computational complexity was reduced significantly. While the number of configurations still grew exponentially as $`\lambda ^{W_{\mathrm{max}}}`$ the value of $`\lambda `$ was reduced from $`\lambda =3`$ to $`\lambda 2`$ with the improved algorithm. Furthermore, for any $`W`$ we know that contributions will start at $`4W`$ since the smallest polygons have to span a $`W\times W`$ rectangle, so for each configuration we need only retain $`4(W_{\mathrm{max}}W)+2`$ terms of the generating functions while in the original algorithm contributions started at $`2W+2`$ because the polygons were required to span only the length-wise direction.

The generalization to enumeration of punctured polygons is obtained by noting that a closed loop is formed whenever we connect a 1-edge to a 2-edge immediately above. If these two edges were the only loop ends in the intersection we would have formed a valid polygon. In other cases we need to ensure that the resulting polygon is a valid punctured polygon. That is, we must ensure that the separate SAP just formed will be completely enclosed within a larger polygon. So we wish to avoid forming separable polygons, as shown in the upper right panel of figure 5, and ‘holes within holes’. As it turns out the rule for a valid ‘12’ closure is simply that we can connect the two loop ends provided there is an odd number of loop ends below the loop being closed. To see this consider that as we go through a configuration we note that each time we pass a loop end we go from the outside of the polygon to the inside and visa versa. So in a configuration all lattice cells between the first and second loop ends will lie inside the finished polygon, lattice cells between the second and third loop end will lie outside the polygon, and so on. Thus we see that by closing a loop which has an odd number of loop ends below it we are closing off a part of the lattice which will lie inside the finished polygon. In particular we see that we avoid the situation shown in the upper right panel of figure 5 with graphs containing disconnected pieces one over another. Likewise we avoid creating graphs with disconnected pieces sitting side by side.

We also avoid forming holes within holes because closing a loop around the hole would be prohibited since there would be an even number of loop ends below the ‘1’-edge (except of course when forming a completed punctured polygon). Let us look at the possible edge-configurations around a puncture in some detail. First look at the configuration {$`\mathrm{}1122\mathrm{}`$} (where the $`\mathrm{}`$ are any edge configurations with an even number of edges making the total configuration of edges a valid intersection). The outer ‘12’-edges can’t be connected (the number of edges below the ‘1’ is even) so we have to connect two other edges on either side of the hole diminishing the number of edges by two and possibly changing the edge labels on either side of the hole (in which case we end up with one of the subsequent cases). In the top panel of figure 6 we show how the simplest interesting case {11112222} leads to only valid punctured polygons. Secondly, in the case {$`\mathrm{}1121\mathrm{}`$}, we can connect the two ‘1’-edges forming a partial loop enclosing the hole. In doing so the matching ‘2’-edge of the second ‘1’-edge is changed to a ‘1’-edge because it now is the new lower edge of the larger loop formed by connecting the two original ‘1’-edges. The important thing to notice is that the number of edges below this new ‘1’-edge is even so we cannot connect it to a ‘2’-edge immediately above and we therefore do not form a hole within a hole. In the middle panel of figure 6 we show the simplest interesting case {1111212222} and demonstrate that it leads to valid punctured polygons. Thirdly, the configuration $`\{\mathrm{}2121\mathrm{}`$} is not very interesting since connecting the ‘21’-edges in front of the hole does not result in a loop partially enclosing the hole. Finally, the configuration {$`\mathrm{}2122\mathrm{}`$} is obviously merely the mirror image of the second case. Note that this covers all cases since more complicated configurations merely would correspond to more convoluted loop structures.

In this work we use the generalisation of the algorithm of  to count the number of punctured SAPs with up to 3 holes. Obviously, the smallest hole we can make in a SAP has perimeter 4 so the number of punctured SAPs is obtained correctly up to perimeter $`4W_{\mathrm{max}}+2+4k`$.

The algorithm used for the enumeration of punctured SAPs by area is a simple variation of the algorithm described above. The encoding of configurations along the intersection and the transformation of these configurations as the intersection is moved remain the same. The only change is that the weights are different. In order to count the enclosed area we proceed as follows: A unit of area may be added as the kink is moved to a new lattice cell (in figure 3 this is the cell in which the dotted lines meet). Whether or not a unit of area is added is determined by whether or not this lattice cell is inside or outside the polygon. But we already know from the arguments given above that a lattice cell is inside the polygon if the number of loop ends below the cell is odd (note that any loop end along the vertical edge cutting the horizontal part of the kink is included in the count). So in the case shown in figure 3 the lattice cell to which the kink is moved lies outside the polygon (there are two loop ends below the kink) and no unit of area is added. Note that a unit of area may be added for any given output configuration. In this case the area generating function is obtained correctly to $`2W_{\mathrm{max}}+3k`$. The factor $`3k`$ arise since it takes 3 lattice cells to completely surround the simplest puncture, which is just a single cell.

As far as we are aware, series for punctured polygons have not been derived previously. In such circumstances it is even more important than usual to undertake careful tests to ensure that the series are correct. To this end, a second algorithm was implemented, by one of the authors, to independently evaluate the series coefficients. The complete agreement we obtained between the two data sets reassures us as to the correctness of our results.

### 2.2 Enumeration of punctured staircase polygons

The enumeration of punctured staircase polygons is much simpler. In fact as we shall demonstrate it is a problem for which the computational complexity grows only as a polynomial in the number of terms. As stated in the Introduction we can think of staircase polygons as consisting of two non-intersecting directed walks on the square lattice. Punctured staircase polygons naturally will have more than two walks, in fact up to $`2k+2`$ walkers can be present in any given column. Due to the restrictions on staircase polygons it follows that a punctured staircase polygon is formed by requiring that the walks be directed and that any two walkers starting at one point join each other later without intersecting other walks. This can be encoded in a transfer matrix calculation as follows: Again we count the number of punctured polygons in rectangles as before and draw a line through the lattice intersecting the $`W+2`$ edges. In this case we need only specify whether or not each edge is part of the polygon or not, so each configuration can be represented by an ordered set of edge states $`\{n_i\}`$, where

$$n_i=\{\begin{array}{cc}\hfill 0& \text{empty edge},\hfill \\ \hfill 1& \text{part of loop closed to the left}.\hfill \end{array}$$

This uniquely specifies the configuration because the first occupied edge is connected to the last occupied edge and any edges in between are paired, e.g., the second and third (fourth and fifth and so on) edges form a loop to the left and has to be connected to each other later on. The rules for updating the configurations are as follows: From the local input state ‘00’ we can always get the output ‘00’, and the output ‘11’ provided there is an odd number of edges below the kink and also that the edge directly below the kink is empty. The rules for the output ‘11’ ensure that the new pair of walkers lie within the enclosing staircase polygon but not inside an internal staircase polygon (thus preventing holes within holes), and that the lower walk does not intersect other walks. From the local input state ‘01’ and ‘10’ we can always produce the output ‘01’, and the output ‘10’ provided either that the edge directly below the kink is empty or there is an even number of edges below the kink. These rules ensure that the walkers do not intersect other walks except when we close a valid loop. Finally from the local input state ‘11’ we always get the output state ‘00’. The weights associated with these updates are obtained in the same way as for the SAP enumeration, whether the enumeration is done by perimeter or area.

As in the previous case we calculate the number of punctured staircase polygons spanning the rectangles. Adding the contributions from all rectangles of width $`WW_{\mathrm{max}}`$ and length $`WL2W_{\mathrm{max}}W+1`$ the number of punctured staircase polygons is obtained correctly up to perimeter $`4W_{\mathrm{max}}+2+4k`$. Note that for fixed $`k`$ the maximum number of configurations $`N_C`$ grows as a polynomial in $`W_{\mathrm{max}}`$

$$N_C=\underset{j=0}{\overset{k}{}}\left(\begin{array}{cc}W_{\mathrm{max}}+1& \\ 2j+2& \end{array}\right).$$

So in this case the algorithm is of polynomial complexity.

## 3 Expected behaviour

As mentioned in the introduction, the problem of enumerating SAP by area has been extensively studied . It has been shown that $`\kappa ^{(k)}=\kappa ^{(0)},`$ and that if the exponents exist, $`\beta _k=\beta _0+k`$ . That is to say, for $`k`$ finite, the connective constant for $`k`$-punctured discs, by area, is the same as that for unpunctured discs, while the critical exponent increases by 1 for each puncture. We repeat that the punctures are disjoint.

As far as we are aware, there has been no previous work on the problem of the perimeter generating function of punctured SAP. We first give a proof that a 1-punctured SAP on the square lattice has the same connective constant as an unpunctured SAP, and indicate how that proof can be generalised to $`k`$-punctured SAP.

Before stating and proving the relevant theorem, we require certain preliminary results. For an (unpunctured) SAP of perimeter $`2n`$ one knows,  eqn. (7.101), that

$$\mathrm{exp}[b\sqrt{n}]\mu ^{2n}/4np_{2n}^{(0)}\mu ^{2n},$$
(3)
where $`b`$ is a constant. Further, it is obvious that a polygon of perimeter $`2n`$ has maximum area $`n^2/4,`$ which occurs when the shape is an $`n/2\times n/2`$ square. Thus the area of a polygon of perimeter $`2n,`$ denoted $`𝒜_{2n},`$ satisfies $`𝒜_{2n}n^2/4.`$

Now consider 1-punctured polygons of total perimeter $`2n,`$ with inner perimeter $`2m`$ and outer perimeter $`2(nm)`$. It is possible for the inner polygon to be Hamiltonian, in which case its perimeter is greater than that of the surrounding polygon. Indeed, a square polygon of side $`2n+1,`$ and hence of perimeter $`8n+4`$ can contain an internal polygon of perimeter as large as $`4n^2.`$ Hence the semi-perimeter of the inner polygon, $`m,`$ can range from a minimum value of 2 to a maximum value of $`n2\sqrt{2n}+2<n\sqrt{n}`$ for $`n>1.`$

With these preliminaries, we can now state and prove the following theorem:
Theorem: $`lim_n\mathrm{}\frac{1}{2n}\mathrm{log}p_{2n}^{(1)}=\mathrm{log}\mu ,`$ where $`\mu `$ is the same constant as appears in the corresponding limit for unpunctured SAP.
Proof: 1-punctured polygons of total perimeter $`2n`$ are constructed by placing polygons $`P`$ of perimeter $`2m`$ inside polygons $`Q`$ of perimeter $`2n2m.`$ Let $`w(P,Q)`$ denote the number of ways of placing polygon $`P`$ inside polygon $`Q,`$ (which of course depends on both $`P`$ and $`Q.`$) Then

$$p_{2n}^{(1)}=\underset{m=2}{\overset{n\sqrt{n}}{}}\underset{P}{}\underset{Q}{}w(P,Q).$$

We bound this summand by the product of three factors. The first two factors are the number of polygons of perimeter $`2m`$ and $`2n2m`$ respectively. The third is the number of ways the smaller polygon can be placed inside the surrounding polygon and is clearly less than or equal to the area of the surrounding polygon. Explicitly,

$$p_{2n}^{(1)}\underset{m=2}{\overset{n\sqrt{n}}{}}(nm)^2\mu ^{2m}\mu ^{2(nm)}/4\mu ^{2n}/4\underset{m=2}{\overset{n}{}}(nm)^2n^3\mu ^{2n}.$$

To obtain a lower bound, consider a $`3\times 3`$ square polygon with a unit square hole at its centre. This unique realisation of $`p_{16}^{(1)}`$ can be uniquely concatenated with each unpunctured polygon by joining them at a solitary specified edge, and then deleting that edge. For each unpunctured polygon we take the set of left-most vertical edges and choose the bottom edge from this set, and we choose the right-most, top-most vertical edge of $`p_{16}^{(1)}`$ (a similar operation is shown for polyominoes in figure A1). The concatenation operation then gives, for each of the $`p_{2n14}`$ unpunctured polygons, a unique member of the set of 1-punctured polygons in which the puncture is a single cell. Thus we have $`p_{2n14}^{(0)}p_{2n}^{(1)}.`$ From the above equations we thus obtain

$$\mathrm{exp}[b\sqrt{n7}]\mu ^{2n14}/4(n7)p_{2n14}^{(0)}p_{2n}^{(1)}n^3\mu ^{2n}$$

for $`n>8.`$ The theorem then follows immediately on taking logarithms, dividing through by $`2n`$ and taking the limit as $`n\mathrm{}.`$

This proof can clearly be extended to two-punctured polygons, then to three-punctured polygons etc., by concatenating unpunctured polygons with minimal two-punctured, three-punctured etc. polygons. In the appendix we prove that $`k`$-punctured polyominoes have the same growth constant as unpunctured SAP by area.

We have been unable to prove a result analogous to the result for the exponent of punctured SAP by area, but give an argument that depends on certain assumptions that are generally accepted, though not proved. In the case of staircase polygons however our assumptions have been proved, and so our result will be rigorously true.

The key results we need are those obtained in  to the effect that the mean area of SAP of perimeter $`2n`$ is proportional to $`n^{1.5}.`$ This is true both for SAP and for staircase polygons, and in the latter case it has been proved. More precisely, we need the following result. There exists constants $`D_1`$ and $`D_2`$ such that the mean area $`\overline{A}_{2n}`$ of polygons of perimeter $`2n`$ satisfies

$$D_1n^{\frac{3}{2}}\overline{A}_{2n}=\frac{_QA_Q}{p_{2n}}D_2n^{\frac{3}{2}},$$

where the sum is taken over all $`p_{2n}`$ polygons $`Q`$ of perimeter $`2n.`$

Further, there exists constants $`C_1`$ and $`C_2`$ such that the number of polygons of perimeter $`2n`$ satisfies

$$C_1\mu ^{2n}n^{\frac{5}{2}}p_{2n}^{(0)}C_2\mu ^{2n}n^{\frac{5}{2}}.$$

For staircase polygons, $`\mu =2`$ and the exponent is $`\frac{3}{2}`$ instead of $`\frac{5}{2}`$ in the above equation.

As above, let $`Q`$ denote a polygon of perimeter $`2n2m,`$ with area $`A_Q.`$ The number of ways of placing a given polygon $`P`$ of perimeter $`2m`$ inside $`Q`$ is clearly less than $`A_Q.`$

Thus the number of ways of placing $`P`$ inside any polygon of perimeter $`2n2m`$ is less than

$$\underset{Q}{}A_Q=p_{2n2m}^{(0)}\overline{A}_{2n2m},$$
(4)
where the sum is over all polygons $`Q`$ of perimeter $`2n2m.`$

Hence the number of ways of placing all of the $`p_{2m}^{(0)}`$ polygons of perimeter $`2m`$ inside any polygon of perimeter $`2n2m`$ is less than $`p_{2m}^{(0)}p_{2n2m}^{(0)}\overline{A}_{2n2m},`$ and so

$$p_{2n}^{(1)}<\underset{m=2}{\overset{n\sqrt{n}}{}}p_{2m}^{(0)}p_{2n2m}^{(0)}\overline{A}_{2n2m}<D_2C_2^2\mu ^{2n}2^5S_{2n},$$
(5)
where

$$S_{2n}=\underset{m=2}{\overset{n\sqrt{n}}{}}m^{\frac{5}{2}}/(nm).$$

The last sum may be evaluated in a variety of ways. Using Maple we find that

$$S_{2n}=0.341../n+O(n^{\frac{3}{2}}).$$

Thus we obtain the bound

$$p_{2n}^{(1)}<const.\mu ^{2n}/n.$$
(6)
For staircase polygons the analogous calculation is slightly simpler, and we find

$$S_{2n}1+\zeta (3/2)+O(n^{\frac{1}{2}}),$$

where $`\zeta (z)`$ is the Riemann zeta function.

To obtain a lower bound, we restrict the inner polygon to be of perimeter 4, that is, a unit square. The number of polygons punctured by a unit square clearly provides a lower bound to the number of 1-punctured polygons. A unit square can be placed anywhere in a polygon $`Q`$ of perimeter $`2n4`$ except on a boundary site. There are $`2n4`$ boundary sites. The mean area of $`Q`$ is greater than $`D_1(n2)^{\frac{3}{2}},`$ so we obtain the bound

$$p_{2n}^{(1)}>p_{2n4}(D_1(n2)^{\frac{3}{2}}2n+4)>C_1\frac{\mu ^{2n4}}{(2n4)^{\frac{5}{2}}}(D_1(n2)^{\frac{3}{2}}2n)>const.\mu ^{2n}/n.$$
(7)
Combining the two bounds gives the result

$$E_1\mu ^{2n}/np_{2n}^{(1)}E_2\mu ^{2n}/n,$$

where $`E_1<E_2`$ are constants. Accepting the usual asymptotic form that is expected for such models, we conclude that

$$p_{2n}^{(1)}const.\mu ^{2n}/n.$$
(8)
(For staircase polygons the analogous result is $`p_{2n}^{(1)}const.4^n.`$) Since the number of unpunctured polygons grows like $`\mu ^nn^{5/2},`$ we see that the exponent is predicted to increase by $`\frac{3}{2}`$ as the result of a single puncture, while in the case of the area generating function, the exponent is found to increase only by 1. We show in the next section that this prediction is borne out by our numerical calculations.

For punctured staircase polygons, a similar conclusion holds. That is, the connective constant is unaltered at $`\mu _{stair}=2,`$ but the exponent increases by $`1.5`$ over its unpunctured counterpart when enumerating 1-punctured staircase polygons by perimeter. In fact, we have been able to calculate the generating function for staircase polygons with a single puncture of perimeter 4, and also with a single puncture of perimeter 6. The generating functions for these special cases are given below, and are precisely in accordance with the more general results given above.

We also note that Cardy  recently considered the problem of the number of punctured SAP with $`k`$ concentric, mutually self-avoiding SAPs surrounding a fixed point of the dual lattice. Thus for $`k=1`$ this corresponds to 1-punctured SAPs surrounding a fixed point. In that case Cardy finds  for the number of such configurations $`b_{2n}^{(1)}`$ that $`b_{2n}^{(1)}=\mu ^{2n}\frac{\mathrm{ln}n}{64\pi ^2n}.`$ The above calculation may be repeated for Cardy’s problem. All that is required is to add a factor $`\overline{A}_m`$ to the summand in eqn. 5, as the surrounded point can be anywhere inside the inner polygon. Thus we must multiply by the mean area of that polygon. The sum defining $`S_{2n}`$ then becomes $`S_{2n}=_{m=2}^{n\sqrt{n}}\frac{1}{m(nm)}\mathrm{log}n/n,`$ in agreement with Cardy’s result. (Our lower bound for this problem is too weak, but with more effort could be improved.)

## 4 Analysis of the series

All the series we have investigated are characterised by coefficients that grow exponentially, with sub-dominant term given by a critical exponent. The generic generating function behaviour is $`G(z)=_ng_nz^nA(z)(1\sigma z)^\xi ,`$ and hence the coefficients of the generating function $`g_n=[z^n]G(z)A(1/\sigma )/\mathrm{\Gamma }(\xi )\sigma ^nn^{\xi 1}.`$ Generally speaking the existence of the growth constant $`\sigma `$ has been proved, but except for exactly solvable models, such as staircase polygons, the existence of the critical exponent $`\xi `$ has only been conjectured, though its existence has never been doubted. The radius of convergence of the generating function is usually given by the critical point, which is at $`z=1/\sigma .`$

We principally used two methods to analyse all the series studied in this paper. Firstly, to obtain the singularity structure of the generating function we used the numerical method of differential approximants . In particular, we used this method to estimate the growth constant $`\sigma `$ and the critical exponent $`\xi .`$ We were invariably able to conjecture an exact value for $`\xi ,`$ which was always an integer or half-integer for all the problems we investigated. Imposing this conjectured exponent permitted a refinement of the estimate of the growth constant — providing so-called biased estimates.

Once the exact value of the exponent was conjectured, and the growth constant accurately estimated, we turned our attention to the “fine structure” of the asymptotic form of the coefficients, by fitting the coefficients to the assumed form $`g_n=[z^n]G(z)\sigma ^nn^{\xi 1}_{i0}c_i/n^{f(i)}.`$ If there is no non-analytic correction term, then $`f(i)=i,`$ while a square-root correction term means $`f(i)=i/2.`$ For all the series studied, only these two situations were encountered.

In all cases, our procedure is to assume a particular form for $`f(i),`$ and see how it fits the data. With the very long series we now have at our disposal, it is usually easy to see if the wrong assumption has been made — the sequence of amplitude estimates $`c_i`$ either diverges to infinity or converges to zero. Once the correct assumption is made, convergence is usually rapid and obvious. A detailed demonstration of the method can be found in .

As an example of the sort of results we obtained, we show sequences of estimates of the coefficients of the perimeter generating function of unpunctured SAP in Table 2. In that case we conjectured that $`f(i)=i.`$ Because of the large amount of tabular data generated by the method, we have not given this level of detail for the many series investigated here. We show only the results for two series. For the others, we just give our assessment of the apparent convergence of the sequences $`c_i,`$ and the estimated value of the limits. As the equations involved are linear, the method is easy to implement, and interested readers can readily generate the relevant data themselves. Some subtleties nevertheless exist. For example, for punctured staircase polygons, the perimeter generating function has two singularities on the circle of convergence, and so both must be taken into account. We discuss this in more detail in the relevant section below.

As for the first stage of the analysis, the method of differential approximants, we proceeded as follows: Estimates of the critical point and critical exponent were obtained by averaging values obtained from first order $`[L/N;M]`$ and second order $`[L/N;M;K]`$ inhomogeneous differential approximants. For each order $`L`$ of the inhomogeneous polynomial we averaged over those approximants to the series which used at least the first 80% - 90% of the terms of the series, and used approximants such that the difference between $`N`$, $`M`$, and $`K`$ didn’t exceed 2. These are therefore “diagonal” approximants. Some approximants were excluded from the averages because the estimates were obviously spurious. The error quoted for these estimates reflects the spread (basically one standard deviation) among the approximants. Note that these error bounds should not be viewed as a measure of the true error as they cannot include possible systematic sources of error. However systematic error can also be taken into account in favourable situations, as for example in the case of SAP enumerated by perimeter . Again, in the interests of space, we present only our results, and not the intermediate detail from which our estimates were made. An example in full detail for one of the series investigated in this study can be found in . We turn now to the analysis of all the series.

### 4.1 Staircase polygons

For (unpunctured) staircase polygons, the multi-variable width, height and area generating function is known . As usual, we denote

$$(a)_n=\underset{i=0}{\overset{n1}{}}(1aq^i).$$

Further, denoting the first two $`q`$-Bessel functions as:

$$J_0(x,y,q)=\underset{n0}{}\frac{(1)^nx^nq^{\left(\genfrac{}{}{0pt}{}{n+1}{2}\right)}}{(q)_n(yq)_n},$$

and

$$J_1(x,y,q)=\underset{n1}{}\frac{(1)^{n1}x^nq^{\left(\genfrac{}{}{0pt}{}{n+1}{2}\right)}}{(q)_{n1}(yq)_{n1}(1yq^n)},$$

the perimeter and area generating function is simply

$$P(x,y,q)=y\frac{J_1(x,y,q)}{J_0(x,y,q)},$$

where $`x`$ $`(y)`$ is the variable conjugate to the horizontal (vertical) semi-perimeter, while $`q`$ counts the area. No analogous results are known for punctured staircase polygons, though the calculation for SAP  by area can be carried over mutatis mutandis to prove that the area generating function for a punctured staircase polygon with $`k`$ holes has the same radius of convergence as the area generating function for unpunctured staircase polygons. Further, the critical exponent increases by 1 for each puncture, and unlike the case for SAP, we can not only prove the existence of a critical exponent for unpunctured staircase polygons, but we know its value. Hence we know the leading term in the asymptotic expansion of the generating function by area for $`k`$-punctured staircase polygons.

For the expected behaviour of the perimeter generating function of punctured staircase polygons, the arguments of the preceding section apply directly. The radius of convergence, and hence the connective constant remains unchanged, and the argument given in the preceding section suggests that the critical exponent should increase by 1.5 for each puncture. The results of our analysis, presented below, bear this out.

#### 4.1.1 Area generating function

For unpunctured staircase polygons, the area generating function is given by

$$A(q)=\underset{n1}{}a_n^{(0)}q^n=\frac{J_1(1,1,q)}{J_0(1,1,q)}.$$

By inspection, this has poles at the zeros of $`J_0(1,1,q).`$ The nearest zero is at  $`1/q=\eta =2.30913859330,`$ and there is a simple pole at that point. The next zero is well separated (at $`1/q=\lambda =1.4435..`$) and so the asymptotic form of the generating function is

$$A(q)D/(1\eta q)+E/(1\lambda q)+\mathrm{},$$

and hence

$$a_n^{(0)}=[q^n]A(q)\eta ^n(c_0+O((\lambda /\eta ))^n).$$

Our analysis bears this out, and we estimate $`c_0=0.12881579.`$

For 1-punctured discs, our analysis, based on more than 100 series coefficients, convincingly suggests the following asymptotic form:

$$A^{(1)}(q)D^{(1)}/(1\eta q)^2+E^{(1)}/(1\eta q)^{1.5}+F^{(1)}/(1\eta q)+\mathrm{},$$

and hence

$$a_n^{(1)}=[q^n]A^{(1)}(q)\eta ^nn\underset{i0}{}c_i/n^{i/2}.$$

The sequences of amplitude estimates, assuming this asymptotic form, are shown in Table 3. The apparent convergence of the amplitude estimates is, as explained above, our source of evidence for this asymptotic form.

For 2-punctured discs, a similar analysis, based on some 86 series coefficients, convincingly suggests the following asymptotic form:

$$A^{(2)}(q)D^{(2)}/(1\eta q)^3+E^{(2)}/(1\eta q)^{2.5}+\mathrm{},$$

and hence

$$a_n^{(2)}=[q^n]A^{(2)}(q)\eta ^nn^2\underset{i0}{}a_i/n^{i/2}.$$

For 3-punctured discs our analysis was based on an 89 term series. We found that the above pattern persists, so that the generating function has the following asymptotic form:

$$A^{(3)}(q)D^{(3)}/(1\eta q)^4+E^{(3)}/(1\eta q)^{3.5}+\mathrm{},$$

and hence

$$a_n^{(3)}=[q^n]A^{(3)}(q)\eta ^nn^3\underset{i0}{}a_i/n^{i/2}.$$

Estimates of the various amplitudes defined above are shown in Table 4.

#### 4.1.2 Perimeter generating function

For unpunctured staircase polygons, the perimeter generating function (ignoring the distinction between height and width) is given by

$$P(x)=\frac{12x\sqrt{14x}}{2},$$

which is, apart from suppression of the first term, the generating function for Catalan numbers. Hence

$$p_{2n}^{(0)}=[x^n]P(x)=\frac{1}{n}\left(\genfrac{}{}{0pt}{}{2n2}{n1}\right)4^n/n^{\frac{3}{2}}\underset{i0}{}c_i/n^i.$$

The values of $`c_i`$ follow immediately from the exact solution. They are given in Table 4.

For 1-punctured discs, our analysis, again based on more than 100 series coefficients, was more equivocal than that of the 1-punctured area generating function. The method of differential approximants clearly identified a singularity at the known critical point, $`x_c=1/4,`$ but almost all approximants had a double root, implying a confluent singularity. The leading exponent was estimated to be $`1,`$ implying a pole in the generating function, (and hence immediately lending support to our conjectured change in the critical exponent of 3/2 as a result of puncturing), but we were unable, from this method, to identify the confluent exponent. Further, a second singularity was identified, of the form const.$`(1+x/x_c)^{6.5}.`$ At this point we wish to remark on the close similarity between a recently solved model of polygons, the three-choice polygon model  and 1-punctured staircase polygons. 1-punctured staircase polygons can be thought of as being constructed from two three-choice polygons, with common edges deleted. This geometric similarity is borne out by the fact that the two models have identical singularity distributions, with even the exponents being the same at all non-physical singularities.

For 1-punctured discs then, the perimeter generating function is expected to be of the following asymptotic form:

$$P^{(1)}(x)=\underset{n}{}p_{2n}^{(1)}x^nB^{(1)}(x)+C^{(1)}(x)(14x)^1+D^{(1)}(x)(1+4x)^{6.5},$$

and hence

$$p_n^{(1)}=[x^n]P^{(1)}(x)4^n\underset{i0}{}c_i/n^{f(i)}+(4)^nn^{7.5}\underset{i0}{}d_i/n^i,$$

where $`B^{(1)}(x),`$ $`C^{(1)}(x)`$ and $`D^{(1)}(x)`$ are assumed regular in the disc $`|4x|1.`$

Assuming $`f(i)=i,`$ which implies only analytic correction-to-scaling terms, gave unsatisfactory results. Notably, we observed that the estimates of the amplitude $`c_1`$ in the above asymptotic form were steadily increasing, suggesting that the assumed form did not properly account for the correction-to-scaling terms. With $`f(i)=i/2`$ the amplitude estimates were much more stable. This then implies that the generating function in fact behaves as

$$P^{(1)}(x)B^{(1)}(x)+C^{(1)}(x)(14x)^1+E^{(1)}(x)(14x)^{\frac{1}{2}}+D^{(1)}(x)(1+4x)^{6.5}.$$
(9)
The amplitude estimates, assuming this asymptotic form, are shown in Table 4. The apparent convergence of the amplitude estimates is, as explained above, our source of evidence for this asymptotic form. We did not tabulate the non-physical amplitudes, as they are of little interest to our investigation. However they need to be included to stabilise estimates of the physical amplitudes. We mention in passing that $`d_00.14.`$

We also generated series for staircase polygons with two and three punctures, with perimeter 150 and 134 steps respectively. Our differential approximant analysis lent support to our expectation that the critical exponent increases by $`1.5`$ per puncture. This turned out to be true for the non-physical singularity at $`x=\frac{1}{4}`$ as well. A similar analysis to that for the one-punctured polygons strongly supported the analogous asymptotic forms,

$$P^{(k)}(x)=\underset{n}{}p_{2n}^{(k)}x^nB^{(k)}(x)+C^{(k)}(x)(14x)^{0.51.5k}+D^{(1)}(x)(1+4x)^{81.5k}$$
(10)
for $`k>0.`$ Hence the asymptotic form of the coefficients is conjectured to be

$$p_{2n}^{(k)}4^nn^{\frac{3(k1)}{2}}\underset{i0}{}c_i^{(k)}/n^{i/2}+(4)^nn^{1.5k9}\underset{i0}{}d_i^{(k)}/n^i,$$

where the amplitude estimates $`c_i^{(k)}`$ are given in Table 4.

We have also determined the exact generating function for staircase polygons punctured by a single hole of perimeter 4, and also the generating function for staircase polygons punctured by a single hole of perimeter 6. We obtained these generating functions by generating the coefficients using the algorithm discussed, and then searching for an underlying differential equation. As a result we find the perimeter generating functions $`P_4^s(x)`$ and $`P_6^s(x)`$ for 1-punctured staircase polygons with a hole of perimeter 4 and 6 respectively.

$$P_4^s(x)=\frac{2x^416x^3+20x^28x+1}{2(14x)}\frac{16x+10x^24x^3}{2\sqrt{14x}}$$
(11)
and

$`P_6^s(x)`$ $`=`$ $`{\displaystyle \frac{126x+228x^2906x^3+1709x^41378x^5+322x^6}{2(14x)^{\frac{5}{2}}}}`$ (12)
$`{\displaystyle \frac{32x^6404x^5+815x^4586x^3+182x^224x+1}{2(14x)^2}}.`$
Note that both these exact solutions display the confluent square root correction that we have found in our numerical investigations in the more general case. In the next subsection we analyse the analogous generating function for SAP.

### 4.2 Self-avoiding polygons

For unpunctured SAP, the perimeter generating function was recently extended  to 90 step polygons, and the asymptotics clearly identified. The polygon generating function is defined to be

$$P^{(0)}(x)=\underset{n}{}p_{2n}^{(0)}x^nA^{(0)}(x)+B^{(0)}(x)(1\mu ^2x)^{\frac{3}{2}},$$
(13)
where the functions $`A^{(0)}`$ and $`B^{(0)}`$ are believed to be regular in the vicinity of $`x_c=1/\mu ^2.`$ We estimated $`\mu =2.63815853034(10)`$. From this equation follows the asymptotic form of the coefficients,

$$p_{2n}^{(0)}\mu ^{2n}n^{\frac{5}{2}}[c_1+c_2/n+c_3/n^2+c_4/n^3+\mathrm{}].$$
(14)
We show in Table 2 the sequence of estimates of $`c_i,`$ and in Table 4 our estimates of the amplitudes, being the limits of the sequences $`\{c_i\}.`$

The area generating function was first studied in , where the first 20 terms were given, and the asymptotic form estimated to be

$$A^{(0)}(q)=\underset{n}{}a_nq^nD(q)+E(q)\mathrm{log}(1\kappa q),$$
(15)
where $`\kappa 3.97087,`$ and the logarithm in the above equation was understood to include the possibility of some power of a logarithm other than unity. (Though our analysis below implies that this is not the case.) In the present work we extend the series to 42 terms.

#### 4.2.1 Area generating function

Using our greatly extended 42 term series, our analysis of the unpunctured disc area generating function was carried out by standard methods. We used differential approximants  and found unbiased critical point and critical exponent estimates. The unbiased exponent estimate had absolute value less than $`10^5,`$ totally supporting our view that it is exactly zero. Assuming this, a biased estimate of the critical point is possible, and in this way we estimate $`\kappa =3.97094397(9).`$

We then proceeded to seek the asymptotic form of the coefficients by writing

$$a_n=[q^n]A^{(0)}(q)\kappa ^n/n\underset{i0}{}c_i/n^{f(i)}.$$

Our numerical results were well converged, demonstrating very convincingly that $`f(i)=i.`$ This is the asymptotic form consistent with a pure logarithmic singularity, not raised to any power. Estimates of the amplitudes are given in Table 4.

For 1-punctured discs, our analysis, based on the series known to $`q^{45}`$, convincingly suggests the following asymptotic form:

$$A^{(1)}(q)D^{(1)}(q)+E^{(1)}(q)/(1\kappa q),$$

and hence

$$a_n^{(1)}=[q^n]A^{(1)}(q)\kappa ^n\underset{i0}{}c_i/n^i.$$

Attempts to fit to alternative forms, corresponding to a confluent logarithm or a non-analytic correction-to-scaling term were unsuccessful, adding to our confidence that the above form is correct. Estimates of the amplitudes are given in Table 4.

For twice-punctured discs, our analysis, also based on the series known to $`q^{48}`$ convincingly suggests the following asymptotic form:

$$A^{(2)}(q)D^{(2)}(q)+E^{(2)}(q)/(1\kappa q)^2,$$

and hence

$$a_n^{(2)}=[q^n]A^{(2)}(q)\kappa ^nn\underset{i0}{}c_i/n^i.$$

The conjectured asymptotic form for $`k`$ punctured polygons, by area, is thus

$$A^{(k)}(q)D^{(k)}(q)+E^{(k)}(q)/(1\kappa q)^k.$$

#### 4.2.2 Perimeter generating function

The results for unpunctured polygons are fully discussed at the beginning of subsection 4.2. As we found with staircase polygons, the generating function for punctured discs by perimeter was a more challenging numerical analysis problem than either its unpunctured counterpart, or its area counterpart.

We found that the method of differential approximants was not particularly satisfactory. Given that we needed some 100 terms to successfully analyse the (presumably simpler) problem of punctured staircase polygons by perimeter, it is not surprising that for punctured SAP, for which we have 33 non-zero coefficients (corresponding to perimeters up to 84 steps), the method was not satisfactory. However, it did indicate the presence of a confluent singularity. As we found a confluent square-root singularity for staircase polygons, it is hardly surprising that a confluent singularity is detected for the punctured SAP generating function. In fact an exponent shift of around 1.7 was seen, compared to the expected value 1.5. We attribute this to the “short” series, coupled with the well known deleterious effect of confluent terms in such an analysis. Nevertheless, subsequent analysis of the asymptotic form of the coefficients, assuming an exponent shift of 1.5, together with a square-root confluent term, as found for punctured staircase polygons, gave satisfactory results.

We denote the generating function for $`k`$-punctured SAP, by perimeter, as

$$P^{(k)}(x)=\underset{n}{}p_{2n}^{(k)}x^nB^{(k)}(x)+C^{(k)}(x)(1\mu ^2x)^{1.51.5k},$$
(16)
where the exponent is conjectured. For 1-punctured SAP, the vanishing of the exponent implies a logarithmic singularity. We fitted the coefficients to the asymptotic form appropriate to $`\mathrm{log}(1\mu ^2x),`$ so that the asymptotic form of the coefficients just involves decreasing integer powers of $`n.`$ We then found that the estimates of the leading amplitude were monotonically increasing, which implies that the asymptotic form is wrong — too weak. Including a confluent square root singularity, as was found for punctured staircase polygons, stabilised the estimates. Accordingly, we conjecture that the asymptotic form is dominated by a logarithmic singularity, with a sub-dominant square root singularity, so that

$$p_{2n}^{(1)}=[x^n]P^{(1)}(x)\mu ^{2n}/n\underset{i0}{}c_i/n^{\frac{i}{2}}$$
(17)
Estimates of the amplitudes $`c_i`$ are given in Table 4.

For twice punctured discs, a similar analysis suggested that the asymptotic form of the generating function is

$$P^{(2)}(x)=\underset{n}{}p_{2n}^{(2)}x^nB^{(2)}(x)+C^{(2)}(x)(1\mu ^2x)^{1.5},$$
(18)
again with evidence of a square root confluent term. As for 1-punctured SAP, we give estimates of the amplitudes $`c_i`$ defined by

$$p_{2n}^{(2)}\mu ^{2n}n^{\frac{1}{2}}\underset{i0}{}c_i/n^{\frac{i}{2}}$$

in Table 4.

For thrice punctured discs, a similar analysis suggested that the asymptotic form of the generating function is

$$P^{(3)}(x)=\underset{n}{}p_{2n}^{(3)}x^nB^{(3)}(x)+C^{(3)}(x)(1\mu ^2x)^3,$$
(19)
again with evidence of a square root confluent term. As for 1-punctured SAP, we give estimates of the amplitudes $`c_i`$ defined by

$$p_{2n}^{(3)}\mu ^{2n}n^2\underset{i0}{}c_i/n^{\frac{i}{2}}$$

in Table 4.

The conjectured asymptotic form for $`k`$ punctured polygons, by perimeter, is then

$$P^{(k)}(x)=\underset{n}{}p_{2n}^{(k)}x^nB^{(k)}(x)+C^{(k)}(x)(1\mu ^2x)^{1.51.5k}$$
(20)
where for $`k>0`$ we find strong evidence for a square root correction-to-scaling term.

### 4.3 Polyominoes

The problem of polyominoes has a long and interesting history, and has been well discussed in the popular scientific literature .

The enumeration of square lattice polyominoes to 24 steps  was given in 1981, extended to 25 steps in 1995  and currently stands at 28 steps .

We have analysed the latest series by the method of differential approximants, and find the generating function behaves as

$$𝒫(y)=\underset{n}{}a_ny^nG(y)+H(y)\mathrm{log}(1\tau y),$$
(21)
where $`a_n`$ is the number of polyominoes of area $`n`$. In  the estimate $`\tau =4.06265(5)`$ was given. The extra terms now available allow us to make the refined estimate $`\tau =4.062591(9).`$ Analysis of the asymptotic form of the coefficients is totally consistent with a simple logarithm in the generating function. Thus

$$a_n=[y^n]𝒫(y)\tau ^n\underset{i0}{}c_i/n^{i+1}.$$

The amplitude estimates are given in Table 4. The leading amplitude is in complete agreement with, but 3 orders of magnitude more accurate than that given in , while the order of the leading term — $`O(1/n)`$ — was predicted by physical arguments back in 1981 . In that work, the logarithmic singularity in the generating function of strongly embedded site-animals was obtained. As discussed in the introduction, these are isomorphic to polyominoes.

As well as extending the polyomino series, Oliveira e Silva has enumerated $`k`$-punctured polyominoes  for $`k6.`$ Clearly, $`k=0`$ polyominoes are just SAP, and as we have seen in the previous section, these grow as $`\kappa ^n`$ where $`\kappa =3.9709..<\tau .`$ The arguments in  can be modified and applied to show that a $`k`$-punctured polyomino has the same growth constant $`\lambda `$ as its unpunctured counterpart, and this is done in the Appendix.

Thus the situation is that, for any finite number of punctures,

$$a_n^{(k)}=[y^n]𝒫^{(k)}(y)\kappa ^n,$$

but that

$$a_n=[y^n]𝒫(y)=\underset{k0}{}a_n^{(k)}\tau ^n.$$

We have analysed the series for $`k`$-punctured polyominoes, $`k=0,1,2,`$ and find the asymptotic form of the coefficients to be

$$a_n^{(k)}=[y^n]𝒫^{(k)}(y)\kappa ^nn^{k1}\underset{i0}{}c_i^{(k)}/n^i,$$

corresponding to a generating function for $`k`$-punctured polyominoes having a $`k^{th}`$ order pole, viz:

$$𝒫^{(k)}(y)=\underset{n}{}a_n^{(k)}y^nH^{(k)}(y)(1\kappa y)^k,$$
(22)
where $`k=0`$ is to be interpreted as a logarithm. Thus just as for punctured polygons, it is found that the exponent of $`k`$-punctured polyominoes increases by 1 for each puncture. This is also proved in the Appendix.

The coefficients $`c_i^{(k)}`$ are given in Table 4.

## 5 Conclusion

We have investigated the effect of punctures on SAP and staircase polygons enumerated both by area and perimeter.

In order to do this, we have developed a new algorithm, exponentially faster than direct counting, whereby we have radically extended a number of series. This extension was necessary in order to probe some rather subtle numerical behaviour.

We found that, in every case, a finite number of punctures does not change the exponential growth factor associated with the unpunctured counterpart of the punctured object being enumerated.

This latter conclusion was also proved for polyominoes. The effect of punctures was also investigated numerically for polyominoes.

Writing $`b_n^{(k)}`$ for the $`n^{th}`$ coefficient in the generating function for some $`k`$punctured object, so that

$$b_n^{(k)}(\omega ^{(k)})^nn^{\theta (k)},$$

we found $`\omega ^{(k)}=\omega ^{(0)}`$ in all cases. We found, further, that $`\theta (k)=\theta (0)+k`$ if enumerating any of the objects we have considered by area. This can be proved, though subject, in some cases, to the existence of the exponent in question. Subject to certain unproved assumptions we also showed that $`\theta (k)=\theta (0)+3k/2`$ if enumerating by perimeter.

We have, for the first time, obtained good numerical estimates of the sub-dominant terms for a range of problems, thus identifying both the nature of the generating function and any correction-to-scaling terms.

We have also obtained an exact solution for the generating functions of staircase polygons, enumerated by perimeter, punctured by a single hole of perimeter 4 and of perimeter 6. These exact solutions provide additional support for our numerically based conjectures of the correction-to-scaling exponent in the general case.

A more accurate estimate of the growth constant for SAP enumerated by area has been given, complementing our earlier work on the perimeter growth constants .

## E-mail or WWW retrieval of series

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

## Acknowledgements

We have derived great benefit from discussions of aspects of this problem with Mireille Bousquet-Mélou, John Cardy, Aleks Owczarek, Buks van Rensburg, and Stuart Whittington. We are particularly grateful to Mireille Bousquet-Mélou for her careful reading of the manuscript, to Buks van Rensburg for clarifying aspects of his earlier work, and to John Cardy for querying some earlier, incorrect results. IJ and AJG gratefully acknowledge financial support from the Australian Research Council. LHW would like to thank The University of Melbourne for the Melbourne Research Scholarship.

## Appendix. Growth constants and exponents for $`k`$-punctured polyominoes

In this appendix we show that the growth constants are the same for all $`k`$-punctured polyominoes. Our method of proof is based on that of van Rensburg and Whittington , making the necessary changes for the polyomino problem, and discussing in detail certain special cases. That is to say, if $`\kappa `$ denotes the connective constant for SAP enumerated by area, then this is the growth constant for $`k`$-punctured polyominoes for any finite $`k`$. Further, if the usual asymptotic form for the number of $`k`$-punctured polyominoes is assumed, $`s_n^{(k)}C_kn^{\varphi _k}\kappa ^n`$, and $`\varphi _0`$ exists, then $`\varphi _k=\varphi _0k.`$

### A.1 Operations and mappings on punctured polyominoes

Let the set of all $`k`$-punctured polyominoes with $`n`$ cells be denoted by $`\mathrm{\Phi }_n^{(k)}`$, and the set of all polyominoes of $`n`$ cells be $`\mathrm{\Phi }_n`$. Then, for each $`n`$, $`\mathrm{\Phi }_n=_k\mathrm{\Phi }_n^{(k)}`$. Throughout, let $`s_n^{(k)}`$ denote the cardinality of $`\mathrm{\Phi }_n^{(k)}`$ and let $`s_n`$ denote the cardinality of $`\mathrm{\Phi }_n`$. Hence $`s_n=_ks_n^{(k)}`$.

Following  we now define operations on punctured polyominoes, including concatenation, drilling and surgery, which are needed in subsequent proofs. Concatenation allows us to change the size of polyominoes, while drilling and surgery are concerned with changing the number of holes in a polyomino.

#### A.1.1 Concatenation

The concatenation mapping defined here is similar to that in . Consider the bounding rectangle $`R(P)`$ of any polyomino $`P\mathrm{\Phi }_n^{(k)}`$. Define the top (bottom) edge of $`P`$ to be the top (bottom) edge along the east (west) side of $`R(P)`$.

Now, the concatenation of two polyominoes $`P\mathrm{\Phi }_n^{(h)}`$ and $`Q\mathrm{\Phi }_m^{(k)}`$ is defined by joining $`P`$ and $`Q`$ while superimposing the top edge of $`P`$ and the bottom edge of $`Q`$. The result is an $`(h+k)`$-punctured polyomino with $`m+n`$ cells, see figure A1. Hence we have a map

$$T:\mathrm{\Phi }_n^{(h)}\times \mathrm{\Phi }_m^{(k)}\mathrm{\Phi }_{m+n}^{(h+k)}$$
(A1)
###### Lemma A1

For all non-negative values of $`h`$ and $`k`$, $`s_n^{(h)}s_m^{(k)}s_{m+n}^{(h+k)}`$ and in particular, $`s_{n+1}^{(h)}s_n^{(h)}`$

Proof: For the mapping $`T`$ defined above, every pair of polyominoes in the domain can be concatenated to form a larger polyomino in the codomain. Conversely, every such polyomino in the codomain can be uniquely broken up into the original ones. However, there are some polyominoes in the codomain which cannot be formed by concatenating two smaller polyominoes. An example is a $`2\times 2`$ polyomino. Hence we get the first part of lemma as

$$|\mathrm{\Phi }_n^{(h)}\times \mathrm{\Phi }_m^{(k)}||\mathrm{\Phi }_{m+n}^{(h+k)}|.$$

Putting $`k=0`$ and $`m=1`$ and noting that $`|\mathrm{\Phi }_1^{(0)}|=1`$, we get the second part. $`\mathrm{}`$

#### A.1.2 Drilling

Simply creating a hole inside an unpunctured polyomino by removing some interior cells does not allow us to drill certain classes of polyominoes. For instance, a polyomino composed of a single linear sequence of cells cannot be ‘drilled’ since removing any cell will either disconnect the polyomino or shorten the sequence.

The following definition of drilling differs somewhat from that in , though the underlying idea is the same.

First, we drill one hole. Let $`P_0\mathrm{\Phi }_n^{(0)}`$ be the polyomino that we are to puncture. Cover $`P_0`$ by a grid system, with each grid square of size $`b\times b`$ cells for any $`b5`$ (the minimum size will be justified later). Say $`b=5.`$ Pick any grid square $`G`$ that covers at least one cell of $`P_0`$ and drill a hole there, as detailed below.

Step 1 Remove all cells within the selected grid square $`G`$. If, after this step, we are left with a $`1`$-punctured polyomino, we are done. If not, go to step 2.

Step 2 Check the corners of $`G`$: If there are 2 disconnected components touching each other only at the corner, connect them by adding a cell at the appropriate corner of $`G`$.

Step 3 Put a $`1`$-punctured polyomino with 8 cells at the centre of $`G`$.

Step 4 Reconnect the disconnected components outside $`G`$ to the $`1`$-punctured polyomino by adding linear sequences of cells (non-unique). These steps are illustrated in figure A2.

Note that the minimum value of $`b`$ is 5 because if we have anything less than that, we might create extra holes unexpectedly as the $`1`$-punctured polyomino in Step 3 must touch the boundary of the grid square. An example in figure A3 will illustrate this. In this example the original polyomino has no holes and the drilled polyomino has two holes, whereas with a grid of size 5 (or more), the number of holes in the drilled polyomino is only one.

After this operation, the maximum number of cells that could be removed is $`b^2=25`$ (finishing at step 1). The maximum number of cells that could be added is $`b^22`$ (this occurs when there was only 1 cell in the square before drilling and we end up with a 1-punctured polyomino with $`b^21`$ cells after the operation). So, depending on each instance of the drilling operation, we obtain a resulting polyomino $`P_1\mathrm{\Phi }_j^{(1)}`$ where $`j`$ could be anything from $`nb^2`$ to $`n+(b^22)`$. The drilling operation thus defines a map

$$D:\mathrm{\Phi }_n^{(0)}\underset{j=nb^2}{\overset{n+(b^22)}{}}\mathrm{\Phi }_j^{(1)}.$$
(A2)
###### Theorem A1

There exists a real constant $`C`$ such that for all $`b5`$,

$$s_n^{(0)}Cs_{n+(b^22)}^{(1)}.$$
(A3)
Proof: Consider the intermediate (possibly disconnected) polyomino after Step 1, call it $`P^i`$. There are many possible initial polyominoes which gives the same $`P^i`$, i.e., all polyominoes with the same configuration outside $`G`$. They all map to “almost” the same resulting polyomino (the uncertainty implied in “almost” comes from the non-uniqueness of reconnecting in Step 4 which will be discussed later). Therefore, the mapping from domain to codomain is $`M^{}`$-to-one where $`M^{}`$ is bounded by the number of ways that at most $`b^2`$ cells can be connected within $`G`$. Let the bound be $`M`$.

On the other hand, Step 4 of the drilling process is not unique. There might be more than one way to connect those disconnected components to the punctured polyomino. But since $`G`$ is finite, the number of ways of reconnection is bounded above. Let that bound be $`C^{}`$.

Together with the drilling mapping $`D`$ defined above, we can write,

$$s_n^{(0)}C^{}M\underset{j=nb^2}{\overset{n+(b^22)}{}}s_j^{(1)}.$$

By the increasing monotonicity of $`s_n^{(h)}`$ over $`n`$ (proven in lemma (A1)),

$$s_n^{(0)}(2b^21)C^{}Ms_{n+(b^22)}^{(1)}.$$

Set $`C=(2b^21)C^{}M`$ and the result follows. $`\mathrm{}`$

Now consider the drilling of $`h`$ holes. Denote by $`c`$ the ceiling of $`c`$, the smallest integer greater than or equal to $`c`$. Now, we could choose $`h`$ drilling locations from at least $`n/b^2`$ grid squares. To see this, consider first a polyomino with $`n`$ cells, where $`nb^2`$. The situation where we have the least number of drilling sites is when $`n`$ cells fall exactly within 1 grid square. In this way, we only have 1 possible grid square where we could drill holes. Similarly, if we have $`b^2n+w`$ cells, where $`n1`$ and $`0w<b^2`$, the minimal number of grid squares is when $`b^2n`$ cells fall into exactly $`n`$ grid squares and the other $`w`$ cells falls into another single grid. Then we have $`n+1`$ possible grid squares to drill holes. Therefore we have at least $`n/b^2`$ grid squares where we could drill holes.

Letting $`H`$ be a set of grid squares, let $`\mathrm{\Phi }_n^{[H]}`$ denote the set of all polyominoes with $`n`$ cells and $`|H|`$ holes where there is a hole in each grid square of $`H`$. Let $`s_n^{[H]}`$ denote its cardinality. One property of this set is that, for any $`n_1,n_2Z^+`$ such that $`n_1n_2`$,

$$s_{n_1}^{[H]}s_{n_2}^{[H]}.$$

###### Theorem A2

There exists a real constant $`K`$ such that for all $`b5`$,

$$\left(\genfrac{}{}{0pt}{}{n/b^2}{h}\right)s_n^{(0)}Ks_{n+(b^22)h}^{(h)},bn/b^2.$$
(A4)
Proof: Place the polyomino in the grid system. Let $`A^{}`$ be the set of all grid squares where we could drill holes. From previous results, we know there are at least $`n/b^2`$ elements in $`A^{}`$. Truncate the set $`A^{}`$ with only the first $`n/b^2`$ elements and call this set of drilling sites $`A`$.

Pick a subset $`HA`$ such that $`|H|=h`$ $`(n/b^2)`$ and drill holes one by one in each grid square in $`H`$, leading to a series of mappings:

$$M:\mathrm{\Phi }_n^{(0)}\mathrm{\Phi }_{m_1}^{(1)}\mathrm{\Phi }_{m_2}^{(2)}\mathrm{}\mathrm{\Phi }_{m_h}^{[H]}$$

where $`m_i`$’s are appropriate constants depending on each instance of operation and $`m_in+(b^22)h`$, $`i`$. So, from this mapping $`M`$, we have

$$s_n^{(0)}C^hs_{m_h}^{[H]}C^hs_{n+(b^22)h}^{[H]}.$$
(A5)
Hence

$`{\displaystyle \underset{HA}{}}s_n^{(0)}`$ $``$ $`{\displaystyle \underset{HA}{}}C^hs_{n+(b^22)h}^{[H]}`$
$`\left({\displaystyle \genfrac{}{}{0pt}{}{n/b^2}{h}}\right)s_n^{(0)}`$ $``$ $`C^h{\displaystyle \underset{HA}{}}s_{n+(b^22)h}^{[H]}`$
$`\left({\displaystyle \genfrac{}{}{0pt}{}{n/b^2}{h}}\right)s_n^{(0)}`$ $``$ $`Ks_{n+(b^22)h}^{(h)}.`$
The last line arises since $`_{HA}\mathrm{\Phi }_n^{[H]}\mathrm{\Phi }_n^{(|H|)}\mathrm{\Phi }_n^{(h)}`$ and put $`K=C^h`$. $`\mathrm{}`$

#### A.1.3 Surgery

The surgery operation removes a linear sequence of cells inside the polyomino and concatenates the sequence to the external boundary of the polyomino. Our objective is to join two holes thereby reducing the number of holes by one.

We divide our domain into three classes:

(1) the set of polyominoes which have at least one hole with size one cell; denote this set by $`\dot{\mathrm{\Phi }}_n^{(h)}`$ with cardinality $`\dot{s}_n^{(h)},`$ see figure A4(a);

(2) the set of polyominoes which are not in (1), and have three holes touching each other at corners around a single cell; denote this set by $`\ddot{\mathrm{\Phi }}_n^{(h)}`$ with cardinality $`\ddot{s}_n^{(h)}`$, see figure A4(b);

(3) the set of polyominoes which are not in (1) nor (2) (the general case), denote this set by $`\stackrel{~}{\mathrm{\Phi }}_n^{(h)}=\mathrm{\Phi }_n^{(h)}\backslash (\dot{\mathrm{\Phi }}_n^{(h)}\ddot{\mathrm{\Phi }}_n^{(h)})`$ with cardinality $`\stackrel{~}{s}_n^{(h)}`$.

##### The general case:

First, let’s look at polyominoes in $`\stackrel{~}{\mathrm{\Phi }}_n^{(h)}`$. Consider a polyomino $`\alpha _n^{(h)}\stackrel{~}{\mathrm{\Phi }}_n^{(h)}`$. Let $`L`$ be the set of all loop-free sequences<sup>1</sup><sup>1</sup>1A loop-free sequence is a sequence of cells where each successive pair of cells are joined and no cell appears more than once in the whole sequence of cells in $`\alpha _n^{(h)}`$, one end of which must touch the boundary of one hole and the other end must touch the boundary of another hole (or the exterior boundary). Define the length of the sequence to be the number of cells in the sequence. Let the set $`Z`$ be the set of sequences in $`L`$ that has minimum length.

Now pick one sequence $`zZ`$ and cut it out. This step will join two holes together (or join one hole and the exterior together). Next, using the definition in section A.1.1, concatenate $`\alpha _n^{(h)}czDz`$ ($`D`$ is a $`3\times 3`$ polyomino block), and the resulting polyomino has one less hole and nine more cells than the original one. So we have a mapping

$$S_3:\stackrel{~}{\mathrm{\Phi }}_n^{(h)}\mathrm{\Phi }_{n+9}^{(h1)}$$

This mapping is at most $`n`$-to-one because for any polyomino in the codomain, we could find at most $`n`$ locations to connect $`z`$ back. As a result, we find that,

$$\stackrel{~}{s}_n^{(h)}ns_{n+9}^{(h1)}.$$
(A6)
##### Special cases:

There are some subclasses of case (1) where we cannot apply the general surgery operation. One example is shown in figure A5 . In this example, the minimum length of the loop-free sequence connecting two boundaries is two, but when we try to remove any such minimum sequence, we will disconnect the polyomino.

To deal with case (1) polyominoes, we simply fill such a hole with a single cell. This satisfies our objective of reducing the number of holes by one. So we find a mapping

$$S_1:\dot{\mathrm{\Phi }}_n^{(h)}\mathrm{\Phi }_{n+1}^{(h1)}.$$

The mapping $`S_1`$ is $`n^{}`$ to 1 where $`n^{}`$ is less than $`n`$. To see this, consider a polyomino in the codomain. To map back to the domain, we can choose any interior cell to remove and the number of choices is obviously less than the total number of cells in the polyomino, so $`n^{}n`$. Therefore,

$$\dot{s}_n^{(h)}ns_n^{(h1)}.$$
(A7)
Finally, consider polyominoes in $`\ddot{\mathrm{\Phi }}_n^{(h)}`$. The problem with these polyominoes is that when we try to remove the cell between the 3 holes, we will inevitably join 3 holes together instead of joining 2. One way to get around this is to deliberately choose another sequence (which is also a cell in this class) in the set $`Z`$. In particular, we choose to remove the cell $`z^{}Z`$ such that $`z^{}`$ is the top cell (in lexicographic ordering) of $`Z`$ which is not the problem cell. Similarly to the general case, we find a mapping

$$S_2:\ddot{\mathrm{\Phi }}_n^{(h)}\mathrm{\Phi }_{n+9}^{(h1)}$$

and hence

$$\ddot{s}_n^{(h)}ns_{n+9}^{(h1)}.$$
(A8)
Adding (A6), (A7), (A8), we have

$$s_n^{(h)}=\dot{s}_n^{(h)}+\ddot{s}_n^{(h)}+\stackrel{~}{s}_n^{(h)}ns_{n+9}^{(h1)}+ns_{n+9}^{(h1)}+ns_n^{(h1)}$$

and hence:

###### Theorem A3

The number of polyominoes of area $`n`$ with $`h`$ holes is bounded above by $`3n`$ times the number of polyominoes of area (n+9) with one less hole. That is,

$$s_n^{(h)}3ns_{n+9}^{(h1)}$$
(A9)
### A.2 Growth constants (by area)

The following theorem proves the existence and equality of all the growth constants for $`k`$-punctured polyominoes for all finite $`k`$.

###### Theorem A4

There exists a constant $`\beta _0`$ such that for all $`h0`$

$$\underset{n\mathrm{}}{lim}n^1\mathrm{log}(s_n^{(h)})=\mathrm{log}(\beta _0).$$
(A10)
Proof: We use induction. First, $`\beta _0`$ exists . Assume $`\beta _h=lim_n\mathrm{}n^1\mathrm{log}(s_n^{(h)})`$ exists. From the results of the concatenation and surgery operations, we have

$$s_{nm}^{(h)}s_m^{(1)}s_n^{(h+1)}3ns_{n+9}^{(h)}.$$

Choose some value $`m`$ such that $`0<s_m^{(1)}<\mathrm{}`$, for example $`m=8`$. Take the logarithm, divide by $`n`$ and take the limit $`n\mathrm{}`$. This gives

$$\mathrm{log}(\beta _h)\mathrm{log}(\beta _{h+1})\mathrm{log}(\beta _h)$$

and hence $`\beta _h=\beta _{h+1}`$. Iterating from $`\beta _0`$ gives $`\beta _0=\beta _h`$ for all finite $`h`$. $`\mathrm{}`$

### A.3 Relationships between critical exponents

The following theorem establishes the relationships between critical exponents should they exist.

###### Theorem A5

Assume for all $`h`$ $`s_n^{(h)}C_hn^{\varphi _h}\beta _0^n`$ where $`C_h`$ is a $`h`$-dependent constant. Then

$$\varphi _h=\varphi _0h.$$
(A11)
Proof: From the results of the drilling operation, we have

$$C^h\left(\genfrac{}{}{0pt}{}{n/b^2}{h}\right)s_{n(b^21)}^{(0)}s_n^{(h)}.$$
(A12)
Since

$`\left({\displaystyle \genfrac{}{}{0pt}{}{n/b^2}{h}}\right)`$ $``$ $`{\displaystyle \frac{1}{h!}}n/b^2^h,`$
substituting the assumed asymptotic form $`s_n^{(h)}C_hn^{\varphi _h}\beta _h^n`$ into (A12), dividing by $`\beta _0^n`$, taking logarithms, dividing by $`\mathrm{log}(n)`$, letting $`n\mathrm{}`$ and using the above result, we get

$$h\varphi _0\varphi _h.$$
(A13)
Next, from the result of the surgery operation, we have

$`s_n^{(h)}3ns_{n+9}^{(h1)}`$ $``$ $`(3n)3(n+9)s_{n+18}^{(h2)}\mathrm{}`$
$``$ $`(3n)3(n+9)\mathrm{}3(n+9h)s_{n+9(h1)}^{(0)}[3(n+9h)]^hs_{n+9h}^{(0)}.`$
Again by substituting the assumed asymptotic form, dividing by $`\beta _0^{(n)}`$, taking logarithms, dividing by $`\mathrm{log}(n)`$ and letting $`n\mathrm{}`$, we get

$$\varphi _hh\varphi _0.$$
(A14)
Hence $`\varphi _h=\varphi _0h.`$ $`\mathrm{}`$