text:
# Gas Dynamics in Central Parts of Galaxies

## 1 Introduction

Nuclear starbursts and Active Galactic Nuclei (AGN) have to be fed by mass transport into the active regions. Due to dissipation, gas inflow is easiest to trigger. There are two obvious dynamical mechanisms causing gas inflow in galaxies: interactions and asymmetries in galactic potentials. The first one involves violent and transient phenomena with strong gas inflow towards the centres of the interacting galaxies seen in numerical simulations (see e.g. Mihos & Hernquist 1996). Closer analysis of the modeled merging sequence shows that the inflow process has two stages: first, tidal perturbations destabilize stellar discs, and then, each disturbed stellar disc exerts torques on the gas, and triggers inflow. Thus even in systems apparently as order-less as merging galaxies, gas dynamics is governed by the perturbed potentials of each individual galaxy.

Here I review the dependence of gas dynamics on the gravitational potential in isolated galaxies, where it can be seen more clearly, since these galaxies remain roughly unchanged on timescales corresponding to the formation of steady gas flow patterns. Inflow in isolated galaxies can occur steadily, which may in principle provide continuous feeding of the central activity. However, the presence and the positions of dynamical resonances determine the mode of gas flow in the nuclear region: it may be quite complex, as the recent high-resolution observations of galactic centres indicate. Structured dust lanes, nested bars, as well as nuclear rings and spirals have been observed (Regan & Mulchaey 1999, Martini & Pogge 1999, Pogge & Martini 2002, Erwin & Sparke 2002).

First, I briefly review global gas flow patterns in the most common asymmetry in disc galaxies, i.e. in the presence of a stellar bar (§2), and then, in subsections of §3, I explore various modes of gas flow in the nuclear region. I also attempt to determine which modes may trigger significant inflow capable of feeding the AGN. Placing a Massive Black Hole (MBH) at the galaxy centre may add, remove or shift resonances, therefore such a black hole may regulate gas flow around itself. Finally, I make a few remarks on numerical methods (§4).

## 2 Gas dynamics in a barred galaxy outside the nuclear ring

Galactic discs commonly develop stellar bars, and the majority of observed galaxies is barred (e.g. Eskridge et al. 2000). The most established observational notion about gas morphology in barred galaxies is the existence of two straight or curved dust lanes in the main body of the stellar bar, which are symmetric around the galaxy centre, and slightly tilted to the bar’s major axis (see e.g. NGC 1097, NGC 1300, NGC 4303, and NGC 6951 in the Hubble Atlas of Galaxies, Sandage 1961). Strong discontinuities in gas velocity, indicative of shocks, have been measured along the dust lanes (Regan et al. 1997). This morphology and kinematics have been seen in numerical models of gas flow in a fixed potential of the bar since the work of Sørensen et al. (1976). It has been explained in terms of stable periodic orbits by Athanassoula (1992).

There are two major families of stable prograde periodic orbits inside the bar: the $`x_1`$ orbits parallel to the bar, and the $`x_2`$ orbits orthogonal to it (Contopoulos & Papayannopoulos 1980). The $`x_2`$ orbits are present only in galaxies with sufficient central mass concentrations and moderate bar pattern speeds — in linear analysis this condition corresponds to the existence of an Inner Lindblad Resonance (ILR, see Binney & Tremaine 1987, and §3). By definition, the $`x_2`$ orbits exist only inside this resonance. Representative orbits of each family are plotted in Figure 1.

In the outer parts of the bar, gas tends to follow the $`x_1`$ orbits. Further inwards, the $`x_1`$ orbits develop cusps and loops: stars can follow such orbits, but gas on self-intersecting orbits gets shocked. A robust feature in the majority of numerical models is a pair of principal shocks that develop on the leading side of the bar (marked PS in Fig.1). There, torques transfer angular momentum from the gas to the bar, and the gas falls towards the centre. Numerical models indicate considerable gas inflow in the principal shocks: up to 1 M$`_{}`$yr<sup>-1</sup> for typical gas densities in the disc (Athanassoula 1992; Regan et al. 1997). The shocked gas usually does not flow straight to the galaxy centre: it tends to settle on the $`x_2`$ orbits, which have lower value of the Jacobi integral $`E_J`$ (corresponding to the energy in the rotating frame). Consequently, the straight shocks are slightly tilted to the bar’s major axis. If there were no $`x_2`$ orbits, straight shocks would develop on the bar’s major axis, possibly extending directly to the centre of the galaxy (Athanassoula 1992, Piner et al. 1995). Only a couple of galaxies with such central (on-axis) straight shocks have been observed (Athanassoula 1992), which suggests that most barred galaxies possess the ILR. The gas settling on $`x_2`$ orbits, well inside the ILR (see Fig.1), often accumulates in nuclear rings, where high density and no shear make favorable conditions for star formation (Piner et al. 1995). Star-forming nuclear rings have been observed in many galaxies (e.g. NGC 1512, Maoz et al. 2001; NGC 4314, Benedict et al. 2002). Nuclear rings may also form in a potential lacking an ILR, given that this potential is almost axisymmetric at the nucleus (Wada, priv. comm.).

## 3 Gas dynamics around resonances inside the nuclear ring

At resonances, the orbital frequency (Doppler-shifted to the reference frame of the bar) is commensurate with the frequency of radial oscillations. Thus a particle (or a gas cloud) at the resonance is subject to a monotonic force, which moves it away from there. The influence of the ILR is the strongest: it creates the principal shocks described in §2. For a bar rotating with a pattern speed $`\mathrm{\Omega }_B`$, the ILR occurs where the orbital frequency $`\mathrm{\Omega }`$ and the frequency of radial oscillations $`\kappa `$ are related by $`\mathrm{\Omega }\mathrm{\Omega }_B=\kappa /2`$. Since $`\mathrm{\Omega }_B=`$const, it is convenient to rewrite this condition as $`\mathrm{\Omega }\kappa /2=\mathrm{\Omega }_B`$, where $`\mathrm{\Omega }\kappa /2`$ is called the frequency curve.

### 3.1 The Triple ILR scenario

The frequency curve always declines towards zero at large radii. It goes to $`+\mathrm{}`$ at the galactic centre if a central MBH is present, but even without a central MBH it can reach high values (larger than $`\mathrm{\Omega }_B`$) at the nucleus, if the galaxy mass is sufficiently centrally concentrated. Thus, in such galaxies one should expect at least one ILR. However, if a galaxy has a central MBH, but a low central mass concentration otherwise, the resonance curve may not be monotonic, and 3 ILRs can occur (Fig.2). They are called (from the innermost out) the nuclear ILR, the inner ILR and the outer ILR.

The nuclear ILR is essentially related to the presence of the central MBH. Fukuda et al. (1998, 2000) showed that gas flow around this resonance closely resembles the flow around the outer ILR described in §2. This gas morphology may prove essential in finding central MBHs in galaxies: as can be seen in Figure 2, the nuclear ILR is located at the radius a few times larger than the MBH’s radius of influence (minimum of the rotation curve). Therefore we can learn about the central MBH from the gas morphology at radii larger than the ones useful in the conventional methods of stellar kinematics. On the other hand, shocks around the nuclear ILR strongly disturb gas kinematics; in this case, the assumption of a circularly rotating disc, common when MBH mass is derived from kinematics of ionized gas (e.g. Macchetto et al. 1997), may no longer be valid.

Because of the small size of the ring inside the nuclear ILR, the gas accumulating there may reach very high densities, so that self-gravity in gas will be the dominating force. Fukuda et al. (2000) investigated the evolution of such a self-gravitating ring, and showed that strong inflow induced by clumps in self-gravitating gas will occur. This is one of the viable mechanisms of feeding the AGN; here the central MBH enhances inflow onto itself due to the resonance that it generates. Another possible scenario involves efficient star formation in the ring, which leads to a nuclear starburst driven by the MBH.

Gas dynamics at the inner ILR is considerably different from that at the nuclear ILR and at the outer ILR. The orbital transition is different there as well: at the outer ILR, gas moves from the $`x_1`$ orbits with higher Jacobi integral $`E_J`$ to the nearly-circular $`x_2`$ ones with lower $`E_J`$. Inside the inner ILR, gas has to populate the $`x_1`$ orbits, which may be highly elongated at the centre (see Fig.1). The transition from the flow along the $`x_2`$ orbits to the one following the $`x_1`$ family creates a leading spiral at the inner ILR (Wada 1994, Knapen et al. 1995). The dynamics of this spiral is markedly different from the trailing spirals that form at the nuclear and outer ILRs: it reflects a converging flow rather than a shock, with the hot gas upstream from the arm (Wada & Koda 2001).

### 3.2 Central rotation curves and mass profiles in galactic centres

The triple ILR scenario predicts reversal of the spiral structure in the central parts of galaxies. Going inwards from the trailing spiral arms and the principal shocks in the bar, one should encounter the leading spiral at the inner ILR, and then a trailing one again in the innermost regions. This is not a commonly observed structure (§3.3), which undermines the viability of the frequency curve in Figure 2. Furthermore, fits to surface brightness profiles of disc galaxy nuclei show a very good agreement with power-law profile (Balcells et al. 2001). For a power-law luminosity profile, the rotation curve is also a power law (with positive exponent), as is the frequency curve (negative exponent here), which means that the latter rises monotonously inwards, and therefore has one ILR only.

The frequency curve is a derivative of the rotation curve, and is therefore more sensitive to the variation of the rotation curve’s shape than to its value. Consider two rotation curves (Fig.3), both flat at large radii, but one a power-law in the inner part, while the other one has a linear inner rise. Although they both may fit a given set of observed data equally well, they generate diametrically different frequency curves. The rotation curve with an inner power-law results in an inner power-law frequency curve, which rises monotonically inwards (innermost radii in left panel of Fig.3). This implies one ILR only. On the other hand, a linearly rising rotation curve indicates solid body rotation, for which the frequency curve is identically equal to zero (innermost radii in right panel of Fig.3). This forces an additional ILR on the system, and now both outer ILR and inner ILR are present. Therefore, a linear rise of the central rotation curve, although a good approximation to the data, may not properly reflect the mass distribution in the inner parts of galaxies, or the resonances that it induces. Special care has to be taken when the presence and the position of resonances are deduced from the observed rotation curve. Resolution effects, like beam smearing, give the impression of a linear rise in the central parts, which unavoidably implies the presence of an inner ILR.

### 3.3 Gas density waves within a single ILR

Recent high-resolution maps of galactic centres reveal intricate dust structures, which are often organized in a spiral pattern (Regan & Mulchaey 1999, Pogge & Martini 2002). In disc galaxies, these nuclear spirals, often winding by more than a $`2\pi `$ angle, seem to connect smoothly to the outer spiral arms or to the principal dust lanes in the bar. This makes an indirect argument against the presence of an inner ILR, since theory and modeling indicate the reversal of the spiral there (§3.1). On the other hand, the standard model of gas flow in the bar (§2) does not account for nuclear spirals either: gas settles there in the nuclear ring. Models by Athanassoula (1992) show curling of the inner parts of the principal shock, but the resolution is too low there to follow this feature further inwards.

For a given potential, gas flow depends sensitively on one parameter only: the velocity dispersion in the gas, which for an isothermal gas is represented by its temperature (Englmaier & Gerhard 1997; Patsis & Athanassoula 2000). Gas dynamics in the nuclear region changes remarkably with the sound speed in the isothermal gas ($`c_S`$): cold gas ($`c_S5`$ km/s) settles on the nuclear ring, as described in §2, while hot gas ($`c_S>10`$ km/s) forms a nuclear spiral (Maciejewski 1998). In the evolutionary sequence of gas flow in the bar, the principal shocks form first as trailing spirals, and then straighten as the bar grows in strength, with their inner parts curling inwards (Fig.4). In hot gas, the curling shock from one side of the bar never meets its counterpart from the opposite side, and both of them propagate towards the centre as an $`m`$=2 spiral mode. The convergence of the flow in such created nuclear spiral is high (large negative $`div𝐯`$), which is indicative of a shock. Significant gas inflow to the galactic centre may occur if the spiral shock proves to be efficient in transferring angular momentum from the gas to the bar. On the other hand, the principal shock in cold gas ceases to propagate inwards. Gas emerging from its innermost parts remains on a relatively open trajectory (Fig.4), and hits the shock on the opposite side of the bar. Hydrodynamical interaction of the flow with the shock results here in gas settling into the nuclear ring.

A good grasp of the dynamics of nuclear spirals can be gained when one considers them as density waves generated by gravitational torques in a non-selfgravitating gas (Englmaier & Shlosman 2000). In the linear approximation, the dispersion relation for the $`m`$=2 mode allows for wave propagation inside the corotation only when $`\mathrm{\Omega }\kappa /2>\mathrm{\Omega }_B`$, i.e. between the outer ILR and the inner ILR (if present) or the galactic centre. Since nuclear spirals are observed in disc galaxies down to the HST resolution limit, the inner ILR must reside inside that radius in these galaxies, if it is present at all.

