text:
# Three-Dimensional Hydrodynamic Simulations of Accretion Tori in Kerr Spacetimes

## 1 Introduction

More than three decades have passed since the first suggestions that black hole accretion is a ubiquitous astrophysical phenomenon (e.g. Lynden-Bell 1969), and there is now substantial and compelling observational evidence for it. The range of permissible masses for compact, unseen companions puts many binary systems squarely in the black hole domain (e.g. Cyg X-1, A0620-00; van Paradijs & McClintock 1995), while stellar kinematical observations of galactic cores (reviewed by Tremaine 1997) imply large amounts of mass ($`10^610^9M_{}`$) concentrated in small (parsec-sized) volumes. Indirect indications of accretion include copious nonthermal radiation from active galactic nuclei, relativistic expansion, and jets. Iron line observations (e.g., Tanaka et al. 1995) show how, in principle, observations can probe directly the physics of the strong-field region, distinguishing black hole accretion from accretion onto white dwarfs or neutron stars.

With the emergence of a unified model to account for the wide variety of observational characteristics of both X-ray binaries and active galactic nuclei (AGN), it becomes increasingly important to support theory and observation with detailed numerical simulations of the physics of black hole accretion. Distinctive and astrophysically interesting effects are expected from accretion flows in the Kerr metric, and modeling will require development of new fully general relativistic (GR) simulation codes. This becomes especially relevant in light of results that suggest that most galactic core black holes are rotating (Elvis, Risaliti, & Zamorani 2002).

Most of the simulations to date of relativistic black hole accretion have been for hydrodynamics alone. Extensive work was done in the 1980s on the numerical simulation of hydrodynamics around black holes (Hawley, Smarr & Wilson 1984a, 1984b; hereafter HSWa and HSWb). This work was based on numerical techniques developed through the pioneering efforts of Wilson (1972). The HSWb GR hydrodynamics code was written in spherical polar Boyer-Lindquist coordinates, and used operator-split finite-differencing on a staggered grid, with piecewise linear representations for the fundamental variables. The GR hydrodynamic studies using this code (HSWb; Hawley & Smarr 1986; Hawley 1986) were two-dimensional (assuming axisymmetry) and dealt with the prompt infall of rotating gas toward a central black hole, the role of the centrifugal barrier, centrifugally-driven outflows, and the formation of pressure supported thick disks or tori. Since then, there have been other axisymmetric GR hydrodynamic simulations using a similar numerical approach (Yokosawa 1995; Igumenshchev & Beloborodov 1997).

It is now recognized that magnetic fields play an essential role in the outward transport of angular momentum in accretion disks through the action of the magnetorotational instability of Balbus and Hawley (1991). This implies the need for a general relativistic magnetohydrodynamics (MHD) simulation code. Here again, Wilson did pioneering work, carrying out two-dimensional simulations of magnetized accretion over two and a half decades ago (Wilson 1975; 1977; 1978). The task was arguably beyond the computers then available, however. Only recently have research groups returned to full GR MHD simulations (e.g., Koide, Shibata, & Kudoh 1999; Komissarov 1999).

Rather than employing full relativity, the most recent three-dimensional global black hole MHD accretion simulations (e.g., Hawley & Balbus 2002) use a pseudo-Newtonian potential that can emulate certain important characteristics of the Schwarzschild metric. Computing a three dimensional MHD accretion flow in a full Kerr metric is a more ambitious undertaking. It is our aim to develop a new three-dimensional GR MHD code to enable such global accretion disk studies. This paper represents a first step, specifically the initial development and application of the hydrodynamic portion of a Kerr metric GR code.

The culmination of the earlier GR hydrodynamics effort was a fully three-dimensional numerical investigation (Hawley 1991; hereafter H91) of the Papaloizou-Pringle instability (Papaploizou & Pringle 1984) for accretion tori in a Schwarzschild metric. In this paper we extend those simulations to Kerr (rotating) black holes, and consider the effect of black hole angular momentum on prograde and retrograde torus orbits. The Papaloizou-Pringle instability (PPI) remains an interesting topic of research, and it provides a nontrivial application of full three-dimensional general-relativistic hydrodynamics.

The plan of this paper is as follows. In §2 we outline the equations and the general numerical procedures. In §3 we discuss the properties of equilibrium gas tori in the Kerr metric, and the properties of the Papaloizou-Pringle instability to which these tori are vulnerable. The results of a series of numerical simulations of tori in the Kerr metric are presented in §4, and these results are discussed in §5.

## 2 Equations and Numerical Methods

We wish to study the evolution of a fluid in the background spacetime of a Kerr (rotating) black hole. We adopt Boyer-Lindquist coordinates, $`(t,r,\theta ,\varphi )`$, for which the line element has the form,

$$ds^2=g_{tt}dt^2+2g_{t\varphi }dtd\varphi +g_{rr}dr^2+g_{\theta \theta }d\theta ^2+g_{\varphi \varphi }d\varphi ^2.$$
(1)
In keeping with Misner, Thorne, & Wheeler (1973), we use the metric signature $`(,+,+,+)`$, along with geometrodynamic units where $`G=c=1`$; the black hole mass is unity, $`M=1`$. The determinant of the 4-metric is $`g`$, and $`\sqrt{g}=\alpha \sqrt{\gamma }`$ where $`\alpha `$ is the lapse function, $`\alpha =1/\sqrt{g^{tt}}`$, and $`\gamma `$ is the determinant of the spatial $`3`$-metric.

For a relativistic test fluid described by a density $`\rho `$, specific internal energy $`ϵ`$, and $`4`$-velocity $`U^\mu `$, we define the transport velocity $`V^\mu `$ as follows:

$$V^\mu =\frac{U^\mu }{U^t},$$
(2)
where $`U^t=W/\alpha `$, and $`W`$ is the gravitational redshift factor,

$$W=\frac{1}{(1V^\mu V_\mu )^{1/2}}.$$
(3)
We also define the momentum,

$$S_\mu =\rho hWU_\mu ,$$
(4)
and auxiliary density and energy functions $`D=\rho W`$ and $`E=Dϵ`$.

Using these variables the conservation laws can be rewritten into a form suitable for finite differencing. The equation for mass conservation is written

$$_tD+\frac{1}{\sqrt{\gamma }}_j(D\sqrt{\gamma }V^j)=0.$$
(5)
Conservation of the fluid energy-momentum tensor $`_\mu T^{\mu \nu }=0`$ yields momentum conservation equations,

$$_tS_i+\frac{1}{\sqrt{\gamma }}_j(S_i\sqrt{\gamma }V^j)+\alpha _jP+\frac{1}{2}\frac{S_\mu S_\nu }{S^t}_jg^{\mu \nu }=0,$$
(6)
and an internal energy conservation equation

$$_tE+\frac{1}{\sqrt{\gamma }}_j(E\sqrt{\gamma }V^j)+P_tW+\frac{P}{\sqrt{\gamma }}_j(W\sqrt{\gamma }V^j)=0.$$
(7)
Spatial indices are indicated by roman characters $`i,j=1,2,3`$. We assume an ideal gas equation of state $`P=\rho ϵ(\mathrm{\Gamma }1)`$, where $`\mathrm{\Gamma }`$ is the adiabatic exponent. This is the same system of equations as described in HSWa.

The GR hydrodynamics code evolves time-explicit, operator-split, finite difference forms of equations (5)—(7). The algorithm is a three-dimensional generalization of the solver described in HSWb, with some additional modifications. For example, velocity renormalization, $`U^\mu U_\mu =1`$, which is invoked in the source and transport steps following the solution of the discretized momentum equation (6) is now implemented using the algebraically equivalent condition $`S^\mu S_\mu =(D+\mathrm{\Gamma }E)^2`$ for improved numerical stability. Several improvements are possible because of the increase in computer power since 1984. Analytically determined metric terms and metric derivatives are now calculated and stored at all needed grid locations rather than averaged. Array operations have also been streamlined wherever possible using Fortran 90 array syntax. The code uses message passing parallelism with a form of domain decomposition, where the global grid is partitioned into subgrids, with each subgrid assigned to a processor. Data on each subgrid is evolved independently during the source and transport phases of a timestep and data on subgrid boundaries is exchanged at the end of each phase through message-passing calls. This results in a highly scalable code that exhibits good speedup over the full range of practically realizable subgrids.

The present simulations have been designed as a follow-on to the results in H91, but unlike the code used there, we do not need to impose equatorial symmetry. The calculations performed on a $`64\times 32\times 64`$ $`(r,\theta ,\varphi )`$-grid in H91 correspond to a $`64\times 64\times 64`$ grid here. The ergosphere can be part of the computational domain, but at the event horizon, $`r_H=M+\sqrt{M^2a^2}`$, some of the Boyer-Lindquist metric terms are singular, so the inner edge of the radial grid must lie at some point outside the horizon, continuing out to some selected outer grid boundary. Although the horizon itself can never be part of the computational grid, tests indicate that it is possible to come arbitrarily close to the horizon (e.g., $`r_{in}=2.001M`$ for the Schwarzschild case). But this comes at a steep price: the narrow region where the metric terms diverge rapidly must be covered by a large number of grid zones (in the extreme Kerr limit, the number of zones inside the ergosphere could be as high as the number of zones that lie outside it). In addition to consuming large amounts of memory in a three dimensional simulation, these innermost grid zones can significantly reduce the time step size, greatly increasing the total computational time. The specific location of the inner edge of the grid is determined by a balance between the physical phenomena of interest in a simulation and the memory and performance demands. For the hydrodynamic torus simulations considered here, the region near the horizon is of secondary interest, so the inner grid boundary, $`r_{min}`$, was set at or just inside the static limit on the equator, $`r_{min}r_{static}|_{\theta =\pi /2}=2M`$. In the runs reported here, this ranged from $`r_{in}=1.90M`$ for the smaller grids to $`r_{in}=2.05M`$ for the larger grids. The polar grid extends almost to the rotational axes at $`\theta =0`$ and $`\pi `$, but the axes themselves are not included. As in non-relativistic codes, the polar axes are problematic due to the shrinking of zone volumes, and the collapse of the inner radial zone face at $`r=0`$. These problems are further compounded by the fact that, here again, some metric components are singular. Torus stability studies do not require the axes, and their omission is not important to the present simulations. Regularized operators will be introduced in future work. The radial and polar grids use logarithmic scaling for all runs in order to concentrate zones near the horizon and the equator. The azimuthal angle $`\varphi `$ spans the full range from $`0`$ to $`2\pi `$ with equally-spaced grid zones and periodic boundary conditions. Previous studies of the Papaloizou-Pringle instability show that using the full $`\varphi `$ range is important as the most important unstable modes have azimuthal wavenumber $`m=1`$.

