text:
# The Optical Emission from Gamma-Ray Quasars

## 1 Introduction

Active galactic nuclei (AGN) are powerful emitters in all wavebands. Emission processes from different parts of the spectrum in many cases appear to be linked. For example, the electrons in the jet which produce the synchrotron radiation may at the same time upscatter photons into the X-ray and gamma-ray regime via the inverse Compton (IC) process (see for example the analysis by Bloom & Marscher 1996). In order to unravel the contribution from each emission process, the shape of each distinct emission component needs to be determined, and the associated physical parameters measured.

The IC emission is produced by seed photons (internal, from the jet itself, or external, from the accretion disc or the emission-line clouds) which are upscattered by the electrons (or other relativistic particles) in the jet. The energy of the photons increases by a factor of $`\gamma ^2`$, where $`\gamma `$ is the Lorentz factor of the relativistic particle. Consequently, there are two crucial factors that determine the nature of the IC emission: the maximum energy of the jet electrons, which determines the turnover frequency of the synchrotron component; and the typical energy of the seed photons, which is governed by the relative strengths of the thermal and non-thermal components. The optical regime is where these components are both energetically important, and so a good understanding of the emission at optical frequencies is crucial for accurate modelling of the high-energy emission.

Before COMPTEL was shut down in June 2000, this telescope detected nine blazars in the range 0.75 – 30 MeV which were also detected by EGRET in the gamma-ray regime $``$ 100 MeV (Schönfelder et al. 2000; Collmar 2001). For many of the sources in this sample, the data coverage in the optical/near infrared (NIR) region is very poor. To better understand this important part of the spectral energy distribution (SED), six blazars from the COMPTEL sample, chosen subject to visibility during an observation period from 27 August to 2 September 2001, were observed in the photometric bands $`B`$, $`V`$, $`R`$, $`I`$, $`J`$, $`H`$, and $`K_n`$. Basic information for this sample is shown in Table 1. We note that the redshift for PKS 1622$``$297 ($`z=0.815`$) is taken from the Parkes Catalogue (Wright & Otrupcek 1990, hereafter PKSCAT90), wherein no reference is given. No other determination of its redshift has been made in the literature, and, while we use this value herein, it should be regarded as potentially suspect.

Detailed models of the multiwavelength spectra of several of the sources in this sample have been published in the literature. For example, Mukherjee et al. (1999) describe a model with both synchrotron and IC emission components fitted to nine epochs of monitoring data for the source PKS 0528$`+`$134. While considerable effort has been expended in fitting the high-energy spectrum (X-rays and above) of this source, the information in the low-energy part of the spectrum, from optical to radio energies, has largely been ignored. Rather than modelling one source in detail, by considering a modest sample of six sources more general conclusions can be drawn about the nature of the emission processes.

As synchrotron emission causes linear polarisation, optical polarisation measurements can be used to check the predictions from the modelling. Optical polarimetry of most of the sources in this sample is available from the literature. It will be shown that the sources with high optical polarisation are generally distinguished by their radio to optical spectral shapes, and are consistent with having a synchrotron component extending to optical wavelengths.

Section 2 describes the photometric optical and NIR data for the six sources. The SEDs in each case are adequately modelled by a power law in Section 3. Finally, in Section 4 we examine the implications of our observations for multiwavelength modelling of gamma-ray sources, as well as the demographics of the sample as a whole.

## 2 Observations and Data

### 2.1 Observations

The six blazars were observed in optical bands $`B`$, $`V`$, $`R`$, and $`I`$ on the night of 29 August 2001 with the ANU 40 in telescope at Siding Spring Observatory (SSO). The exposure time in each band was 300 seconds. The seeing was typically $`2.5{}_{}{}^{\prime \prime }3^{\prime \prime }`$, but conditions were generally photometric (although the weather did force some breaks in the observing, and the loss of the $`B`$ band data for PKS 1622$``$297). Cross-checking of field stars from Raiteri et al. (1998) enables verification of the photometric conditions. Observations of standard stars from Graham (1982) enabled calibration to an absolute magnitude/flux scale.

Observations in the NIR bands $`J`$, $`H`$, and $`K_n`$ were made on the nights of 31 August and 1 September 2001 with the CASPIR (Cryogenic Array Spectrometer/Imager—see McGregor et al. 1994) instrument on the ANU 2.3 m telescope, SSO. Each image is a median-combination of four (or, in the case of the fainter sources, 20) dithered images of 60 seconds exposure time each (these individual exposures were: 2 cycles of 30 seconds for $`J`$; 6 cycles of 10 seconds in $`H`$; and 12 cycles of 5 seconds in $`K_n`$).

Poor weather on the first night allowed us to observe only PKS 1622$``$297 and two standard stars, but the second night had good conditions, allowing us to observe the remaining sources. The seeing was better than for the optical (typically $`1.5{}_{}{}^{\prime \prime }2^{\prime \prime }`$) and conditions were generally photometric. Individual frames where the background increased due to the presence of cloud were discarded prior to co-addition. Observations of IRIS standard stars from the Carter system (Carter & Meadows 1995) were used for calibration.

### 2.2 Data

The data reduction was performed in a standard way using IRAF software. Optical images were bias- and overscan-subtracted, and then flat-fielded. The flat-field images were generated from a combination of night sky flats (created from all object images such that any sources are excluded) and flats made from exposures of a uniformly cloudy sky. NIR images were bias- and dark-subtracted, and then flat-fielded using flats derived from the difference of exposures of the dome with and without illumination. This removes telescope emission and improves photometric accuracy. The sky emission was then subtracted from the individual images using a median of all dither positions (excluding those with high background due to cloud), and then the individual dither positions were added together to produce the final image.

The photometry was performed using the IRAF daophot package. The resulting optical and NIR magnitudes are shown in Tables 2 and 3 respectively. The errors quoted are $`1\sigma `$ errors resulting from the sum in quadrature of the random photometric errors, as calculated by the daophot package, and the systematic errors that result from the spread in zero points for the various standard star observations over the course of the night. The photometry includes corrections for the varying atmospheric extinction due to the different airmasses of the sources. For two sources (PKS 0506$``$612 and PKS 1622$``$297), subtraction of nearby stars was required to provide accurate background subtraction for the aperture photometry (although in no case were the images of the star and quasar overlapping). Figure 1 plots the observed photometric SED for each source as a function of observed wavelength.

The magnitudes were converted into fluxes to enable modelling and plotting of the data. The zero magnitude fluxes for the optical bands were taken from Bessell, Castelli, & Plez (1998), while the CASPIR zero magnitudes are obtained from a 11200 K blackbody normalised to $`F_\lambda (555\mathrm{n}\mathrm{m})=3.44\times 10^{12}`$ W cm-2$`\mu \mathrm{m}^1`$ (Bersanelli, Bouchet, & Falomo 1991). The values of the zero points are shown in Table 4. Note that although the observations used the $`K_n`$ filter, they were calibrated according to the $`K`$ magnitudes of the IRIS standards without a colour-correction term, and so have been normalised using the $`K`$ band zero point.

### 2.3 Galactic Extinction

Values for the Galactic extinction, taken from the maps of Schlegel, Finkbeiner, & Davis (1998) are included in Table 1. For two objects, PKS 0528$`+`$134 and PKS 1622$``$297, the extinction is large due to their location close to the Galactic plane. We note that Zhang et al. (1994) found evidence for excess absorption over the Galactic value for PKS 0528$`+`$134 from X-ray observations, which would imply a larger $`E(BV)`$ value than that from Schlegel et al. (1998) (and hence a larger correction to the observed magnitudes). For our analysis, however, we just use the Galactic value.

The values of E$`(BV)`$ were transformed to an absorption $`A_\lambda `$ for each photometric band using the conversion factors given in Schlegel et al. (1998). These absorption values were then subtracted from the observed magnitudes to give the extinction-corrected magnitudes, and both sets of data (raw and corrected) are included in Tables 2 and 3. The photometric extinction-corrected SEDs for each object have also been plotted in Figure 1.

## 3 Model Fitting

Each of the extinction-corrected SEDs has been fitted by a power law, given by $`f_\lambda (\lambda )\lambda ^\alpha `$, using a $`\chi ^2`$-minimisation method. The fit is deemed ‘acceptable’ if the reduced $`\chi ^2`$ (i.e. $`\chi ^2/\nu `$, where $`\nu `$ is the number of degrees of freedom in the fitted model) is below the 99% confidence level. For the power-law model fitted to seven data points, $`\nu =5`$ and the acceptance threshold is $`\chi ^2=13.08`$. A fit with a $`\chi ^2`$ value greater than this is rejected (with 99% confidence). The resulting fits are shown in Figure 1, and the fitted parameter values—the spectral index $`\alpha `$ and the reduced chi-squared value $`\chi ^2/\nu `$—are shown in Table 5.

We find that the fitted power laws can be separated into two groups. Three sources—PKS 0506$``$612, PKS 0528$`+`$134, and PKS 1622$``$297—are quite blue (although PKS 0506$``$612 is a poor fit, as discussed below, and should probably not be included in this group). These slopes are the same as those found at the blue end of the distribution of fitted slopes for the Parkes Half-Jansky Flat-spectrum Sample (PHFS) quasars in Whiting, Webster, & Francis (2001). A natural interpretation of such blue slopes is that they are due to emission from the accretion disc. This would be expected for high-redshift sources such as PKS 0528$`+`$134, as the observed optical emission is probing the rest-frame UV.

The remaining three sources have much redder optical slopes. Again comparing to the PHFS sources in Whiting et al. (2001), we find that these slopes are in the middle of the distribution, corresponding to sources with synchrotron-dominated spectra.

One source had a poor fit: PKS 0506$``$612. There is a large offset between the optical and NIR points, giving what appears to be substantial spectral curvature. However, it is likely that the source has faded in the three days between the two observations—if the spectrum is a smooth continuum then a reduction in flux of $``$ 1 mag is implied. This also implies the intrinsic spectrum is a red power law of slope $`1.3`$. A similar, albeit much weaker effect is seen in PKS 0528$`+`$134, although this could also be due to the $`I`$ band data having excess emission over the power law continuum. Contamination by the Mg ii emission line (which would appear at $`\lambda _{\mathrm{obs}}=0.86\mu \mathrm{m}`$) is a likely origin of this excess.

We can use data from the literature to support the interpretation of the different optical slopes. Synchrotron radiation is intrinsically highly polarised, and so a significant net polarisation from a source is a good indication that the emission is dominated by synchrotron emission. Impey & Tapia (1990) measured polarisation for five of the sources (excluding PKS 1622$``$297), and their values are also shown in Figure 5. It can be seen that the three sources with red optical slopes are the sources that have been observed to exhibit high optical polarisations, indicating a likely synchrotron origin for the emission.

## 4 Discussion

### 4.1 Demographics

All the sources we have studied here are radio-loud, flat-spectrum quasars (FSRQs). We would like to know how typical these objects are compared to other FSRQs, or whether they are exceptional in some way that would explain their gamma-ray emission. We examine here two properties—the redshift and the radio luminosity—and compare it to a sample of flat-spectrum radio-loud objects. For this comparison sample we use the PHFS (Drinkwater et al. 1997), a sample of southern, radio-bright, flat-spectrum objects. We also include in this discussion the COMPTEL sources that we did not observe: PKS 1222$`+`$216 ($`z=0.435`$), PKS 1226$`+`$023 (3C 273, $`z=0.158`$), and PKS 1253$``$055 (3C 279, $`z=0.536`$). Note that the latter two sources are also in the PHFS.

The cumulative redshift distributions for both samples are plotted in Figure 2a. Note that we have excluded from the PHFS distribution the 33 objects with no known redshift. These objects include a few BL Lac objects, which have no emission line from which to measure a redshift, but are mostly optically faint, meaning no suitable spectrum has been obtained. The COMPTEL sample tends to be biased towards lower redshift values, with only one source at high redshift $`z1`$ (PKS 0528$`+`$134 at $`z=2.07`$, Hunter et al. 1993). This is likely to be a consequence of the sensitivity of the COMPTEL detector, with sources at suitably high redshift not being bright enough at gamma rays for detection. We do note, however, that a Kolmogorov-Smirnov (K-S) test shows that the distributions are different only at the 90% confidence level (the probability that they are the same is 9.9%).

Figure 2b shows the cumulative distributions of the radio luminosity for both samples. This is calculated from the 2.7 GHz fluxes from PKSCAT90 (assuming $`H_0=75`$ km s-1 Mpc-1 and $`q_0=0.5`$). We see that the COMPTEL sources tend to have, on average, higher radio luminosities than the PHFS sources (there is a difference in average luminosity of 0.56 dex, or a factor of $`3.6`$). Again, however, a K-S test shows this is not very significant (the probability that the distributions are different is 12.6%). We can, however, discuss the trends shown in the graphs.