Although some nuclear spirals in strongly barred galaxies may be well out of the linear regime, the density-wave mechanism applies to a weak spiral inside the nuclear ring seen in the models by Maciejewski et al. (2002). It is also extremely efficient in very weak bars, where gas dynamics closely follows the linear theory. Its predictions are the same for any weak nonaxisymmetry of the potential, including lopsidedness, because low harmonics of the mass distribution play a role there. Maciejewski (2001) studied gas flow in the gravitational potential with the quadruple moment of the bar 10 times lower, and the bar’s axis ratio half that of the standard model presented in §2. Such a bar is likely to be too weak to be detected observationally, and it does not generate principal shocks of any shape. Nevertheless, a clear spiral density enhancement is seen inside the ILR. In accordance with the density-wave theory, it extends all the way to the galactic centre if there is no inner ILR, or ceases at the inner ILR if it is present. Unlike nuclear spirals in strong bars, which unwind outwards rapidly in order to match the principal shock in the bar, nuclear spirals in very weak bars follow the linear mode longer, and wind up to 3 times around the centre (6$`\pi `$ angle). Thus tightly wound nuclear spirals may be observationally associated with galaxies where no bar has been detected.

### 3.4 Sorting out distinct scenarios of gas flow in nested bars

Galaxies possessing nested bars (a secondary bar inside the main one) have been noticed already almost 30 years ago (de Vaucouleurs 1974), but only recent studies (Laine et al. 2002, Erwin & Sparke 2002) showed them to be quite a common phenomenon: up to 30% of barred galaxies contain nested stellar bars. The random relative orientation of the bars suggests that they rotate independently. For a while nested bars became a favourite mechanism to fuel the nucleus: if gas flow in the secondary bar is analogous to the one in the main bar, then the secondary bar may force gas deeper into the potential well. Meanwhile, the pioneering work of Shlosman et al. (1989) examined fueling mechanism that involves instability in the self-gravitating nuclear gaseous disc. As the gas inflow in the large bar stagnates at the nuclear ring or disc, gas accumulates there, and can become gravitationally unstable. By analogy to instability of the galactic stellar disc, Shlosman et al. proposed that a gaseous inner bar forms inside the large bar (gaseous bar scenario).

The gaseous bar scenario may be related to the observed nested stellar bars, when star formation occurs in the dense gaseous inner bar, or when the inner gaseous bar drags stars of the nuclear component (Combes 1994). Neither of these links is likely though, since the gaseous bar scenario relies on global instability in self-gravitating gas, while cold massive nuclear discs are more likely to become unstable locally, and form clumps. This is clearly seen in simulations by Heller & Shlosman (1994): although the nuclear gaseous disc first collapses into a bar-like shape, clumps develop very quickly, so that no bar is seen a few Myr after its formation. This time is too short to form stars or to alter stellar dynamics, and the star-forming clumps do not follow bar-like orbits either. Therefore nested stellar bars are unlikely to be produced by the gaseous bar scenario. However, this scenario may lead to an episodic gas inflow: one galaxy where it may operate is Circinus (Maiolino et al. 2000) — one of the three closest AGN. A strong gas inflow on a scale $``$50 pc is seen on one side of the nucleus, but it is not accompanied by any outflow. It is likely that the unstable gaseous disc developed an $`m=1`$ instability mode there, and is feeding the nucleus in a short episode.

In order to estimate the dynamical significance of the gaseous component Friedli & Martinet (1993) and Combes (1994) looked for formation of nested bars in numerical simulations involving stellar and gaseous particles. They were able to create secondary stellar bars, which decouple from the large bars after the central mass concentration increases substantially as a result of gas inflow in the large bar. These inner stellar bars can then drag gas with them. In this mechanism, the role of gas is to increase central mass concentration, and to prevent heating of the central stellar component. No gravitational instability in the gaseous disc is involved, contrary to the mechanism proposed by Shlosman et al. (1989). Note however that these simulations use numerical methods suffering from excessive gas inflow to the centre (see §4), and therefore central mass concentrations may be unrealistically high, which undermines their conclusion. Also, the size ratio of the bars produced in this mechanism (Set III in Friedli & Martinet 1993) is 2 — much lower than the observed values. It indicates that these models may miss the physical mechanisms operating in the observed systems. A more realistic size ratio (3.8) is produced in Set II of Friedli & Martinet, but there the gas only prevents stars from heating, and the initial conditions are particularly chosen so that two bars are formed.

Since inner bars are often observed in near infrared, they are most likely made out of old stellar populations, and are long-lived. They are often devoid of gas, hence it is likely that purely stellar-dynamical mechanisms forming decoupled nested bars are in operation there (Rautiainen & Salo 1999). Maciejewski & Sparke (2000) explored limitations which are imposed on nested stellar bars by their orbital structure. Since orbits in a potential of two independently rotating bars (a class of pulsating potentials) do not close in general, they generalized the concept of the orbit to pulsating potentials, and found sets of particles populating closed curves that support bars in their motion. If resonant coupling is assumed to minimize the number of chaotic zones (Sygnet et al. 1988), then curves hosting sets of particles that support the inner bar are round and end well within its corotation, so that no shocks in gas are induced. Maciejewski et al. (2002) confirmed this with hydrodynamical simulations: gas flow in the inner stellar bar, which is in resonant coupling with the large bar, does not develop straight shocks, and does not increase the inflow. It instead organizes itself in an ellipsoidal pattern around the inner bar, like the one seen in color maps of NGC 3081 (Regan & Mulchaey 1999). Nevertheless, nested bars may not follow the resonant coupling assumed in this scenario, and other modes of gas flow in the inner bar are possible, including straight dust lanes (Maciejewski 2002).

## 4 Remarks about numerical methods

Numerical methods used in modeling of gas flow in galaxies can generally be divided into 3 classes: one treats gas as a set of ballistic clouds with some rules for energy and momentum exchange (sometimes called ’sticky particles’), the second solves hydrodynamical equations involving the SPH algorithm (Smooth Particle Hydrodynamics), while the third solves these equations directly on the grid. The first method, and to some extent the second, involves a number of free parameters in gas description, which can be adjusted to describe particular features of the ISM. On the other hand, the third method may involve algorithms (like the Riemann solver) which limit the set of free parameters to purely numerical ones, and which attempt to give the best approximation to the dynamics of non-viscous gas. This leads to unique hydrodynamical solutions, but it remains unclear how the real ISM differs from non-viscous fluid.

Consequences of numerical viscosity are particularly significant in studies of gas inflow. Bulk viscosity in SPH-based codes is necessary to stabilize post-shock oscillations, and it is quite large. One can compare the evolution of density within the inner 300 pc (inside the nuclear ring) in the grid-based model (Piner et al. 1995) to that in the SPH model (Patsis & Athanassoula 2000) for the same potential: in the first model, the density remains roughly constant, while in the second one it increases by more than a factor of 50 after 8 pattern rotations. Such a central mass concentration may be not realistic, and its strong dynamical influence, postulated by Friedli & Martinet (1993) and others, may not be relevant to present-day galaxies. Also the importance of self-gravity in gas may need to be scaled down. Heller & Shlosman (1994) measure mass of gas within 500 pc at some stage of their SPH simulations to be half of the dynamical mass. This is significantly more than the observed values (see e.g. Sakamoto et al. 1999). They also measure the mass of the central gas clump to be $`3\times 10^9`$M$`_{}`$, when usually no more than $`10^9`$M$`_{}`$of gas is seen in galaxy centres. In the SPH simulations by Berentzen et al. (1998) a considerable amount of gas quickly falls onto the ’central accreting object’, so that mass in gas accumulated there amounts to 1.6% of the total galaxy mass. This is an order of magnitude larger than the currently measured central MBH masses (Kormendy & Gebhardt 2001).

All 3 classes of methods discussed above usually represent the ISM as an isothermal fluid. Although this is a valid approximation (Cowie 1980), which has produced many significant results so far, one will eventually need a more realistic description of the ISM. Immense progress has been made here in the recent years by Wada & Norman (1999, 2001), who developed a grid-based code which allows for cooling and heating of the gas, as well as for self-gravity. As a result, a multiphase medium forms, with lognormal density distribution over 4 orders of magnitude. Shear by differential rotation seems to be sufficient to create flocculent nuclear spirals, and dynamics of this medium in presence of a weak bar has been studied recently (Wada & Koda 2001).

## 5 Conclusions

Asymmetries in the galactic potential can alter gas dynamics, generate inflow, and consequently speed up galaxy evolution. Gas dynamics in centres of galaxies strongly depends on the presence of resonances created by asymmetries in the potential. The presence and the number of resonances is regulated by the central mass distribution. In particular, the MBH in the centre of galaxy can create and/or move the resonances, and thus regulate gas flow around itself. Grand design nuclear spirals observed in centres of galaxies are consistent with the presence of one ILR only. Some of them may be sites of shocks that increase inflow to the galactic centre.

Two fundamentally different gas flows have been labeled as ’gas flow in nested bars’, which often causes confusion. One mode of the flow, originally proposed by Shlosman et al. (1989), refers to gaseous inner bars, which form in unstable gaseous nuclear discs. In this scenario the bar is the flow: it involves rapid radial gas motion on timescales of the order of dynamical time. It may be able to feed the AGN, but modeling work was not able yet to find connection between it and the nested stellar bars that we observe in galaxies. The second scenario refers to the flow in nested stellar bars. Stable orbits supporting such systems enable their long lifetime, and at the same time they restrict gas flow to follow elliptical trajectories (Maciejewski & Sparke 2000). In this scenario nested bars do not increase inflow to the galactic centre.

Acknowledgments. I would like to thank the organizers of the Galactic Dynamics Workshop at the JENAM 2002 Symposium for their invitation to give this review. Insightful comments by Lia Athanassoula, Eric Emsellem, Peter Erwin, Roberto Maiolino, Alessandro Marconi, Eva Schinnerer, Peter Teuben and Keiichi Wada have been invaluable to this review.


format: text

text:
# THE OPTICAL TO X-RAY SPECTRUM OF THE BROAD-LINE ULTRASOFT AGN RE J2248-511

## INTRODUCTION

Studies of ROSAT PSPC slopes in AGN have concluded that sources with steep soft X-ray (0.1-2.4 keV) continuum slopes in addition to broad optical emission lines (H$`\beta `$ FWHM$``$3000 km s-1) are simply not found in nature (Boller, Brandt & Fink 1996; Wandel & Boller 1998). The steep soft X-ray slopes characteristic of Narrow-line Seyfert 1 galaxies (NLS1s; whose H$`\beta `$ FWHM$``$2000 km s-1) are thought to be due to a small mass black hole with a high accretion rate (eg Pounds, Done & Osborne 1995) or perhaps due to a broad line region (BLR) which forms relatively far out from the central black hole (Puchnarewicz et al. 1992).

RE J2248-511 is a radio-quiet EUV-selected AGN ($`z`$ = 0.101) discovered by the ROSAT Wide Field Camera (Pounds et al. 1993), and classified as a normal broad-line Seyfert 1 galaxy (Mason et al. 1995) with H$`\beta `$ FWHM$``$3000 km s-1. A subsequent observation with the ROSAT PSPC (Puchnarewicz et al. 1995) found that the soft X-ray spectrum was very soft and was best-fitted with a broken power law where $`\mathrm{\Gamma }`$ = 4.13 at energies below 0.25 keV (and $`\mathrm{\Gamma }`$ = 2.62 at higher energies). Thus it differs from a NLS1 in that it has high velocity BLR clouds, yet a very steep, ultrasoft X-ray spectrum. This puts it in a category comprising, to our knowledge, only 3 AGN with both these characteristics (RE J2248-511, 1H 0419-577 and RXJ0437.4-4711).

RE J2248-511 was also observed by ASCA in 1997 (Breeveld et al. 2001). These observations showed that the soft excess (below 2.0 keV) had completely disappeared and the entire 0.5-10 keV spectrum could be represented by a single power law. Large amplitude variability in the soft X-ray region has also been observed in another of the broad-line ultrasoft AGN 1H 0419-577 (Guainazzi et al. 1998, Turner et al. 1999, Page et al. 2002). This two-state behaviour has led to comparisons with Galactic black holes (GBHs). Page et al. (2001) show that the 0.2-10 keV spectrum of 1H 0419-577 closely resembles that of two Comptonised components, although these cannot be well constrained, and that the variability between XMM-Newton and ASCA is explained by a change in optical depth of the Comptonising material. Comptonisation has also been proposed as a possible source for the X-ray emission of RE J2248-511 (Breeveld et al. 2001) which we will investigate in this paper.

Here, we present the results of a recent multi-wavelength monitoring campaign consisting of observations from XMM-Newton and the South African Astronomical Observatory (SAAO), to examine the long-term variability of this source and the relationship between the X-ray emitting components and the response of the BLR.

## XMM-NEWTON OBSERVATIONS

RE J2248-511 was observed by the EPIC (Strüder et al. 2001; Turner et al. 2001) and RGS (den Herder et al. 2001) instruments on-board XMM-Newton (Jansen et al. 2001) on 2000 October 26 (17.6 ks pn), 2001 April 21 ($``$8.5 ks) and 2001 October 31 (15.4 ks pn). The Optical Monitor was switched off since a nearby bright star saturates the field. All EPIC instruments were in small window mode through the medium filter. About 4.5 ks of the 2001 April observation was contaminated by an extremely high background flare which makes the source difficult to distinguish above the noise. For this reason we do not use these data here. We will concentrate on the EPIC pn data for this analysis since the MOS cameras were in small window free running mode which has not yet been fully calibrated, and the pn data have the higher net count rate.