The time step $`\mathrm{\Delta }t`$ is determined by the extremal light-crossing time for a grid zone, as described in HSWb. This time step size remains fixed for the duration of the simulation since it is a purely geometric result.

The initial state for a simulation is generated from the thick disk equilibrium solution (described in §3) with given input parameters. This initial configuration is then perturbed with small, random enthalpy fluctuations with a maximum amplitude of 1%. These perturbations are the seeds from which the full PPI will develop.

## 3 Tori and the Kerr Metric

The initial stationary, axisymmetric analytic torus solution is presented in HSWa. For convenience, the essential results are repeated here.

The initial state for the simulations is a stationary, axisymmetric solution ($`_t=_\varphi =0`$) to equations (5)—(7) with no internal poloidal motion, i.e. $`U^r=U^\theta =0`$. The combined centrifugal and gravitational accelerations (which together make up an effective potential) are balanced by pressure gradients, keeping the disk in equilibrium. To develop this solution, define the specific angular momentum ($`l`$) and angular velocity ($`\mathrm{\Omega }`$), as

$$U_\mu =U_t(1,0,0,l)$$
(8)
$$U^\mu =U^t(1,0,0,\mathrm{\Omega }).$$
(9)
Applying the orthogonality condition $`U^\mu U_\mu =1`$ leads directly to the following expression for $`U_t`$

$$U_t=\left(g^{tt}2lg^{t\varphi }+l^2g^{\varphi \varphi }\right)^{1/2}$$
(10)
which is related to the binding energy, $`e_{bind}=U_t`$. Here, we note that equation (94a) of HSWa should be extended to allow a more general disk outer boundary, which is specified by the binding energy at the surface $`e_{surf}=U_{t_{lim}}`$. The more general form of equation (94a) is

$$ϵ=\frac{1}{\mathrm{\Gamma }}\left(\frac{U_{t_{lim}}}{U_t}1\right).$$
(11)
For a constant entropy adiabatic gas the pressure is given by $`P=\rho ϵ(\mathrm{\Gamma }1)=K\rho ^\mathrm{\Gamma }`$, and density is given by $`\rho =\left[ϵ(\mathrm{\Gamma }1)/K\right]^{1/(\mathrm{\Gamma }1)}`$. These relations completely specify the initial equilibrium torus.

A particular constant angular momentum thick disk solution is specified by choosing the angular momentum $`l`$, the binding energy at the surface of the torus (as above), and the entropy parameter $`K`$. In keeping with HSWb, $`K=0.01`$ is fixed for all simulations, so only the specific angular momentum $`l`$, surface binding energy $`e_{surf}`$, and the Kerr parameter $`a`$ are varied.

### 3.1 The Papaloizou-Pringle Instability

The stability of tori with constant specific angular momentum $`l`$ has been of interest since the work of Papaloizou and Pringle (1984), who established that these tori are unstable to non-axisymmetric global modes. It has since been demonstrated that although other angular momentum distributions are also unstable, the constant-$`l`$ tori are the most susceptible to the PPI. Subsequent work, notably Narayan, Goldreich, and Goodman (1987), consolidated the central features of the instability into the following picture. The global unstable modes have a co-rotation radius within the torus; the co-rotation is located in a narrow region where waves cannot propagate; this region separates inner and outer regions where wave propagation is possible; waves can tunnel through the corotation zone and interact with waves in the other region; and the transmitted modes are amplified only if there is a feedback mechanism, usually in the form of a reflecting boundary at the inner and/or outer edge of the torus.

However, linear stability analysis can go only so far, and numerical work is required to probe the non-linear effects that help determine the final amplitude of the global modes, and hence the cumulative effect of the instability on the torus. Hawley (1987) carried out a relativistic numerical study of the evolution of the PPI into the fully non-linear regime in 2D height-integrated disks. Although this paper verified that the azimuthally averaged rotation law ($`\mathrm{\Omega }r^q`$) exhibits the expected effects of angular momentum redistribution ($`q_{sat}\sqrt{3}`$), it also showed that the distribution of matter was not azimuthally symmetric, but instead took the form of counter-rotating epicyclic vortices, or “planets”, with $`m`$ planets emerging from the growth of a mode of order $`m`$. The stability of thick tori in general relativity was studied by Blaes and Hawley (1988) who evolved unstable tori in a two-dimensional $`(r,\varphi )`$ Schwarzshild metric. This work led to the full three dimensional simulations in H91. The present work extends these calculations to include the effects of nonzero black hole angular momentum. Although we do not repeat the linear analysis of Blaes and Hawley (1988) with a non-zero Kerr parameter $`a`$, we know from previous work that mode growth depends very sensitively on the initial torus parameters. Black hole rotation simply adds additional complexity.

A point of particular interest has been the effect of accretion on the evolution of the PPI. Blaes (1987) computed growth rates for two-dimensional tori in the Schwarzschild metric and found that they went to zero for models with even a small net accretion flow past the location of the marginally-stable orbit. It was argued that the loss of the inner reflecting boundary through the development of an accretion flow suppresses mode growth. This result was subsequently confirmed analytically by Gat and Livio (1992) for a Newtonian potential. In the three-dimensional case, things are not so clear cut. H91 argues that since the accretion flow is confined to the equatorial plane, there is still a reflecting boundary above and below the plane. However, Dwarkadas & Balbus (1996) argued that the stabilizing effect of an accretion flow is due more to the dynamics at the corotation point than to the absence of reflection.

### 3.2 Disk Properties and the Kerr Parameter

One of the goals of H91 was to study the effect of the proximity of the cusp in the relativistic potential to the inner edge of the torus, and how accretion into the black hole affects mode growth. Thus the majority of the H91 runs dealt with tori with angular momentum above the marginally bound value, that is the smallest specific angular momentum for which a particle orbit has an inner turning point. In H91, an initial torus configuration was specified by choosing the angular momentum and the surface binding energy. Here, the structural properties of the initial torus are also influenced by the Kerr parameter. We therefore need to understand the interdependence of $`l`$, $`U_{t_{lim}}`$, and $`a`$ in order to specify initial configurations for the simulations.

To help guide the choice of $`l`$, it is useful to consider the quantity $`l_{mb}`$ that gives the specific angular momentum for a test particle on a marginally bound orbit. Particles on retrograde orbits require larger absolute values of the specific angular momentum to remain out of the hole than do particles on prograde orbits. This can be seen in Figure 1, which plots the dependence of $`l_{mb}`$ on the Kerr parameter $`a`$; the relation $`l_{mb}(a)`$ can be found, for instance, in Frolov and Novikov (1998). The quantity $`l_{mb}`$ also helps distinguish families of tori that have stable orbits from those that do not. Tori with $`l=l_{mb}`$ are referred to as marginal tori from hereon.

The PPI affects slender, intermediate, and wide tori in different ways, so the choice of torus dimensions is important. Given a value for $`l`$, the choice of surface binding energy, $`U_{t_{lim}}`$, sets the overall size of the disk. Figure 2 illustrates how different choices of $`U_{t_{lim}}`$ determine the width of the torus, here shown for $`l=4.5`$ and a range of choices for $`a`$. This illustration reinforces the notion that the torus is created by “filling” a local minimum in the energy surface; it also shows that an inner bound for a torus may not exist as $`a`$ is made increasingly negative, and the reason for this is best seen in the next figure.

The role of the Kerr parameter $`a`$ on the disk configuration is illustrated by Figure 3, which shows equipotentials for disks with $`l=4.5`$, and the same values of $`a`$ as Figure 2. The equipotentials are defined by the Boyer energy parameter $`\mathrm{\Phi }=\mathrm{log}(U_t)`$. This figure presents polar “slices” through the axisymmetric disk. The set of level contours is the same in each plot, so the effect of $`a`$ is readily apparent: the equipotentials shrink in the prograde case, and expand in the retrograde case, a result consistent with what was shown in Figure 2 for the equatorial plane. So, for a given choice of surface binding energy, a retrograde torus would be larger, and a prograde torus would be more compact, than the Schwarzschild case.

In H91 the value of $`l=4.5`$ was chosen since it corresponded to a bound torus in the Schwarzschild metric. Figures 1 and 2 hint that this choice of $`l`$ may not correspond to bound tori for all values of $`a`$ down to $`a=1`$. Figure 3 further reinforces this observation, and also provides a clearer picture: many of the contours for the $`a=0.8`$ case are open onto the black hole, and the viable range of parameters for which a bound torus can be constructed is growing smaller. For choices of $`a`$ below $`a=0.8`$ virtually all contours are open onto the hole, so an equilibrium initial state cannot be constructed.

The transition from closed to open contours can best be seen in Figure 4, which shows a plot of equipotentials for $`l`$ at and just below the marginally bound value for a disk orbiting a Kerr hole with $`a=0.5`$ (the point in Figure 1 where the $`l_{mb}`$ curve passes through $`l=4.5`$). The effect near the cusp is striking: a small reduction in $`l`$ below the marginally bound value opens a pathway to the black hole in the equatorial plane, so that a disk whose structural parameters give it an inner bound near $`r=4.95M`$ would not be expected to remain in equilibrium, and any small departure from equilibrium would trigger accretion. The opening of the potential cusp can also be accomplished by choosing $`a<0.5`$ for an $`l=4.5`$ torus, which explains the shape of the contours in the left panel of Figure 3.