Even though, as we see from the redshift distribution, the COMPTEL sample lacks the high-redshift sources (which one would expect to be more luminous than low-redshift sources in a flux-limited sample), the average radio luminosity of this sample is still greater than that of the PHFS. This implies that the sample selected at gamma rays lacks sources of low radio luminosity. This is probably an indication that a certain radio luminosity is required for a source to be active at high-energy gamma rays.

### 4.2 Multiwavelength SEDs

In Figure 3 we plot the multiwavelength (radio to gamma-ray) SEDs constructed from published data, as well as the new data presented herein. Two epochs of radio observation are shown. One set of data (indicated by the asterisks) is from PKSCAT90, while the second (indicated by crosses) is from Kovalev et al. (1999). The latter does not include PKS 0208$``$512 and PKS 0506$``$612. The optical/NIR data (from this paper) have been corrected for Galactic extinction. The X-ray data come from Brinkmann, Yuan, & Siebert (1997). The ‘bowtie’ represents the twin constraints given by the errors on the total 0.1–2.4 keV flux and the spectral slope. (Note that the constraints on the slope for PKS 1622$``$297 are rather poor.) The COMPTEL data, taken from Schönfelder et al. (2000), are indicated by the four points between $`10^{20}`$ Hz and $`10^{22}`$ Hz. The horizontal error bars indicate the energy range of each bin, while the vertical error bars indicate the flux error. Upper limits ($`2\sigma `$) are indicated by arrows. The final point is the EGRET datum from the Third EGRET Catalog (Hartman et al. 1999), using their approximation to generate a 400 MeV flux density by multiplying the catalogued flux by 1.7. Note that the observations presented here are not simultaneous.

The SEDs provide an interesting glimpse of the range of different spectral shapes present in radio-loud quasars. The optically-red sources, particularly PKS 0208$``$512 and PKS 2251$`+`$158, show optical–NIR SEDs that can quite easily be connected to the radio data, as would be expected if both regimes are dominated by synchrotron emission from the relativistic jet.

For other sources, however, it is less obvious what, if any, connection exists between these two regimes. For instance, the $`K`$ band flux of PKS 1622$``$297 is approximately equal (on a $`\nu F_\nu `$ scale) to the 20 GHz flux from Kovalev et al. (1999), while the spectral index in the two regimes is approximately the same. This, however, should not be surprising if we attribute the optical emission to the accretion disc, as it will then not be the same component as the radio emission. In such a picture, the synchrotron emission from the jet has turned over somewhere in the infrared, leaving the accretion disc emission to dominate the optical. Note that the flux level at radio frequencies is approximately equal for all sources, whereas the redder sources tend to be brighter in the optical/NIR, indicating the presence of an additional emission component.

Clearly, while the non-thermal processes (i.e. synchrotron and IC) are undoubtedly important for the radio and high-energy emission, the optical regime is not always going to be dominated by them. This was shown in the detailed study of the optical/NIR properties of radio-loud quasars by Whiting et al. (2001), which found that a large fraction of sources showed no evidence for non-thermal synchrotron emission at optical wavelengths. The observations of strong emission-line spectra and low polarisation in the optical for the blue objects in this sample also back this up.

This raises an important issue when it comes to modelling the multiwavelength SEDs of such objects. Detailed models do exist in the literature, and are used with good effect to model the high-energy emission (X-rays through gamma rays). See, for instance, Mukherjee et al. (1999) for modelling of the multiwavelength spectra of PKS 0528$`+`$134. However, many of these models do not take the lower-energy (i.e. radio through optical) data fully into account. As an example, we refer the reader to Figures 8 and 9 of Mukherjee et al. (1999), where the model flux in the optical goes in the opposite direction to the data. This seems to be due to the fitted models not fully taking into account the presence of significant emission from the accretion disc.

This is an important issue, since if the thermal accretion disc emission is dominating the optical, the synchrotron emission in that frequency range is much less than previously assumed. This will then affect the derived values for the Lorentz factors, which affect the generation of the gamma rays via the IC process.

The photon density at the jet of the seed photons is another key parameter for calculating the IC luminosity. For a given observed total optical flux, this parameter will be different for synchrotron-dominated and thermally-dominated spectra, and its spectral shape will differ as well. This will alter the predicted shape of the IC emission.

It can be seen, therefore, that a good understanding of the makeup of the optical emission is crucial for accurate modelling of the high-energy emission. This is something that needs to be considered for future multiwavelength campaigns involving high-energy observations and subsequent modelling, particularly with the advent of new gamma-ray telescopes such as INTEGRAL and GLAST.

## 5 Conclusions

We have observed six COMPTEL-detected quasars at optical and NIR wavelengths, and we present photometry in seven wavebands, corrected for Galactic extinction.

We find large differences in the optical properties between the sources. Two of the sources have intrinsically blue optical SEDs, and are likely to be dominated by thermal emission from an accretion disc. A further three have much redder optical slopes, and are more likely to be dominated in the optical by non-thermal synchrotron emission from the relativistic jet.

This identification of accretion disc emission in some gamma-ray sources has important implications for the sources’ broad-band modelling. It implies both that the thermal emission is stronger and the non-thermal emission weaker in the optical regime than previously considered, indicating that broad-band models used to explain the high-energy emission will need to be revisited, particularly as far as the optical regime is concerned.

## Acknowledgments

The authors would like to thank A. Melatos for some very helpful discussions. We also want to thank the MSSSO TAC for granting us time for the observations at the 40 in and 2.3 m telescopes at the SSO.

This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

## References

Bersanelli, M., Bouchet, P., & Falomo, R. 1991, A&A, 252, 854

Bessell, M.S., Castelli, F., & Plez, B. 1998, A&A, 333, 231

Brinkmann, W., Yuan, W., & Siebert, J. 1997, A&A, 319, 413

Bloom, S.D., & Marscher, A.P. 1996, ApJ, 461, 657

Carter, B.S., & Meadows, V.S. 1995, MNRAS, 276, 734

Collmar, W. 2001, in High Energy Gamma-Ray Astronomy, AIP Proceedings Volume 558, ed. F.A. Aharonian $`\&`$ H.J. Völk (Melville, NY: American Institute of Physics), 656

Drinkwater, M.J., et al. 1997, MNRAS, 284, 85

Graham, J.A. 1982, PASP, 94, 244

Hartman, R.C., et al. 1999, ApJS, 123, 79

Hunter, S.D., et al. 1993, ApJ, 409, 134

Impey, C.D., & Tapia, S. 1990, ApJ, 354, 124

Kovalev, Y.Y., Nizhelsky, N.A., Kovalev, Yu.A., Berlin, A.B., Zhekanis, G.V., Mingaliev, M.G., & Bogdantsov, A.V. 1999, A&AS, 139, 545

McGregor, P., et al. 1994, in Infrared Astronomy with Arrays: The Next Generation, ed. I.S. McLean (Dordrecht: Kluwer), 299

Mukherjee, R., et al. 1999, ApJ, 527, 132

Raiteri, C.M., Villata, M., Lanteri, L., Cavallone, M., & Sobrito, G. 1998, A&AS, 130, 495

Schlegel, D.J., Finkbeiner, D.P., & Davis, M. 1998, ApJ, 500, 525

Schönfelder, V., et al. 2000, A&AS, 143, 145

Whiting, M.T., Webster, R.L., & Francis, P.J. 2001, MNRAS, 323, 718

Wright, A., & Otrupcek, R. 1990, Parkes Catalogue (Australian Telescope National Facility) (PKSCAT90)

Zhang, Y.F., Marscher, A.P., Aller, H.D., Aller, M.F., Terasranta, H., & Valtaoja, E. 1994, ApJ, 432, 91


format: text

text:
# The enumeration of simple permutations

## 1 Introduction and definitions

The permutation $`2647513`$ maps the interval $`\mathrm{2..5}`$ onto the interval $`\mathrm{4..7}`$. In other words, it has a *segment* (set of consecutive positions) whose values form a *range* (set of consecutive values). Such a segment is called a *block* of the permutation. Every permutation has singleton blocks, together with the block $`1..n`$. If these are the only blocks the permutation is called *simple*. For example, $`58317462`$ is simple and the simple permutations of length up to $`5`$ are as follows:.

| Length | Simple permutations |
| --- | --- |
| $`1`$ | 1 |
| $`2`$ | 12, 21 |
| $`3`$ | None |
| $`4`$ | 2413, 3142 |
| $`5`$ | 24153, 25314, 31524, 35142, 41352, 42513 |

Simple permutations have recently had important applications in the study of pattern closed classes of permutations .

Let $`s_n`$ denote the number of simple permutations of length $`n`$. We shall be concerned with properties of the sequence $`(s_n)`$. Consider the ordinary generating functions:

$`F(x)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k!x^k;`$
$`S(x)`$ $`=`$ $`{\displaystyle \underset{k=4}{\overset{\mathrm{}}{}}}s_kx^k.`$
We start $`S(x)`$ from $`x^4`$ because simple permutations of length $`1`$ and $`2`$ need special treatment. Later in this section we will see that the coefficients of $`S`$ differ from those of $`F^1`$ (functional inverse, not reciprocal) alternately by $`2`$ and $`2`$. The coefficients of $`F^1(x)`$ were considered by Comtet \[4, p. 171\] without any combinatorial interpretation. The sequence of absolute values of these coefficients appears as sequence A059372 of , and the first few terms are:

$$1,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}4},\mathrm{\hspace{0.17em}4},\mathrm{\hspace{0.17em}48},\mathrm{\hspace{0.17em}336},\mathrm{\hspace{0.17em}2928},\mathrm{\hspace{0.17em}28144},\mathrm{\hspace{0.17em}298528},\mathrm{\hspace{0.17em}3454432},\mathrm{\hspace{0.17em}43286528}.$$

So we shall see that the numbers $`s_n`$ are:

$$1,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}6},\mathrm{\hspace{0.17em}46},\mathrm{\hspace{0.17em}338},\mathrm{\hspace{0.17em}2926},\mathrm{\hspace{0.17em}28146},\mathrm{\hspace{0.17em}298526},\mathrm{\hspace{0.17em}3454434},\mathrm{\hspace{0.17em}43286526}.$$

In section 2 we shall prove that $`(s_n)`$ is not P-recursive (it cannot be defined by a linear recurrence with polynomial coefficients). In section 3 we derive the asymptotic behaviour of $`s_n`$ (the main term is $`n!/\mathrm{e}^2`$) and section 4 gives various congruences satisfied by the numbers $`s_n`$.

In the remainder of this section we derive a structure theorem that shows how arbitrary permutations are built from simple ones, and read off from it equations satisfied by generating functions. We begin with some terminology and notation that will be used throughout.

A block decomposition of a permutation $`\sigma `$ is a partition of $`\sigma `$ into blocks. Of course, if $`\sigma `$ is simple there will only be the two trivial block decompositions. An example of a non-trivial decomposition is $`\sigma =67183524`$ with blocks $`(67)(1)(8)(3524)`$.

Given a block decomposition of $`\sigma `$, its *pattern* is the permutation defined by the relative order of the blocks. In the example above, the pattern of the block decomposition $`(67)(1)(8)(3524)`$ is $`3142`$. We may think of the permutation $`67183524`$ as being constructed from the permutation $`3142`$ by inflating each of the elements into a block, in this case the blocks $`12`$, $`1`$, $`1`$, and $`2413`$ (we view each block as a permutation in its own right). We write:

$$67183524=(3142)[12,1,1,2413].$$

This example is further illustrated in Figure 1. The inflation procedure is an instance of the wreath product for permutations .

A permutation which cannot be written in the form $`(12)[\alpha ,\beta ]`$ is called plus indecomposable, and one which cannot be written in the form $`(21)[\alpha ,\beta ]`$ is called minus indecomposable. Let $`i_n`$ denote the number of plus indecomposable permutations of length $`n`$. The number of minus indecomposable permutations of length $`n`$ is also $`i_n`$ as is easily seen by considering the bijection on permutations of length $`n`$ which sends $`\pi `$ to $`\pi ^{}`$ where $`\pi ^{}(t)=n+1\pi (t)`$.

###### Theorem 1

For every non-singleton permutation $`\pi `$ there exists a unique simple non-singleton permutation $`\sigma `$, and permutations $`\alpha _1,\alpha _2,\mathrm{},\alpha _k`$ such that

$$\pi =\sigma [\alpha _1,\alpha _2,\mathrm{},\alpha _k].$$

Moreover, if $`\sigma 12,21`$ then $`\alpha _1,\alpha _2,\mathrm{},\alpha _k`$ are also uniquely determined. If $`\sigma =12`$ (respectively $`21`$) then $`\alpha _1`$ and $`\alpha _2`$ are uniquely determined subject to the additional condition that $`\alpha _1`$ be plus (respectively minus) indecomposable.

The caveat added for the case where $`\sigma =12`$ (or 21) is necessary, as is easily seen by considering $`\pi =123`$. This can be decomposed as $`(12)[1,12]`$ or as $`(12)[12,1]`$. However, only the former decomposition has a plus indecomposable first part.