The raw data from all observations were processed with the XMM-Newton SAS (Science Analysis Software) version 5.3. We extracted a circular region around the source with radius 41.3 arcsec. Background subtraction was done using a nearby source-free region 4x larger in radius than the source extraction region. We summed the spectra from the 2000 October and the 2001 October pn observations to increase signal to noise, having first confirmed that no significant variability occurs either during or between the two observations. The latest EPIC pn response file for small window medium filter single events has been used. The resulting spectrum is grouped such that a minimum of 20 counts fall in each bin so $`\chi ^2`$ statistics apply. The mean net count rate is 5.71$`\pm `$0.02 counts s-1.

## X-RAY SPECTRAL ANALYSIS

Spectral analysis was done using the X-ray spectral fitting package XSPEC, using the $`\chi ^2`$ minimization technique. We have restricted the considered energy range to 0.3-10.0 keV as there may be uncertainties in the instrument calibration at energies outside this range. We find that a simple power-law cannot fit the spectrum over this range largely due to the presence of a strong soft X-ray excess (Figure 1).  Fig. 1. The summed EPIC pn spectrum is shown against a power-law model with galactic absorption.   Fig. 2. The unfolded EPIC pn spectrum with the best-fitting model of 2 blackbodies and a power-law plus galactic absorption.

The pn spectrum is well fitted with a model consisting of two blackbodies dominating at low energies with temperatures kT$`{}_{1}{}^{}=0.09\pm 0.01`$ keV and kT$`{}_{2}{}^{}=0.21\pm 0.03`$ keV, plus a power-law with $`\mathrm{\Gamma }`$ = 1.8$`\pm `$0.08 ($`\chi ^2`$/dof = 458/468; Figure 2). The power-law slope is consistent with the overall 2-10 keV spectrum of Seyferts in general, which have indices typically $`\mathrm{\Gamma }`$1.8-2.0. The X-ray flux in the 0.3-10.0 keV energy range is 9.83$`\times `$10-12$`\text{erg}\text{cm}^2\text{s}^1`$, corresponding to a luminosity of 4.74 $`\times `$1044 erg s-1 (assuming H0 = 50.0 and q0 = 0.0). Upper limits on the equivalent widths of an Iron line at 6.4 keV (in the rest frame of the source) are 60 eV for a narrow line ($`\sigma `$ = 0.01) and 50 eV for a broad line ($`\sigma `$ = 0.3). Inclusion of the narrow line gives an improvement according to the F-test of F statistic value = 2.14 and probability = 9.7$`\times `$10-2. The Galactic column density was fixed at 1.4$`\times `$1020 cm-2 (Puchnarewicz et al. 1995) in all fits. Adding an intrinsic absorption component gives no significant improvement to the fit.

### Physical Models

We tried a simple power-law-with-reflection model (Magdziarz & Zdziarski 1995), where incoming photons may be reflected off ionised material such as an accretion disk. This gives a good fit only if the X-ray flux is reflected with a reflection fraction of $``$15 and $`\chi ^2`$/dof = 721/659.

Comptonisation has often been suggested as a source of both the soft X-ray and hard X-ray spectrum so we investigated a number of such possibilties using the comptt model within XSPEC (Titarchuk 1994). The best fitting model ( $`\chi ^2`$/dof = 672/656; Figure 3) is a hot optically thin component (kT1 = 83 keV, $`\tau _1`$ = 0.8) irradiated by a 30 eV photon temperature, alongside a cooler optically thick component (kT2 = 3.5 keV, $`\tau _2`$ = 2.8) irradiated by a 60 eV temperature medium, both upscattering the incident photons. However, while the parameters may appear physically realistic they cannot be reasonably constrained.

We also note that a multi-temperature disk blackbody plus Comptonisation can also fit the spectrum well with $`\chi ^2`$/dof = 713/654, but again the parameters are unconstrained.

Fig. 3. The unfolded EPIC pn spectrum and double Comptonisation model.

Fig. 4. The XMM-Newton EPIC pn data in the 0.7-10 keV range, together with the power-law model which best fits the ASCA GIS2 and GIS3 data. The lower panel shows the ratio of the fit to the data, clearly showing excess emission in the pn continuum below 1 keV.

## COMPARISON WITH ASCA AND ROSAT

Comparing the data from all three missions (see also Puchnarewicz et al. 1995 and Breeveld et al. 2001) in the overlapping energy range of 0.7-2.0 keV by simply fitting a power-law (Table 1), we find that in the XMM-Newton observations RE J2248-511 is in a much softer state than at the time of the ASCA observation ($`\mathrm{\Gamma }`$ = 2.46 in XMM-Newton but 2.05 in ASCA ; for ASCA , the SIS response is uncertain below 1 keV, so only the GIS data were used). The power-law slope measured with the ROSAT PSPC is softer still ($`\mathrm{\Gamma }`$ = 2.7), although it should be borne in mind that there are reports that the PSPC may over-estimate $`\mathrm{\Gamma }`$ by as much as 0.4 in which case the XMM-Newton and ROSAT PSPC data would not be inconsistent. We also note that the error on the PSPC slope is large enough to be consistent with either the XMM-Newton or ASCA slopes over this range.

Thus, the soft X-ray component first observed with ROSAT which subsequently disappeared in the ASCA observation (Breeveld et al. 2001), now seems to have re-appeared in XMM-Newton (Figure 4). The 2-10 keV spectrum has changed relatively little (in flux and spectral form) between the ASCA and XMM-Newton observations. From the two-blackbody-plus-power-law fit to the XMM-Newton EPIC pn data (kT$`{}_{1}{}^{}=0.09\pm 0.01`$ keV, kT$`{}_{2}{}^{}=0.21\pm 0.03`$ keV, $`\mathrm{\Gamma }=1.8\pm 0.08`$), we calculate a flux in the soft component (ie in the two blackbodies; 0.3-2.0 keV) of 2.49$`\times `$10-12$`\text{erg}\text{cm}^2\text{s}^1`$.

## OPTICAL AND IR OBSERVATIONS

RE J2248-511 was observed on 2000 October 19-24 using the South African Astronomical Observatory’s 1.9m Radcliffe telescope. Three observations were made with the grating spectrograph and SITe CCD one night apart and measurements of the IR magnitude in the J,H,K bands were made on a single night with the infrared photometer. Spectra were taken using both a wide (700 $`\mu `$m) and a narrow (300 $`\mu `$m) slit with overlapping red and blue gratings ($`\lambda _{central}`$ = 4600,7800Å). Arc spectra were taken before and after every target and every standard spectrum. The spectra have been extracted using the Image Reduction and Analysis Facility (IRAF).

## RESULTS FROM THE OPTICAL SPECTROPHOTOMETRY

Seeing, humidity and cloud cover were extremely variable during the observing period and consequently only one of the nights we consider to be photometric. The spectra were not taken at the parallactic angle, but there appears to have been little or no loss of blue light since the spectra become stronger towards the blue end.

### The Broad Line Region

Both H$`\alpha `$ and H$`\beta `$ are best fitted with three gaussians with broad components of FWHM approximately 3700Å. The ratio of H$`\alpha `$ to H$`\beta `$ line flux, excluding the narrow-line component, is H$`\alpha `$/H$`\beta `$ = 2.67. The line parameters we measure are consistent with those measured in the spectrum during 1991 (Mason et al. 1995; H$`\alpha `$ was not observed). The H$`\beta `$ measurements are in approximate agreement with the 1992 measured values (Grupe et al. 1999).

Fig. 5. Spectral energy distribution from optical to X-ray for RE J2248-511 including all available optical, IR and UV data and best-fitting X-ray models spanning 1990 to 2001.

Fig. 6. The SED for October 2000 and two Comptonised disk models with $`\dot{M}=0.8M_{Edd}`$ and $`\mathrm{cos}(i)=0.8`$. The model represented by the dashed line is for black hole mass $`M=10^8\text{M}_{}`$; the dot-dashed line is for $`M=10^{7.9}\text{M}_{}`$.

### Optical Variability and the Big Blue Bump

In Figure 5 all available data for RE J2248-511 is plotted to show the various flux levels that have been observed. We find that the slope of the optical spectrum in the SED plot lies in between the slopes measured in 1991 (Mason et al. 1995) and 1992 (Grupe et al. 1998). Since the soft X-ray data from 2000/01 can be fit with a multi-temperature disk blackbody model, we investigated the possibility of the entire optical to soft X-ray flux originating in an accretion disk. We use Comptonised accretion disk models developed by Czerny & Elvis (1987) and fitted these to the spectral energy distribution for October 2000. Two model fits are shown with the data in Figure 6. These models can fit the optical slope reasonably well and are not inconsistent with the soft X-ray data. These data require models with masses of order $`M=10^8\text{M}_{}`$, high accretion rates of $`\dot{M}0.8M_{Edd}`$ and high disk inclination angles of $`\mathrm{cos}(i)0.8`$.

## DISCUSSION

Of the physical models for the X-ray spectrum tested here we are unable to constrain the key parameters in any case. The power-law-plus-reflection model is able to fit the data only with a high reflection fraction. This multiple reflection scenario could conceivably occur if there were a large number of clouds above or around the X-ray source which could reflect the emission but not transmit it. Alternatively, the accretion disk could be of varying height, like ripples on a pond, each ripple reflecting the incoming X-rays (Ballantyne, Ross & Fabian 2001). However, we do not detect strong Iron emission or the neutral Iron edge at 7.1 keV which are characteristic of strong reflection, except where the reflecting material is highly ionised. This Seyfert lies at the boundary between the quasar and Seyfert galaxy luminosities. Seyfert galaxy X-ray spectra typically show disk features including Iron lines and a reflection component which are less common in quasars. This is thought to be because the disk material becomes fully ionised and does not produce emission lines (eg. Matt et al. 1993). Though we do not see low ionisation Iron K$`\alpha `$ emission in these data, the upper limits on its equivalent width are typical of lines already observed in many Seyfert galaxies. Therefore reflection is unlikely to be making a strong contribution to the X-ray continuum flux in this source.

The double Comptonisation model requires both optically thick and optically thin material, each being illuminated by different temperature photons. A change in optical depth of the material could explain the change between the flat ASCA spectrum and the very soft EPIC pn spectrum. This scenario has also been suggested for the Seyfert 1 galaxy 1H 0419-577 (Page et al. 2001) and is discussed for this source in Breeveld et al. (2001). This model, however, could not account for any correlation with the optical flux.

Currently there is insufficient data to distinguish between many of the models and a 2-blackbodies-plus-power-law model remains the best representation of this data.
From the variability time-scale we estimate the size of the region emitting the flux of the soft X-ray excess to be $``$3.4 l.y. or $``$1 pc, which only locates it to at or within the BLR.

That both the 2000 optical and soft X-ray slopes are steep and both wavebands show great variability suggests a link between them. The optical to soft X-ray spectral energy distribution for 2000 is not inconsistent with Comptonised accretion disk models. Contrary to the proposed NLS1 scenario, the data are best fitted with intermediate to high black hole masses, whilst still favouring the high accretion rates suggested for NLS1 galaxies. Therefore, it is not necessary for a black hole to have a low mass for the formation of an ultrasoft X-ray excess. The relatively high velocity of the BLR may, however, be linked with the black hole mass, perhaps as a consequence of the increased gravitational potential. The mass and accretion rate we estimate from this modelling implies an inner disk effective temperature of $``$25 eV assuming a standard disk model, lower than that fitted to the lower-temperature soft X-ray component. However, the temperature measured in the X-ray blackbody models is the colour temperature which may be much higher than the effective temperature if the spectrum is modified by electron scattering (eg Titarchuk & Shrader 2002).

## CONCLUSIONS

This XMM-Newton observation has shown yet more dramatic activity in the soft X-ray spectrum of RE J2248-511, whilst the hard X-ray power law remains constant. We find that both Comptonisation and reflection models, while appearing to fit the X-ray spectrum well, cannot be constrained, and a model consisting of 2 blackbodies and a power-law remains the best representation of these data. The 2000 optical to soft X-ray spectral energy distribution may be consistent with a face-on Comptonised accretion disk accreting at a rate approaching the Eddington limit onto a 10$`{}_{}{}^{8}\text{M}_{}^{}`$ black hole. Observations at UV wavelengths are needed to tie down the parameters of any such model, and more simultaneous multi-wavelength data are needed to confirm the correlation between optical and soft X-ray variability and the presence of a single Big Blue Bump.

## ACKNOWLEDGEMENTS

This work is based on observations obtained with XMM-Newton , an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). IRAF is distributed by the National Optical Astronomy Observatories, which are operated by AURA, Inc., under cooperative agreement with the National Science Foundation. We would like to thank Aneta Siemiginowska for use of the Comptonised accretion disk models. RLCS acknowledges support from a PPARC studentship.

email address of R. L. C. Starling rlcs@mssl.ucl.ac.uk
Manuscript received 10/12/02; revised 11/02/03; accepted 13/02/03


text:
# On the Generalized Borel Transform and its Application to the Statistical Mechanics of Macromolecules.

## I Introduction

Consider a group of $`n`$ small molecules connected in a sequential manner such that each small molecule has only two nearest neighbours with whom it forms a chemical bond except for the small molecules at the end of the chain that have only one connected nearest neighbour. This chain is called a polymer. In this particular case, the chain ends do not form a chemical bond thus, the polymer is said to have a linear topology (linear polymer). If the ends were to form a chemical bond, then the polymer is said to be a ring (cyclic) polymer. Another well known topology is the star topology where $`m`$ linear polymers grow from the same point and they are connected at this point. This defines an $`m`$-arm star polymer. Combinations of these three topologies or new topologies define more complex macromolecules.

Another characteristic of polymers is that they can be made of two or more different kinds of small molecules.