Finally the analytic thick disk solution is characterized by a number of derived parameters. One is the location of the pressure maximum, $`r_{P_{max}}`$, which is obtained by root-finding using the condition $`_rϵ(r,\pi /2)=0`$ since the pressure and energy maxima are coincident. The angular frequency at the pressure maximum, $`\mathrm{\Omega }_{P_{max}}=(g^{t\varphi }lg^{\varphi \varphi })/(g^{tt}lg^{t\varphi })`$, and the orbital period at the pressure maximum, $`T_{orb}=2\pi /\mathrm{\Omega }_{P_{max}}`$. Time in the simulations will be reported in terms of this orbital period. The boundaries of the disk in the equatorial plane are $`r_{in}`$ and $`r_{out}`$. Since $`r_{in}`$ is an input parameter (through the choice of $`U_{t_{lim}}`$), a simple plot of $`\rho (r,\pi /2,\varphi )`$ helps locate $`r_{out}`$, which is in turn used to specify the required outer boundary of the radial grid (usually set at about $`r_{max}2r_{out}`$).

### 3.3 Torus Diagnostics

The Kerr GR hydrodynamics code is used to study the evolution of the Papaloizou-Pringle instability in thick tori. In keeping with H91, the evolution of the instability is characterized by several reduced quantities extracted from the simulation data. These include the maximum density enhancement at mode saturation, $`\delta \rho /\rho `$, the Fourier power in density for the azimuthal wavenumbers, $`m=1`$ and $`2`$, $`q_{sat}`$, the angular velocity distibution parameter at mode saturation, and the mass flux in the equatorial plane (to help shed light on the relationship between accretion and the growth of the PPI). In all models, the calculations are done $`20`$ times per orbit at the initial pressure maximum. The only exception is the mass flux, which is computed once per orbit.

The density enhancement, $`\delta \rho /\rho `$, is obtained by finding the maximum density at in the data arrays, and computing $`\delta \rho =\rho _{max}|_{t=t_{sat}}\rho _{max}|_{t=0}`$.

We extract the $`m=1`$ and $`m=2`$ Fourier modes by computing azimuthal averages using the numerical equivalents of

$`\mathrm{}k_m(r)`$ $`=`$ $`{\displaystyle _0^{2\pi }}\rho (r,\pi /2,\varphi )\mathrm{cos}(m\varphi )𝑑\varphi `$ (12)
$`\mathrm{}k_m(r)`$ $`=`$ $`{\displaystyle _0^{2\pi }}\rho (r,\pi /2,\varphi )\mathrm{sin}(m\varphi )𝑑\varphi .`$ (13)
The mode power is then

$$f_m=\frac{1}{r_{out}r_{in}}_{r_{in}}^{r_{out}}\mathrm{log}_e\left((\mathrm{}k_m(r))^2+(\mathrm{}k_m(r))^2\right)𝑑r.$$
(14)
A linear fit is performed to the time-sequenced data to extract a mode growth rate. For this calculation $`r_{in}`$ and $`r_{out}`$ are the initial values for the inner and outer edges of the torus in the equatorial plane.

In H91, the parameter $`q_{sat}`$ was obtained by a radial power law fit to the azimuthally averaged angular velocity, $`\overline{\mathrm{\Omega }}(r)r^q`$, at mode saturation. This parameter was used to characterize deviations from a purely Keplerian profile. Since H91 dealt with Schwarzschild black holes, the angular velocity for the equilibrium fat disk had a simple form, $`\mathrm{\Omega }=U^\varphi /U^tg_{tt}r^2`$, i.e. the usual Keplerian profile multiplied by the redshift factor. It was therefore straightforward to compare the final disk profile against the Keplerian case. With Kerr black holes, the equilibrium fat disk has a more complicated radial dependence, $`\mathrm{\Omega }=(g^{t\varphi }lg^{\varphi \varphi })/(g^{tt}lg^{t\varphi })`$. However, since the aim is to measure a change in the equatorial angular velocity profile, we adapt the procedure. As with H91, we obtain the azimuthally averaged angular velocity using the numerical equivalent of

$$\overline{\mathrm{\Omega }}(r)=\frac{1}{2\pi }_0^{2\pi }V^\varphi (r,\pi /2,\varphi )𝑑\varphi .$$
(15)
A power-law fit $`\overline{\mathrm{\Omega }}(r)r^q`$ is obtained from the slope of a log-log plot of $`\overline{\mathrm{\Omega }}(r)`$ for $`r(r_{in},r_{out})`$ at the time step corresponding to mode saturation. We also extract the initial value $`q_0`$ in an analogous manner, and report the change $`\delta q=q_{sat}q_0`$ as a measure of the redistribution of angular velocity.

The azimuthally-averaged mass flux in the equatorial plane is computed using the numerical equivalent of

$$\dot{M}_{eq}(r)=\frac{\sqrt{\gamma (r,\pi /2)}}{2\pi }_0^{2\pi }D(r,\pi /2,\varphi )V^\varphi (r,\pi /2,\varphi )𝑑\varphi .$$
(16)
In the figures below $`\dot{M}_{eq}(r)`$ is plotted at a radius $`r`$ lying just inside the initial inner edge of the torus $`r_{in}`$ and is used to establish the presence of a flow of matter towards the black hole.

## 4 Results

Before proceeding with the simulations in the Kerr metric, test runs were carried out with $`a=0`$ to establish that this new version of the GR hydrodynamics solver was able to reproduce the original results in H91. These tests were successful and are not elaborated upon here.

### 4.1 Parameters for Simulations

The parameters for the 3D runs in H91, summarized in Table 1, fall into two broad categories that are distinguished by specific angular momentum. The $`l=4.5`$ models, A3, B3, and C3, are tori with innermost radii (in the equatorial plane) of $`r_{in}=11.0M`$, $`9.5M`$, and $`8.0M`$ respectively; these represent slender, intermediate, and wide tori. The marginal model, E3, is a wide torus with $`r_{in}=4.4M`$.

In adapting these models for non-zero Kerr parameters, we are confronted with the change in torus width with $`a`$, and also with the restricted range over which the $`l=4.5`$ torus remains bound. To emphasize the role of $`a`$, we chose the extreme Kerr limits wherever possible. Six models were selected, and are listed in Table 2. We use the same naming convention as in H91; models are named according to the choice of the pair of parameters $`(l,r_{in})`$, here augmented by the Kerr parameter $`a`$. Model A3p represents a prograde $`l=4.5`$ torus in the extreme Kerr limit; the slender A3 torus has expanded to intermediate size with this change. Similarly, prograde model B3p has expanded into a wide torus. Model B3r represents the “last” bound retrograde $`l=4.5`$ torus, for $`a=0.8`$; this model is also well below the marginally bound value for $`a=0.8`$ (it could be termed a sub-marginal torus). The intermediate B3 torus has shrunk considerably with this change in $`a`$. The E3 models are more loosely based on the H91 E3 original; this was done to maintain their nature as marginal tori, and allowing the values of $`r_{in}`$ to depart from the H91 value. Model E3p is a marginal prograde torus with $`a=0.5`$; the wide marginal E3 torus has shrunk considerably with this choice of parameters. Model E3r is an intermediate, retrograde torus in the extreme Kerr limit. One additional model, X3p, was introduced to complete the data set with a prograde marginal torus that is related to E3p, but whose inner radius lies closer to the static limit.

Table 3 is a summary of the simulation results. The growth rate of the first and second Fourier modes is given (in units of $`\mathrm{\Omega }`$), as is the time at which mode saturation occurred (in units of $`T_{orb}`$), the change in power-law exponent for the best-fit to the azimuthally-averaged angular velocity at mode saturation, and the density enhancement at mode saturation. Mode saturation is defined as the time at which the first Fourier mode reaches its maximum amplitude. For reference the corresponding results from H91 are given in Table 4.

### 4.2 Slender Torus - Model B3r

Model B3r is a slender torus rotating in the retrograde sense around a Kerr black hole with $`a=0.8`$. The torus has an initial maximum density of $`\rho _{max_0}=0.031`$ at $`r=12.8M`$. Figure 5(a) shows that at mode saturation a planet has formed, with a density maximum at $`r11M`$ and a fractional density enhancement of $`0.26`$. The planet has a crescent shape with some evidence of a tightly-wrapped outward spiral wave in the lower left quadrant of panel (a). A gray-scale enhanced view of the inner region near the static limit clearly shows that a strong inflow has developed that contacts the static limit in the upper left quadrant of panel (b); there is also a hint of frame dragging in the reversal of the sense of the flow near the static limit (in this figure, as in all subsequent ones, the black hole rotates in the counter-clockwise sense, and the retrograde torus here rotates in the opposite sense). There is also a evidence of a weaker inbound spiral of matter in the right half of panel (b). The growth of the PPI modes has two distinct stages, as seen in Figure 5(c). The $`m=2`$ mode grows quickly at the beginning, then levels off at $`2T_{orb}`$ and remains roughly constant until $`7T_{orb}`$ before increasing again, paralleling the growth of the $`m=1`$ mode. The $`m=1`$ mode growth does not become apparent until approximately $`3T_{orb}`$, at which point it exhibits strong linear growth to saturation at $`9T_{orb}`$. Up to $`4T_{orb}`$, the $`m=2`$ mode is dominant, but lags the $`m=1`$ mode in the later stages of the simulation. The development of the accretion flow is clarified in Figure 5(d), which shows the mass influx in the equatorial plane inside the inner edge of the torus. There is no significant accretion until $`8T_{orb}`$. The accretion rate reaches a maximum shortly after mode saturation. This feature, as will be seen in subsequent runs, is associated with all solutions that yield planets.

The effect of frame dragging is illustrated in Figure 6, which is an equatorial slice at $`10.02T_{orb}`$, one orbit after saturation. The left panel of the figure is a density plot with 15 contours, chosen to highlight the stream of high-density material flowing from the inner edge of the grid into the hole. The right panel shows momentum direction vectors inside the flow stream outlined by the overlaid contour. The expected reversal in azimuthal flow direction due to frame dragging begins at $`r4M`$. The direction becomes progressively more prograde as the flow approaches the the inner edge of the grid, which lies at $`2.05M`$ (just outside the static limit). This reversal, from retrograde to prograde, produces the dog-leg pattern seen in the density plot as the flow approaches the static limit. A similar pattern is visible for model E3r (see below).

### 4.3 Intermediate Torus - Model A3p