Proof: We first of all suppose that $`\pi `$ has two distinct maximal proper blocks $`A`$ and $`B`$ that have a non-empty intersection. Then, as the union of intersecting segments is a segment and the union of intersecting ranges is a range, $`AB`$ is a block. Because of the maximality, $`AB=[n]`$. But it is also clear that $`A`$ cannot be an interior segment of $`[n]`$ nor can it define an interior range. In other words we have

$$\pi =\sigma [\alpha ,\beta ]$$

where $`\sigma =12`$ or $`\sigma =21`$. These two possibilities are obviously mutually exclusive. In either case consider all decompositions of $`\pi `$ as $`\sigma [\gamma ,\delta ]`$. The intersection of their $`\gamma `$ parts is also the $`\gamma `$ part of a decomposition of this type. So there is a unique such decomposition with smallest $`\gamma `$ part. Clearly, this part is plus indecomposable in the case $`\sigma =12`$ and minus indecomposable if $`\sigma =21`$.

We next suppose that every pair of distinct maximal blocks has empty intersection. Obviously, then the maximal blocks form a block decomposition of $`\pi `$ and this decomposition must be coarser than every other block decomposition of $`\pi `$. It follows that this decomposition is the only one whose pattern $`\sigma `$ is simple and so we obtain the unique representation claimed for $`\pi `$.

We shall shortly see that this theorem gives relations between the following three generating functions:

$`F(x)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k!x^k;`$
$`I(x)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}i_kx^k;`$
$`S(x)`$ $`=`$ $`{\displaystyle \underset{k=4}{\overset{\mathrm{}}{}}}s_kx^k.`$
Note that our generating functions are all taken to have zero constant term. This slightly unconventional choice turns out to be algebraically convenient at several points.

From Theorem 1 it is easy to see that there is a one to one correspondence between the collection of all permutations with length at least 2 and the collection of sequences:

$$(\sigma ,\alpha _1,\alpha _2,\mathrm{},\alpha _k).$$

Here $`\sigma `$ may be any simple permutation of length $`k2`$, and if $`\sigma 12,21`$ then $`\alpha _1`$ through $`\alpha _k`$ are arbitrary permutations, while if $`\sigma =12`$ (respectively 21), $`\alpha _1`$ is plus-indecomposable (respectively minus indecomposable) and $`\alpha _2`$ is arbitrary.

This correspondence, together with the earlier observation that the numbers of plus and minus indecomposable permutations of length $`n`$ are the same, translates naturally into the following equation:

$$F(x)=x+2I(x)F(x)+(SF)(x).$$
(1)
However, since a plus indecomposable permutation cannot correspond to a sequence beginning with $`12`$, while all other sequences do represent plus indecomposables, it is also clear from the correspondence that

$$I(x)=x+I(x)F(x)+(SF)(x).$$

Solving this latter equation for $`I`$, and then substituting in equation (1) before solving for $`SF`$ gives:

$$(SF)(x)=\frac{F(x)F(x)^2}{1+F(x)}x.$$

Now letting $`t=F(x)`$ we obtain:

$$S(t)=t\frac{2t^2}{1+t}F^1(t).$$
(2)
We can also obtain an equation for the ordinary generating function of plus indecomposable permutations through the observation that every permutation decomposes into a sequence of plus indecomposable permutations so

$$F(x)=\frac{I(x)}{1I(x)}$$

or equivalently

$$I(x)=\frac{F(x)}{1+F(x)}.$$
(3)
Denoting the coefficient of $`t^n`$ in $`F^1(t)`$ by $`\mathrm{Com}_n`$ (in reference to Comtet who initiated the consideration of this sequence) we obtain directly from equation (2) the simple relationship that for $`n4`$:

$$s_n=\mathrm{Com}_n+(1)^{n+1}2.$$

## 2 Non P-recursiveness

A sequence of numbers $`(a_n)`$ is called P-recursive if it satisfies a linear recurrence with polynomial coefficients. A power series is called D-finite if it satisfies a linear differential equation with polynomial coefficients. A sequence $`(a_n)`$ is P-recursive if and only if its ordinary generating function $`A(x)=_na_nx^n`$ is D-finite. More information on D-finiteness and P-recursiveness can be found in Stanley \[13, Chapter 6\]. If $`a_n=n!`$ then $`a_nna_{n1}=0`$, and thus the sequence $`(n!)`$ is P-recursive. We show that on the other hand neither sequence $`(i_n)`$ nor $`(s_n)`$ is P-recursive. By (2), instead of the latter sequence we can work with $`(\mathrm{Com}_n)`$.

###### Proposition 2

The power series $`I(x)`$ and $`C(x)=F^1(x)=_{k=1}^{\mathrm{}}\mathrm{Com}_kx^k`$ satisfy the differential equations

$`I^{}`$ $`=`$ $`x^2I^2+(x^2+x^1)Ix^1;`$
$`C^{}`$ $`=`$ $`{\displaystyle \frac{C^2}{x(1+x)C}}.`$
Proof: It follows from the recurrence for $`n!`$ that $`F(x)`$ satisfies $`x+xF+x^2F^{}=F`$. Thus $`F^{}=((1x)Fx)/x^2`$. Combining this with $`F=I/(1I)`$ we obtain the differential equation for $`I(x)`$. Similarly, $`C^{}=1/F^{}(C)=C^2/((1C)xC)`$ which is the differential equation for $`C(x)`$.

Klazar  used the following method to show that a sequence $`(a_n)`$ is not P-recursive. Suppose that the ordinary generating function $`A(x)`$ is non-analytic and satisfies a first order differential equation $`A^{}=R(x,A)`$ where $`R`$ is some expression. Differentiating this relationship and replacing $`A^{}`$ by $`R(x,A)`$, the derivatives of $`A`$ are expressed as $`A^{(k)}=R_k(x,A)`$; $`R_0(x,A)=A`$ and $`R_1(x,A)=R(x,A)`$. Substituting $`R_k(x,A)`$ in the equation of D-finiteness

$$b_0A+b_1A^{}+b_2A^{\prime \prime }+\mathrm{}+b_sA^{(s)}=0,$$

where $`s1`$, $`b_i𝐂(x)`$ and $`b_s0`$, we get a non-differential equation $`_{k=0}^sb_kR_k(x,A)=0`$. If $`R`$ is such that the expressions $`R_0,R_1,R_2,\mathrm{}`$ are (i) analytic or even algebraic and (ii) linearly independent over $`𝐂(x)`$, we have a nontrivial analytic equation for $`A`$. This implies that $`A`$ is analytic (see Klazar’s paper  for more details) which is a contradiction. So $`A`$ cannot be D-finite and the sequence of its coefficients cannot be P-recursive.

To state the result of  precisely, we remind the reader that a power series $`R(x,y)𝐂[[x,y]]`$ is analytic if it absolutely converges in a neighborhood of the origin and that $`R(x,y)𝐂((x,y))`$ is an analytic Laurent series if, for some positive integer $`k`$, $`(xy)^kR(x,y)𝐂[[x,y]]`$ is analytic. Theorem 1 of  says that if $`A𝐂[[x]]`$ is non-analytic, $`R(x,y)𝐂((x,y))`$ is analytic, $`A^{}=R(x,A)`$, and $`R`$ contains at least one monomial $`ax^iy^j`$, $`a0`$, with $`j<0`$, then $`A`$ is not D-finite. This result applies directly neither to $`I(x)`$ nor $`C(x)`$ (see Proposition 2) because in the case of $`I(x)`$ the last condition on $`R`$ is not satisfied and in the case of $`C(x)`$ the right hand side $`R`$ even cannot be expanded as a Laurent series.

However, the substitution $`x(1+x)C(x)=\theta (x)`$ transforms the second differential equation of Proposition 2 into

$$\theta ^{}=\frac{x^2}{1+x}\frac{1}{\theta }+\frac{1+2x}{1+x}.$$

Now all conditions are satisfied ($`F(x)`$ is clearly non-analytic which implies that $`C(x)`$ and $`\theta (x)`$ are non-analytic) and thus $`\theta (x)`$ is not D-finite by Theorem 1 of . The dependence of $`C(x)`$ and $`S(x)`$ on $`\theta (x)`$ and the fact that D-finite power series form a $`𝐂(x)`$-algebra (\[13, Theorem 6.4.9\]) shows that neither $`C(x)`$ nor $`S(x)`$ is D-finite.

In order to deal with the case of $`I(x)`$, we use this opportunity to complement Theorem 1 of  in which $`R𝐂((x,y))`$ by the following theorem which treats the case $`R𝐂(x,y)`$. Neither of the theorems subsumes the other because not every rational function in $`x`$ and $`y`$ can be represented by an element of $`𝐂((x,y))`$ (as we have seen) and, of course, not every Laurent series sums up to a rational function. However, the next theorem seems to be more useful because in both examples in  and both examples here the right hand side $`R(x,y)`$ is, in fact, a rational function.

###### Theorem 3

Let $`P,Q𝐂[x,y]`$ be two nonzero coprime polynomials and $`A𝐂[[x]]`$ be a non-analytic power series which satisfies the differential equation

$$A^{}=\frac{P(x,A)}{Q(x,A)}.$$

If $`\mathrm{deg}_yQ=0`$ and $`\mathrm{deg}_yP1`$ then $`A`$ is, trivially, D-finite. In all remaining cases $`A`$ is not D-finite.

Proof: The first claim is clear. If $`\mathrm{deg}_yQ=0`$ and $`r=\mathrm{deg}_yP2`$ then $`A^{}=a_0+a_1A+\mathrm{}+a_rA^r`$ where $`a_i𝐂(x)`$, $`r2`$, and $`a_r0`$. Differentiation by $`x`$ gives

$$A^{(k)}=R_k(x,A)=a_{0,k}+a_{1,k}A+\mathrm{}+a_{krk+1,k}A^{krk+1}$$

where $`a_{i,j}𝐂(x)`$ and

$$a_{krk+1,k}=r(2r1)(3r2)\mathrm{}((k1)rk+2)a_r^k0.$$

Thus $`R_k(x,y)𝐂(x)[y]`$ have $`y`$-degrees $`krk+1`$, $`k=0,1,2,\mathrm{}`$, which is for $`r2`$ a strictly increasing sequence. Therefore $`R_0,R_1,R_2,\mathrm{}`$ are linearly independent over $`𝐂(x)`$ and, by the above discussion, $`A`$ is not $`D`$-finite.

In the remaining case $`\mathrm{deg}_yQ1`$. Differentiation of $`A^{}=R(x,A)=P(x,A)/Q(x,A)`$ by $`x`$ gives $`A^{(k)}=R_k(x,A)`$ where $`R_k(x,y)𝐂(x,y)`$. For example,

$`R_2`$ $`=`$ $`{\displaystyle \frac{(P_x+P_yR_1)QP(Q_x+Q_yR_1)}{Q^2}}`$
$`=`$ $`{\displaystyle \frac{P_xQPQ_x}{Q^2}}+{\displaystyle \frac{P(P_yQPQ_y)}{Q^3}}.`$
Let $`\alpha `$, $`Q(x,\alpha )=0`$, be a pole of $`R_1(x,y)`$ of order $`\mathrm{ord}_\alpha (R_1)=\mathrm{ord}_\alpha (P/Q)=\mathrm{ord}_\alpha (Q)=l1`$. We have $`\mathrm{ord}_\alpha ((P_xQPQ_x)Q^2)2l`$ and $`\mathrm{ord}_\alpha (P(P_yQPQ_y)Q^3)=3l+\mathrm{ord}_\alpha (P_yQPQ_y)=2l+1`$ since $`\mathrm{ord}_\alpha (P)=0`$, $`\mathrm{ord}_\alpha (P_yQ)l`$, and $`\mathrm{ord}_\alpha (PQ_y)=l+1`$. So $`\mathrm{ord}_\alpha (R_2)=2l+1`$. In general, the same argument shows that $`\mathrm{ord}_\alpha (R_{k+1})=2\mathrm{ord}_\alpha (R_k)+1`$. Hence $`\mathrm{ord}_\alpha (R_k)=2^{k1}l+2^{k1}1`$, $`k=1,2,\mathrm{}`$. This is a strictly increasing sequence and we conclude again, since $`R_0,R_1,R_2,\mathrm{}`$ are linearly independent over $`𝐂(x)`$, that $`A`$ is not $`D`$-finite.

Proposition 2 and Theorem 3 show that $`I(x)`$ is not $`D`$-finite and we can summarize the results of this section in the following corollary.

###### Corollary 4

The sequences $`(i_n)`$, $`(\mathrm{Com}_n)`$, and $`(s_n)`$ are not P-recursive.

## 3 Asymptotics

We turn now to the computation of an asymptotic expansion for the numbers $`s_n`$. We will prove that:

###### Theorem 5

$$s_n=\frac{n!}{\mathrm{e}^2}\left(1\frac{4}{n}+\frac{2}{n(n1)}+O(n^3)\right).$$