We conclude this section by speculating about other possible models for which the GBT might provide some useful insight. But, before we proceed, we should note that extensions of the present calculation to other polymer topologies like rings and stars are straightforward. For example, in the case of rings the cyclic topology of the polymer generates an extra constraint in the polymer propagator of the form

$$\delta \left(\underset{j=1}{\overset{n}{}}𝐑_j\right),$$
(1)
which can be easily handled by the GBT.

A generalization of Eq. (23) that takes into account more complex systems like copolymers and the effects of external fields has the mathematical expression given by

$$P(𝐑,n,\left\{p_j^\alpha \right\})=d\left\{𝐑_k\right\}\underset{j=1}{\overset{n}{}}\left(\underset{\alpha =1}{\overset{s}{}}\left(p_j^\alpha \tau _j^\alpha \left(𝐑_j\right)\right)\right)\delta \left(\underset{j=1}{\overset{n}{}}𝐑_j𝐑\right)\mathrm{exp}\left(\omega \left(\left\{𝐑_j\right\}\right)\right),$$
(2)
where $`s`$ is the total number of chemical species forming the copolymer, $`p_j^\alpha `$ is the probability of finding the $`j`$-th segment in the $`\alpha `$-th chemical species and $`\omega \left(\left\{𝐑_j\right\}\right)`$ has the mathematical form

$$\omega \left(\left\{𝐑_j\right\}\right)=\underset{j=1}{\overset{n}{}}\eta \left(𝐑_j\right),$$
(3)
$`\eta \left(𝐑\right)`$ being any function. In particular, the effects of external vectorial $`\left(\eta \left(𝐑\right)=𝐅𝐑\right)`$ and quadrupolar $`\left(\eta \left(𝐑\right)=Q_{ij}R_iR_j\right)`$ fields can be studied exactly. In general, the statistical properties of the model defined by Eq. (2) can be computed exactly using GBT.

In recent years, new experimental tools have been developed to study the behavior of single macromolecules. These tools have given researchers the ability to manipulate polymers individually. Perhaps, the force-elongation measurements done in different biopolymers using Atomic Force Microscopes and Optical Tweezers are the best known examplesMolecules . This new experimental capability has proved that our most advanced models are still incomplete and cannot describe modern experimental data adequately. One example of this statement is the force-elongation curve of the polyprotein made of the I27 module of human cardiac titinMolecules\_2 . This curve shows discontinuities that occur when one module of the polyprotein unravels. The data is generally fitted to the WCM. But, this model cannot describe structures like $`\alpha `$-helices or $`\beta `$-sheets present in polyproteins.

The solution of the WCM is mathematically involved and the addition of secondary structure like $`\alpha `$-helices does not simplify the model. On the contrary, the mathematical complexity of this “extended” WCM is expected to be even more involved. As a consequence of this, we have started to explore new mathematical tools that might lead to alternative approaches to the models of polymers. One tool that has proved to be potentially useful is the Generalized Borel Transform (GBT) which we present in this paper. This technique was borrowed from Quantum Mechanics (QM) and Quantum Field Theories (QFT) where it has been successfully applied to solve different problems such as heavy quark-antiquark interaction in QFT tecnica2  and singular scattering potentials in QM tecnica3 . In addition, this mathematical tool allows one to obtain analytical solutions of parametric integrals like Laplace/Mellin transforms for all range of its parameters. Therefore, it is extremely useful to study non-perturbative regimes. In this paper, we show how to use this computational tool in the field of statistical mechanics of single polymers by computing exactly the polymer propagator of the RFM of linear macromolecules. In the companion paper, we generalize these ideas to the statistical mechanics of linear semiflexible polymer chains, as described by the WCM, and present new results for this model.

The RFM is the simplest attempt to describe the statistical properties of flexible polymer chains. Thus, it is the first model where the GBT should be applied. In this model, polymers are assumed to be made of n segments of length l (=Kuhn length) which are connected in a sequential manner. These segments are free to rotate around any axis as long as the move does not break the polymer chain. Figure 1 depicts the RFM.

This paper is organized as follows. In section II we show how is possible to calculate Fourier sin/cos transforms using the Generalized Borel Transform and we present a brief summary of the mathematical aspects of this technique. Then in section III we solve exactly a particular Fourier sin transform, namely the polymer propagator (normalized distribution function) of the RFM. Finally, in section IV we present the conclusions and discuss some other models in polymer physics where this new application of the GBT might lead to exact solutions.

## II The Generalized Borel Transform

In the previous section we have shown that the Fourier Sine transform plays a very important role in the evaluation of the statistical properties of models for single macromolecules. Therefore, let us start the description of the GBT by writing the general expression of a Fourier Sine transform of a function $`H(k,a)`$

$$Q(R,a)=_0^{\mathrm{}}\mathrm{sin}\left(Rk\right)H(k,a)𝑑k.$$
(4)
Furthermore, consider the Laplace transform of the same function which is given by the definition

$$S(g,a)=_0^{\mathrm{}}H(x,a)\mathrm{exp}\left(gx\right)𝑑xg>0.$$
(5)
Then, we observe that we can obtain the Fourier Sine transform, $`Q(R,a)`$, from the Laplace transfom, $`S(g,a)`$, as the analytic continuation of $`S(g,a)`$ to the complex plane as follows

$$Q(R,a)=Im\left\{S\left(g=iR,a\right)\right\}.$$
(6)
Consequently, we will focus on the evaluation of Laplace transforms. For this purpose, we will employ the GBT technique described hereafter.

The main goal of the GBT is to obtain analytical solutions of parametric integrals of the Mellin/Laplace typetecnica1 ; tecnica2 ; tecnica3  for all the range of values of the parameters. Therefore, this technique is very useful to study non-perturbative regimes. The basic idea of method consists of introducing two auxiliary functions, $`S(g,a,n)`$ and $`B_\lambda (s,a,n)`$(the Generalized Borel Transform). These functions depend on auxiliary parameters called $`n`$ and $`\lambda `$. These parameters have no physical meaning and are introduced for the sole purpose of helping in the computation of an explicit mathematical expression for $`S(g,a)`$ for all values of the true parameters $`g`$ and $`a`$.

Let us start with the mathematical definition of $`S(g,a,n)`$ which is the following

$$S(g,a,n)=_0^{\mathrm{}}x^nH(x,a)\mathrm{exp}(gx)dx.g>0$$
(7)
We have explicitly extracted a factor $`x^n`$ from the function to be transformed. This integral is related to the Laplace transform,Eq. (5), by the following relationship

$$S(g,a,n)=()^n\frac{^n}{g^n}S(g,a),$$
(8)
which can be inverted to give

$$S(g,a)=()^n\underset{n}{\underset{}{𝑑g\mathrm{}𝑑g}}S(g,a,n)+\underset{p=0}{\overset{n1}{}}c_p(a,n)g^p.$$
(9)
The finite sum comes from the indefinite integrations. Note that all the coefficients vanish whenever the Laplace transform, Eq. (5), fulfills the following asymptotic condition

$$\underset{g\mathrm{}}{lim}S(g,a)=0.$$
(10)
In addition, the expression given by Eq. (9), is valid for any value of the parameter $`n,`$ in particular for $`n1`$ where the GBT provides an approximate analytical expression for $`S(g,a,n)`$ as we describe below.

Let us define the Generalized Borel Transform of $`S(g,a,n)`$ as follows

$$B_\lambda (s,a,n)_0^{\mathrm{}}\mathrm{exp}\left[s/\eta \left(g\right)\right]\left[\frac{1}{\lambda \eta \left(g\right)}+1\right]^{\lambda s}\frac{S(g,a,n)}{\left[\eta \left(g\right)\right]^2}\frac{\eta \left(g\right)}{g}𝑑g,Re\left(s\right)<0$$
(11)
where $`\lambda `$ is any real, positive and non-zero number, and $`\eta `$ is defined as follows $`1/\eta \lambda \left(\mathrm{exp}\left(g/\lambda \right)1\right)`$. Then, $`B_\lambda (s,a,n)`$ is an analytic function for real values of $`s`$ less than zero. Moreover, the analytic continuation to the other half of the complex plane generates an analytic function with a cut on the positive real axis.

In order to invert the transform defined by Eq.(11), we note that the change of variables $`u\left(g\right)=1/\eta \lambda \mathrm{ln}\left[1+1/\lambda \eta \right]`$ transforms the integral, Eq.(11), into a Laplace transform

$$B_\lambda (s,a,n)_0^{\mathrm{}}\mathrm{exp}\left[su\right]L_\lambda (S,a,n,u)𝑑u,Re\left(s\right)<0$$
(12)
where $`L_\lambda (S,a,n,u)`$ depends on $`S(g,a,n)`$. Consequently, the inverse Laplace transform of Eq.(12) provides a procedure for the evaluation of $`S(g,a,n)`$ by integrating $`B_\lambda (s,a,n)`$ on the imaginary axis or over the discontinuity of $`B_\lambda (s,a,n)`$ on the cut. After a change of variables we can write $`S(g,a,n)`$ as follows

$$S(g,a,n)=2\lambda ^2\left(1\mathrm{exp}\left(g/\lambda \right)\right)_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\mathrm{exp}\left[G(w,t,g,\lambda ,a,n)\right]𝑑w𝑑t.$$
(13)
$`G(w,t,g,\lambda ,a,n)`$ is given by the following expression (for more details see Ref. tecnica2 )

$$\begin{array}{c}G(w,t,g,\lambda ,a,n)=s\left(t\right)u\left(g\right)+t\mathrm{ln}\left\{\mathrm{\Gamma }\left[\lambda \left(s\left(t\right)+x\left(w\right)\right)\right]\right\}\\ +\left\{\lambda \left[s\left(t\right)+x\left(w\right)\right]1\right\}\mathrm{ln}\left(\lambda s\left(t\right)\right)\lambda s\left(t\right)+w+\mathrm{ln}\left[x\left(w\right)^nH\left(x\left(w\right)\right)\right],\end{array}$$
(14)
where $`s\left(t\right)=\lambda \mathrm{exp}\left(t\right)`$ and $`x\left(w\right)=\mathrm{exp}\left(w\right)`$.

Note that the expression given by Eq. (13) is valid for any non-zero, real and positive value of the parameter $`\lambda .`$ However, the resulting expression for $`S(g,a,n)`$ does not depend on $`\lambda `$ explicitly. Thus, each value of the parameter $`\lambda `$ defines a particular Borel transform. Consequently, we can choose the value of this parameter in such a way that it allows us to solve Eq. (13).

The dominant contribution to the double integral is obtained using steepest descentnorman ; copson  in the variables $`t`$ and $`w`$. For this purpose, we first compute the expressions of the saddle point $`t_o(g,a,n)`$ and $`w_o(g,a,n)`$ in the limit $`\lambda 1.`$The results are the following

$$t_o(g,a,n)=\mathrm{ln}\left[\frac{x_o^2(g,a,n)}{f(x_o(g,a,n),a,n)}\right],w_o(g,a,n)=\mathrm{ln}\left[x_o(g,a,n)\right],$$
(15)
where $`x_o(g,a,n)`$ is the real and positive solution of the implicit equation coming from the extremes of the function $`G(w,t,g,\lambda ,a,n)`$ in the asymptotic limit in $`\lambda `$. Explicitly, the equation is

$$x_o^2g^2=f(x_o,a,n)\left[f(x_o,a,n)+1\right],$$
(16)
where

$$f(x_o,a,n)1+n+x_o\frac{d\mathrm{ln}\left[H(x_o,a)\right]}{dx_o}.$$
(17)
Afterward, we check the positivity conditionjefrey  (the Hessian of $`G(w,t,g,\lambda ,a,n)`$ at the saddle point must be positive). Let us call the Hessian $`D(x_o,a,n)`$. Its mathematical expression is

$$D(x_o,a,n)x_o\frac{df(x_o,a,n)}{dx_o}\left[1/2+f(x_o,a,n)\right]+f(x_o,a,n)\left[1+f(x_o,a,n)\right].$$
(18)
Observe that in the range of the parameters where $`f(x_o,a,n)1`$, which is fulfilled when $`n1`$, we can keep the second order term in the expansion of $`G(w,t,g,\lambda ,a,n)`$ around the saddle point. Then, we can approximate the double integral in Eq. (13) as follows

$$S_{Aprox}(g,a,n)=4\pi \frac{\lambda ^2\left(1\mathrm{exp}\left(g/\lambda \right)\right)}{\sqrt{D[x_o,a,n]}}\mathrm{exp}\left[G(w_o,t_o,g,\lambda ,a,n)\right].$$
(19)
In the limit $`\lambda \mathrm{}`$ we obtain the following approximate expression for $`S(g,a,n)`$

$$S_{Aprox}(g,a,n)=\sqrt{2\pi }e^{1/2}\frac{\sqrt{f[x_o,a,n]+1}}{\sqrt{D[x_o,a,n]}}\left[x_o\right]^{n+1}H[x_o,a]\mathrm{exp}\left[f[x_o,a,n]\right].$$
(20)
Note that the expression given by Eq. (20) is valid for functions $`H(x,a)`$ that fulfill the following general conditions. Firstly, the relationship given by Eq. (16) must be biunivocal. Secondly, $`D(x_o,a,n)`$ must be positive at $`x_o`$. Thirdly, $`f(x_o,a,n)`$ must be larger than one. In particular, this condition is fulfilled when $`n1.`$ These conditions provide the range of values of the parameters where the approximate solution, Eq. (20), is valid.