Model A3p is an intermediate torus rotating in the prograde sense around a Kerr black hole with $`a=1.0`$. The torus has an initial maximum density of $`\rho _{max_0}=0.23`$ at $`r=17.2M`$. Figure 7(a) shows that at mode saturation a planet has formed, with a density maximum at $`r11M`$ and a fractional density enhancement of $`0.57`$. The planet has a narrow crescent shape with an outward spiral wave visible in the upper left quadrant of panel (a). A gray-scale highly enhanced ($`\times 50`$) view of the inner region near the static limit clearly shows that a weak inspiral of material is underway at saturation; however, this is by far the weakest of the inflows in the set of planet-forming simulations. The growth of the PPI modes is strong from the outset, as seen in Figure 7(c). The $`m=1`$ and $`m=2`$ modes track one another, with the $`m=1`$ mode growing at a slightly greater rate. The $`m=1`$ mode is also less erratic than the $`m=2`$ mode. Both modes exhibit strong linear growth up to saturation at $`10T_{orb}`$. A small amount of inflow is present at saturation, as can be seen in Figure 7(d), which shows the mass influx in the equatorial plane inside the inner edge of the torus. The inflow of matter is quiescent until $`8T_{orb}`$, and reaches a maximum at mode saturation. The apparent reversal of the flow after $`10T_{orb}`$, combined with the lack of evidence for substantial inflow in panel (b), suggests that the planet is settling down to a new equilibrium state, and that this apparent outward flow of matter is simply indicative of the reorganization of matter within the disk.

### 4.4 Wide Torus - Model B3p

Model B3p is a wide torus rotating in the prograde sense around a Kerr black hole with $`a=1.0`$. The torus has an initial maximum density of $`\rho _{max_0}=0.72`$ at $`r=17.2M`$. Although the disk extends to $`r90M`$, most of the matter is concentrated inside $`r30M`$. Figure 8(a) shows that at mode saturation a planet has formed, with a density maximum at $`r11M`$ and a fractional density enhancement of $`0.35`$. The planet has a tightly-wrapped crescent shape with a noticeable outward spiral wave in the upper right quadrant of panel (a). A gray-scale enhanced view of the inner region near the static limit clearly shows that an inflow has developed that contacts the static limit in the lower right quadrant of panel (b). There is also a hint of frame dragging of the flow near the static limit where dense material appears to smear out into a ring. A weaker inbound spiral of matter in the upper half of panel (b) is also feeding the inner ring of matter. The growth of the PPI modes is strong from the outset, as seen in Figure 8(c). The $`m=1`$ and $`m=2`$ modes track one another, with the $`m=1`$ mode growing a a slightly greater rate. The $`m=1`$ mode is also less erratic than the $`m=2`$ mode. Both modes exhibit strong linear growth up to saturation at $`13T_{orb}`$; the saturation peak is very broad, making a precise determination of $`t_{sat}`$ difficult. As with the previous two models, inflow is triggered at saturation, as can be seen in Figure 8(d), which shows the mass influx in the equatorial plane inside the inner edge of the torus. The flow of matter is quiescent until $`11T_{orb}`$, and reaches a maximum at $`13T_{orb}`$, again corresponding with mode saturation.

### 4.5 Marginal Torus - Model E3p

Model E3p is a marginal torus rotating in the prograde sense around a Kerr black hole with $`a=0.5`$. The torus has an initial maximum density of $`\rho _{max_0}=0.86`$ at $`r=6.8M`$. Figure 9(a) shows that at mode saturation an elongated planet has formed, with a density maximum at $`r6M`$ and a fractional density decrease of $`0.01`$. A gray-scale view ($`\times 2`$ enhancement) of the inner region near the static limit clearly shows that an inflow has developed that contacts the static limit in the upper right quadrant of panel (b). Figure 9(c) shows that the $`m=1`$ mode undergoes steady growth from the outset, with some variability in the first $`10`$ orbits. The $`m=2`$ mode tracks the $`m=1`$ mode until $`12T_{orb}`$, where it flattens out until $`25T_{orb}`$, then grows again until the $`m=1`$ mode saturates at $`33T_{orb}`$. Perhaps the most striking feature of this model can be seen in the mass influx curve, Figure 9(d), again plotted for a point in the equatorial plane inside the inner edge of the disk. The transition between quiescence and accretion is less sharply defined than in other models, and a strong inflow begins after $`25T_{orb}`$, roughly corresponding to the point where the $`m=2`$ mode begins to grow again. Compared to the other models, the amount of accretion in this model is remarkable; panels (a) and (b) clearly show very prominent in-bound spirals of matter.

### 4.6 Marginal Torus - Model E3r

Model E3r is a marginal torus rotating in the retrograde sense around a Kerr black hole with $`a=1.0`$. The torus has an initial maximum density of $`\rho _{max_0}=0.43`$ at $`r=15.7M`$. Figure 10(a) shows that at mode saturation the disk remains largely unperturbed, with an insignificant fractional density decrease of $`0.01`$. Figure 10(b) shows a gray-scale enhanced view of the inner region near the static limit clearly shows that an inflow has developed that contacts the static limit in the upper right quadrant of panel (b); it must be noted that this inflow is extremely weak ($`\times 100`$ enhancement) yet it bears some similarity to the other retrograde model, B3r, in that the stream of matter exhibits a change in flow direction near the static limit. The growth of the PPI modes is weak, with the $`m=1`$ mode levelling off quickly after a brief initial growth, and the $`m=2`$ mode showing little activity, as seen in Figure 10(c). The $`m=1`$ mode reaches a very broad, weak maximum at $`20T_{orb}`$, making a precise determination of $`t_{sat}`$ difficult. Unlike the previous models, a noticeable amount of inflow is present from the outset as can be seen in Figure 10(d), which shows the mass influx in the equatorial plane inside the inner edge of the torus. Since this model is marginal, the early onset of an inflow is to be expected. In comparing panels (c) and (d), we see that the inflow is established before any significant mode growth has occured; inflow levels off to a more or less constant value in the first 5 orbits; the mode growth also levels off in this period. This would suggest that early accretion is an effective mechanism for inhibiting mode growth. This observation is strengthened by the last model to be considered, X3p.

### 4.7 Marginal Torus - Model X3p

Model X3p is a marginal torus rotating in the prograde sense around a Kerr black hole with $`a=0.5`$. The torus has an initial maximum density of $`\rho _{max_0}=1.58`$ at $`r=6.4M`$. Figure 11(a) shows that at mode saturation the disk remains largely unperturbed, with an insignificant fractional density decrease of $`0.01`$. Figure 11(b) shows a gray-scale enhanced view of the inner region near the static limit clearly shows that a highly symmetric, weak inflow has developed. The growth of the PPI modes is weak, with both modes tracking one another, as seen in Figure 11(c). Although Table 3 would suggest that these modes have a vigorous growth rate, their maximum amplitude is comparable to that of model E3r, and the outcome is very similar to that model, which is interesting since model X3p is much more closely related, in terms of its structural parameters, to model E3p. The $`m=1`$ mode reaches a very broad, weak maximum at $`5.5T_{orb}`$. As with model E3r, a noticeable amount of inflow is present from the outset as can be seen in Figure 11(d), which shows the mass influx in the equatorial plane inside the inner edge of the torus. As with E3r, this model is marginal, and the early onset of an inflow is to be expected. In comparing panels (c) and (d), we see that the inflow is established before any significant mode growth has occured; inflow levels off to a more or less constant value in the first 4 orbits; the mode growth also levels off in this period. As with model E3r, early accretion appears to be an effective mechanism for inhibiting mode growth.

## 5 Conclusions

The numerical simulations presented in this paper were undertaken to test the hydrodynamic portions of a three-dimensional general relativistic magnetohydrodynamic solver that is being developed by the authors. We studied the growth of the Papaloizou-Pringle instability in tori orbiting Kerr black holes. As was noted in H91, the parameter space spanned by tori in the Schwarzschild spacetime is a rich one, and the detailed behaviour of the unstable modes is very much model-dependent. This is certainly still true now that we have added black role angular momentum to the mix. However, several general features do emerge from the study of tori in the Kerr spacetime.

Previous analytic and numerical work has established that the rate of growth of unstable modes decreases as the width of the torus increases. Our results for models A3p, B3p, and B3r show that this remains the case for tori orbiting Kerr black holes.

As in H91, we see a decrease in the density enhancement with torus width in the bound tori. However, only models A3p and B3p are bound tori, with $`l>l_{mb}`$, (model B3r is sub-marginal, and the other three models are marginal). Furthermore, as was noted in H91, the detailed behaviour of mode growth depends sensitively on the choice of initial model, and mode growth dictates the redistribution of matter and hence the magnitude of the density enhancements.

Two-dimensional studies have shown that accretion suppresses mode growth due to the loss of the inner reflecting boundary, which is a prerequisite for mode amplification. This conclusion was supported by numerical and analytic work in previous papers, as discussed in section 3.1. However, the interplay between accretion and mode growth in three-dimensional simulations is less clear. The argument has been made, in H91, that since inflow is primarily concentrated on the equatorial plane, the loss of a reflecting boundary is localized to this region, and mode growth can still occur given the right circumstances in off-equator regions of the disk. Models B3r, A3p, and B3p show that mode growth in the absence of early inflow produces planets, while still generating appreciable transient inflows as a consequence of the redistribution of matter accompanying planet formation. In the converse case, we have seen that early inflow (marginal models E3r and X3p) effectively inhibits mode growth; the establishement of a steady inflow coincides in these two models with the capping of mode growth at a very low level, insufficient to disturb the disk. Model E3p would seem to straddle these two extremes in that it exhibits both early inflow and the development of unstable modes. Model E3p, like models E3r and X3p, is a marginal torus. E3p and X3p are very similar in their initial conditions (see Table 2), although model E3p has an initial inner edge that is slightly farther out that model X3p. In spite of these structural similarities, model E3p develops a significant PPI mode while E3r and X3p do not. As marginal tori, all of these three models show an early inflow of matter. However, the mass influx rate is less for model E3p (about $`0.01`$ in units normalized to the peak density) than it is for the other marginal models (between $`0.04`$ and $`0.08`$). Model E3p appears to be a “transition” model, where the presence of (weak) early inflow is insufficient to prevent mode growth. Mode growth eventually begins, and progresses much more slowly than in the other models that yield planets; the growth rate of the $`m=1`$ mode for model E3p is one third the rate of the next slowest mode growth, for model B3p. Blaes (1987) found that although the PPI modes can be stabilized by accretion, there is a finite transition that permits both accretion and mode growth. E3p appears to be such a case. These results for 3D tori in the Kerr metric reinforce the idea that accretion through the inner boundary at the equatorial plane suppresses the growth of the PPI, and that this suppression is as effective in three dimensions as it is in two.