Our methods are such that, in principle, higher order terms could be obtained as a matter of brute force computation. In order to carry out this expansion we will first consider permutations which may not be simple, but whose non-trivial blocks all have length greater than some fixed value $`m`$. We will apply inclusion-exclusion arguments (dressed in the form of generating functions ), an argument which allows us to reduce the number of terms considered, and a bootstrapping approach.

The case $`m=2`$, was already considered by Kaplansky . Permutations of this type are those in which no two elements consecutive in position are also consecutive in value (in either order). These were called irreducible permutations by Atkinson and Stitt , but there is no standard terminology in the field. Indeed the permutations that we have referred to as plus and minus indecomposable have also been called irreducible in other contexts.

An amusing equivalent form for the case $`m=2`$ is that the number of such permutations is also the number of ways of placing $`n`$ mutually non-attacking krooks on an $`n\times n`$ chessboard. A krook is a piece which can move either like a king, or a rook in chess. Kaplansky’s expansion is:

$$\frac{n!}{\mathrm{e}^2}\left(1\frac{2}{n(n1)}+O(n^3)\right).$$

In fact he derives asymptotic forms for the number of permutations containing exactly $`r`$ blocks of length 2 for any $`r`$. Our methods parallel his, and could also be used to derive such detailed information.

The decomposition provided by Theorem 1 of a permutation into its maximal proper blocks represents a top down view of how non-simple permutations are constructed from simple ones. There is a corresponding bottom-up view that focuses on minimal blocks, put together in an arbitrary order. By a minimal block in $`\pi `$ we mean a non-singleton block in $`\pi `$ minimal with respect to inclusion. Note that the pattern of each minimal block is that of a simple permutation. Any permutation can be decomposed into minimal blocks and singletons, e.g., $`3524716=(3524)(7)(1)(6)`$. However, this decomposition is not unique, for two essentially different reasons. The first one is that decompositions $`\pi =\sigma [\alpha _1,\alpha _2,\mathrm{},\alpha _k]`$, where $`\sigma `$ is arbitrary and $`\alpha _i`$ are simple, are not unique because it may be possible to coalesce singletons into simple blocks, or vice versa. Thus besides $`3524716=2413[2413,1,1,1]`$ we also have $`3524716=3524716[1,1,1,1,1,1,1]`$. The second problem is that we require any two minimal blocks to be disjoint. While this is necessarily true whenever either of them has length more than 2, two minimal blocks of length 2 may intersect, as in $`123`$. Thus we consider decompositions $`\pi =\sigma [\alpha _1,\alpha _2,\mathrm{},\alpha _k]`$ where $`\sigma `$ is arbitrary and each $`\alpha _i`$ is either 1, a simple permutation of length at least 4, or the identity permutation of length at least 2 or its reverse. We refer to blocks of the latter type as clusters in $`\pi `$.

By using clusters we have solved the second problem but the non-uniqueness remains and, moreover, we have introduced another source of it: consecutive (reversed) identical permutations may coalesce into longer (reversed) identical permutations, as in $`345612=21[1234,12]=231[12,12,12]`$. To remedy the non-uniqueness we introduce the notion of marking a permutation. A marked permutation $`(\pi ,M)`$ consists of a permutation $`\pi `$ and a collection $`M`$ of minimal blocks of $`\pi `$. A marked cluster in $`(\pi ,M)`$ is a maximal chain of marked overlapping minimal blocks of length 2 (a marked cluster may be a proper subset of a maximal cluster). Let $`_1`$ denote the set of all simple permutations of length at least 4 and $`_2`$ denote the set of all identical permutations of length at least 2 and their reversals. Marking makes our decomposition unique:

###### Theorem 6

Let $`X`$ be the set of all marked permutations $`(\pi ,M)`$ and $`Y`$ be the set of all sequences $`(\sigma ;\alpha _1,\alpha _2,\mathrm{},\alpha _k)`$ where $`\sigma `$ is any permutation of length $`k1`$ and $`\alpha _i\{1\}_1_2`$. There is a bijection between the sets $`X`$ and $`Y`$ such that if $`(\pi ,M)(\sigma ;\alpha _1,\alpha _2,\mathrm{},\alpha _k)`$, where $`r`$ of the $`\alpha _i`$ belong to $`_1`$ and $`s`$ of them to $`_2`$, then

$$\pi =\sigma [\alpha _1,\alpha _2,\mathrm{},\alpha _k]$$

and $`|M|=r+ls`$ where $`l`$ is the total length of the $`\alpha _i`$ belonging to $`_2`$.

Proof: Given a marked permutation, collapse its marked minimal blocks of length at least 4 and its marked clusters into singletons. This gives the permutation $`\sigma `$. If the $`i`$-th term of $`\sigma `$ was not obtained by collapse then $`\alpha _i=1`$, otherwise $`\alpha _i`$ equals to the corresponding element of $`_1_2`$. Since each $`\alpha _i_1`$ contributes 1 to $`|M|`$ and each $`\alpha _i_2`$ of length $`m`$ contributes $`m1`$, we have $`|M|=r+ls`$. It is clear that $`\pi =\sigma [\alpha _1,\alpha _2,\mathrm{},\alpha _k]`$ and that $`(\pi ,M)`$ can be uniquely recovered from $`(\sigma ;\alpha _1,\alpha _2,\mathrm{},\alpha _k)`$.

Now suppose $`m`$ to be some fixed value (we will later make choices of $`m`$ suitable for our purposes, but will always assume that $`m2`$ since smaller values of $`m`$ are trivial). Each permutation $`\pi `$ has an associated collection $`B_m(\pi )`$ consisting of the minimal blocks of $`\pi `$ whose length is less than or equal to $`m`$. So, if $`\pi `$ is simple and of length greater than $`m`$, $`B_m(\pi )`$ is empty, while for $`\pi =5672413`$, $`B_2(\pi )=\{56,67\}`$, and $`B_4(\pi )=\{56,67,2413\}`$. An $`m`$-marking of $`\pi `$ is simply a subset of $`B_m(\pi )`$. We consider the generating function:

$$F_m(x,v)=\underset{\pi }{}x^{|\pi |}\underset{MB_m(\pi )}{}v^{|M|}=\underset{\pi }{}x^{|\pi |}(1+v)^{|B_m(\pi )|}.$$

Then of course $`F_m(x,1)`$ is the ordinary generating function for permutations all of whose non-singleton blocks have length greater than $`m`$.

We remark that $`F_m(x,t1)`$ is the generating function where the coefficient of $`x^nt^k`$ is precisely the number of permutations of length $`n`$ with $`k`$ minimal blocks of length less than or equal to $`m`$.

Let

$$S_m(x)=\underset{j=4}{\overset{m}{}}s_jx^j.$$

We apply the bijection of Theorem 6 to marked permutations which contain no marked minimal blocks of length more than $`m`$. It follows that the generating function of the corresponding permutations $`\alpha \{1\}_1_2`$, in which $`x`$ counts the length and $`v`$ the contribution to $`|M|`$, is

$$x+vS_m(x)+\frac{2vx^2}{1vx}.$$

So:

$$F_m(x,v)=\underset{k1}{}k!\left(x+\frac{2vx^2}{1vx}+vS_m(x)\right)^k$$

from which it follows that:

$$f_m(x):=F_m(x,1)=\underset{k1}{}k!\left(x\frac{2x^2}{1+x}S_m(x)\right)^k.$$
(4)
Before using this equation to derive asymptotic information about $`s_n`$ we digress briefly to show how it can be used to obtain an alternative derivation of (2). Instead of using $`S_m(x)`$ in (4), use $`S(x)`$. This gives us $`f_{\mathrm{}}(x)`$, an ordinary generating function for permutations having no minimal block. The only such permutation is $`1`$ so $`f_{\mathrm{}}(x)=x`$. That is:

$$x=F(x\frac{2x^2}{1+x}S(x))$$

which yields (2) after applying $`F^1`$ to both sides.

Now recall that $`f_m(x)`$ is the generating function for permutations all of whose blocks have length greater than $`m`$. In order to make use of these generating functions in the asymptotic computation of $`s_n`$ we must determine a suitable value of $`m`$ so that $`f_m`$ provides useful information about $`s_n`$. To that end the following lemma is useful.

###### Lemma 7

If $`p_{n,k}`$ denotes the number of permutations of length $`n`$ which contain a minimal block of length $`k`$ then for any fixed positive integer $`c`$:

$$\underset{k=c+2}{\overset{nc}{}}\frac{p_{n,k}}{n!}=O(n^c).$$

Proof: First observe that

$$p_{n,k}s_k(nk+1)(nk+1)!$$

since the right hand side counts the number of ways to choose the structure of a block of length $`k`$, to choose its minimal element, and to arrange it with other elements, so it overcounts permutations with more than one such block.

The estimate given then follows directly by using the fact that $`s_kk!`$. Only the two extreme terms in the sum can have magnitude as large as $`O(n^c)`$, and the remaining terms have magnitude $`O(n^{c1})`$. Since there are fewer than $`n`$ terms, the result follows.

So when seeking an asymptotic expansion of $`s_n/n!`$ with an error term of $`O(n^{c1})`$ we may count instead the permutations which contain no blocks of length less than or equal to $`c+2`$, or greater than or equal to $`nc`$. In particular, as a direct consequence of the result quoted above due to Kaplansky  we obtain:

###### Observation 8

$$\frac{s_n}{n!}=\frac{1}{\mathrm{e}^2}+O(n^1).$$

An alternative proof of this result follows from a more general theorem of Bender and Richmond  which provides the first order asymptotics of a class of series which include the inverse series of $`F(x)`$.

We will set as our goal to obtain the asymptotics of $`s_n/n!`$ with error term $`O(n^3)`$. However, the technique we use is completely general, and could be applied, at the expense of a great deal of tedious computation, to any fixed error bound of this type. By the remarks above, we may ignore minimal block sizes between $`5`$ and $`n3`$ inclusive. We first consider $`f_4(x)`$ which enumerates permutations having no minimal blocks of size less than or equal to 4. Recall that:

$$f_4(x)=\underset{k1}{}k!\left(x\frac{2x^2}{1+x}2x^4\right)^k.$$

So:

$`{\displaystyle \frac{1}{n!}}[t^n]f_4(t)`$ $`=`$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}k![t^n]\left(t{\displaystyle \frac{2t^2}{1+t}}2t^4\right)^k`$ (5)
$`=`$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}k![t^{nk}]\left(1{\displaystyle \frac{2t}{1+t}}2t^3\right)^k`$
$`=`$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \underset{l=0}{\overset{n}{}}}(nl)![t^l]\left(1{\displaystyle \frac{2t}{1+t}}2t^3\right)^{nl}`$
$`=`$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \underset{l=0}{\overset{n}{}}}(nl)!{\displaystyle \underset{i=0}{\overset{l}{}}}(2)^i\left({\displaystyle \genfrac{}{}{0pt}{}{nl}{i}}\right)[t^l]\left({\displaystyle \frac{t}{1+t}}+t^3\right)^i`$
$`=`$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \underset{l=0}{\overset{n}{}}}(nl)!{\displaystyle \underset{i=0}{\overset{l}{}}}(2)^i\left({\displaystyle \genfrac{}{}{0pt}{}{nl}{i}}\right)[t^{li}]\left({\displaystyle \frac{1}{1+t}}+t^2\right)^i.`$
Consider now any fixed value of $`l`$ in equation (5). In order to obtain terms whose order in $`n`$ is $`n^2`$ or more, we need only consider the values $`l2il`$. Despite the fact that we sum over values of $`l`$ running from $`0`$ through $`n`$, we may safely ignore the other terms. As we shall see in computing the three significant terms the summation over $`l`$ does not affect the order of the terms.

So, the three terms that we need to consider are:

$$\begin{array}{c}\frac{(nl)!}{n!}(2)^l\left(\genfrac{}{}{0pt}{}{nl}{l}\right)+\hfill \\ \frac{(nl)!}{n!}\left((2)^{l1}\left(\genfrac{}{}{0pt}{}{nl}{l1}\right)(l+1)\right)+\hfill \\ \frac{(nl)!}{n!}\left((2)^{l2}\left(\genfrac{}{}{0pt}{}{nl}{l2}\right)\left((l+2)(l+1)/2+l2\right)\right).\hfill \end{array}$$
(6)
Each of these terms will be converted to the form:

$$\frac{(2)^l}{l!}\left(\text{an asymptotic expansion in }n\right).$$

Since the first two and the first part of the third, are the same as those arising in the $`m=2`$ case, we can make use of their known form, that is, use the asymptotics from Kaplansky’s result, leaving only the term

$`{\displaystyle \frac{(nl)!(2)^{l2}(l2)}{n!}}\left({\displaystyle \genfrac{}{}{0pt}{}{nl}{l2}}\right)`$ $`=`$ $`{\displaystyle \frac{(2)^l}{l!}}\left({\displaystyle \frac{l(l1)(l2)}{4}}{\displaystyle \frac{(nl)!(nl)!}{n!(n2l+2)!}}\right)`$
$`=`$ $`{\displaystyle \frac{(2)^l}{l!}}\left({\displaystyle \frac{l(l1)(l2)}{4n(n1)}}+O(n^3)\right)`$
Summing this expression over $`l`$ gives $`2\mathrm{e}^2/n(n1)+O(n^3)`$.

Now we combine this additional term with Kaplansky’s results to give the asymptotic expansion of $`[t^n]f_4(t)`$ through three terms as:

$$[t^n]f_4(t)=\frac{n!}{\mathrm{e}^2}\left(1\frac{4}{n(n1)}+O(n^3)\right).$$

Finally we use this in establishing the second order asymptotics of $`s_n`$. From Observation 8 applied to $`s_{n1}`$ we obtain:

$$s_{n1}=\frac{n!}{\mathrm{e}^2}\left(\frac{1}{n}+O(n^2)\right).$$

Furthermore, the number of permutations of length $`n`$ containing a simple block of length $`n1`$ is precisely $`4s_{n1}`$. Since, in computing the $`1/n`$ term in the expansion of $`s_n`$ we can ignore contributions arising from blocks of length $`n2`$, and since the events of having a simple block of length from 2 to 4, and having a simple block of length $`(n1)`$ are disjoint:

$`s_n`$ $`=`$ $`[t^n]f_4(t)4s_{n1}+O(n^2n!);`$
$`=`$ $`{\displaystyle \frac{n!}{\mathrm{e}^2}}\left(1{\displaystyle \frac{4}{n}}+O(n^2)\right).`$
We apply this bootstrap approach once more to get the second order behaviour. We now know that:

$`s_{n1}`$ $`=`$ $`{\displaystyle \frac{n!}{\mathrm{e}^2}}\left({\displaystyle \frac{1}{n}}{\displaystyle \frac{4}{n(n1)}}+O(n^3)\right)`$
$`s_{n2}`$ $`=`$ $`{\displaystyle \frac{n!}{\mathrm{e}^2}}\left({\displaystyle \frac{1}{n(n1)}}+O(n^3)\right).`$
Furthermore there are $`18s_{n2}`$ permutations of length $`n`$ containing a simple block of length $`n2`$. However, of these $`8s_{n2}`$ also contain a simple block of length $`2`$. So:

$`s_n`$ $`=`$ $`[t^n]f_4(t)4s_{n1}10s_{n2}+O(n^3n!)`$
$`=`$ $`{\displaystyle \frac{n!}{\mathrm{e}^2}}\left(1{\displaystyle \frac{4}{n}}+{\displaystyle \frac{2}{n(n1)}}+O(n^3)\right),`$
as we claimed at the beginning of this section.

Finally, in this section we note that the asymptotic estimate of $`s_n`$ is, as might be expected, a poor approximation. For example, $`s_{20}=264111424634864638`$ and our asymptotic estimate has a relative error of about $`3.89\times 10^3`$.

## 4 Congruences

In this section we derive congruence properties of the numbers $`\mathrm{Com}_n`$ for the moduli $`2^a`$ and $`3`$ (from which follow similar congruences for $`s_n`$). Our main tool is the following result that follows immediately from the Lagrange inversion formula.

###### Lemma 9

$$n\mathrm{Com}_n=[x^{n1}]\left(\underset{k0}{}(1)^k(2!x+3!x^2+\mathrm{})^k\right)^n.$$

For a prime $`p`$, let $`\mathrm{ord}_p(n)`$ denote the largest integer $`m`$ such that $`p^m`$ divides $`n`$. As the following table shows, $`\mathrm{ord}_2(\mathrm{Com}_n)`$ is unexpectedly large:

| $`n`$ | $`1`$ | $`2`$ | $`3`$ | $`4`$ | $`5`$ | $`6`$ | $`7`$ | $`8`$ | $`9`$ | $`10`$ | $`11`$ | $`12`$ | $`13`$ | $`14`$ | $`15`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`\mathrm{ord}_2(\mathrm{Com}_n)`$ | $`0`$ | $`1`$ | $`1`$ | $`2`$ | $`2`$ | $`4`$ | $`4`$ | $`4`$ | $`4`$ | $`5`$ | $`5`$ | $`15`$ | $`13`$ | $`12`$ | $`12`$ |

| $`16`$ | $`17`$ | $`18`$ | $`19`$ | $`20`$ | $`21`$ | $`22`$ | $`23`$ | $`24`$ | $`25`$ | $`26`$ | $`27`$ | $`28`$ | $`29`$ | $`30`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`8`$ | $`8`$ | $`9`$ | $`9`$ | $`10`$ | $`10`$ | $`12`$ | $`12`$ | $`14`$ | $`14`$ | $`15`$ | $`15`$ | $`17`$ | $`17`$ | $`22`$ |

In Theorem 11 we give a lower bound on $`\mathrm{ord}_2(\mathrm{Com}_n)`$ which is tight for infinitely many $`n`$ and we completely characterize the values of $`n`$ for which the equality is attained.

For convenience we note the following result that follows directly from the well-known formula

$$\mathrm{ord}_p(m!)=\frac{m}{p}+\frac{m}{p^2}+\mathrm{}$$

###### Lemma 10

For all $`m`$, $`\mathrm{ord}_2((m+1)!)\frac{m}{2}`$ where equality holds if and only if $`m=1`$ or $`2`$. Also, $`\mathrm{ord}_3(m!)m1`$ for all $`m`$.

###### Theorem 11

Let $`m=n/2`$. Then

$$\mathrm{ord}_2(\mathrm{Com}_n)\frac{n1}{2}.$$

Equality holds if and only if $`\left(\genfrac{}{}{0pt}{}{3m}{m}\right)`$ is odd and this happens if and only if the binary expansion of $`m`$ has no two consecutive unit digits.

Proof: Let the numbers $`b_k`$, $`k0`$, be defined by

$$\underset{k0}{}b_kx^k=\underset{k0}{}(1)^k(2!x+3!x^2+\mathrm{})^k.$$

Thus $`b_0=1`$ and for $`k1`$,

$$b_k=\underset{\stackrel{c_1,c_2,\mathrm{},c_s1}{c_1+c_2+\mathrm{}+c_s=k}}{}(1)^s(c_1+1)!(c_2+1)!\mathrm{}(c_s+1)!.$$

By Lemma 9,

$$n\mathrm{Com}_n=\underset{\stackrel{k_1,k_2,\mathrm{},k_n0}{k_1+k_2+\mathrm{}+k_n=n1}}{}b_{k_1}b_{k_2}\mathrm{}b_{k_n}.$$

By Lemma 10, $`\mathrm{ord}_2((c+1)!)c/2`$ for all $`c`$. Hence, for all $`k`$ and $`n`$,

$$\mathrm{ord}_2(b_k)\frac{k}{2}\text{ and }\mathrm{ord}_2(n\mathrm{Com}_n)\frac{n1}{2}.$$

In particular, for odd $`n`$ we have $`\mathrm{ord}_2(\mathrm{Com}_n)=\mathrm{ord}_2(n\mathrm{Com}_n)(n1)/2`$.

To obtain the more exact result of the theorem we need the following better estimates for $`\mathrm{ord}_2(b_k)`$:

$$\mathrm{ord}_2(b_k)\{\begin{array}{cc}=k/2\hfill & \text{for even }k\text{;}\hfill \\ =(k+1)/2\hfill & \text{for }k1\text{ mod }4\text{;}\hfill \\ >(k+1)/2\hfill & \text{for }k3\text{ mod }4\text{.}\hfill \end{array}$$

To prove them we look more closely at the sum for $`b_k`$. Suppose first that $`k`$ is even. Then the sum has exactly one summand with $`\mathrm{ord}_2`$ equal to $`k/2`$, namely that with $`c_1=c_2=\mathrm{}=c_{k/2}=2`$ (by Lemma 10, $`\mathrm{ord}_2((c+1)!)=c/2`$ only if $`c=2`$), and the other summands have $`\mathrm{ord}_2`$ bigger than $`k/2`$. Hence $`\mathrm{ord}_2(b_k)=k/2`$. Now suppose that $`k`$ is odd. Then each summand has an odd number of odd $`c_i`$’s. The summands $`t`$ with three and more odd $`c_i`$’s satisfy $`\mathrm{ord}_2(t)(k+3)/2`$ (each odd $`c_i`$ contributes $`1/2`$ to $`k/2`$). The same is true if $`t`$ has only one odd $`c_i`$ but that $`c_i`$ is not $`1`$ (by Lemma 10, $`\mathrm{ord}_2((c+1)!)(c+3)/2`$ for odd $`c>1`$), or if some even $`c_i`$ is not $`2`$ (Lemma 10). The remaining summands $`t`$, in which $`c_i=2`$ with multiplicity $`(k1)/2`$ and once $`c_i=1`$, satisfy $`\mathrm{ord}_2(t)=(k+1)/2`$. We see that, for odd $`k`$, $`\mathrm{ord}_2(b_k)=(k+1)/2`$ if and only if the number of the remaining summands is odd. This number equals $`(k1)/2+1=(k+1)/2`$. So $`\mathrm{ord}_2(b_k)=(k+1)/2`$ if and only if $`k1`$ mod $`4`$.

Let $`n=2m+1`$ be odd. If $`s`$ is a summand of the above sum for $`n\mathrm{Com}_n`$, then $`\mathrm{ord}_2(s)=(n1)/2`$ if and only if all $`k_i`$ in $`s`$ are even; other summands $`t`$ have $`\mathrm{ord}_2(t)>(n1)/2`$. It follows that $`\mathrm{ord}_2(\mathrm{Com}_n)=(n1)/2`$ if and only if the number of the former summands $`s`$ is odd. This number equals

$$[x^{n1}]\left(\underset{r0}{}x^{2r}\right)^n=[x^{n1}]\frac{1}{(1x^2)^n}=[x^{n1}]\underset{r0}{}\left(\genfrac{}{}{0pt}{}{n+r1}{r}\right)x^{2r}=\left(\genfrac{}{}{0pt}{}{3m}{m}\right).$$

Let $`n=2m`$ be even. We know that $`\mathrm{ord}_2(b_k)=k/2`$ for even $`k`$ and $`\mathrm{ord}_2(b_k)(k+1)/2`$ for odd $`k`$. In the sum for $`n\mathrm{Com}_n`$, every composition $`k_1+k_2+\mathrm{}+k_n=n1`$ of $`n1`$ has an odd number of odd parts. For any $`t`$-tuple $`l_1,l_2,\mathrm{},l_t`$, where $`t`$ and all $`l_i`$ are odd and $`l_1+\mathrm{}+l_tn1`$, we let $`S(l_1,l_2,\mathrm{},l_t)`$ denote the sum of those $`b_{k_1}b_{k_2}\mathrm{}b_{k_n}`$ with $`k_1+k_2+\mathrm{}+k_n=n1`$ in which $`k_i=l_i`$, $`1it`$, and $`k_i`$ is even for $`i>t`$. It follows that

$$n\mathrm{Com}_n=\left(\genfrac{}{}{0pt}{}{n}{t}\right)S(l_1,l_2,\mathrm{},l_t)$$

where we sum over all mentioned $`t`$-tuples $`l_1,l_2,\mathrm{},l_t`$. By the properties of $`\mathrm{ord}_2`$ and of the numbers $`b_k`$, $`\mathrm{ord}_2(S(l_1,l_2,\mathrm{},l_t))(n+t1)/2`$. Also, for odd $`t`$ we have $`\mathrm{ord}_2(\left(\genfrac{}{}{0pt}{}{n}{t}\right))=\mathrm{ord}_2(\frac{n}{t}\left(\genfrac{}{}{0pt}{}{n1}{t1}\right))=\mathrm{ord}_2(n)\mathrm{ord}_2(t)+\mathrm{ord}_2(\left(\genfrac{}{}{0pt}{}{n1}{t1}\right))\mathrm{ord}_2(n)`$, and $`\mathrm{ord}_2(\left(\genfrac{}{}{0pt}{}{n}{1}\right))=\mathrm{ord}_2(n)`$. It follows that $`\mathrm{ord}_2(\mathrm{Com}_n)n/2`$ and, moreover, $`\mathrm{ord}_2(\mathrm{Com}_n)=n/2`$ if and only if

$$\mathrm{ord}_2\left(\underset{ln,l\mathrm{odd}}{}S(l)\right)=n/2.$$

In the last sum still many summands have $`\mathrm{ord}_2`$ bigger than $`n/2`$: if $`l3`$ mod $`4`$ then $`\mathrm{ord}(S(l))>n/2`$. On the other hand, if $`l1`$ mod $`4`$ then each summand $`b_lb_{k_2}\mathrm{}b_{k_n}`$ in $`S(l)`$ has $`\mathrm{ord}_2(b_lb_{k_2}\mathrm{}b_{k_n})=n/2`$. We conclude that $`\mathrm{ord}_2(\mathrm{Com}_n)=n/2`$ if and only if the number $`c(n)`$ of compositions of $`n1`$ into $`n`$ parts, where the first part is $`1`$ mod $`4`$ and the remaining $`n1`$ parts are even (zero parts are allowed), is odd. We have

$`c(n)`$ $`=`$ $`[x^{n1}]{\displaystyle \frac{x}{1x^4}}{\displaystyle \frac{1}{(1x^2)^{n1}}}=[x^{n1}]{\displaystyle \frac{x}{1+x^2}}{\displaystyle \frac{1}{(1x^2)^n}}`$
$``$ $`[x^{n1}]{\displaystyle \frac{x}{1x^2}}{\displaystyle \frac{1}{(1x^2)^n}}=[x^{n1}]{\displaystyle \frac{x}{(1x^2)^{n+1}}}\mathrm{mod}2`$
$`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{3m1}{m1}}\right){\displaystyle \frac{3m}{m}}\left({\displaystyle \genfrac{}{}{0pt}{}{3m1}{m1}}\right)\mathrm{mod}2`$
$`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{3m}{m}}\right).`$
It was noted by Kummer , see also Singmaster , that $`\mathrm{ord}_p(\left(\genfrac{}{}{0pt}{}{a+b}{b}\right))`$ is equal to the number of carries required when adding $`a`$ and $`b`$ in the $`p`$-ary notation. Applying this for $`p=2`$, $`a=m`$, and $`b=2m`$, we get the stated criterion.

###### Corollary 12

For all $`n3`$,

$$s_n\{\begin{array}{ccc}\hfill 2& \mathrm{mod}2^{(n1)/2}\hfill & \text{ for odd }n;\hfill \\ \hfill 2& \mathrm{mod}2^{n/2}\hfill & \text{ for even }n\text{.}\hfill \end{array}$$

Let

$$C_n=\frac{1}{n+1}\left(\genfrac{}{}{0pt}{}{2n}{n}\right)$$

be the $`n`$th Catalan number.

###### Proposition 13

For all $`n`$, $`\mathrm{Com}_nC_{n1}\mathrm{mod}3`$.

Proof: We have, for every non-negative integer $`k`$,

$$(2!x+3!x^2+\mathrm{})^k=(2x)^k+3a_k(x)$$

with $`a_k(x)𝐙[[x]]`$. Thus

$`{\displaystyle \underset{k0}{}}(1)^k(2!x+3!x^2+\mathrm{})^k`$ $`=`$ $`{\displaystyle \frac{1}{1+2x}}+3{\displaystyle \underset{k0}{}}(1)^ka_k(x)`$
$`=`$ $`{\displaystyle \frac{1}{1+2x}}+3b(x)`$
with $`b(x)𝐙[[x]]`$. Let $`m=\mathrm{ord}_3(n)`$. Since $`\mathrm{ord}_3(k!)k1`$ for every $`k`$ (Lemma 10), we have

$$\mathrm{ord}_3\left(3^k\left(\genfrac{}{}{0pt}{}{n}{k}\right)\right)m+1\text{for }k=1,2,\mathrm{},n\text{.}$$

By Lemma 9,

$`n\mathrm{Com}_n=[x^{n1}]\left({\displaystyle \frac{1}{1+2x}}+3b(x)\right)^n`$ $``$ $`[x^{n1}]{\displaystyle \frac{1}{(1+2x)^n}}\mathrm{mod}3^{m+1}`$
$`=`$ $`(2)^{n1}\left({\displaystyle \genfrac{}{}{0pt}{}{2n2}{n1}}\right).`$
Canceling in the last congruence the common factor $`3^m`$, we get

$$\frac{n}{3^m}\mathrm{Com}_n\frac{(2)^{n1}}{3^m}\left(\genfrac{}{}{0pt}{}{2n2}{n1}\right)\frac{1}{3^m}\left(\genfrac{}{}{0pt}{}{2n2}{n1}\right)\mathrm{mod}3.$$

Since $`n/3^m0\mathrm{mod}3`$, we can divide by it and get

$$\mathrm{Com}_n\frac{1}{n}\left(\genfrac{}{}{0pt}{}{2n2}{n1}\right)\mathrm{mod}3.$$

###### Corollary 14

For all $`n>2`$,

$$s_nC_{n1}+(1)^n\mathrm{mod}3.$$

## 5 Concluding remarks

The simplicity property for permutations does not seem to have been studied until very recently . We have begun the study of the numbers $`s_n`$ by showing that they are not P-recursive, giving the first few terms of their asymptotic expansion, and showing that they satisfy some unexpected congruence properties.

These results suggest a number of natural continuations. Although, in principle, we could obtain more terms of the asymptotic expansion the entire expansion remains elusive, and computing it seems to be rather a difficult problem. On the other hand we have some computational evidence to suggest that the sequence $`\mathrm{Com}_n`$ has additional congruence properties, particularly with respect to odd primes.

We suggest also some algorithmic problems that are natural counterparts to the enumerative results:

* How can one efficiently generate simple permutations in lexicographic order?
* Is it possible to generate simple permutations uniformly at random in worst-case linear time per permutation?
* How efficiently can one recognise a simple permutation?

With regards to the final question, there is a natural dynamic programming algorithm that achieves the task in $`O(n^2)`$ time; so the issue is whether one can do better.


text:
# Magnetism in the Hubbard model-An improved treatment

## I Introduction

While much progress has been made in recent years in the theoretical treatment of magnetism in the transition metals and similar materials, a complete treatment of the many aspects of these familiar but complicated systems is still lacking. The electrons retain many of their energy band properties as is reflected by their metallic behaviour, but the magnetic properties reveal that strong correlations exist between the electrons on the individual atoms, and exchange occurs between them.

The simplest treatment of metallic ferromagnets is due to Stoner E.C.Stoner (1938) who assumed that the electron interaction gave rise to a molecular field which predicted ferromagnetism in a simple band model if the ratio of the interaction to the band width (or more precisely the density of states of the Fermi surface) was adequately large. But beyond the mean field the model gives no correlation between the magnetic electrons. More complex band structures also predict antiferromagnetism in some cases D.R.Penn (1966).

In recent years much more sophisticated treatments of the band structure using established methods such as those based on the spin dependent local density approximation K.Terakura et al. (1984); A.Svane and O.Gunnarsson (1990); V.I.Anisimov et al. (1997); A.I.Lichtenstein and M.I.Katsnelson (1998) have been able to make much better predictions of the Fermi surface and related magnetic properties, but again a full treatment of the atomic electron correlations is lacking. Hybrid treatments such as that of Gutzwiller M.C.Gutzwiller (1963); F.Gebhard (1990); J.Bünemann et al. (1998) which start from a mixture of correlated and uncorrelated states can give an improved treatment of the magnetic properties at the expense of some of the metallic ones.

An alternative approach to this problem has been through the study of simple band models with electron interactions of which that due to Hubbard J.Hubbard (1963) is the simplest which appears to contain the essence of the problem. He restricts the electron-electron interaction to occur on single atomic sites which in the case of a single band only occurs between electrons of opposite spin because of the exclusion principle. Hubbard, in a series of papers J.Hubbard (1963, 1964), pioneered different approximate treatments to this model Hamiltonian. The simplest replaces the electron-electron interaction by its average so that an electron of one spin sees a field proportional to the density of the other spin and the system has simple energy bands separated by an energy proportional to the magnetisation and hence is completely equivalent to Stoner’s mean field approximation. In a different vocabulary it can be regarded as a virtual crystal approximation (VCA) where the electrons move in energy bands derived from an average potential.

In order to give a better treatment of the electron-electron scattering Hubbard J.Hubbard (1964) pioneered another approximation where it was assumed that electrons of one spin were frozen on their sites and the electrons of the other spin moved through the random potential so created. He proceeded to solve this problem by determining an effective potential (which was complex and energy dependent) which took account of the average scattering from this array in an analogy with an alloy. This method has been independently used in a wide variety of alloy problems and is usually referred to as the coherent potential approximation (CPA) M.Lax (1951); D.W.Taylor (1967); P.Soven (1967); R.J.Elliott et al. (1974). It is clear that this approximation over-estimates the electron-electron scattering and in particular causes the energy bands of each spin character to split if the interaction is large enough. As a result the tendency to ferromagnetism is under-estimated and it has been demonstrated J.Schneider and V.Drchal (1975); Harris and V.Lange (1967) that no magnetic order is predicted except in the special case of a half filled band.

Several attempts have been made to improve on these approximate treatments, for example using a dynamical CPA Y.Kakehashi (1992, 2002) based on the integration of this approximation into functional integral techniques. There has been a number of other attempts to modify the alloy analogy W.Nolting and W.Borgiel (1989); T.Hermann and W.Nolting (1996); M.Potthoff et al. (1998); W.Borgiel et al. (2001). In this work we propose a simple physical model which is intermediate between the straightforward VCA and CPA treatments. In it we assume that the random scattering centres experienced by an electron of one type in the alloy analogy exist only for a certain time $`\tau `$ which is proportional to a typical hopping time between sites and hence inversely to the band width. At later times the electrons see a completely average potential reflecting the movement of the scattering centres. Thus the propagators are calculated in the CPA approximation for the short time interval $`\tau `$ and have to be matched to the VCA propagators which hold for longer times. Although this is a straightforward procedure in principle certain further approximations are necessary to make it tractable. The properties of the system are evaluated in this so-called $`\tau `$-CPA. The single particle properties such as the density of states show, as expected, a mixture of characteristics associated with the VCA and CPA. Ferromagnetism is predicted to occur for a range of parameters which are somewhat more restrictive than those in the simple Stoner model.
This paper is organised as follows: Sec.II introduces the model used. The VCA and CPA methods are outlined in Sec.III. The $`\tau `$-CPA, mixing the VCA and CPA, is constructed in Sec.IV. Numerical results for the density of states and the magnetisation using a simple cubic band are presented and discussed in Sec.V.

## II Model Hamiltonian

The Hubbard model is currently widely considered to be the most concise model that captures the essence of the interactions in an electronic system. It allows the hopping of the electrons from site to site to compete with the Coulomb repulsion while including the Pauli principle. We will consider the simplest case of the Hubbard model, a single-band with nearest-neighbour hopping $`t_{ij}`$ and on-site repulsion :

$$H=\underset{H_0}{\underset{}{\underset{i,j,\sigma }{}t_{ij}c_{i\sigma }^+c_{j\sigma }}}+U\underset{i}{}c_i^+c_ic_i^+c_i$$
(1)
$`c_{i\sigma }^+`$ and $`c_{i\sigma }`$ are the usual creation and destruction operators for electrons in the Wannier states. The free part of the Hamiltonian $`H_0`$ is evaluated in the Bloch states $`H_0=_kϵ_kc_{k\sigma }^+c_{k\sigma }`$ using a tight-binding approximation for the band structure $`ϵ_k`$. We will concentrate in this paper on the simple cubic band of half band width $`w=1`$. As for the interaction, it is repulsive and only felt when one electron of spin $`\sigma `$ meets another of spin $`\sigma `$ at the same site. This model simulates relatively well a real material possessing narrow energy bands, like the transition metals and their alloys. We will work on a three dimensional lattice of $`N=L^3`$ sites, at $`T=0`$. We use the framework of the canonical ensemble and fix the number $`N_e`$ of electrons per atom that distribute themselves at random over the lattice sites. The Fermi statistics takes care that no site is occupied by two electrons of the same spin while the interaction favours singly occupied sites. The propagation of an electron of spin $`\sigma `$ in the lattice can be described by the retarded Green’s function

$`G_{ij}^\sigma (t,t^{})`$ $`=i\mathrm{\hspace{0.17em}2}\pi \theta (tt^{})\{c_{i\sigma }(t),c_{j\sigma }^+(t^{})\}_+`$
$`=:c_{i\sigma }(t);c_{j\sigma }^+(t^{})`$ (2)
where $`\mathrm{}`$ has been set to $`1`$. The following Fourier transform (where no particular origin is specified)

$$G^\sigma (E)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}d(tt^{})G^\sigma (tt^{})e^{iE(tt^{})}$$
(3)
will often be used, and $`t^{}`$ will be set to $`0`$.
We are interested in the properties of this system averaged over all possible configurations, for which the Green’s function $`G_k^\sigma (E)`$ can be defined by

$$G_{ij}^\sigma (E)=\frac{1}{N}\underset{k}{}G_k^\sigma (E)e^{ik(r_ir_j)}$$
(4)
The density of states $`\rho ^\sigma `$ is obtained from the single-particle Green’s function with

$$\rho ^\sigma (E)=\frac{1}{\pi }𝙸𝚖F^\sigma (E)$$
(5)
where $`F^\sigma (E)`$ is the trace of the Green’s function. It can be obtained as

$$F^\sigma (E)=\frac{1}{N}\underset{k}{}G_k^\sigma (E)$$
(6)
The relative proportions of the number of up-spins and down-spins must be determined self-consistently together with the chemical potential $`\mu `$. Defining $`N^{}`$ and $`N^{}`$ respectively as the average number of spins up and down per atom, the equilibrium conditions read

$$N^\sigma =_{\mathrm{}}^\mu 𝑑E\rho ^\sigma (E)$$
(7)
with $`\sigma \{,\}`$

## III VCA and CPA

The virtual-crystal approximation and the coherent-potential approximation have been both applied with various success to both the Hubbard model and the somewhat related problem of the random binary alloy $`A_{1c}B_c`$. It is worth recalling briefly the development of the VCA and the CPA in the binary alloy before turning to the Hubbard model.

### III.1 VCA and CPA in the binary alloy

The analogy between the Hubbard model and a binary alloy was noted by Hubbard in the last paper of the series J.Hubbard (1964). We consider a particle propagating on a lattice where each site is occupied at random either by a host atom $`A`$ or an impurity atom $`B`$. The concentration of impurities is $`c`$. The model is represented by $`H_{alloy}=H_o+V`$, where $`V`$ is a quadratic operator of the form $`V=_m|mϵ_mm|`$. The atomic potentials can take two values, either $`ϵ_A`$ or $`ϵ_B`$, depending on whether the particle is on an $`A`$ site or a $`B`$ site. Setting $`ϵ_A`$ to $`0`$ and $`ϵ_B`$ to $`U`$, we find that this is very similar to the Hubbard model (hence the alloy analogy), for which the spin $`\sigma `$ propagating in the lattice interacts with energy $`U`$ only on sites where a $`\sigma `$ spin is present. The analogy would be exact if the spins of opposite direction could be considered as frozen in the Hubbard model.
The VCA and CPA are well established methods B.Velický et al. (1968); R.J.Elliott et al. (1974) and will only be outlined here.
If $`P`$ describes the propagator in the perfect crystal defined by $`H_0`$, the perturbation expansion in the whole Hamiltonian $`H_{alloy}`$ reads

$$G=P+PVP+PVPVP+\mathrm{}=P+PVG$$
(8)
The problem of evaluating the configurational average $`G=P+PVG`$ is solved in the VCA by the random phase approximation decoupling $`VGVG`$. With $`ϵ_m=cU`$, the propagator for the binary alloy in the VCA is

$$G_k^{VCA}(E)=\frac{1}{Eϵ_kcU}$$
(9)
While the VCA is only a valid approximation for small perturbations $`V`$, the CPA is a particularly successful method that interpolates between the limits of strong and weak disorder and interaction. The CPA assumes the existence of an effective medium described by the Green’s function $`G_e`$ and self-energy $`\mathrm{\Sigma }`$ that obey the Dyson equation

$$G_e=P+P\mathrm{\Sigma }G_e$$
(10)
The true Green’s function $`G`$ can be developed in the perturbation expansion (8), which, together with (10), lead to a new expansion of $`G`$ in terms of the medium

$$G=G_e+G_e\left(V\mathrm{\Sigma }\right)G$$
(11)
A particle moving in the effective medium sees the impurities as embedded in a uniform potential. As shown in Refs.B.Velický et al., 1968; R.J.Elliott et al., 1974, it is also possible to expand the true $`G`$ in terms of $`G_e`$ and a scattering $`T`$ matrix so that

$$G=G_e+G_eTG_e$$
(12)
The $`T`$ matrix can be expressed as a function of the individual atomic scattering matrices $`T_m`$

$$T=\underset{m}{}T_m+\underset{nm}{}T_nG_eT_m+\underset{nml}{}T_nG_eT_mG_eT_l+\mathrm{}$$
(13)
Projected on the the Wannier states, the scattering on the $`m`$ atom can be found to be B.Velický et al. (1968); R.J.Elliott et al. (1974)

$$t_m=\frac{ϵ_m\mathrm{\Sigma }}{1(ϵ_m\mathrm{\Sigma })F}$$
(14)
$`F`$ is the diagonal element of $`G_e`$, $`F=m\left|G_e\right|m`$ (see 6). A central step of the CPA is to demand that the medium propagator $`G_e`$ corresponds to the true propagator averaged over the disorder: $`G^{CPA}:=GG_e`$. Consequently, from (12), the total scattering of impurities is set to zero: $`T=0`$. $`T`$ can then be decoupled by a random phase approximation on the single $`T_m`$ matrices. We note that this is a much better approximation than the VCA, since the repeated scattering by the same site is forbidden, as (13) shows. In this form, the CPA neglects the scattering from clusters of atoms, though some extensions of the theory do include this effect.
The CPA condition becomes therefore

$$t_m=0$$
(15)
This is the central equation of the CPA. The self-energy is calculated such that the scattering generated by any single impurity is zero on average. The CPA approximation for the one particle Green’s function is then given by

$$G_k^{CPA}(E)=\frac{1}{Eϵ_k\mathrm{\Sigma }(E)}$$
(16)
with the self-energy $`\mathrm{\Sigma }`$ determined self-consistently with $`F`$ from (15) and (14). From these two equations and for the binary allow with $`ϵ_A=0`$, $`ϵ_B=U`$, the relations between $`F`$ and $`\mathrm{\Sigma }`$ can be expressed as

$$\mathrm{\Sigma }=\frac{cU}{1(U\mathrm{\Sigma })F}$$
(17)
where $`F`$ is the diagonal element of $`G_k^{CPA}`$. While the VCA self-energy is a real number that merely shifts the energy levels by $`cU`$, the CPA $`\mathrm{\Sigma }`$ is an energy-dependant, complex number which accounts for the damping of the quasiparticle states by the impurities.

### III.2 VCA and CPA in the Hubbard model

The VCA for the Green’s function (II) in the Hubbard model (1) can be obtained in two different ways. The first approach consists in applying the Hartree-Fock approximation to the Hubbard model and effectively reducing it to a Stoner model

$$H_{Stoner}=H_0+U\underset{i}{}\left(N^{}\widehat{n}_i+N^{}\widehat{n}_i\right)$$
(18)
The usual number operator notation $`\widehat{n}_{i\sigma }=c_{i\sigma }^+c_{i\sigma }`$ has been used. The Stoner model (18) leads directly to the set of Green’s function analogous to (9)

$$G_k^{\sigma VCA}(E)=\frac{1}{Eϵ_kN^\sigma U}$$
(19)
The up-spins see the (frozen) down-spins as impurities of concentration $`N^{}`$ with whom they interact with energy $`U`$, and vice-versa for the down-spins.
The second way of obtaining the virtual-crystal approximation is to decouple the equation of motion for the Green’s function (II) in the Hubbard model

$`i{\displaystyle \frac{d}{dt}}c_{i\sigma }(t),c_{j\sigma }^+`$ $`=2\pi \delta (t)\delta _{ij}+{\displaystyle \underset{l}{}}t_{il}c_{l\sigma }(t),c_{j\sigma }^+`$
$`+U\widehat{n}_{i,\sigma }(t)c_{i\sigma }(t),c_{j\sigma }^+`$ (20)
If $`\widehat{n}_{i,\sigma }(t)`$ is taken out of the last Green’s function in (III.2) and averaged separately, $`\widehat{n}_{i,\sigma }(t)N^\sigma `$, the equation of motion is decoupled and the solution is given again by (19). If, on the other hand, $`\widehat{n}_{i,\sigma }(t)`$ is considered at $`t=0`$ and takes the value $`0`$ or $`1`$ at random, applying thus the alloy analogy, the CPA result (16,17) can be retrieved.

$$G_k^{\sigma CPA}(E)=\frac{1}{Eϵ_k\mathrm{\Sigma }^\sigma (E)}$$

$$\mathrm{\Sigma }^\sigma (E)=\frac{N^\sigma U}{1[U\mathrm{\Sigma }^\sigma (E)]F^\sigma (E)}$$
(21)
This is however not a straightforward procedure. Hubbard developed an elaborate scheme based on the alloy analogy in the last paper of the series J.Hubbard (1964), calling it the “scattering correction”. It was noted later by Velický et al that the scattering correction Green’s function is the same as the CPA result for the alloy model.

## IV The $`\tau `$-CPA

It has been known for some time now that the CPA does not allow for ferromagnetic solutions Harris and V.Lange (1967); J.Schneider and V.Drchal (1975). On the other hand, the VCA treatment of the Hubbard model leads to an artificially strong ferromagnetic phase. We propose to combine these two different approximations in order to keep the advantages of the CPA while allowing for spontaneous magnetic ordering to occur. As was seen in the previous section, the VCA and CPA take two very different views when considering the motion of a particle in the lattice. The CPA sees the scatterers as fixed at the positions given at $`t=0`$, while the VCA considers the scattering events only as averaged over a time $`t\mathrm{}`$. In other words, the CPA is valid for the Hubbard model only at $`t=0`$ and the VCA at $`t\mathrm{}`$. A better physical picture should emerge if the CPA treatment is applied for a short time only, when it is still reasonable to consider the scatterers as fixed. This time is typically the time needed for the particle to travel between sites. Thereafter we can assume that the VCA average provides a reasonably good treatment.

Let $`G_k^{VCA}(t)`$ and $`G_k^{CPA}(t)`$ be respectively the time Fourier transform of the VCA and CPA Green’s functions for $`t0`$. The central assumption is that up to a time $`\tau >0`$ and $`\tau 1/w`$ the system is described by the CPA Green’s function, whereas beyond $`\tau `$ it is described by the VCA Green’s function. A new function $`g_k^{\tau CPA}(t)`$ is defined as

$$g_k^{\tau CPA}(t)=i2\pi \theta (\tau t)G_k^{CPA}(t)\alpha i2\pi \theta (t\tau )G_k^{VCA}(t)$$
(22)
At $`t=0`$, $`g_k^{\tau CPA}`$ is equal to the CPA Green’s function and is thus properly normalised. A factor $`\alpha `$, the “matching factor”, has been introduced to ensure the continuity of the wave function at time $`\tau `$.
The retarded Green’s functions for the VCA and CPA are given respectively by $`G_k^{VCA}(E)`$ (9) and $`G_k^{CPA}(E)`$ (16). The Fourier transform reads:

$$G_k(t)=_{\mathrm{}}^{\mathrm{}}𝑑Ee^{iEt}G_k(E)$$
(23)
It is straight-forward to obtain $`G_k^{VCA}(t)`$ by contour integration in the lower half of the complex plane. There is one pole at $`z_0=ϵ_k+cUi\delta `$, so that

$$G_k^{VCA}(t)=i2\pi \theta (t)e^{i(ϵ_k+cU)t}$$
(24)
It is more difficult to do a similar calculation for the CPA. The problem lies in the nature of the CPA poles. It has been noted by Velický et al B.Velický et al. (1968) that the two complex CPA poles cannot in all generality be interpreted as quasiparticle energy. A calculation for the CPA similar to the one made for the VCA does not lead to a correct result. In particular, the correct inverse Fourier transform is not retrieved. In order to keep the denominator of the energy Green’s function as $`Eϵ_k\mathrm{\Sigma }(E)`$, we assume instead, by similarity with the VCA:

$$G_k^{CPA}(t)=i2\pi \theta (t)e^{i(ϵ_k+\mathrm{\Sigma })t}$$
(25)
At this stage the fact that $`\mathrm{\Sigma }`$ is an energy-dependant quantity is ignored, and $`\mathrm{\Sigma }`$ is treated as a fixed quantity rather than a function. By doing so, Fourier-transforming (25) gives the expected result (16). The substitution $`\mathrm{\Sigma }\mathrm{\Sigma }(E)`$ must be made thereafter. The whole $`g_k^{\tau CPA}(t)`$ is still normalised at $`t=0`$, but because of the dependence of $`\mathrm{\Sigma }`$ on the energy, this does not guarantee the normalisation of $`g_k^{\tau CPA}(E)`$. The normalisation will have to be corrected numerically.
The integrals for the Fourier transform on the partial intervals give, without taking into account the matching factor $`\alpha `$,

$$𝒢_k^{VCA}(E)=\frac{1}{2\pi }_0^\tau 𝑑te^{iEt}G_k^{VCA}(t)=\frac{e^{i\left(Eϵ_kcU\right)\tau }}{Eϵ_kcU}$$
(26)
$$𝒢_k^{CPA}(E)=\frac{1}{2\pi }_\tau ^{\mathrm{}}𝑑te^{iEt}G_k^{CPA}(t)=\frac{1e^{i\left(Eϵ_k\mathrm{\Sigma }\right)\tau }}{Eϵ_k\mathrm{\Sigma }}$$
(27)
The $`𝒢_k^{VCA}(E)`$ and $`𝒢_k^{CPA}(E)`$ are the standard VCA and CPA Green’s functions, weighted by a complex $`\tau `$-dependant quantity.
The factor $`\alpha `$ that forces the CPA and VCA Green’s functions to match at $`t=\tau `$ is found to be

$$\alpha =e^{i(\mathrm{\Sigma }cU)\tau }$$
(28)
This value for $`\alpha `$ can be reintroduced into (22). The result for $`g_k^{\tau CPA}(E)`$ is then obtained from the Fourier transform of (22):

$`g_k^{\tau CPA}(E)`$ $`={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}𝑑tG_k^{\tau CPA}(t)e^{iEt}`$
$`={\displaystyle \frac{1e^{i(Eϵ_k\mathrm{\Sigma })\tau }}{Eϵ_k\mathrm{\Sigma }}}+{\displaystyle \frac{e^{i(Eϵ_k\mathrm{\Sigma })\tau }}{Eϵ_kcU}}`$ (29)
$`\mathrm{\Sigma }`$ is now read as $`\mathrm{\Sigma }(E)`$, as calculated by the pure CPA equation (17). The expression (25) now reproduces exactly the CPA behaviour in the limit $`\tau \mathrm{}`$. The VCA approximation is retrieved in the limit $`\tau 0`$.

We recall that the normalisation of the Green’s function has been lost by choosing to ignore the energy dependence of $`\mathrm{\Sigma }`$. It is therefore necessary to calculate numerically for each $`k`$ the factor $`\beta _k`$, so that $`\beta _k^1g_k^{\tau CPA}(E)`$ is normalised.
The result (IV) is however not giving a physical density of states: $`\rho (E)`$ changes sign along the $`E`$ axis. By stopping the Green’s function at a time $`\tau `$ we have in fact introduced a non-physical oscillatory term coming from the factor weighing VCA and CPA

$$e^{i(\omega ϵ_k\mathrm{\Sigma })\tau }=e^{i(\omega ϵ_k𝚁𝚎\mathrm{\Sigma })\tau }e^{𝙸𝚖\mathrm{\Sigma }\tau }$$
(30)
The oscillatory term destroys the analytical properties of the CPA and VCA. The imaginary part of $`N^1_k\pi ^1g_k^{\tau CPA}(E+i\delta )`$ does not give the physical density of states, as the branch cut no longer runs along the x-axis. It is thus necessary to smooth out the Green’s function by putting the oscillatory term to unity. We obtain finally

$$G_k^{\tau CPA}(E)=\frac{1}{\beta _k}\left(\frac{1e^{𝙸𝚖\mathrm{\Sigma }(E)\tau }}{Eϵ_k\mathrm{\Sigma }(E)}+\frac{e^{𝙸𝚖\mathrm{\Sigma }(E)\tau }}{Eϵ_kcU}\right)$$
(31)
This function now possesses all the properties required. It interpolates between a pure VCA Green’s function at $`\tau 0`$ and the CPA Green’s function at $`\tau \mathrm{}`$. We observe that (31) is decomposed into a CPA and VCA contribution weighted by a real $`\tau `$-dependant factor. Although $`\tau `$ is not strictly determined by the problem, a reasonable value for it is given on physical ground by the band width.

## V Numerical results and discussion

In this section, the $`\tau `$-CPA Green’s function (31) is used to obtain the density of states for both species of spins in the Hubbard model. Various values of the impurity concentration $`c`$ ($`c`$ corresponds to either $`N^{}`$ or $`N^{}`$), of the interaction $`U`$ and of $`\tau `$ are considered. Calculations are made for a simple cubic band with $`w=1`$. Another one-particle property is then evaluated, the magnetisation $`m`$ versus the interaction strength $`U`$, for various range of parameters.

### V.1 Density of states in the $`\tau `$-CPA

The density of states is obtained from the Green’s function (31) using (5,6). Fig.1 shows the density of states obtained for a strong disorder and interaction strength, $`U=6.0`$ and $`c=0.4`$.

$`\tau `$ is varied from $`0`$ (pure VCA) to $`\tau \mathrm{}`$ (pure CPA). At $`\tau =0`$ the spectrum exhibits the single VCA band, centred at $`E=cU`$. As $`\tau `$ increases, the two CPA sub-bands appear, at $`E=0`$ for the main band and $`E=U`$ for the impurity band. Weight is transfered from the VCA band to both the CPA sub-bands, affecting all the three sub-band heights, but not their general shape. At $`\tau \mathrm{}`$ the VCA band has disappeared and only the two CPA sub-bands survive, with their original, pure CPA respective weight. In this graph the $`\tau `$-CPA merges with the pure CPA at $`\tau 600`$.
More dramatic effects on the band shapes occur at smaller $`U`$, when the bands mix. Fig.2, for $`U=0.8`$ and $`c=0.3`$, depicts a case for which $`U`$ and $`c`$ are sufficiently small so that the pure CPA sub-bands are not split. As $`\tau `$ increases, the height of the VCA band decreases while that of the main CPA sub-band increases. The general picture produced in that case is that for $`\tau >0`$, the degenerate states near the centre of the VCA band start to scatter, while the edges of the band are still untouched. For larger $`\tau `$, the states continue to diffuse to the impurity band from the centre to the right edge. At the same time, states of lower energy are being occupied, re-centring the created sub-band around $`0`$. Cases at intermediate $`U`$ are not so smooth, as can be seen in Figs.3(a) and 3(b), where $`U=2.0`$ and $`c=0.2`$. In this figure the density of states of the $`\tau `$-CPA becomes equal to that of the CPA for $`\tau 250`$. The mixing of the bands produces irregular features as the diminishing VCA meets the increasing main CPA band at larger $`\tau `$, Fig.3(b). This is due to the fact that a large imaginary part at an energy $`E`$ favours the CPA behaviour compared to that of the VCA. The CPA self-energy has an imaginary part only within the CPA bands. When the VCA band is centred within the main CPA band but is sufficiently distinct from it, its flat top can be sharpened if close to the top of the CPA band, producing the small peaks in Fig.3(b). Such drops in the VCA part of the density of states also occur for the parameters of Fig.2. However under these conditions the VCA band is only slightly shifted from the CPA band, and the sharp peaks are compensated by the CPA part.

### V.2 Magnetisation in the $`\tau `$-CPA

The $`\tau `$-CPA offers a treatment that leads to density of states that interpolate between the pure VCA and the pure CPA. As will be shown in this section, the method also interpolates the value of the magnetisation versus the interaction strength $`U`$ between the VCA ground state (ordered ground state for sufficiently large $`U`$) and the CPA ground state (no ordered ground state at any $`U`$). We only consider in this paper the competition between paramagnetic and ferromagnetic ground states.
In practice, the equilibrium conditions for a certain $`U`$ and a given, fixed number of spins $`N_e`$, are found by iteration. For example, a value for $`N^{}`$ is given as an input. This allows the determination of $`G^{}`$ (31) and the density of states from (5,6). $`N^{}`$, deduced from $`N_e=N^{}+N^{}`$, determines the chemical potential $`\mu `$ from (7) with $`\sigma =`$. $`N^{}`$ also generates $`G^{}`$, which together with the chemical potential, determine a new value for $`N^{}`$. The procedure is repeated until $`N^{}`$ converges. The paramagnetic solution $`N^{}=N^{}`$ is always a solution. If another solution is found, its stability must be tested against the paramagnetic case by measuring the total energy of both solutions. The total energy $`E`$ of the system is found by summing up the separate energies of the two spins J.Schneider and V.Drchal (1975); W.Nolting and W.Borgiel (1989)

$$E=\frac{1}{2}_{\mathrm{}}^\mu 𝑑\omega \underset{k,\sigma }{}(\omega +ϵ_k)A_\sigma (k,\omega )$$
(32)
where $`A_\sigma (k,\omega )`$ is the spectral density defined by $`A_\sigma (k,\omega )=\pi ^1𝙸𝚖G_k^\sigma (\omega )`$. Diagrams of magnetisation $`m=N^{}N^{}`$ versus interaction strength $`U`$ can be calculated for any value of $`N_e`$. In Fig.4 the magnetisation is depicted for a range of values of $`U`$.

$`\tau `$ is fixed to $`1`$ and the curves calculated for various filling. The first observation is that at no value of $`N_e`$ ($`N_e<1`$) is there saturated ferromagnetism. This is due to the residual CPA part around $`E=0`$ that is never totally empty. We can also see that at $`N_e0.4`$ the maximum of the magnetisation is not at large $`U`$: $`m`$ increases at around $`U8`$ to reach its peak at $`U2`$. Keeping the filling constant ($`N_e=0.4`$ and $`N_e=0.7`$ respectively), while varying $`\tau `$ give the plots of Figs.5 and 6. The magnetisation is interpolated between its VCA value at $`\tau =0`$ and the CPA value at $`\tau =\mathrm{}`$. At $`\tau =0`$ and for the fillings considered, the system is fully polarised at large enough $`U`$. The magnetisation has already dropped to zero for any $`U`$ at $`\tau =10`$ for $`N_e=0.4`$, and at $`\tau =50`$ for $`N_e=0.7`$. The increase in $`m`$ before the drop in magnetisation is particularly acute for $`N_e=0.4`$ and $`\tau =2.0`$, and is an artifact of the model caused by the position of the chemical potential in the VCA and CPA sub-bands. Fig.7 shows the respective positions of the CPA and the VCA sub-bands and the chemical potential for $`N_e=0.4`$, compared with $`N_e=0.7`$ at $`\tau =2.0`$.

At large $`U`$, the density of states of the up-spins (the minority spins) exhibits three well separated sub-bands for both fillings $`N_e=0.4`$ and $`0.7`$. The lower CPA sub-band is completely full. As $`U`$ decreases, the VCA band gets closer to the lower CPA sub-band, in both the up- and down-spin DOS. At first, only the down-spin bands are affected, meaning that the chemical potential travels down the energy scale while the proportion of up-spins and down-spins remains more or less constant. This is the situation illustrated in the two uppermost plots of Fig.7. When $`U`$ is further decreased, the chemical potential will eventually reach the edge of the CPA lower sub-band. At a filling of $`N_e=0.4`$, this is happening at a fairly large $`U`$, $`U8.0`$, so that the VCA band is still well separated (plot in the middle of the left-hand side, Fig.7). $`N`$ will thus decrease until the VCA band merges with the CPA band (bottom plot on the left, Fig.7). For a higher filling, the two sub-bands merge before or at the same time as the chemical potential gets into the CPA band. The density of states at the chemical potential thus rapidly increases for the up-spins, and the magnetisation decreases accordingly as the up-spin band matches the down-spin band.
The increase in $`m`$ at intermediate $`U`$, though also occurring at other small fillings, becomes dramatic around $`N_e0.4`$ and $`\tau 2.0`$. It is at this order of magnitude of $`\tau `$ that the system, at $`N_e0.4`$ , exhibits a balanced CPA-VCA behaviour. At much smaller $`\tau `$, the relative weight of the CPA band becomes very small compared to the VCA part. The chemical potential falls within the CPA sub-band of the up-spin over a large range of $`U`$ before the up-spin VCA band reaches it. At large $`\tau `$ the system becomes more CPA in character with a small VCA part responsible for the low magnetisation.

## VI Conclusion

A new treatment of the Hubbard model is proposed that combines the physical pictures and respective advantages of the CPA and the VCA. It is based on the fact that the CPA on the Hubbard model considers the scatterers as fixed at their positions at $`t=0`$, whereas the VCA performs a time-average of the system. The one-particle Green’s function in the so-called $`\tau `$-CPA is built from the one-particle time Green’s functions of the CPA and VCA. The CPA Green’s function is applied up to a time $`\tau `$, thereafter the VCA Green’s function, multiplied by a factor for continuity, is used. This time $`\tau `$ is taken to be of the order of magnitude of the time the spins travel between sites, that is $`\tau `$ has the value of the inverse of the band width. The resulting Green’s function, which must be first smoothed out and normalised, interpolates between the VCA propagator at $`\tau =0`$ and the CPA propagator at $`\tau =\mathrm{}`$. The density of states at finite $`\tau `$ exhibits both the CPA sub-bands (or one CPA band, if $`U`$ is small) and a VCA band. The bands are weighted according to the value of $`\tau `$. The CPA sub-bands are fixed at their positions at $`E=0`$ and $`E=U`$, while the VCA band, centred at $`E=cU`$, moves along with the concentration $`c`$. Equilibrium conditions with a polarised ground state are therefore found that are more stringent than the Stoner criterion for finite $`\tau `$.
Artifacts of the $`\tau `$-CPA mixture include some non-regular features in the density of states for a range of parameters (intermediate $`U`$ and concentration), and peculiar behaviour of the $`m`$ versus $`U`$ curves at large $`\tau `$ and small concentration. The $`\tau `$-CPA is however a successful method for the Hubbard model that is based on the better physical picture given by the combination of the CPA and VCA. It allows the treatment of the scattering by the CPA while also predicting ordered ground state for large enough $`U`$. In a next publication, it will be shown that the $`\tau `$-CPA can also be used to calculate two-particle Green’s functions. In particular the dynamical susceptibility in the Hubbard model can be evaluated using this method.

###### Acknowledgements.

AU acknowledges the support of the Swiss National Science Foundation, the Berrow scholarship trust of Lincoln College and the ORS award scheme.