Finally, we replace the expression given by Eq. (20) into Eq. (9) to obtain an approximate analytical expression for the Laplace transform $`S(g,a)`$. In particular, in the limit $`n\mathrm{}`$, we obtain the following analytical solution for $`S(g,a)`$

$$S(g,a)=\begin{array}{c}lim\\ n\mathrm{}\end{array}()^n\underset{n}{\underset{}{𝑑g\mathrm{}𝑑g}}S_{Aprox}(g,a,n).$$
(21)
One particular case of this result is the one where $`H(x,a)`$ does not contribute to the saddle point. This is the case when $`f(x_o,a,n)`$ can be approximated by $`1+n`$ (the derivative of $`H(x_o,a)`$ is negligible). Then, the saddle point solution is $`x_o(g,a,n)\left(n+3/2\right)/g`$ and the expression of $`S_{Aprox}(g,a,n)`$ is

$$S_{Aprox}(g,a,n)\frac{\mathrm{\Gamma }\left(n+1\right)}{g^{n+1}}H[x_o(g,a,n),a]n1.$$
(22)
Another important property of the expression given by Eq. (21) is that, in the limit $`n\mathrm{}`$, the approximate solution, Eq. (20), becomes an exact solution for Eq. (7). Thus, as long as the $`n`$ indefinite integrals are calculated without approximations, then Eq. (21) is an exact solution for Eq. (5).

In summary, the procedure to use the GBT to compute Fourier Sine transforms is the following. Firstly, one has to solve the implicit equation, Eq. (16), for $`n1`$ to obtain the mathematical expression of $`x_o(g,a,n)`$. Replacing this expression into Eq. (20) and doing the $`n`$ indefinite integrals in Eq. (21), we get the exact/approximate expression for $`S(g,a)`$. Finally, one has to compute the analytic continuation of $`S(g,a)`$, Eq. (6), to get the exact/approximate solution of the Fourier Sine transform, Eq. (4).

In the next section we apply this technique to solve exactly the statistical mechanics of flexible macromolecules.

## III Application

Let us start by analyzing the polymer propagator predicted by the Random-Flight Model which is given by

$$P(𝐑,n)=d\left\{𝐑_k\right\}\underset{j=1}{\overset{n}{}}\tau \left(𝐑_j\right)\delta \left(\underset{j=1}{\overset{n}{}}𝐑_j𝐑\right),$$
(23)
where $`𝐑_j`$ is the bond vector between the $`(j1)`$-th and $`j`$-th beads, $`n`$ is the total number of segments per polymer chain, $`𝐑`$ is the end-to-end vector and $`\tau \left(𝐑_j\right)`$ is given by the formula

$$\tau \left(𝐑_j\right)=\frac{\delta \left(\left|𝐑_j\right|l\right)}{4\pi l^2}.$$
(24)
Using the Fourier representation of the delta function, we obtain

$$\begin{array}{c}P(𝐑,n)=\frac{d^3k\mathrm{exp}\left(i𝐑𝐤\right)}{\left(2\pi \right)^3\left(4\pi l^2\right)^n}\left[d\left\{𝐑_k\right\}\underset{j=1}{\overset{n}{}}\delta \left(\left|𝐑_j\right|l\right)\mathrm{exp}\left(i\underset{j=1}{\overset{n}{}}𝐑_j𝐤\right)\right]\\ =\frac{d^3k}{\left(2\pi \right)^3}\mathrm{exp}\left(i𝐑𝐤\right)K(𝐤,n,l),\end{array}$$
(25)
where the characteristic function, $`K(𝐤,n,l)`$, is

$$K(𝐤,n,l)=\left(\frac{\mathrm{sin}\left(\left|𝐤\right|l\right)}{\left|𝐤\right|l}\right)^n.$$
(26)
The evaluation of the angular integrals in Eq. (25) is straightforward. After rescaling $`𝐑`$ and $`𝐤`$ with the Kuhn length, $`l`$, we obtain the final expression for the polymer propagator

$$P(R,n)=\frac{2}{\left(2\pi \right)^2R}_0^{\mathrm{}}𝑑k\left[\mathrm{sin}\left(kR\right)\left(\frac{\mathrm{sin}\left(k\right)}{k}\right)^nk\right],$$
(27)
where $`k=\left|𝐤\right|`$ and $`R=\left|𝐑\right|`$.

This integral representation of the polymer propagator is a Fourier Sine transform and can be solve exactly using GBT. Then, our first step consists of expressing the polymer propagator, Eq. (27), in terms of a Laplace Transform. For this purpose we define the function

$$S(b,n)_0^{\mathrm{}}𝑑w\left[\mathrm{exp}\left(wb\right)\left(\frac{\mathrm{sin}\left(w\right)}{w}\right)^nw\right].$$
(28)
from where we recover the expression of the polymer propagator, Eq. (27), as the analytic continuation of the function $`S(b,n)`$ to the complex plane

$$P(R,n)=\frac{2}{\left(2\pi \right)^2R}Im\left\{S\left(b=iR,n\right)\right\}.$$
(29)
Let us now rewrite Eq. (28) as follows

$$S(b,n)=\frac{^n}{c^n}\left\{_0^{\mathrm{}}𝑑w\left[w\mathrm{exp}\left(wb\right)\mathrm{exp}\left(c\frac{\mathrm{sin}\left(w\right)}{w}\right)\right]\right\}_{c=0}=\frac{^n}{c^n}\left\{GA(b,c)\right\}_{c=0},$$
(30)
where

$$GA(b,c)_0^{\mathrm{}}𝑑ww\mathrm{exp}\left(wb\right)H(w,c),$$
(31)
and

$$H(w,c)\mathrm{exp}\left[c\frac{\mathrm{sin}\left(w\right)}{w}\right].$$
(32)
Then, the integral expressed by Eq. (31) satisfies all the requirements of the GBT technique tecnica3 . Consequently, we evaluate this integral in the following way

$$GA(b,c)=\begin{array}{c}lim\\ N\mathrm{}\end{array}\begin{array}{c}()^N\underset{}{𝑑b\mathrm{}\mathrm{}\mathrm{}𝑑b}GA_N,\\ N\end{array}$$
(33)
where we have defined

$$GA_N(b,c)_0^{\mathrm{}}𝑑w\left[w^{N+1}\mathrm{exp}\left(wb\right)H(w,c)\right].$$
(34)
In the asymptotic limit of $`N\mathrm{}`$, the GBT provides an analytical solution for Eq. (34). Following the technique, we first solve the following implicit equation for $`w_o`$, Eq. (16),

$$\left\{N+1+w_o\frac{}{w_o}\mathrm{ln}\left[H(w_o,c)\right]\right\}\left\{N+2+w_o\frac{}{w_o}\mathrm{ln}\left[H(w_o,c)\right]\right\}=w_o^2b^2,$$
(35)
whose asymptotic solution is

$$w_o\frac{N+5/2}{b}N1.$$
(36)
Replacing this expression for $`w_o`$ in the expression provided by the GBT, Eq. (20), we obtain

$$GA_N(b,c)\frac{\mathrm{\Gamma }\left(N+2\right)}{b^{N+2}}H(\frac{N+5/2}{b},c)N1.$$
(37)
Furthermore, we replace Eq. (37) into Eq. (33) and the resulting expression into Eq. (30), then we obtain

$`S(b,n)`$ $`=\begin{array}{c}lim\\ c0\end{array}{\displaystyle \frac{^n}{c^n}}\{\begin{array}{c}lim\\ N\mathrm{}\end{array}\begin{array}{c}()^N\underset{}{𝑑b\mathrm{}\mathrm{}\mathrm{}𝑑b}\frac{\mathrm{\Gamma }\left(N+2\right)}{b^{N+2}}H(\frac{N+5/2}{b},c)\hfill \\ N\hfill \end{array}\}.`$ (44)
We now proceed to exchange the order of the operators, first we evaluate the $`n`$-th derivative of the function $`H`$ with respect to $`c`$ and, afterward, we take the limit of $`c0`$ . As a result, we obtain

$$S(b,n)=\begin{array}{c}lim\\ N\mathrm{}\end{array}()^N𝑑b\mathrm{}𝑑b\frac{\mathrm{\Gamma }\left(N+2\right)}{b^{N+2}}\left(\frac{\mathrm{sin}\left(\frac{N}{b}\right)}{\left(\frac{N}{b}\right)}\right)^n.$$
(45)
Next, we solve the $`N`$ integrations. Using properties of the function $`\mathrm{sin}\left(x\right)`$ we can write $`S(b,n)`$ for any odd number of segments as follows

$$S(b,n)=\frac{1}{2^{n1}}\underset{k=0}{\overset{\frac{n1}{2}}{}}()^{\frac{n1}{2}+k}\left(\begin{array}{c}n\\ k\end{array}\right)M(N,n,k,b),$$
(46)
where

$$M(N,n,k,b)\begin{array}{c}lim\\ N\mathrm{}\end{array}()^NIm\underset{r=0}{\overset{\mathrm{}}{}}\frac{\left(i\left(n2k\right)\right)^rN^{rn}}{r!}𝑑b\mathrm{}𝑑b\frac{\mathrm{\Gamma }\left(N+2\right)}{b^{N+2n+r}}.$$
(47)
We note that the only powers of $`b`$ in Eq. (47) that fulfill the asymptotic behavior of the function $`S(b,n)`$, Eq. (10), are those that satisfy the condition $`r\left(n1\right)`$. Consequently, the $`N`$ indefinite integrations are exactly doable, the result is

$$𝑑b\mathrm{}𝑑b\frac{1}{b^{N+2n+r}}=\frac{\mathrm{\Gamma }\left(2+rn\right)}{\mathrm{\Gamma }\left(N+2n+r\right)}\frac{()^N}{b^{2+rn}}.$$
(48)
Replacing Eq. (48) into Eq. (47) and introducing the dummy variable $`r=x+n1`$ we can write

$$\begin{array}{c}M(N,n,k,b)Im\frac{1}{b}\left(i\left(n2k\right)\right)^{n1}\underset{x=0}{\overset{\mathrm{}}{}}\left(\frac{i\left(n2k\right)}{b}\right)^x\frac{\mathrm{\Gamma }\left(x+1\right)}{\mathrm{\Gamma }\left(x+n\right)}\\ \times \begin{array}{c}lim\\ N\mathrm{}\end{array}\frac{N^{x1}\mathrm{\Gamma }\left(N+2\right)}{\mathrm{\Gamma }\left(N+x+1\right)},\end{array}$$
(49)
which, after using the asymptotic properties of the Gamma function abramovich

$$\begin{array}{c}lim\\ N\mathrm{}\end{array}\frac{N^{x1}\mathrm{\Gamma }\left(N+2\right)}{\mathrm{\Gamma }\left(N+x+1\right)}=1,$$
(50)
becomes

$$M(n,k,b)=\frac{1}{b}Im\underset{x=0}{\overset{\mathrm{}}{}}\left(i\left(n2k\right)\right)^{n1}\left(\frac{i\left(n2k\right)}{b}\right)^x\frac{\mathrm{\Gamma }\left(x+1\right)}{\mathrm{\Gamma }\left(x+n\right)}.$$
(51)
The sum over $`x`$ is doable, the result gives the following expression for $`M(n,k,b)`$

$$M(n,k,b)=\frac{1}{b}Im\left[\left(i\left(n2k\right)\right)^{n1}FD(n,k,b)\right],$$
(52)
where we have defined

$$\begin{array}{c}FD(n,k,b)\frac{\mathrm{\Gamma }\left(\frac{1}{2}\right)_3F_2([1,1,\frac{1}{2}],[\frac{n+1}{2},\frac{n}{2}],\frac{\left(n2k\right)^2}{b^2})}{\sqrt{\pi }\mathrm{\Gamma }\left(n\right)}\\ +\frac{i\left(n2k\right)}{b}\frac{{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{b^2})}{\mathrm{\Gamma }\left(n+1\right)}.\end{array}$$
(53)
$`{}_{3}{}^{}F_{2}^{}([,,],[,],x)`$ is the Generalized Hypergeometric function hyper . From Eq. (52) we can see that the imaginary part affects only the function $`FD(n,k,b)`$. Thus, we obtain the final expression for $`S(b,n)`$

$$S(b,n)=\underset{k=0}{\overset{\frac{n1}{2}}{}}()^k\left(\begin{array}{c}n\\ k\end{array}\right)\frac{\left(n2k\right)^n{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{b^2})}{b^22^{n1}\mathrm{\Gamma }\left(n+1\right)}.$$
(54)
The last step to obtain the analytical expression of the polymer propagator consists of inserting Eq. (54) into Eq. (29) and computing the analytic continuation of the resulting expression to the complex plane through the substitution $`b=iR`$. After doing these computations, we arrived at the following expression for the polymer propagator

$$\begin{array}{c}P(R,n)=\frac{1}{2^n\pi ^2R^3}\underset{k=0}{\overset{\frac{n1}{2}}{}}()^{k+1}\left(\begin{array}{c}n\\ k\end{array}\right)\left(n2k\right)^n\frac{1}{\mathrm{\Gamma }\left(n+1\right)}\\ \times Im\left\{{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{R^2})\right\}.\end{array}$$
(55)
This expression can be simplified even further if we use the well known analytical properties of the Hypergeometric function $`{}_{3}{}^{}F_{2}^{}\left(z\right)`$ hyper  which is an analytic function for values of $`\left|z\right|<1`$ and its continuation to the rest of complex plane generates one cut on the positive real axis starting at $`Re\left(z\right)=1`$. This implies that only values of $`\frac{\left(n2k\right)^2}{R^2}1`$ will contribute to the imaginary part of $`{}_{3}{}^{}F_{2}^{}\left(z\right)`$. Consequently, this condition reduces the number of terms in the $`k`$-sum such that the last term of Eq. (55) is $`k=\left[\frac{nR}{2}\right]`$.