Perhaps most interestingly, we have seen evidence of frame-dragging in those models that develop an accretion flow. Models B3r and E3r show evidence of a retrograde inflow changing its winding sense as it enters the innermost regions near the static limit, and model B3p shows a prograde accretion flow getting smeared into a ring in the vicinity of the static limit. As has been suspected from the earliest days of numerical simulations of GR hydrodynamics, it is in the deepest regions of the black hole potential well that we find the most interesting, and perhaps most astrophysically distinct consequences of general relativity.

This work was supported by NSF grant AST-0070979 and NASA grant NAG5-9266. The simulations were carried out on the Origin 2000 system at NCSA, and the Bluehorizon system of NPACI.


format: text

text:
# 1 Introduction

## 1 Introduction

Let $`\pi S_n`$ be a permutation of $`\{1,2,\mathrm{},n\}`$ written in one-line notation. Let $`\alpha S_m`$. We say that $`\pi `$ contains the pattern $`\alpha `$ if there exist indices $`i_1,i_2,\mathrm{},i_m`$ such that $`\pi _{i_1}\pi _{i_2}\mathrm{}\pi _{i_m}`$ is equivalent to $`\alpha `$, where we define equivalence as follows. Define $`\overline{\pi }_{i_j}=|\{x:\pi _{i_x}\pi _{i_j},1xm\}|`$. If $`\alpha =\overline{\pi }_{i_1}\overline{\pi }_{i_2}\mathrm{}\overline{\pi }_{i_m}`$ then we say that $`\alpha `$ and $`\pi _{i_1}\pi _{i_2}\mathrm{}\pi _{i_m}`$ are equivalent.

Define $`S_n^k(T)`$, $`TS_m`$, to be the set of permutations with $`k`$ fixed points which avoid all pattern in $`T`$.

In recent paper \[RSZ\], the authors began the study of $`S_n^k(T)`$ for $`TS_3`$, $`|T|=1`$. In this paper we enumerate $`S_n^k(T)`$ for all $`TS_3`$ with $`|T|2`$. Thus, all that remains to enumerate are $`S_n^k(123)`$ and $`S_n^k(231)`$ for all $`0kn`$.

We start by giving two bijections which preserve the parameter “number of fixed points.” Let $`I:S_nS_n`$ be the group theoretical inverse. Let $`R:S_nS_n`$ be the reversal bijection, i.e. $`R(\pi _1\pi _2\mathrm{}\pi _n)=\pi _n\pi _{n1}\mathrm{}\pi _1`$, and let $`C:S_nS_n`$ be the complement bijection, i.e. $`C(\pi _1\pi _2\mathrm{}\pi _n)=(n+1\pi _1)(n+1\pi _2)\mathrm{}(n+1\pi _n)`$.

It is easy to see that the bijections $`I`$ and $`RC`$ preserve the number of fixed points in a permutation. This reduces the number of cases that we have to consider.

In the sequel, we will be using the following notation. Let $`I_S`$ be the characteristic function, i.e. $`I_S=1`$ if $`S`$ is true and $`I_S=0`$ if $`S`$ is false. We also define $`s_n^j(T)=0`$ for $`j<0`$ or $`j>n`$ for any $`TS_3`$.

## 2 The Case $`|T|=2`$

Using $`I`$ and $`RC`$ we see that we have the following cases.

(1) $`\overline{\{123,321\}}=\{\{123,321\}\}`$

(2) $`\overline{\{123,132\}}=\{\{123,132\},\{123,213\}\}`$

(3) $`\overline{\{123,231\}}=\{\{123,231\},\{123,312\}\}`$

(4) $`\overline{\{132,213\}}=\{\{132,213\}\}`$

(5) $`\overline{\{132,231\}}=\{\{132,231\},\{132,312\}\}`$

(6) $`\overline{\{132,321\}}=\{\{132,321\},\{213,321\}\}`$

(7) $`\overline{\{213,231\}}=\{\{213,231\},\{213,312\}\}`$

(8) $`\overline{\{231,312\}}=\{\{231,312\}\}`$

(9) $`\overline{\{231,321\}}=\{\{231,321\},\{312,321\}\}`$

###### Theorem 2.1

$`\{s_n^0(123,321)\}_{n0}=1,0,1,2,4,0,0,`$, $`\{s_n^1(123,321)\}_{n0}=0,1,0,2,0,0,`$, $`\{s_n^2(123,321)\}_{n0}=0,0,1,0,0`$, and $`s_n^k(123,321)=0`$ for all $`3kn`$.

Proof. Obvious. $`\mathrm{}`$

###### Theorem 2.2

For $`n1`$, $`s_n^k(123,132)=0`$ for $`3kn`$,