The explicit evaluation of $`Im\left\{{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{R^2})\right\}`$ can be found in Appendix A. The final expression is

$$\begin{array}{c}Im\left\{{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{R^2})\right\}\\ =\frac{R^2\pi }{2\left(n2k\right)^n}\frac{\mathrm{\Gamma }\left(n+1\right)}{\mathrm{\Gamma }\left(n1\right)}\left[n2kR\right]^{n2}.\end{array}n2$$
(56)
Finally, we replace Eq. (56) into Eq.(55) to obtain the exact expression for the polymer propagator

$$P(R,n)=\frac{1}{2^{n+1}\pi R}\underset{k=0}{\overset{\left[\frac{nR}{2}\right]}{}}()^k\left(\begin{array}{c}n\\ k\end{array}\right)\frac{\left[n2kR\right]^{n2}}{\mathrm{\Gamma }\left(n1\right)}.$$
(57)
Equation (57) is valid for odd number of segments but, it is extended to polymers with any number of segments larger than two via analytic continuation.

Therefore, we have obtained the well-known Wang ; Hs  exact analytical expression for the polymer propagator of flexible chains, Eq. (27), with any number of segments, $`n`$, and any end-to-end distance, $`R`$.

Observe that the result given by Eq. (57) can be used to describe the statistical properties of polymers with other topologies. For example, consider the case of a flexible $`m`$-arm star polymer as shown in Fig. . Since the polymer is flexible, then each arm behaves independently from the other ones except for the fact that all of them start at the origin. Thus, the probability of finding the end of the $`j`$-th arm in the shell of radius $`R_j`$ with thickness $`dR_j`$ centered at the origin is

$$4\pi R_j^2P(R_j,n_j)dR_j,$$
(58)
where $`n_j`$ is the number of segments in the $`j`$-th arm. If we consider all the arms, then the probability of finding the end of the first, second, etc. arms in the shells of radii $`R_1`$,$`R_2`$, etc. with thicknesses $`dR_1`$, $`dR_2`$, etc. centered at the origin is

$$\left(4\pi \right)^m\underset{j=1}{\overset{m}{}}R_j^2P(R_j,n_j)dR_j.$$
(59)
Other probability distributions for star polymers can also be computed easily.

Another example is the case of ring (cyclic) polymers. Figure shows this topology. From this figure and following the steps presented in this paper for linear polymers it can be proved that the probability of finding any pair of segments separated by a distance $`R`$ should be proportional to the product of two propagators of the form given by Eq. (57)

$$P_{Ring}(R,s,Ls)P(R,s)P(R,ns),$$
(60)
where $`n`$ is the total number of segments in the ring and $`s`$ is the number of segments (along the contour of the polymer chain) between the two chosen segments.

The aformentioned two examples clearly show that the results obtained for linear polymers using the GBT can be used for polymers with other topologies thus, increasing the range of models that are mathematically tractable with the GBT.

## IV conclusions

In this paper we have described a new mathematical method called the Generalized Borel Transform and applied it to compute some statistical properties (polymer propagator) of models of flexible polymers. Specifically, we showed how the GBT is constructed and how to use it to compute Mellin/Laplace transforms. Moreover, some mathematical properties were presented. The application of this technique to the statistical mechanics of single flexible polymers led to the exact solution for the polymer propagator of linear polymers. The propagator obtained turned out to be a finite sum of polynomials valid for any end-to-end distance, $`R`$, and number of segments, $`n`$. Furthermore, this result was used to construct distribution functions for two other topologies, rings and stars.

The exact computation of the polymer propagator of the RFM is a straightforward calculation that requires simple mathematics when the GBT is used. This mathematical simplicity of the GBT makes it a potentially very useful computational tool for more complex models of single polymer chains because it does not add any complexity to the physics of the starting model.

Equation (57) together with its extensions to stars and rings, Eqs. (59) and (60), and the discussion presented in the introduction show that the GBT can solve exactly a wide range of models for polymers. However, more advanced models of single polymer chains like the Wormlike Chain Model or helical polymers where the bond vectors are correlated with each other through potential interactions are not exactly soluble with the GBT at present. A generalization of the GBT to multidimensional integrals is required to address these models.

## ACKNOWLEDGMENTS

We acknowledge the National Science Foundation, Grant # CHE-0132278 (CAREER), the Ohio Board of Regents Action Fund, Proposal # R566 and The University of Akron for financial support.

## Appendix A evaluation of the Imaginary part of the Hypergeometric function

In this appendix, we calculate the expression $`Im\left\{{}_{3}{}^{}F_{2}^{}\left(z\right)\right\}`$. For this purpose, we use the following integral representation of the Hypergeometric function tablarusa

$$\begin{array}{c}{}_{3}{}^{}F_{2}^{}([\nu ,\frac{\lambda }{2},\frac{\lambda +1}{2}],[\frac{\lambda +\mu }{2},\frac{\lambda +\mu +1}{2}],\frac{q^2}{w^2})=\left[w^{2\gamma }q^{\lambda +\mu 1}B(\lambda ,\mu )\right]^1\\ \times _0^qx^{\lambda 1}[qx]^{\mu 1}[x^2+w^2]^\nu dx,\lambda ,\mu >0,Re\left(\frac{q}{w}\right)>0\end{array}$$
(61)
where $`B(\lambda ,\mu )`$ is the Beta function tablarusa .

We now assign the values $`\nu =1,\lambda =2,\mu =n1,q=n2k,w=b`$ to the parameters in Eq. (61) to obtain

$$\begin{array}{c}{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{b^2})=\frac{1}{b^2\left(n2k\right)^nB(2,n1)}\\ \times _0^{n2k}x[n2kx]^{n2}[x^2+b^2]^1dx.\end{array}$$
(62)
This integral representation is valid only for $`n2`$. Therefore, when we take the analytic continuation to the complex plane as before $`(b=iR)`$, we can express the imaginary part of the Hypergeometric function as follows

$$\begin{array}{c}Im\left\{{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{R^2})\right\}=\frac{R^2}{\left(n2k\right)^nB(2,n1)}\\ \times Im_0^{n2k}x\left[n2kx\right]^{n2}\left[x^2R^2\right]^1𝑑x.\end{array}$$
(63)
Thus, we need to evaluate

$$LIm_0^{n2k}x\left[n2kx\right]^{n2}\left[xR\right]^1\left[x+R\right]^1𝑑x.$$
(64)
After analyzing the analytical behavior of the integrand, we concluded that we can exchange the operations of integration and imaginary part to obtain

$$L=_0^{n2k}x\left[n2kx\right]^{n2}\left[x+R\right]^1Im\left\{\left[xR\right]^1\right\}𝑑x.$$
(65)
Thus, we have to compute

$$LS=Im\left\{\frac{1}{\left(xR\right)}\right\},$$
(66)
first and, afterward, we have to solve the integral given by Eq. (65).

The analytical behavior of the function $`\left(xR\right)^1`$ is well known. It is an analytic function for $`\left|x\right|>R`$ but, its analytic continuation to the complex plane generates a cut on the real axis in the range $`R<Re\left(x\right)<R`$. This cut generates its imaginary part which is

$$Im\left\{\frac{1}{\left(xR\right)}\right\}=\pi \delta \left(xR\right).$$
(67)
Thus, replacing the expression given by Eq. (67) into Eq. (65) and making the change of variables $`y=xR`$, we obtain

$$L=\pi _R^{n2kR}F_k(y,n,R)\delta \left(y\right)𝑑y,$$
(68)
where

$$F_k(y,n,R)\left(y+R\right)\left[n2kRy\right]^{n2}\left[y+2R\right]^1.$$
(69)
The result of the integration gives

$$L=\frac{\pi }{2}\left[n2kR\right]^{n2}.$$
(70)
Finally, we replace Eq. (70) into Eq. (63) to obtain the final expression

$$\begin{array}{c}Im\left\{{}_{3}{}^{}F_{2}^{}([1,1,\frac{3}{2}],[\frac{n+1}{2},\frac{n+2}{2}],\frac{\left(n2k\right)^2}{R^2})\right\}\\ =\frac{\pi }{2}\frac{R^2}{\left(n2k\right)^nB(2,n1)}\left[n2kR\right]^{n2}.\end{array}$$
(71)
### List of Figures

* Linear, ring and 4-arm star topologies for flexible polymers. $`𝐑`$ indicates the end-to-end vector for the linear topology and the relative position of two segments in the case of a ring polymer. $`𝐑_𝐣`$ indicates the position of the end of the $`j`$-th arm in the star topology.
* The Random Flight Model of polymer chains. $`l`$ is the Kuhn length and $`\theta `$ is the bond angle.


text:
# Peano-like bounds for some Newton-Cotes Formulae

## 1 Introduction

In this paper we present an error analysis for some Newton-Cotes quadrature formulae. We consider Simpson’s rule, 3/8 Simpson rule and Boole’s rule. A similar error analysis for Simpson’s rule have been investigated more recently (, , , , ) with the view of obtaining bounds on the quadrature rule in terms of a variety of norms involving, at most, the first derivative. It is well known that if the mapping $`f`$ is neither four times differentiable nor is the fourth derivative $`f^{(4)}`$ bounded, then we cannot apply the classical Simpson’s quadrature formula, which, actually is one of the most used quadrature formulas in practical applications. Thus, the above mentioned analysis is important as well as the analysis presented here.

The current work brings results for the above mentioned Newton-Cotes quadrature rules giving explicit error bounds and using results from the modern theory of inequalities. The used inequalities are known in the literature as inequalities of Ostrowski-Grüss type. The error bounds are expressed in terms of second derivatives. As we have already mentioned for Simpson’s rule, the general approach used in the past involves the assumption of bounded derivatives of degree higher than two. We also mention that the obtained results can be derived using the Peano kernel theorem. In any case, these bounds are generally, but not always, better that the usual Peano error bounds (see Remarks 5, 7 and 9).

Here we do not consider composite quadrature rules since they can be formed in the usual way. However, the analysis presented here allows the determination of the partition required that would assure the accuracy the result would be within a prescribed error tolerance.

In Section 2 we established some auxiliary results. We use these results in further sections. In Section 3 we consider Simpson’s rule. In Section 4 we consider 3/8 Simpson rule and in Section 5 we consider Boole’s rule.

## 2 Preliminary results

###### Lemma 1

Let $`IR`$ be an open interval and $`a,bI,`$ $`a<b.`$ Let $`f:IR`$ be a twice differentiable function and let $`x[a,b]`$ be a fixed element. Then we have

$`f(x)(ba)(x{\displaystyle \frac{a+b}{2}})\left[f(b)f(a)\right]{\displaystyle \underset{a}{\overset{b}{}}}f(t)𝑑t`$
$`=`$ $`{\displaystyle \frac{1}{ba}}{\displaystyle \underset{a}{\overset{b}{}}}{\displaystyle \underset{a}{\overset{b}{}}}p(x,t)p(t,s)f^{\prime \prime }(s)𝑑s𝑑t,`$
where

$$p(x,t)=\{\begin{array}{c}ta,t[a,x]\\ tb,t(x,b]\end{array}.$$
(2)
Proof. Integrating by parts, we have

$$\underset{a}{\overset{b}{}}p(x,t)f^{}(t)𝑑t=f(x)(ba)\underset{a}{\overset{b}{}}f(t)𝑑t,$$

i. e.

$$f(x)=\frac{1}{ba}\underset{a}{\overset{b}{}}f(t)𝑑t+\frac{1}{ba}\underset{a}{\overset{b}{}}p(x,t)f^{}(t)𝑑t.$$

If we substitute $`ff^{}`$ in the above relation, then we get

$$f^{}(t)=\frac{1}{ba}\left[f(b)f(a)\right]+\frac{1}{ba}\underset{a}{\overset{b}{}}p(t,s)f^{\prime \prime }(s)𝑑s.$$

Thus,

$`{\displaystyle \underset{a}{\overset{b}{}}}p(x,t)f^{}(t)𝑑t`$
$`=`$ $`{\displaystyle \underset{a}{\overset{b}{}}}p(x,t)\left[{\displaystyle \frac{1}{ba}}\left[f(b)f(a)\right]+{\displaystyle \frac{1}{ba}}{\displaystyle \underset{a}{\overset{b}{}}}p(t,s)f^{\prime \prime }(s)𝑑s\right]𝑑t.`$
We also have

$$\underset{a}{\overset{b}{}}p(x,t)𝑑t=(ba)(x\frac{a+b}{2}).$$

From the above relations it follows

$`{\displaystyle \frac{1}{ba}}{\displaystyle \underset{a}{\overset{b}{}}}{\displaystyle \underset{a}{\overset{b}{}}}p(x,t)p(t,s)f^{\prime \prime }(s)𝑑s𝑑t`$
$`=`$ $`f(x)(ba)(x{\displaystyle \frac{a+b}{2}})\left[f(b)f(a)\right]{\displaystyle \underset{a}{\overset{b}{}}}f(t)𝑑t.`$
###### Lemma 2

Let $`p(x,t)`$ be defined by (2). Then we have

$`q(x,s)`$ $`=`$ $`{\displaystyle \underset{a}{\overset{b}{}}}p(x,t)p(t,s)𝑑t`$ (3)
$`=`$ $`\{\begin{array}{c}(ba)(x\frac{a+b}{2})(sa)\frac{ba}{2}(sa)^2,s[a,x]\\ (ba)(x\frac{a+b}{2})(sb)\frac{ba}{2}(sb)^2,s(x,b]\end{array},`$ (6)
where $`x[a,b].`$

Proof. For $`s[a,x]`$ we have

$`q(x,s)`$ $`=`$ $`{\displaystyle \underset{a}{\overset{s}{}}}(ta)(sb)𝑑t+{\displaystyle \underset{s}{\overset{x}{}}}(ta)(sa)𝑑t+{\displaystyle \underset{x}{\overset{b}{}}}(tb)(sa)𝑑t`$
$`=`$ $`(sb){\displaystyle \frac{(sa)^2}{2}}+(sa){\displaystyle \frac{(xa)^2(sa)^2}{2}}(sa){\displaystyle \frac{(xb)^2}{2}}`$
$`=`$ $`(ba)(x{\displaystyle \frac{a+b}{2}})(sa){\displaystyle \frac{ba}{2}}(sa)^2.`$
For $`s(x,b]`$ we have

$`q(x,s)`$ $`=`$ $`{\displaystyle \underset{a}{\overset{x}{}}}(ta)(sb)𝑑t+{\displaystyle \underset{x}{\overset{s}{}}}(tb)(sb)𝑑t+{\displaystyle \underset{s}{\overset{b}{}}}(tb)(sa)𝑑t`$
$`=`$ $`(sa){\displaystyle \frac{(xa)^2}{2}}+(sb){\displaystyle \frac{(xb)^2+(sb)^2}{2}}(sa){\displaystyle \frac{(sb)^2}{2}}`$
$`=`$ $`(ba)(x{\displaystyle \frac{a+b}{2}})(sb){\displaystyle \frac{ba}{2}}(sb)^2.`$
From the above relations we see that (3) holds.

###### Corollary 3

Let the assumptions of Lemma 2 be satisfied. Then we have

$`q(a,s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(ba)^2(sb){\displaystyle \frac{ba}{2}}(sb)^2,`$
$`q(b,s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(ba)^2(sa){\displaystyle \frac{ba}{2}}(sa)^2,`$
$`q({\displaystyle \frac{a+b}{2}},s)`$ $`=`$ $`\{\begin{array}{c}\frac{ba}{2}(sa)^2,s[a,\frac{a+b}{2}]\\ \frac{ba}{2}(sb)^2,s(\frac{a+b}{2},b]\end{array},`$
$`q({\displaystyle \frac{3a+b}{4}},s)`$ $`=`$ $`\{\begin{array}{c}\frac{1}{4}(ba)^2(sa)\frac{ba}{2}(sa)^2,s[a,\frac{3a+b}{4}]\\ \frac{1}{4}(ba)^2(sb)\frac{ba}{2}(sb)^2,s(\frac{3a+b}{4},b]\end{array},`$
$`q({\displaystyle \frac{a+3b}{4}},s)`$ $`=`$ $`\{\begin{array}{c}\frac{1}{4}(ba)^2(sa)\frac{ba}{2}(sa)^2,s[a,\frac{a+3b}{4}]\\ \frac{1}{4}(ba)^2(sb)\frac{ba}{2}(sb)^2,s(\frac{a+3b}{4},b]\end{array},`$
$`q({\displaystyle \frac{2a+b}{3}},s)`$ $`=`$ $`\{\begin{array}{c}\frac{1}{6}(ba)^2(sa)\frac{ba}{2}(sa)^2,s[a,\frac{2a+b}{3}]\\ \frac{1}{6}(ba)^2(sb)\frac{ba}{2}(sb)^2,s(\frac{2a+b}{3},b]\end{array},`$
$`q({\displaystyle \frac{a+2b}{3}},s)`$ $`=`$ $`\{\begin{array}{c}\frac{1}{6}(ba)^2(sa)\frac{ba}{2}(sa)^2,s[a,\frac{a+2b}{3}]\\ \frac{1}{6}(ba)^2(sb)\frac{ba}{2}(sb)^2,s(\frac{a+2b}{3},b]\end{array}.`$
## 3 Simpson’s rule

###### Theorem 4

Let the assumptions of Lemma 1 hold and let $`\gamma ,\mathrm{\Gamma }`$ be real numbers such that $`\gamma f^{\prime \prime }(t)\mathrm{\Gamma },`$ $`t[a,b]`$. Then we have

$$\left|\frac{ba}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right]\underset{a}{\overset{b}{}}f(t)𝑑t\right|\frac{\mathrm{\Gamma }\gamma }{162}(ba)^3.$$
(12)
Proof. For $`x=a`$ the left-hand side in (1) is equal to

$$(ba)\frac{f(b)+f(a)}{2}\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(13)
For $`x=b`$ the left-hand side in (1) is equal to

$$(ba)\frac{f(b)+f(a)}{2}\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(14)
For $`x=\frac{a+b}{2}`$ the left-hand side in (1) is equal to

$$(ba)f(\frac{a+b}{2})\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(15)
If we now multiply (15) by 2 and add (13) or (14) then we get

$$\frac{ba}{2}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right]3\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(16)
The corresponding right-hand side is

$`R(a,b)`$ $`=`$ $`{\displaystyle \frac{1}{ba}}{\displaystyle \underset{a}{\overset{b}{}}}{\displaystyle \underset{a}{\overset{b}{}}}\left[2p({\displaystyle \frac{a+b}{2}},t)+p(a,t)\right]p(t,s)f^{\prime \prime }(s)𝑑t𝑑s`$
$`=`$ $`{\displaystyle \frac{1}{ba}}{\displaystyle \underset{a}{\overset{b}{}}}\left[2q({\displaystyle \frac{a+b}{2}},s)+q(a,s)\right]f^{\prime \prime }(s)𝑑s.`$
Let $`K_1(s)=2q(\frac{a+b}{2},s)+q(a,s)`$. Then we have

$$\underset{a}{\overset{b}{}}K_1(s)𝑑s=0$$

such that

$$\frac{1}{ba}\underset{a}{\overset{b}{}}K_1(s)f^{\prime \prime }(s)𝑑s=\frac{1}{ba}\underset{a}{\overset{b}{}}K_1(s)\left[f^{\prime \prime }(s)\frac{\mathrm{\Gamma }+\gamma }{2}\right]𝑑s$$

and

$`\left|R(a,b)\right|`$ $``$ $`{\displaystyle \frac{1}{ba}}\underset{s[a,b]}{\mathrm{max}}\left|f^{\prime \prime }(s){\displaystyle \frac{\mathrm{\Gamma }+\gamma }{2}}\right|{\displaystyle \underset{a}{\overset{b}{}}}\left|K_1(s)\right|𝑑s`$
$``$ $`{\displaystyle \frac{\mathrm{\Gamma }\gamma }{2(ba)}}{\displaystyle \underset{a}{\overset{b}{}}}\left|K_1(s)\right|𝑑s,`$
since

$$\underset{s[a,b]}{\mathrm{max}}\left|f^{\prime \prime }(s)\frac{\mathrm{\Gamma }+\gamma }{2}\right|\frac{\mathrm{\Gamma }\gamma }{2}.$$
(18)
Hence,

$$\left|R(a,b)\right|\frac{\mathrm{\Gamma }\gamma }{2(ba)}\underset{a}{\overset{b}{}}\left|2q(\frac{a+b}{2},s)+q(a,s)\right|𝑑s.$$
(19)
We now calculate

$`{\displaystyle \underset{a}{\overset{b}{}}}\left|2q({\displaystyle \frac{a+b}{2}},s)+q(a,s)\right|𝑑s`$
$`=`$ $`{\displaystyle \underset{a}{\overset{\frac{a+b}{2}}{}}}\left|2{\displaystyle \frac{ba}{2}}(sa)^2{\displaystyle \frac{(ba)^2}{2}}(sb){\displaystyle \frac{ba}{2}}(sb)^2\right|𝑑s`$
$`+{\displaystyle \underset{\frac{a+b}{2}}{\overset{b}{}}}\left|2{\displaystyle \frac{ba}{2}}(sb)^2{\displaystyle \frac{(ba)^2}{2}}(sb){\displaystyle \frac{ba}{2}}(sb)^2\right|𝑑s.`$
From the equation

$$2\frac{ba}{2}(sa)^2\frac{(ba)^2}{2}(sb)\frac{ba}{2}(sb)^2=0$$

we find the solutions

$$s_1=a,s_2=\frac{2a+b}{3}.$$
(21)
From the equation

$$2\frac{ba}{2}(sb)^2\frac{(ba)^2}{2}(sb)\frac{ba}{2}(sb)^2=0$$

we find the solutions

$$s_3=b,s_4=\frac{a+2b}{3}.$$
(22)
From (3)-(22) we have

$`{\displaystyle \underset{a}{\overset{b}{}}}\left|2q({\displaystyle \frac{a+b}{2}},s)+q(a,s)\right|𝑑s`$
$`=`$ $`{\displaystyle \underset{a}{\overset{\frac{2a+b}{3}}{}}}\left[2{\displaystyle \frac{ba}{2}}(sa)^2{\displaystyle \frac{(ba)^2}{2}}(sb){\displaystyle \frac{ba}{2}}(sb)^2\right]𝑑s`$
$`+{\displaystyle \underset{\frac{2a+b}{3}}{\overset{\frac{a+b}{2}}{}}}\left[2{\displaystyle \frac{ba}{2}}(sa)^2+{\displaystyle \frac{(ba)^2}{2}}(sb)+{\displaystyle \frac{ba}{2}}(sb)^2\right]𝑑s`$
$`+{\displaystyle \underset{\frac{a+b}{2}}{\overset{\frac{a+2b}{3}}{}}}\left[2{\displaystyle \frac{ba}{2}}(sb)^2+{\displaystyle \frac{(ba)^2}{2}}(sb)+{\displaystyle \frac{ba}{2}}(sb)^2\right]𝑑s`$
$`+{\displaystyle \underset{\frac{a+2b}{3}}{\overset{b}{}}}\left[2{\displaystyle \frac{ba}{2}}(sb)^2{\displaystyle \frac{(ba)^2}{2}}(sb){\displaystyle \frac{ba}{2}}(sb)^2\right]𝑑s`$
$`=`$ $`4{\displaystyle \frac{1}{108}}(ba)^4={\displaystyle \frac{1}{27}}(ba)^4.`$
From the above relation and (16)–(19) we get

$$\left|\frac{ba}{2}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right]3\underset{a}{\overset{b}{}}f(t)𝑑t\right|\frac{\mathrm{\Gamma }\gamma }{54}(ba)^3.$$
(23)
From (23) we easily get (12).

###### Remark 5

The usual Peano error bound is

$$\left|\frac{ba}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right]\underset{a}{\overset{b}{}}f(t)𝑑t\right|\frac{(ba)^3}{81}f^{\prime \prime }_{\mathrm{}}.$$
(24)
If we choose

$$\gamma =\underset{s[a,b]}{inf}f^{\prime \prime }(s)\text{}\mathrm{\Gamma }=\underset{s[a,b]}{sup}f^{\prime \prime }(s)$$

then $`\frac{\mathrm{\Gamma }\gamma }{2}f^{\prime \prime }_{\mathrm{}}`$ and it is obvious that (12) is better than (24). In fact, these two bounds are equal if and only if $`\mathrm{\Gamma }=\gamma .`$ This case ($`\mathrm{\Gamma }=\gamma `$) is very rare in practice. Specially, if $`\mathrm{\Gamma }`$ is large and $`\mathrm{\Gamma }\gamma `$ then (12) is much better than (24).

## 4 3/8 Simpson rule

###### Theorem 6

Under the assumptions of Theorem 4 we have

$$\left|\frac{ba}{8}\left[f(a)+3f(\frac{2a+b}{3})+3f(\frac{a+2b}{3})+f(b)\right]\underset{a}{\overset{b}{}}f(t)𝑑t\right|\frac{\mathrm{\Gamma }\gamma }{384}(ba)^3.$$
(25)
Proof. If we substitute $`x=\frac{2a+b}{3}`$ in (1) then we get the corresponding left-hand side

$$\frac{ba}{6}\left[f(b)f(a)\right]+f(\frac{2a+b}{3})(ba)\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(26)
For $`x=\frac{a+2b}{3}`$ we have the corresponding left-hand side

$$\frac{ba}{6}\left[f(b)f(a)\right]+f(\frac{a+2b}{3})(ba)\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(27)
If we now multiply (26) and (27) by 3 and add (13) and (14) then we get the left-hand side of the form

$$(ba)\left[f(a)+3f(\frac{2a+b}{3})+3f(\frac{a+2b}{3})+f(b)\right]8\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(28)
The corresponding right-hand side is

$$R(a,b)=\frac{1}{ba}\underset{a}{\overset{b}{}}\left[q(a,s)+3q(\frac{2a+b}{3},s)+3q(\frac{a+2b}{3},s)+q(b,s)\right]f^{\prime \prime }(s)𝑑s.$$
(29)
Let $`K_2(s)=q(a,s)+3q(\frac{2a+b}{3},s)+3q(\frac{a+2b}{3},s)+q(b,s)`$. Then we have

$$\underset{a}{\overset{b}{}}K_2(s)𝑑s=0$$