$$\begin{array}{cc}s_{2n+i}^2(123,132)=\{\begin{array}{cc}\frac{4^{n1}+2}{3}\hfill & \mathrm{𝑖𝑓}i=0\hfill \\ & \\ 0\hfill & \mathrm{𝑖𝑓}i=1\hfill \end{array},\hfill & s_{2n+i}^1(123,132)=\{\begin{array}{cc}\frac{2(4^{n1}1)}{3}\hfill & \mathrm{𝑖𝑓}i=0\hfill \\ & \\ \frac{4^n+2}{3}\hfill & \mathrm{𝑖𝑓}i=1\hfill \end{array},and\\ & \\ s_{2n+i}^0(123,132)=\{\begin{array}{cc}4^{n1}\hfill & \mathrm{𝑖𝑓}i=0\hfill \\ & \\ \frac{2(4^n1)}{3}\hfill & \mathrm{𝑖𝑓}i=1\hfill \end{array}.\hfill & \end{array}$$

Proof. Let $`\pi S_n(123,132)`$. It is easy to prove by induction that we must have $`\pi =(\pi (1),\pi (2),\mathrm{},\pi (m))`$ where $`\pi (j)=(t_{j1}1,t_{j1}2,\mathrm{},t_j+1,t_{j1})`$ for some $`n=t_0>t_1>t_2>\mathrm{}>t_m=0`$. We call $`\pi (j)`$ a $`j`$-block of $`\pi `$. (See \[Man\] for more details.)

Cleary $`s_n^k(123,132)=0`$ for $`3kn`$ since three fixed points yield a $`123`$ pattern. Hence, we need only consider $`k=0,1,2`$. However, using a result in \[SS\] and the fact that $`s_n^k(123,132)=0`$ for $`3kn`$, we immediately obtain the stated equation for $`s_n^0(123,231)`$ upon proving the formulas for $`k=1,2`$. Hence, we need only consider $`k=1,2`$.

We start with $`k=2`$. Consider the graph of $`\pi `$, which consists of the points $`(i,\pi _i)^2`$. By the symmetries of the graph of $`\pi `$, we see that the number of intersections with the line $`y=x`$ (i.e. the number of fixed points of $`\pi `$) is two if and only if $`m`$ is odd for $`\pi =(\pi (1),\pi (2),\mathrm{},\pi (m))`$ and the intersections are in the $`\frac{m1}{2}`$-block of $`\pi `$. So, the number of elements in the $`\frac{m1}{2}`$-block must be even, as well as the total number of elements of $`\pi `$. Hence, we have $`s_n^2(123,132)=0`$ for $`n`$ odd. For $`n`$ even, let $`2d`$ be the number of elements in the $`\frac{m1}{2}`$-block. Since the other blocks cannot contain a fixed point, we have $`(s_{nd}(123,132))^2`$ permutations for $`d=1,2,\mathrm{},n`$. Using a result in \[SS\], this gives $`1+_{d=1}^{n1}(2^{nd1})^2`$ permutations. Hence, $`s_{2n}^2(123,132)=\frac{4^{n1}+2}{3}`$.

Next, we have $`k=1`$ and again we consider the graph of $`\pi `$.

For $`n`$ odd, we have a situation similar to the $`k=2`$ case. We must have the middle block containing an odd number of elements. Let this middle block have $`2d+1`$ elements. Using a result in \[SS\], this gives $`1+_{d=1}^{n1}(2^{nd1})^2`$ permutations. Hence, $`s_{2n+1}^1(123,132)=\frac{4^{n1}+2}{3}`$.

For $`n`$ even, we have three possibilities.

1. $`\pi (m)`$ contains $`2d_1`$ elements, for some $`2d_1n1`$, with one fixed point, or

2. $`\pi (1)`$ contain $`2d_2`$ of elements, for some $`2d_2n1`$, with one fixed point, or

3. There exists exactly one block with $`2d_3+1`$ elements, for some $`2d_3n1`$, containing one fixed point.

Letting $`N=2n`$, using a result in \[SS\], and summing over these three cases we get

$$\underset{d_1=2}{\overset{n1}{}}2^{N2d_11}+\underset{d_2=2}{\overset{n1}{}}2^{N2d_21}+\underset{d_3=1}{\overset{n1}{}}(2^{nd_31})^2$$

permutations. This gives $`s_{2n}^1(123,132)=\frac{2(4^{n1}1)}{3}`$, thereby completing the proof. $`\mathrm{}`$

###### Theorem 2.3

For $`n2`$, $`s_n^k(123,231)=0`$ for $`3kn`$,

$$\begin{array}{c}s_n^2(123,231)=\{\begin{array}{cc}\frac{n(n6)}{24}+\frac{n}{2}\hfill & \mathrm{𝑖𝑓}n0(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n1)(n+1)}{24}\hfill & \mathrm{𝑖𝑓}n1,5(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n4)(n2)}{24}+\frac{n}{2}\hfill & \mathrm{𝑖𝑓}n2,4(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n3)(n+3)}{24}\hfill & \mathrm{𝑖𝑓}n3(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \end{array},\hfill \\ \\ s_n^1(123,231)=\{\begin{array}{cc}\frac{n(n6)}{12}+6\left(\genfrac{}{}{0pt}{}{\frac{n+6}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n0(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n3)(n1)}{8}+\frac{(n7)(n1)}{12}+6\left(\genfrac{}{}{0pt}{}{\frac{n+5}{6}}{2}\right)+\frac{n+2}{3}\hfill & \mathrm{𝑖𝑓}n1(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{n(n2)}{12}+6\left(\genfrac{}{}{0pt}{}{\frac{n+4}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n2(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n3)(n1)}{8}+\frac{(n5)(n3)}{12}+6\left(\genfrac{}{}{0pt}{}{\frac{n+3}{6}}{2}\right)+\frac{2n+3}{3}\hfill & \mathrm{𝑖𝑓}n3(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n12)(n+2)}{12}+6\left(\genfrac{}{}{0pt}{}{\frac{n+8}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n4(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n3)(n1)}{8}+\frac{(n5)(n3)}{12}+6\left(\genfrac{}{}{0pt}{}{\frac{n+1}{6}}{2}\right)+n\hfill & \mathrm{𝑖𝑓}n5(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \end{array},and\hfill \end{array}$$

and $`s_n^0(123,231)=\left(\genfrac{}{}{0pt}{}{n}{2}\right)+1s_n^1(123,231)s_n^2(123,231)`$.

Proof. Cleary $`s_n^k(123,231)=0`$ for $`3kn`$ since three fixed points yield a $`123`$ pattern. Hence, we need only consider $`k=0,1,2`$. However, using a result in \[SS\] and the fact that $`s_n^k(123,231)=0`$ for $`3kn`$, we immediately obtain the stated equation for $`s_n^0(123,231)`$ upon proving the formulas for $`k=1,2`$. Hence, we need only consider $`k=1,2`$.

Let $`\pi S_n(123,231)`$ with $`\pi _i=n`$. If $`i1`$ then we must have $`\pi =(i1)(i2)\mathrm{}21n(n1)(n2)\mathrm{}i`$ for any $`2in`$. This case contributes $`n1`$ permutations to $`S_n^1(123,231)`$ for $`n`$ odd and $`\frac{n}{2}`$ to $`S_n^2(123,231)`$ for $`n`$ even. For $`i=1`$, we have $`\pi =n(n1)\mathrm{}(nx+1)y(y1)\mathrm{}1(nx)(nx1)\mathrm{}(y+1)`$ for some $`1xn2`$ and $`1ynx1`$, or we have $`\pi =n(n1)\mathrm{}21`$.

We now consider the cases $`k=1,2`$ for $`i=1`$.

We have three potential intervals in which fixed points may reside:

(A) $`n(n1)\mathrm{}(nx+1)`$,

(B) $`y(y1)\mathrm{}1`$, or

(C) $`(nx)(nx1)\mathrm{}(y+1)`$

Hence, we have the following situations in which a fixed point appears.

A. For $`n`$ odd we get one fixed point for any $`\frac{n+1}{2}xn2`$ and $`1ynx1`$, or if $`\pi =n(n1)\mathrm{}1`$. The fixed point occurs at position $`\frac{n+1}{2}`$.

B. For $`x+y`$ odd we get one fixed point for any $`1x\frac{n2}{2}`$ and $`x+1ynx1`$. The fixed point occurs at position $`\frac{x+y+1}{2}`$.

C. For $`n+y`$ odd we get one fixed point for any $`1x\frac{n2}{2}`$ and $`1yn2x1`$. The fixed point occurs at position $`\frac{n+y+1}{2}`$.

We now consider the case $`k=2`$ ($`i=1`$). To have two distinct fixed points we must have both $`B`$ and $`C`$ both occuring. Hence, we get

$$\underset{x=1}{\overset{\frac{n2}{3}}{}}\underset{\genfrac{}{}{0pt}{}{\genfrac{}{}{0pt}{}{y=x+1}{n+yodd}}{x+yodd}}{\overset{n2x1}{}}1=\{\begin{array}{cc}\frac{n(n6)}{24}\hfill & \mathrm{𝑖𝑓}n0(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n1)(n+1)}{24}\hfill & \mathrm{𝑖𝑓}n1,5(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n4)(n2)}{24}\hfill & \mathrm{𝑖𝑓}n2,4(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n3)(n+3)}{24}\hfill & \mathrm{𝑖𝑓}n3(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \end{array}$$

permutations with two fixed points in this case.

We now consider the case $`k=1`$ ($`i=1`$). We may have only A occur, only B occur, or only C occur:

Only A occurs: In this case we require $`n`$ to be odd and we get

$$1+\underset{x=\frac{n+1}{2}}{\overset{n2}{}}\underset{y=1}{\overset{nx1}{}}1=\frac{(n3)(n1)}{8}+1$$

permutations with one fixed point.

Only B occurs: In this case we get

$$\underset{x=1}{\overset{\frac{n1}{3}}{}}\underset{\genfrac{}{}{0pt}{}{y=n2x}{x+yodd}}{\overset{nx1}{}}1+\underset{x=\frac{n1}{3}+1}{\overset{\frac{n2}{2}}{}}\underset{\genfrac{}{}{0pt}{}{y=x}{x+yodd}}{\overset{nx1}{}}1+\underset{\genfrac{}{}{0pt}{}{x=1}{xodd}}{\overset{\frac{n2}{3}}{}}\underset{\genfrac{}{}{0pt}{}{y=x+1}{yeven}}{\overset{n2x1}{}}1=\{\begin{array}{cc}\frac{n(n6)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+6}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n0(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n15)(n1)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+5}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n1(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{n(n2)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+4}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n2(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n9)(n3)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+3}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n3(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n12)(n+2)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+8}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n4(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n5)(n3)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+1}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n5(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \end{array}$$

permutations with one fixed point.

Only C occurs: In this case we get

$$\underset{x=1}{\overset{\frac{n1}{3}}{}}\underset{\genfrac{}{}{0pt}{}{y=1}{n+yodd}}{\overset{x}{}}1+\underset{x=\frac{n1}{3}+1}{\overset{\frac{n2}{2}}{}}\underset{\genfrac{}{}{0pt}{}{\genfrac{}{}{0pt}{}{y=1}{n+yodd}}{x+yeven}}{\overset{n2x1}{}}1+\underset{x=1}{\overset{\frac{n2}{3}}{}}\underset{\genfrac{}{}{0pt}{}{y=x+1}{n+yodd}}{\overset{n2x1}{}}1=\{\begin{array}{cc}\frac{n(n6)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+6}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n0(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n15)(n1)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+5}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n1(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{n(n2)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+4}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n2(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n9)(n3)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+3}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n3(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n12)(n+2)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+8}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n4(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \\ & \\ \frac{(n5)(n3)}{24}+3\left(\genfrac{}{}{0pt}{}{\frac{n+1}{6}}{2}\right)\hfill & \mathrm{𝑖𝑓}n5(\mathrm{𝑚𝑜𝑑}\mathrm{\hspace{0.17em}\hspace{0.17em}6})\hfill \end{array}$$

permutations with one fixed point.

Summing over all cases yields the stated formula. $`\mathrm{}`$

###### Theorem 2.4

For $`n1`$, $`s_n^{n1}(213,132)=0`$, $`s_n^n(213,132)=1`$,

$$s_{2n+i}^0(213,132)=\{\begin{array}{cc}\frac{54^{n1}2}{3}\hfill & \mathrm{𝑖𝑓}i=0\hfill \\ & \\ \frac{2(4^n1)}{3}\hfill & \mathrm{𝑖𝑓}i=1\hfill \end{array},$$

and, for $`1kn1`$,

$$\begin{array}{c}s_{2n+i}^{2k}(213,132)=\{\begin{array}{cc}4^{nk1}\hfill & \mathrm{𝑖𝑓}i=0\hfill \\ & \\ 0\hfill & \mathrm{𝑖𝑓}i=1\hfill \end{array},and\hfill \\ \\ s_{2n+i}^{2k+1}(213,132)=\{\begin{array}{cc}0\hfill & \mathrm{𝑖𝑓}i=0\hfill \\ & \\ 4^{nk1}\hfill & \mathrm{𝑖𝑓}i=1\hfill \end{array}.\hfill \end{array}$$

Proof. The proof is very similar to the proof of Theorem 2.2 and we provide a sketch of this proof.

It is easy to prove by induction that $`\pi =(\pi (1),\pi (2),\mathrm{},\pi (m))`$ where $`\pi (j)=(t_j+1,t_j+2,\mathrm{},t_{j1}1,t_{j1})`$ for some $`0=t_0<t_1<t_2<\mathrm{}<t_m=n`$. We call $`\pi (j)`$ a $`j`$-block of $`\pi `$. (See \[Man\] for more details.)

Recall that the graph of $`\pi `$ consists of the points $`(i,\pi _i)^2`$.

For $`2k`$ fixed points, we see that the middle block of $`\pi =(\pi (1),\mathrm{},\pi (m))`$ must be of length $`2k`$ and must contain all $`2k`$ fixed points and that $`n`$ must be even. Hence, we have $`s_{2n}^{2k}=(s_{nk}(213,132))^2=4^{nk1}`$ and $`s_{2n+1}^{2k}=0`$ for $`1kn1`$.

Similarly, for $`2k+1`$ fixed points, we have $`s_{2n+1}^{2k+1}=(s_{nk}(213,132))^2=4^{nk1}`$ and $`s_{2n}^{2k+1}=0`$ for $`1kn1`$.

Using a result in \[SS\] gives the stated formula for $`s_n^0(213,132)`$. $`\mathrm{}`$

###### Theorem 2.5

For $`n3`$, $`s_n^0(132,231)=\frac{1}{3}(2^{n1}+(1)^n)`$, $`s_n^k(132,231)=\frac{2}{3}(2^{nk}+(1)^{nk+1})`$ for $`1kn2`$, $`s_n^{n1}(132,231)=0`$, and $`s_n^n(132,231)=1`$.

Proof. Let $`\pi S_n(132,231)`$. Clearly, we must have either $`\pi _1=n`$ or $`\pi _n=n`$. For the case $`\pi _n=n`$ we have $`s_{n1}^{k1}(132,231)`$ permutations. Hence, we need only consider $`\pi _1=n`$. Write $`\pi =(n,\pi (1),1,\pi (2))`$. We must have $`\pi (1)`$ be a decreasing sequence and $`\pi (2)`$ be an increasing sequence. Hence, we cannot have a fixed point in $`\pi (2)`$ and we can have at most one fixed point in $`\pi (1)`$. Thus, for $`k2`$ the case $`\pi _1=n`$ is empty giving $`s_n^k(132,231)=s_{n1}^{k1}(132,231)`$ for $`k2`$.

It remains to consider $`k=0,1`$ for $`\pi =(n,\pi (1),1,\pi (2))`$. We consider $`k=0`$ first (the case $`k=1`$ will follow immediately). We may place the elements $`2,3,\mathrm{},n1`$ in $`2^{n2}`$ ways such that $`\pi (1)`$ is decreasing and $`\pi (2)`$ is increasing. To see this, insert the elements $`2,3,\mathrm{},n1`$ in reverse order at $`1`$’s position and slide the element all the way to the left or all the way to right. However, we must subtract those permutations which contain a fixed point, which may only occur in $`\pi (1)`$.

Let $`f`$ be a fixed point and write $`\pi (1)=(\sigma (1),f,\sigma (2))`$. We have $`\left(\genfrac{}{}{0pt}{}{nf1}{f2}\right)`$ choices for the elements in $`\sigma (1)`$. Once these are chosen, there are then $`2^{f2}`$ ways to place the elements $`2,3,\mathrm{},f1`$ into $`\sigma (2)`$ and $`\pi (2)`$ such that $`\sigma (2)`$ is decreasing and $`\pi (2)`$ is increasing. Thus, we have

$$s_n^0(132,231)=2^{n2}\underset{f=2}{\overset{\frac{n}{2}}{}}\left(\genfrac{}{}{0pt}{}{nf1}{f2}\right)2^{f1}.$$

From here, it is easy to check that $`s_n^0(132,231)=s_{n1}^0(132,231)+2s_{n2}^0(132,231)`$, which gives $`s_n^0(132,231)=\frac{1}{3}(2^{n1}+(1)^n)`$ for $`s_1^0(132,231)=0`$ and $`s_2^0(132,231)=1`$.

For the case $`k=1`$, the above argument shows that

$$s_n^1(132,231)=\underset{f=2}{\overset{\frac{n}{2}}{}}\left(\genfrac{}{}{0pt}{}{nf1}{f2}\right)2^{f1}.$$

From here, it is easy to check that $`s_n^1(132,231)=s_{n1}^1(132,231)+2s_{n2}^1(132,231)`$, which gives $`s_n^1(132,231)=\frac{2}{3}(2^{n1}+(1)^n)`$ for $`s_1^1(132,231)=1`$ and $`s_2^1(132,231)=0`$. $`\mathrm{}`$

Remark. $`s_n^0(132,231)=J_{n2}`$ where $`J_n`$ is the $`n^{\mathrm{𝑡ℎ}}`$ Jacobsthal number.

###### Theorem 2.6

For $`n1`$, $`s_n^k(132,321)=nk1`$ for $`0kn1`$ and $`s_n^n(132,321)=1`$.

Proof. Let $`\pi S_n(132,321)`$. It is easy to see that in order to avoid both $`132`$ and $`321`$ we have either $`\pi =(\pi ^{},n)`$ where $`\pi ^{}S_{n1}(132,321)`$ or $`\pi =(j+1)(j+2)\mathrm{}(n1)n12\mathrm{}j`$ for some $`1jn1`$.

Adjusting for the number of fixed points we see that

$$s_n^k(132,321)=s_{n1}^{k1}(132,321)+\underset{j=1}{\overset{n1}{}}I_{k=0}.$$

A straightforward induction on $`k`$ for $`0k<n`$ finishes the proof. $`\mathrm{}`$

As a corollary of Theorem 2.6 we rederive a result given in \[SS\].

###### Corollary 2.7

$`s_n(132,321)=\left(\genfrac{}{}{0pt}{}{n}{2}\right)+1`$

Proof. Summing over $`0kn`$ we obtain the result immediately.

###### Theorem 2.8

For $`n3`$, $`s_n^0(213,231)=\frac{1}{3}(2^{n1}+(1)^n)`$, $`s_n^k(213,231)=\frac{2}{3}(2^{nk}+(1)^{nk+1})`$ for $`1kn2`$, $`s_n^{n1}(213,231)=0`$, and $`s_n^n(213,231)=1`$.

Proof. Using the fact that $`s_n^k(213,231)=s_n^k(213,312)`$ and using the patterns $`213`$ and $`312`$, the proof of this theorem is nearly identical to the proof of Theorem 2.5 and is left to the reader. $`\mathrm{}`$

We now find the need for the following definition.

###### Definition 2.1

Let $`S_1,S_2S_m`$. If $`s_n^k(S_1)=s_n^k(S_2)`$ for all $`0kn`$ for all $`nm`$ we say that $`S_1`$ and $`S_2`$ are Super-Wilf equivalent.

Remark. Using the above definition, Theorems 1.5 and 1.8, and results from \[RSZ\] we see that $`321`$, $`132`$, and $`213`$ are Super-Wilf equivalent, $`231`$ and $`312`$ are Super-Wilf equivalent, and $`\{132,231\}`$, $`\{132,312\}`$, $`\{213,231\}`$ and $`\{213,312\}`$ are Super-Wilf equivalent.

###### Theorem 2.9

For $`n1`$,

$$s_n^k(231,312)=\{\begin{array}{cc}2^{\frac{nk2}{2}}\left(\left(\genfrac{}{}{0pt}{}{\frac{n+k}{2}}{\frac{nk}{2}}\right)+\left(\genfrac{}{}{0pt}{}{\frac{n+k2}{2}}{\frac{nk}{2}}\right)\right)\hfill & \mathrm{𝑓𝑜𝑟}n+k\mathrm{𝑒𝑣𝑒𝑛}\hfill \\ & \\ 0\hfill & \mathrm{𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒}\hfill \end{array}$$

Proof. For $`\pi S_n(231,312)`$ it is easy to see that we must have $`\pi =\pi (1)n(n1)\mathrm{}j`$ with $`\pi (1)S_{j1}^k(231,312)`$ for some $`1jn`$. Hence, we have

$$s_n^k(231,312)=\underset{\genfrac{}{}{0pt}{}{j=1}{njodd}}{\overset{n1}{}}s_{j1}^k(231,312)+\underset{\genfrac{}{}{0pt}{}{j=1}{\genfrac{}{}{0pt}{}{njeven}{orn=j}}}{\overset{n}{}}s_{j1}^{k1}(231,312).$$

From here we see that $`s_n^k(231,312)=2s_{n2}^k(231,312)+s_{n1}^{k1}(231,312)`$. Hence, using initial conditions it is easy to check that the formula given is correct. $`\mathrm{}`$

###### Theorem 2.10

Let $`G_k(x)`$ be the generating function for $`\{s_n^k(231,321)\}_{n0}`$. Then $`G_k(x)=\frac{x^k(1x)^{k+1}}{(1xx^2)^{k+1}}`$. In particular, $`s_n^0(231,321)=F_{n2}`$, for $`n2`$, where $`\{F_n\}_{n0}`$ is the Fibonacci sequence initialized by $`F_0=F_1=1`$.

Proof. Let $`\pi S_n(231,321)`$. It is easy to see that we must have $`\pi =(\pi ^{},\sigma )`$ where $`\pi S_{nj}^{kI_{j=1}}(231,321)`$ and $`\sigma =n(nj+1)(nj+2)\mathrm{}(n1)`$ for some $`1jn`$.

Hence, we have $`s_n^k(231,321)=s_{n1}^{k1}(231,321)+_{j=2}^ns_{nj}^k(231,321)`$ for $`0kn`$. This implies that $`s_n^k(231,321)=s_{n1}^{k1}(231,321)+s_{n2}^k(231,321)+s_{n1}^k(231,321)s_{n2}^{k1}(231,321)`$.

For $`k=0`$ we have $`s_n^0(231,321)=s_{n1}^0(231,321)+s_{n2}^0(231,321)`$, from which it follows that $`G_0(x)=\frac{1x}{1xx^2}`$. For $`k>0`$ we have $`G_k(x)=\frac{x(1x)}{1xx^2}G_{k1}(x)`$. From here all of the results follow. $`\mathrm{}`$

The following corollary rederives a result in \[SS\] .

###### Corollary 2.11

$`_{n0}s_n(231,321)x^n=\frac{1x}{12x}`$, i.e. $`s_n(231,321)=2^{n1}`$ for $`n1`$.

Proof. Summing $`_{k0}G_k(x)`$ gives the desired result. $`\mathrm{}`$

## 3 The Case $`|T|=3`$

Using $`I`$ and $`RC`$ we see that we have the following cases, except that case (1) is grouped together because both $`123`$ and $`321`$ are to be avoided, which is not possible for $`n5`$.

(1) $`\overline{\{123,132,321\}}=\{\{123,132,321\},\{123,213,321\},\{123,231,321\},\{123,312,321\}\}`$

(2) $`\overline{\{123,132,213\}}=\{\{123,132,213\}\}`$

(3) $`\overline{\{123,132,231\}}=\{\{123,132,231\},\{123,132,312\},\{123,213,231\},\{123,213,312\}\}`$

(4) $`\overline{\{123,231,312\}}=\{\{123,231,312\}\}`$

(5) $`\overline{\{132,213,231\}}=\{\{132,213,231\},\{132,213,312\}\}`$

(6) $`\overline{\{132,213,321\}}=\{\{132,213,321\}\}`$

(7) $`\overline{\{132,231,312\}}=\{\{132,231,312\},\{213,231,312\}\}`$

(8) $`\overline{\{132,231,321\}}=\{\{132,231,321\},\{132,312,321\},\{213,231,321\},\{213,312,321\}\}`$

(9) $`\overline{\{231,312,321\}}=\{\{231,312,321\}\}`$

###### Theorem 3.1

Let $`\alpha \{132,213,231,312\}`$. Then

$$\begin{array}{c}\{s_n^0(123,321)\}_{n0}=1,0,1,2,1,0,0,\mathrm{}\mathrm{𝑓𝑜𝑟}\alpha 231,312\hfill \\ \{s_n^0(123,321)\}_{n0}=1,0,1,1,1,0,0,\mathrm{}\mathrm{𝑓𝑜𝑟}\alpha =231\mathrm{𝑜𝑟}\alpha =312\hfill \\ \{s_n^1(123,\alpha ,321)\}_{n0}=0,1,0,2,0,0,\mathrm{}\mathrm{𝑓𝑜𝑟}\alpha 132,213\hfill \\ \{s_n^1(123,\alpha ,321)\}_{n0}=0,1,0,1,0,0,\mathrm{}\mathrm{𝑓𝑜𝑟}\alpha =132\mathrm{𝑜𝑟}\alpha =213\hfill \\ \{s_n^2(123,\alpha ,321)\}_{n0}=0,0,1,0,0\mathrm{},\mathrm{𝑓𝑜𝑟}\mathrm{𝑎𝑛𝑦}\alpha ,\mathrm{𝑎𝑛𝑑}\hfill \end{array}$$

$`s_n^k(123,\alpha ,321)=0`$ for any $`\alpha `$, for all $`3kn`$.

Proof. Obvious. $`\mathrm{}`$

###### Theorem 3.2

For $`n3`$, $`s_n^k(123,132,213)=0`$ for $`3kn`$, $`s_n^2(123,132,213)=F_{\frac{n2}{2}}^2I_{neven}`$, $`s_n^1(123,132,213)=F_{\frac{n1}{2}}^2I_{nodd}`$, and

$$s_n^0(123,132,213)=\{\begin{array}{cc}F_nF_{\frac{n2}{2}}^2\hfill & \mathrm{𝑖𝑓}n\mathrm{𝑖𝑠}\mathrm{𝑒𝑣𝑒𝑛}\hfill \\ & \\ F_nF_{\frac{n1}{2}}^2\hfill & \mathrm{𝑖𝑓}n\mathrm{𝑖𝑠}\mathrm{𝑒𝑣𝑒𝑛}\hfill \end{array},$$

where $`\{F_n\}_{n0}`$ is the Fibonacci sequence intialized by $`F_0=F_1=1`$.

Proof. Clearly, $`s_n^k(123,132,213)=0`$ for $`k3`$ since three fixed points are a $`123`$ pattern. We may then use a result from from \[SS\] to get $`s_n^0(123,132,213)=s_n(123,132,213)s_n^1(123,132,213)s_n^2(123,132,213)=F_ns_n^1(123,132,213)s_n^2(123,132,213)`$. Hence, we only need consider $`k=1,2`$.

We start with $`k=2`$. The fixed points must be adjacent in order to avoid all of $`123,132,213`$. Let the fixed points be $`j`$ and $`j+1`$. Write $`\pi =(\pi (1),j,j+1,\pi (2))`$. Then $`\pi (1)`$ is on the elements $`\{nj+2,\mathrm{},n\}`$. To avoid the $`123`$ pattern, we require $`nj+2j+2`$ and $`nj+1j+1`$. Hence, $`j=\frac{n}{2}`$ (so that $`n`$ must be even). Thus, since it is not possible to have fixed points in $`\pi (1)`$ or $`\pi (2)`$ we have $`s_n^2(123,132,213)=(s_{\frac{n2}{2}}(123,132,213))^2I_{neven}`$, which, using a result from \[SS\], gives the stated formula.

Now consider $`k=1`$. Let $`j`$ be the fixed point. Note that $`j1,n`$ in order to avoid all of $`123,132,213`$. Write $`\pi =(\pi (1),j,\pi (2))`$. In order to avoid $`123`$ and $`213`$, there exists at most one $`x\pi (1)`$ with $`x<j`$. Assume, for a contradiction, that such an $`x`$ exists. Then we must have $`x=\pi _{j1}`$ in order to avoid the $`132`$ pattern. Next, for $`z\pi (2)`$ we must have $`z<j`$ else $`xjz`$ is a $`123`$ pattern. This forces $`x=j1`$ to be a fixed point, a contradiction. Hence, no such $`x`$ exists so that for all $`i\pi (1)`$ we have $`i>j`$. In order to avoid $`123`$ and $`132`$, there exists at most one $`y\pi (2)`$ with $`y>j`$. Assume, for a contradiction, that such a $`y`$ exists. Then we must have $`y=\pi _{j+1}`$ in order to avoid the $`213`$ pattern. Furthermore, for all $`i\pi (1)`$ we must have $`i>j`$ so as not to have a $`213`$ pattern with $`ijy`$. This forces $`y=j+1`$ to be a fixed point, a contradiction. Hence, no such $`y`$ exists. Since we have $`i\pi (1)`$, $`i>j`$ and $`k\pi (2)`$, $`k<j`$, and for $`j`$ to be a fixed point we require $`n`$ to be odd so that $`j=\frac{n+1}{2}`$, we get $`s_n^1(123,132,213)=(s_{\frac{n1}{2}}(123,132,213))^2I_{nodd}`$, which, using a result from \[SS\], gives the stated formula. $`\mathrm{}`$

###### Theorem 3.3

For $`n3`$, $`s_n^0(123,132,231)=\frac{n}{2}`$, $`s_n^1(123,132,231)=\frac{n}{2}+(1)^{n+1}`$, $`s_n^2(123,132,231)=\frac{1}{2}(1+(1)^n)`$, and $`s_n^k(123,132,231)=0`$ for $`3kn`$.

Proof. Let $`\pi S_n(123,132,231)`$. It is easy to see that we must have $`\pi =n(n1)\mathrm{}(nj+1)(nj1)(nj2)\mathrm{}21(nj)`$ for some $`0jn1`$. From here the results follow easily. $`\mathrm{}`$

###### Theorem 3.4

For $`n3`$, $`s_n^0(123,231,312)=s_n^2(123,231,312)=\frac{n}{2}(1I_{nodd})`$,
$`s_n^1(123,231,312)=n(1I_{neven})`$, and $`s_n^k(123,231,312)=0`$ for $`3kn`$.

Proof. Let $`\pi S_n(123,231,312)`$. It is easy to see that we must have $`\pi =j(j1)\mathrm{}1n(n1)\mathrm{}(j+1)`$. From here the results follow easily. $`\mathrm{}`$

###### Theorem 3.5

For $`n3`$, $`s_n^0(132,213,231)=\frac{n}{2}+(\frac{n}{2}+1)I_{neven}`$, $`s_n^1(132,213,231)=\frac{n}{2}I_{nodd}`$, and $`s_n^k(132,213,231)=I_{n=k}`$ for $`2kn`$.

Proof. Let $`\pi S_n(132,213,231)`$. It is easy to see that we must have $`\pi =n(n1)\mathrm{}(nj+1)12\mathrm{}(nj)`$ for some $`1jn`$. From here the results follow easily. $`\mathrm{}`$

###### Theorem 3.6

For $`n3`$, $`s_n^0(132,213,321)=n1`$ and $`s_n^k(132,213,321)=I_{n=k}`$ for $`1kn`$.

Proof. Let $`\pi S_n(132,213,321)`$. Then we must have $`\pi =j(j+1)\mathrm{}n12\mathrm{}(j1)`$ for some $`1jn`$. From here the results follow easily. $`\mathrm{}`$

###### Theorem 3.7

For $`n2`$, $`s_n^0(132,231,312)=\frac{1}{2}(1+(1)^n)`$,

$$s_n^{2k}(132,231,312)=\{\begin{array}{ccc}1+(1)^n\hfill & & \mathrm{𝑖𝑓}n>2k\hfill \\ 1\hfill & & \mathrm{𝑖𝑓}n=2k\hfill \end{array},and$$

$$s_n^{2k1}(132,231,312)=\{\begin{array}{ccc}1+(1)^{n+1}\hfill & & \mathrm{𝑖𝑓}n>2k1\hfill \\ 1\hfill & & \mathrm{𝑖𝑓}n=2k1\hfill \end{array}.$$

Proof. Let $`\pi S_n(132,231,312)`$. It is easy to see that we must have $`\pi =j(j1)\mathrm{}1(j+1)(j+2)\mathrm{}n`$ for some $`1jn`$. From here the results follow easily. $`\mathrm{}`$

###### Theorem 3.8

For $`n3`$,

$$s_n^k(132,231,321)=\{\begin{array}{ccc}1\hfill & & \mathrm{𝑖𝑓}\mathrm{\hspace{0.17em}\hspace{0.17em}0}kn2\hfill \\ 0\hfill & & \mathrm{𝑖𝑓}k=n1\hfill \\ 1\hfill & & \mathrm{𝑖𝑓}k=n\hfill \end{array}.$$

Proof. Let $`\pi S_n(132,231,321)`$. It is easy to see that we must have $`\pi =j12\mathrm{}(j1)(j+1)\mathrm{}n`$ for some $`1jn`$. From here the results follow easily. $`\mathrm{}`$

###### Theorem 3.9

For $`n3`$ and $`0kn`$, $`s_n^k(231,312,321)=\left(\genfrac{}{}{0pt}{}{\frac{n+k}{2}}{k}\right)I_{n+keven}`$.

Proof. Let $`\pi S_n^k(231,312,321)`$. Then $`\pi `$ must be of the form $`(1,\pi ^{})`$ or $`(2,1,\pi ^{\prime \prime })`$ where $`\pi ^{}S_{n1}^{k1}(231,312,321)`$ and $`\pi ^{\prime \prime }S_{n1}^k(231,312,321)`$. This gives us $`s_n^k(231,312,321)=s_{n1}^{k1}(231,312,321)+s_{n2}^k(231,312,321)`$. A straightforward induction on $`n+k`$ finishes the proof. $`\mathrm{}`$

## 4 The Cases $`|T|4`$

The cases $`|T|4`$ are easy and in fact $`s_n^k(T)\{0,1,2\}`$ for all $`TS_3`$, $`|T|4`$, for all $`n1`$ and $`0kn`$.

## References

\[Man\] T. Mansour, Permutations Avoiding a Pattern from $`S_k`$ and at Least Two Patterns from $`S_3`$, Ars Combinatoria 62 (2001),

\[RSZ\] A. Robertson, D. Saracino, and D. Zeilberger, Refined Restricted Permutations, arXiv: math.CO/0203033

\[SS\] R. Simion and F. Schmidt, Restricted Permutations, European Journal of Combinatorics 6 (1985), 383-406.