such that

$$\frac{1}{ba}\underset{a}{\overset{b}{}}K_2(s)f^{\prime \prime }(s)𝑑s=\frac{1}{ba}\underset{a}{\overset{b}{}}K_2(s)\left[f^{\prime \prime }(s)\frac{\mathrm{\Gamma }+\gamma }{2}\right]𝑑s$$

and

$`\left|R(a,b)\right|`$ $``$ $`{\displaystyle \frac{1}{ba}}\underset{s[a,b]}{\mathrm{max}}\left|f^{\prime \prime }(s){\displaystyle \frac{\mathrm{\Gamma }+\gamma }{2}}\right|{\displaystyle \underset{a}{\overset{b}{}}}\left|K_2(s)\right|𝑑s`$
$``$ $`{\displaystyle \frac{\mathrm{\Gamma }\gamma }{2(ba)}}{\displaystyle \underset{a}{\overset{b}{}}}\left|K_2(s)\right|𝑑s,`$
since (18) holds. We now calculate

$`{\displaystyle \underset{a}{\overset{b}{}}}\left|q(a,s)+3q({\displaystyle \frac{2a+b}{3}},s)+3q({\displaystyle \frac{a+2b}{3}},s)+q(b,s)\right|𝑑s`$
$`=`$ $`{\displaystyle \underset{a}{\overset{\frac{2a+b}{3}}{}}}\left|{\displaystyle \frac{1}{2}}(ba)^3{\displaystyle \frac{1}{2}}(ba)(sb)^2{\displaystyle \frac{7}{2}}(ba)(sa)^2\right|𝑑s`$
$`+{\displaystyle \underset{\frac{2a+b}{3}}{\overset{\frac{a+2b}{3}}{}}}|2(ba)([(sa)^2+(sb)^2](ba)^3|ds`$
$`+{\displaystyle \underset{\frac{a+2b}{3}}{\overset{b}{}}}\left|{\displaystyle \frac{1}{2}}(ba)^3{\displaystyle \frac{1}{2}}(ba)(sa)^2{\displaystyle \frac{7}{2}}(ba)(sb)^2\right|𝑑s.`$
From the equation

$$\frac{1}{2}(ba)^3\frac{1}{2}(ba)(sb)^2\frac{7}{2}(ba)(sa)^2=0.$$

we find the solutions

$$s_1=a,s_2=\frac{3a+b}{4}.$$
(31)
From the equation

$$2(ba)([(sa)^2+(sb)^2](ba)^3=0$$

we find the solution

$$s_3=\frac{a+b}{2}.$$
(32)
From the equation

$$\frac{1}{2}(ba)^3\frac{1}{2}(ba)(sa)^2\frac{7}{2}(ba)(sb)^2=0$$

we find the solutions

$$s_4=b,s_5=\frac{a+3b}{4}.$$
(33)
From the above relations we get

$`{\displaystyle \underset{a}{\overset{b}{}}}\left|q(a,s)+3q({\displaystyle \frac{2a+b}{3}},s)+3q({\displaystyle \frac{a+2b}{3}},s)+q(b,s)\right|𝑑s`$
$`=`$ $`{\displaystyle \underset{a}{\overset{\frac{3a+b}{4}}{}}}\left[{\displaystyle \frac{1}{2}}(ba)^3{\displaystyle \frac{1}{2}}(ba)(sb)^2{\displaystyle \frac{7}{2}}(ba)(sa)^2\right]𝑑s`$
$`{\displaystyle \underset{\frac{3a+b}{4}}{\overset{\frac{2a+b}{3}}{}}}\left[{\displaystyle \frac{1}{2}}(ba)^3{\displaystyle \frac{1}{2}}(ba)(sb)^2{\displaystyle \frac{7}{2}}(ba)(sa)^2\right]𝑑s`$
$`+{\displaystyle \underset{\frac{2a+b}{3}}{\overset{\frac{a+2b}{3}}{}}}\left[2(ba)((sb)^2+(sa)^2)(ba)^3\right]𝑑s`$
$`{\displaystyle \underset{\frac{a+2b}{3}}{\overset{\frac{a+3b}{4}}{}}}\left[{\displaystyle \frac{1}{2}}(ba)^3{\displaystyle \frac{1}{2}}(ba)(sa)^2{\displaystyle \frac{7}{2}}(ba)(sb)^2\right]𝑑s`$
$`+{\displaystyle \underset{\frac{a+3b}{4}}{\overset{b}{}}}\left[{\displaystyle \frac{1}{2}}(ba)^3{\displaystyle \frac{1}{2}}(ba)(sa)^2{\displaystyle \frac{7}{2}}(ba)(sb)^2\right]𝑑s`$
$`=`$ $`{\displaystyle \frac{1}{24}}(ba)^4.`$
From the above relation and (28) and (4) we get

$$\left|(ba)\left[f(a)+3f(\frac{2a+b}{3})+3f(\frac{a+2b}{3})+f(b)\right]8\underset{a}{\overset{b}{}}f(t)𝑑t\right|\frac{\mathrm{\Gamma }\gamma }{48}(ba)^3.$$

This completes the proof.

###### Remark 7

The usual Peano error bound is

$$\left|\frac{ba}{8}\left[f(a)+3f(\frac{2a+b}{3})+3f(\frac{a+2b}{3})+f(b)\right]\underset{a}{\overset{b}{}}f(t)𝑑t\right|\frac{f^{\prime \prime }_{\mathrm{}}}{192}(ba)^3.$$

For the reasons given in Remark 5 the estimation obtained in Theorem 6, which is a Peano-like bound, is better than the above Peano bound.

## 5 Boole’s rule

###### Theorem 8

Under the assumptions of Theorem 6 we have

$$\begin{array}{c}|\frac{ba}{90}[7f(a)+32f(\frac{3a+b}{4})+12f(\frac{a+b}{2})+32f(\frac{a+3b}{4})+7f(b)]\\ \underset{a}{\overset{b}{}}f(t)dt|\frac{509}{\mathrm{273\hspace{0.17em}375}}(\mathrm{\Gamma }\gamma )(ba)^3.\end{array}$$
(34)
Proof. We first write left-hand sides of (1) for $`x=\frac{3a+b}{4}`$ and $`x=\frac{a+3b}{4}.`$ We have

$$f(\frac{3a+b}{4})+\frac{ba}{4}\left[f(b)f(a)\right]\underset{a}{\overset{b}{}}f(t)𝑑t$$
(35)
and

$$f(\frac{a+3b}{4})\frac{ba}{4}\left[f(b)f(a)\right]\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(36)
If we now multiply (13) and (14) by 7, (15) by 12, (35) and (36) by 32 and sum the obtained results then we get the left-hand side of the form

$$(ba)\left[7f(a)+32f(\frac{3a+b}{4})+12f(\frac{a+b}{2})+32f(\frac{a+3b}{4})+7f(b)\right]90\underset{a}{\overset{b}{}}f(t)𝑑t.$$
(37)
For the corresponding right-hand side we get

$`R(a,b)`$
$`=`$ $`{\displaystyle \frac{1}{ba}}{\displaystyle \underset{a}{\overset{b}{}}}[7q(a,s)+32q({\displaystyle \frac{3a+b}{4}},s)+12q({\displaystyle \frac{a+b}{2}},s)+32q({\displaystyle \frac{a+3b}{4}},s)`$
$`+7q(b,s)]f^{\prime \prime }(s)ds.`$
Let $`K_3(s)=7q(a,s)+32q(\frac{3a+b}{4},s)+12q(\frac{a+b}{2},s)+32q(\frac{a+3b}{4},s)+7q(b,s)`$. Then we have

$$\underset{a}{\overset{b}{}}K_3(s)𝑑s=0$$

such that

$$\frac{1}{ba}\underset{a}{\overset{b}{}}K_3(s)f^{\prime \prime }(s)𝑑s=\frac{1}{ba}\underset{a}{\overset{b}{}}K_3(s)\left[f^{\prime \prime }(s)\frac{\mathrm{\Gamma }+\gamma }{2}\right]𝑑s$$

and

$`\left|R(a,b)\right|`$ $``$ $`{\displaystyle \frac{1}{ba}}\underset{s[a,b]}{\mathrm{max}}\left|f^{\prime \prime }(s){\displaystyle \frac{\mathrm{\Gamma }+\gamma }{2}}\right|{\displaystyle \underset{a}{\overset{b}{}}}\left|K_3(s)\right|𝑑s`$
$``$ $`{\displaystyle \frac{\mathrm{\Gamma }\gamma }{2(ba)}}{\displaystyle \underset{a}{\overset{b}{}}}\left|K_3(s)\right|𝑑s,`$
since (18) holds. Hence,

$`\left|R(a,b)\right|`$ $``$ $`{\displaystyle \frac{\mathrm{\Gamma }\gamma }{2(ba)}}{\displaystyle \underset{a}{\overset{b}{}}}|7q(a,s)+32q({\displaystyle \frac{3a+b}{4}},s)+12q({\displaystyle \frac{a+b}{2}},s)+32q({\displaystyle \frac{a+3b}{4}},s)`$
$`+7q(b,s)|ds.`$
We now calculate

$`{\displaystyle \underset{a}{\overset{b}{}}}\left|7q(a,s)+32q({\displaystyle \frac{3a+b}{4}},s)+12q({\displaystyle \frac{a+b}{2}},s)+32q({\displaystyle \frac{a+3b}{4}},s)+7q(b,s)\right|𝑑s`$
$`=`$ $`{\displaystyle \underset{a}{\overset{\frac{3a+b}{4}}{}}}\left|{\displaystyle \frac{7(ba)}{2}}\left[(sb)^2+(sa)^2\right]38(ba)(sa)^2+{\displaystyle \frac{7}{2}}(ba)^3\right|𝑑s`$
$`+{\displaystyle \underset{\frac{3a+b}{4}}{\overset{\frac{a+b}{2}}{}}}\left|39(ba)\left[{\displaystyle \frac{(ba)^2}{4}}+(s{\displaystyle \frac{a+b}{2}})^2\right]{\displaystyle \frac{23}{2}}(ba)^3+6(ba)(sa)^2\right|𝑑s`$
$`+{\displaystyle \underset{\frac{a+b}{2}}{\overset{\frac{a+3b}{4}}{}}}\left|39(ba)\left[{\displaystyle \frac{(ba)^2}{4}}+(s{\displaystyle \frac{a+b}{2}})^2\right]{\displaystyle \frac{23}{2}}(ba)^3+6(ba)(sb)^2\right|𝑑s`$
$`+{\displaystyle \underset{\frac{a+3b}{4}}{\overset{b}{}}}\left|{\displaystyle \frac{7(ba)}{2}}\left[(sb)^2+(sa)^2\right]38(ba)(sb)^2+{\displaystyle \frac{7}{2}}(ba)^3\right|𝑑s.`$
Let us denote the integrands of the right-hand side of the above relation by $`Q_1(s)`$, $`Q_2(s)`$, $`Q_3(s)`$ and $`Q_4(s)`$, respectively. These integrands have the next zero points

$$s_0=a,\text{ }s_1=\frac{38a+7b}{45},\text{ }s_2=\frac{2a+b}{3},\text{ }s_3=\frac{a+2b}{3},\text{ }s_4=\frac{7a+38b}{45},\text{ }s_5=b.$$

From the above two relations we get

$`{\displaystyle \underset{a}{\overset{b}{}}}\left|7q(a,s)+32q({\displaystyle \frac{3a+b}{4}},s)+12q({\displaystyle \frac{a+b}{2}},s)+32q({\displaystyle \frac{a+3b}{4}},s)+7q(b,s)\right|𝑑s`$
$`=`$ $`{\displaystyle \underset{a}{\overset{s_1}{}}}Q_1(s)𝑑s{\displaystyle \underset{s_1}{\overset{\frac{3a+b}{4}}{}}}Q_1(s)𝑑s{\displaystyle \underset{\frac{3a+b}{4}}{\overset{s_2}{}}}Q_2(s)𝑑s+{\displaystyle \underset{s_2}{\overset{\frac{a+b}{2}}{}}}Q_2(s)𝑑s`$
$`+{\displaystyle \underset{\frac{a+b}{2}}{\overset{s_3}{}}}Q_3(s)𝑑s{\displaystyle \underset{s_3}{\overset{\frac{a+3b}{4}}{}}}Q_3(s)𝑑s{\displaystyle \underset{\frac{a+3b}{4}}{\overset{s_4}{}}}Q_4(s)𝑑s+{\displaystyle \underset{s_4}{\overset{b}{}}}Q_4(s)𝑑s`$
$`=`$ $`{\displaystyle \frac{2036}{6075}}(ba)^4.`$
Thus, we have

$$\left|R(a,b)\right|\frac{1018}{6075}(\mathrm{\Gamma }\gamma )(ba)^3.$$
(39)
From (37), (5) and (39) we easily get (34).

###### Remark 9

The usual Peano error bound is

$$\begin{array}{c}|\frac{ba}{90}[7f(a)+32f(\frac{3a+b}{4})+12f(\frac{a+b}{2})+32f(\frac{a+3b}{4})+7f(b)]\\ \underset{a}{\overset{b}{}}f(t)dt|\frac{1018}{\mathrm{273\hspace{0.17em}375}}f^{\prime \prime }_{\mathrm{}}(ba)^3.\end{array}$$

For the reasons given in Remark 5 the estimation obtained in Theorem 8, which is a Peano-like bound, is better than the above Peano bound.