format: text

text:
# Chandra OBSERVATIONS OF TWO PULSARS PSR 0628-28 AND PSR 1813-36

## 1 INTRODUCTION

Radio channel is the best way to discover neutron stars, but their characteristics are best revealed via the X-ray observations i.e. whether it’s a cooling neutron star or not. With this in mind, we observed with Chandra two ROSAT All-Sky-Survey (RASS) pulsars, that were claimed to be detected, to confirm their detection and examine their X-ray pulses and spectra. During RASS PSR 1813-36 and PSR 0628-28 had detection likelihood just below the detection threshold, Becker et al. (1992). Therefore both pulsars had follow-up ROSAT pointed observations. PSR 1813-36 was not detected and we also report a non-detection of the source for the Chandra observation. However PSR 0628-28 was detected with $``$60 counts in the ROSAT pointed observation and here we report a positive detection of a point source around the pulsar’s radio position. Due to low countrate spectral analysis did not suffice to distinguish between various pulsar models such as, powerlaw spectrum or a blackbody spectrum. In this article we report on observations acquired for both pulsars, namely PSR 1813-36 and PSR 0628-28 .

## 2 OBSERVATIONS

### 2.1 PSR 0628-28

PSR 0628-28 was observed twice, one on 2001 November 04 and the other on 2002 March 25 for 2000 s and 17000 s, respectively. For both observations photons were collected using the Advanced CCD Imaging Spectrometer (ACIS; Burke et al. (1997)), which has resolution $`\mathrm{\Delta }E/E0.1`$ at 1 keV scaling as 1/$`\sqrt{E}`$ over its 0.2-10 keV range. The pulsar was positioned on the back-illuminated S3 chip of the ACIS-S array, offset by 0.58 from the aim point, where the PSF is undersampled by the 0′′.4920$`\times `$0′′.4920 CCD pixels. Data were collected in the nominal timing mode, with 1.141s exposures between CCD readouts, and in ”FAINT” spectral mode. Data reduction and analysis were performed by CIAO data analysis software package, version 2.2.1. We reprocessed the Level 1 event data to correct the detrimental effects of charge transfer efficiency. The imaging, timing and spectral analysis presented here were done only on the 17 ks observation, the 2.0 ks observation was disregarded due to an interruption by a large solar storm. The background countrate during this solar storm increased by a factor of $``$20. The net source countrate during the first observation period was $`0.0159\pm 0.0050`$ where the error is in the 90% confidence range. Since by taking into account the 2.0 ks observation we gain only $``$27 counts which is not enough to improve our statistics significantly we prefered to disregard these counts, which potentially can be misleading.

The net countrate of the 17 ks observation is $`0.0111\pm 0.0013`$. We also analysed the ROSAT archival data and determined the ROSAT countrate of PSR 0628-28 which after normalizing to Chandra countrate is $`0.0120\pm 0.0025`$. The countrate of the portion of the observation that we used as can be seen is consistent within the errors with both 2.0 ks and ROSAT observations. The measured PSF of PSR 0628-28 is consistent with the ACIS point-source response, hence the ACIS image reveals a pointlike X-ray source at the pulsar position. The Chandra position of PSR 0628-28 is $`\alpha =06^h30^m49^{\prime \prime }.427`$, $`\delta =28^{}34^{}43^{\prime \prime }.60`$ (J2000.0), which considering 0.′′5 rms error and $``$0.′′6 absolute astrometric accuracy of Chandra , is in good agreement with the radio position of $`\alpha `$=06h30m49′′.531, $`\delta =28^{}34^{}43^{\prime \prime }.60`$ (J2000.0), Taylor et al. (1993). We should also note that the Chandra position is also consistent with the ROSAT position of the PSR 0628-28 . We also took into account the proper motion of the pulsar. The Chandra position of the pulsar is $``$0.′′5 away from the expected position at the observation epoch which is extrapolated from the proper motion measurements taken from Taylor et al. (1993).

With this positional coincidence we can firmly say that PSR 0628-28 is one of the handful of X-ray emitting radio-pulsars.

### 2.2 PSR 1813-36

PSR 1813-36 was observed on 2001 October 25 for 30 ks with ACIS-S. Confirming the previous result for the ROSAT pointing observation no source was detected at the pulsar posistion of $`\alpha `$=18h17m05′′.76, $`\delta `$=-361805′′.50 (J2000.0), Taylor et al. (1993). The observation yielded a background countrate per pixel square of 3.4$`\times `$10-6 cps/pixel2, hence for a 50 pixel2 extraction region and a signal to noise ratio of 3 the source should have minimum countrate of 2.3$`\times `$10-4 cps for detectability. Assuming a powerlaw like spectrum this countrate gives an upperlimit for the flux of 1.5$`\times `$10-15 erg cm-2 s-1 , which is in the sensitivity range of ACIS-S (0.2-10 keV), for a detector ontime of 30 ks. In the same range the implied X-ray flux from spin-down energy is 2.1$`\times `$10-15 erg s-1 , which corresponds to a luminosity of $`L_X`$=3.7$`\times `$1030 erg s-1 at an assumed distance of 3.8 kpc.

## 3 RESULTS FOR PSR 0628-28

### 3.1 Timing

To search for pulsations from PSR 0628-28 184 source counts were extracted from a 2′′ radius aperture around the Chandra position of the pulsar. The photon arrival times were barycenter corrected using CIAO axbary program. Then the $`Z_1^2`$ statistic (Buccheri et al. (1983)) was calculated over a narrow range of frequencies at the expected radio frequency ($`f=0.80358910`$ Hz, Taylor et al. (1993)) The resulting periodogram peaked at a value of $`Z_1^2=11.8961`$ at $`f=0.803352585\pm 0.000059`$ Hz referenced to the observation epoch MJD 52,358.9. Which is not consistent with the value obtained for the frequency from extrapolation of the radio ephemeris ($`f`$,$`\dot{f}`$). A wider search including the ACIS sampling frequency of $`f`$=0.8763934 Hz, showed a number of high peaks, indicating that these results are influenced by the beat frequencies between the pulsars rotation and the Chandra ACIS sampling. Continuing to investigate the folding results only around the expected pulsar frequency, we get $`Z_1^2`$=0.7 which is distrubuted like $`\chi ^2`$ with 2 dof, hence the probability of getting this value or higher is 0.713 which is not too remarkable. Looking at the higher harmonics, $`Z_2^2`$=0.92 (the sum of 1st & 2nd harmonics), which has a random probability of 0.92 and $`Z_3^2`$=10.1 which has a random probability of 0.12. The conclusion which can be drawn from this is that the data was not sufficient to claim the existence of pulsed-emission.

### 3.2 Spectral

To fit a spectral model to the pulsar data we extracted the source counts from a circle of radius 2′′ and extracted the background photons from a box of dimensions 131′′$`\times `$135′′, which was taken considerably away from the pulsar where no other sources were present. The extraction radius of 2′′ contains 90% of the events below 2 keV. The countrates for the source and background were 0.0111 cps and 0.214 cps, respectively. The CIAO tool MKARF and MKRMF were used to build the corresponding point-source mirror response matrix. The extracted photons then were binned with 0.5 keV intervals. The resulting spectra were fitted using the generalized fitting engine of CIAO, Sherpa. Two different models were tried, a blackbody and a powerlaw model. For the blackbody fit the $`N_H`$ vs. $`kT`$ error ellipses have been plotted in Figure 1.

The hydrogen column density that we obtained from the dispersion measure is $`N_H=1.07\times 10^{21}`$ cm-2, where we assumed a ratio of $`n_H/n_e=10`$. The blackbody fit had a reduced $`\chi ^2`$ of 0.29. The best fit parameter for this fit was $`kT=0.232`$ keV for the surface temperature of the neutron star and $`N_H=8\times 10^{20}`$ cm-2 for the hydrogen column density (see Table 1). This model gives an absorbed flux of 2.04$`\times `$10-14 erg cm-2 s-1 in the 0.1-10 keV range, with a distance of 2.14 kpc this flux corresponds to a X-ray luminosity of $`L_X`$=1.13$`\times `$1031 erg s-1 , in the same range. The implied source radius is very small (0.17 km); possibly from polar caps.

A powerlaw model also gave an acceptable fit. The reduced $`\chi ^2`$ minimized at a value of 0.32 for the parameters $`\gamma =2.65`$ and $`N_H=1.29\times 10^{21}`$ cm-2. This fit yielded a flux of 4.06$`\times `$10-14 erg cm-2 s-1 , which with a distance of 2.14 kpc corresponds to a X-ray luminosity of $`L_X`$=2.25$`\times `$1031 erg s-1 (see Figure 2).

## 4 SUMMARY & DISCUSSION

We report the discovery of X-ray emission from a long-period, (1.24s) radio-pulsar with a mild rotational energy loss-rate. Expressed in numbers, the measured power-law fit $`L_X`$=2.25$`\times `$1031 erg s-1 which should be compared to the 0.1% of $`\dot{E}`$, a robust formula which gives a good fit to the rotation-powered pulsars (Becker & Trümper (1997)). The 0.1% of $`\dot{E}`$ is $`10^{29}`$erg s-1 . Effectively the observed X-ray Luminosity is 100 times more than that predicted by the magnetospheric model. We should also note that the blackbody fit luminosity is about the same as the powerlaw luminosity. What is the source of this extra luminosity? PSR 0628-28 is an extreme pulsar long-period, small $`\dot{E}`$. The physical process producing the X-rays may include reheating due to vortex creep, accretion from a fall-back disc, magnetic decay like the AXPs or accretion from the ISM. Also the luminosity derived from the blackbody fit ($`L_X`$=1.13$`\times `$1031erg s-1 ) and the characteristic age ($`P/2\dot{P}`$=2.8$`\times `$106 yr) of the pulsar puts it on the standard cooling curve, see Figure 3, although it is on the fast-falling, photon-cooling dominated part.

Deeper observations, to produce sufficient counts with better timing accuracy, are needed in order to resolve the issues stated above. Chandra observations have focused our attention on this extreme pulsar.

## 5 ACKNOWLEDGEMENTS

We thank Patrick O. Slane, Werner Becker, Roberto Mignani for useful discussions. This work has been supported by NASA Grant #02500773.

## 6 REFERENCES

Becker W., J. Trümper, The X-ray luminosity of rotation-powered neutron stars, Astronomy & Astrophysics 326, 682-691, 1997.

Becker W., J. Trümper, H. Ögelman, Search for cooling neutron stars in the ROSAT survey, in Isolated Pulsars, edited by Van Riper K. A., Epstein R. and Ho C., 104, Cambridge University Press, 1992.

Buccheri R., K. Bennett , G. F. Bignami et al., Search for pulsed $`\gamma `$-ray emission from the radio pulsars in the COS-B data, Astronomy & Astrophysics 128, 245-251, 1983.

Burke, B. E., J. Gregory, M. W. Bautz et al., IEEE Trans. Electron Devices, 44, 1633, 1997.

Taylor J. H., R. N. Manchester, A. G. Lyne, Catalog of 558 pulsars, Astrop. J. Suppl. Ser. 88, 529-568, 1993.


text:
# The efficiency of the spiral-in of a black hole to the Galactic centre

## 1 Introduction

The very young objects near the Galactic centre, such as the Quintuplet star cluster (Nagata et al. 1990; Okuda et al. 1990), the Arches cluster (Nagata et al., 1995) and the central star cluster (Tamblyn & Rieke 1993; Krabbe et al. 1995) are of considerable interest. One of the more interesting conundrums is the presence of stars as young as few Myr (Tamblyn & Rieke 1993; Krabbe et al. 1995) within a parsec from the Galactic centre (Gerhard, 2001). In situ formation is problematic, due to the strong tidal field of the Galaxy, which makes this region inhospitable for star formation. One possible solution is provided by Gerhard (2001), who proposes that a $`10^6\mathrm{M}_{}`$ star cluster spirals in to the Galactic centre within a few million years from a distance $`30`$ pc. The infall process is driven by dynamical friction (Chandrasekhar, 1943). A quantitative analysis of this model by McMillan & Portegies Zwart (2002) confirms Gerhard’s results. The main uncertainty in the efficiency of dynamical friction, and therewith the time scale for spiral-in, is hidden in a single parameter, called the Coulomb logarithm $`\mathrm{ln}\mathrm{\Lambda }`$. Accurate determination of this parameter is crucial for understanding this process. Nevertheless, a precise value of $`\mathrm{ln}\mathrm{\Lambda }`$ for the Galactic central region is not available. In this paper we determine $`\mathrm{ln}\mathrm{\Lambda }`$ for the Galactic centre. We focus on the efficiency of the interaction between an intermediate mass black hole (BH hereafter) and the stars in the Galactic central region. In Section 4 we comment on how this approach can be applied to star clusters that sink to the Galactic centre.

In the classical study of Chandrasekhar (1943), dynamical friction is driven by the drag force experienced by a point mass moving through an infinite medium of homogeneous density. The consequences of finiteness and non-homogeneity have been analysed in various works (see Maoz 1993; Milosavljević & Merritt 2001). Just & Peñarrubia (2003) carried out an analytical study of dynamical friction in inhomogeneous systems, leading to a value of the Coulomb logarithm that depends on the infalling object position. Colpi & Pallavicini (1998) developed a general theoretical framework for the interaction of a satellite with a primary galaxy, able to describe dynamical friction in finite, inhomogeneous systems. They applied their theory of linear response to orbital decay of satellites onto a spherical galaxy (Colpi, 1998) and short-lived encounters with a high-speed secondary (Colpi & Pallavicini, 1998). They studied evolution of satellites in isothermal spherical haloes with cores (Colpi et al., 1999), extended in Taffoni et al. (2003), treating satellite finite size and mass loss. Still, the original expression of Chandrasekhar is used to model dynamical friction in many astronomical situations (see Binney & Tremaine 1987, Section 7.1; Hashimoto et al. 2002). The cases we study here are characterised by a point mass, with a very small mass compared to the primary system. Therefore Chandrasekhar’s formulation is appropriate in our cases.

Dynamical friction is important for a large variety of astronomical phenomena, e.g. planet migration (Goldreich & Tremaine 1980; Cionco & Brunini 2002), core collapse in dense star clusters (Portegies Zwart et al., 1999) or mergers in galaxy clusters (Makino 1997; Cora et al. 1997; van den Bosch et al. 1999). The physics of the infall process of a satellite in the parent galaxy is basically the same as in the case of a BH spiralling-in to the Galactic centre. The relevant parameters, however, are quite different in the two cases. For example, an inspiraling galaxy has finite size, whereas a BH is a point mass. Dynamical friction also plays an important role in the evolution of the black hole binary formed after the merging of two galaxies both hosting a BH at their centre (Milosavljević & Merritt, 2001). In this case, dynamical friction is important in the early phase of galaxy merging, when the BHs orbits converge and become bound.

We determine the value of $`\mathrm{ln}\mathrm{\Lambda }`$ for a BH spiralling-in to the Galactic centre by means of self-consistent $`N`$-body simulations. This is by far not an easy task. $`N`$-body models either lack in the number of particles (a direct $`N`$-body code can treat up to about $`10^5`$ particles, compared to $`10^8`$ for the real system) or have to introduce softening (Aarseth, 1963) and approximation of the force calculation (treecode (Barnes & Hut, 1986) or particle-mesh code (Hockney & Eastwood, 1988)). The softening parameter $`ϵ`$ was introduced to limit the strength of the mutual gravitational interaction during close stellar encounters. Without softening, the very high accelerations experienced by the encountering bodies would cause very tiny integration steps, which would result in a substantial freeze of the global system evolution, with consequent dramatic performance degradation. The use of such approximation does not invalidate the numerical results, as long as the simulated system is studied on a time scale shorter than the relaxation time scale (Binney & Tremaine 1987, ch.4, see also discussion in Section 3.3 below). The dynamical friction time scale of the systems we simulate is in all cases shorter than the relaxation time scale, so we can safely use the approximate methods.

Nevertheless, since close encounters have a relevant effect on dynamical friction, decreasing their strength by means of softening also decreases the strength of dynamical friction, i.e. lowers the value of $`\mathrm{ln}\mathrm{\Lambda }`$. The same role of softening is played, in the particle-mesh code, by the grid cell size $`l`$.

Our methodological approach for the present work (see Fig. 1) consists of comparing the “exact” results obtained with the direct method for low particle numbers (up to $`\mathrm{80\hspace{0.17em}000}`$) with the results of the treecode, which are less accurate and influenced by force softening, to understand how the softening $`ϵ`$ influences the results and how they have to be scaled according to the value of $`ϵ`$. Then the results of the treecode are compared to the results of the particle-mesh code, to see how softening (tree) and grid-resolution $`l`$ (particle-mesh) can be compared and scaled. Finally, having the right scaling between the different codes, we will be able to perform high particle number simulations (up to $`410^7`$) with the particle-mesh code to obtain the value of the Coulomb logarithm for the inner Galactic Bulge.

## 2 Methods and model

### 2.1 Direct method

For our direct $`N`$-body calculations we used the kira integrator module of the Starlab software environment111See: http://manybody.org (Portegies Zwart et al., 2001). Conceived and written as an independent alternative to Aarseth’s NBODY4 and NBODY5 (Aarseth, 1985, 1999), the workhorses of collisional $`N`$-body calculations for the past 25 years, kira is a high-order predictor-corrector scheme designed for simulations of collisional stellar systems. This integrator incorporates a Hermite integration scheme (Makino & Aarseth, 1992) and a block time step scheduler (McMillan, 1986) that allows homogeneous treatment of all objects in the system.

While kira is designed to operate efficiently on general-purpose computers, it achieves by far its greatest speed when combined with GRAPE-6 special purpose hardware (Makino et al., 2002)222See: http://www.astrogrape.org. For the work presented here we performed simulations with the GRAPE-6 system at the University of Tokyo with up to $`\mathrm{80\hspace{0.17em}000}`$ particles.

### 2.2 Treecode

The hierarchical treecode (Barnes & Hut, 1986) is widely used for the simulation of collisionless systems. The force on a given particle $`i`$ is computed by considering particle groups of larger and larger size as their distance from $`i`$ increases. Force contributions from such groups are evaluated by using truncated multipole expansions. The grouping is based on a hierarchical tree data structure. Such hierarchical tree is realised by inserting the particles one by one into the initially empty simulation cube. Each time two particles are into the same cube, this is divided into eight ’child’ cubes, whose linear size is one half of its parent’s. This procedure is repeated until each particle finds itself into a different cube. Hierarchically connecting such cubic cells according to their parental relation leads to the hierarchical tree data structure. When force on particle $`i`$ is computed, the tree is traversed looking for cells that satisfy an appropriate acceptance criterion; if the cell is not accepted, then its children are checked. See Salmon & Warren (1994) for a detailed overview on acceptance criteria. By applying this procedure recursively starting from the tree root, i.e. the cell containing the whole system, all the cells satisfying the acceptance criterion are found.

Our treecode simulations were initially performed with both a code written by Jun Makino (Makino, 1991), and with GADGET (Springel et al., 2001). In GADGET each particle is assigned an individual time-step, and at each iteration only those particles having an update time below a certain time are selected for force evaluation. This criterion was originally introduced in the direct $`N`$-body code (cf. Aarseth 1999). This code is parallelized using MPI (Message Passing Interface Forum, 1997). In the parallel version, the geometrical domain is partitioned, and each processor hosts the particles located in the domain partition assigned to it. The computation of forces on the selected $`i`$-particles is performed by scattering the particle data to remote processors. Then partial forces from the particles hosted by the remote processors are computed locally. Finally, calculated forces are received back by the $`i`$-particle host, and added up resulting in the total force on the $`i`$-particles.

### 2.3 Particle-mesh code

To perform calculations using several millions of particles we use a particle-mesh (PM) code named Superbox (Fellhauer et al., 2000). With the particle-mesh technique densities are derived on Cartesian grids. Using a fast Fourier transform algorithm these densities are converted into a grid-based potential. Forces acting on the particles are calculated using these grid-based potentials, making the code nearly collisionless. The current implementation completely neglects two-body relaxation causing it to retain only a small amount of grid-based relaxation (Fellhauer et al., 2000).

The adopted implementation incorporates some differences to standard PM-codes. State of the art PM-codes are using a cloud-in-cell (CIC) scheme to assign the masses of the particles to the grid-cells; which means that the mass of a particle $`i`$ is split up into neighbouring cells according to its distance to the centre of the cell. Forces are then calculated by adding up the same fractions of the forces from all cells to particle $`i`$. In contrast our code uses the “old-fashioned” nearest-grid-point scheme, where the total mass of the particle is assigned to the grid-cell the particle is located in. Forces acting on the particle are then calculated only from the forces acting on this particular cell. To achieve similar precision as CIC, we use space derivatives up to the second order.

To achieve high resolution at the places of interest, we incorporate for every simulated object (e.g. each galaxy and/or star cluster or disc, bulge and halo) two levels of sub-grids which stay focused on the objects of interest while they are moving through the simulated area (see Fig. 2). This provides higher resolution only where it is necessary. Our PM-code is parallelized using MPI.

### 2.4 The theory of the Coulomb logarithm

Dynamical friction affects a mass moving in a background sea of lower mass objects. A practical expression for the strength of the drag force on a point particle with mass $`M_{BH}`$ is (Binney & Tremaine 1987, p. 424):

$$\frac{d𝐯_{BH}}{dt}=4\pi G^2\mathrm{ln}\mathrm{\Lambda }\rho M_{BH}\frac{𝐯_{BH}}{v_{BH}^3}\left[\mathrm{erf}(X)\frac{2X}{\sqrt{\pi }}e^{X^2}\right].$$
(1)
Here $`X=v_{BH}/(\sqrt{2}\sigma )`$, $`\sigma `$ being the Maxwellian velocity dispersion, and $`\rho `$ is the background stellar density.

The classical value of $`\mathrm{\Lambda }`$ is (Binney & Tremaine 1987, p. 423)

$$\mathrm{\Lambda }=\frac{b_{max}v_{typ}^2}{G(M_{BH}+m)}.$$
(2)
Here $`b_{max}`$ is the largest possible impact parameter for an encounter between the massive point particle and a member of the background population, $`v_{typ}`$ is the typical speed of the objects in the background population, and $`m`$ is the mass of each of the background stars. Equation (2) can then be generalised to

$$\mathrm{\Lambda }=\frac{b_{max}}{b_{min}}.$$
(3)
Here $`b_{min}`$ is the distance below which an encountering particle is captured, instead of being scattered by the massive object. It is somewhat smaller than the $`90\mathrm{°}`$ turn-around distance. With the direct $`N`$-body technique, $`\mathrm{\Lambda }`$ can be measured precisely. However, with approximate $`N`$-body methods, such as the treecode or the PM code, we have to take care of the interference of the softening length/cell size with $`b_{min}`$, as discussed in Section 2.5.

McMillan & Portegies Zwart (2002) obtained an analytic expression for the distance $`r(t)`$ of the BH to the Galactic centre, with the assumptions that the BH’s orbits are nearly circular, and the mass profile of the Galaxy is given by a power law:

$$M(R)=AR^\alpha .$$
(4)
They obtained:

$$r(t)=R_0\left[1\frac{\alpha (\alpha +3)}{\alpha +1}\sqrt{\frac{G}{AR_0^{\alpha +3}}}\chi M_{BH}\mathrm{ln}\mathrm{\Lambda }t\right]^{\frac{2}{3+\alpha }},$$
(5)
where

$$\chi =\mathrm{erf}(X)\frac{2X}{\sqrt{\pi }}e^{X^2}\mathrm{and}X=\frac{v_{BH}}{\sqrt{2}\sigma },$$

$`\sigma `$ being the velocity dispersion. In McMillan & Portegies Zwart (2002) the value of $`X`$ in the Galactic centre is also computed, resulting in $`X=\sqrt{2\alpha }`$. Finally, we take $`R_0`$ equal to the half-mass radius $`R_{hm}`$ (see Section 2.6). The best fit of equation (5) on the simulation data gives the value of $`\mathrm{ln}\mathrm{\Lambda }`$ for that simulation. Such values, for all simulation performed, are reported in the last column of Tables 34 and 5.

### 2.5 The role of softening in the determination of the Coulomb logarithm

As already mentioned in the introduction, softening was introduced in numerical stellar dynamics to limit the strength of mutual forces during close stellar encounters, mainly for computational performance purposes. It consists in a modification of the Newton law for the gravity exerted by a particle $`j`$ on a particle $`i`$, as follows:

$$𝐚_i=G\frac{m_j}{(r_{ij}^2+ϵ^2)^{(3/2)}}𝐫_{ij}$$
(6)
where $`𝐫_{ij}=𝐫_j𝐫_i`$, and $`ϵ`$ is the softening parameter. As $`𝐫_{ij}0`$, the presence of $`ϵ`$ causes the force to change from inverse square to elastic, with constant $`Gm_im_j/ϵ^3`$. In this way the strength of the mutual force between encountering particles is no more unbound.

Softening was first used by Aarseth (1963) in a particle-particle (PP) context. Accuracy requirements soon led to a more precise treatment of close encounters and binaries by means of an analytic approach (Kustaanheimo & Stiefel 1965; Aarseth 1972; Mikkola & Aarseth 1990). The softened force in equation (6) remains used in the treecode, where high accuracy in close encounters treatment is not mandatory. Here we will use the softening both in the treecode simulations, where it is necessary, and in the PP code simulations, where it is used to compare the results of the two codes, in order to study the relation between $`ϵ`$ and $`\mathrm{ln}\mathrm{\Lambda }`$.

For the PM code, as described in Section 2.3, force is not computed by using the Newton force, or the softened force in equation (6). Instead, the fact that the gravitational potential on each grid point of the mesh is obtained from a density field defined on the same mesh, leads to an accuracy for the force on each particle limited by the cell size of the grid, $`l`$.

Here, we are concerned with the accuracy of the computation of the encounters experienced by a black hole spiralling-in to the Galactic centre. Since the softening (PP and treecode) and the cell size (PM code) affect this accuracy, we will use $`ϵ`$ and $`l`$ to quantify the accuracy decrement in our simulations. In Section 3.4 we will study quantitatively the dependence of $`\mathrm{ln}\mathrm{\Lambda }`$ on $`ϵ`$ and $`l`$.

The reference value for $`ϵ`$ in the work presented here will be $`ϵ_0=0.003735`$ (units given below in Table 1). This value, according to Athanassoula et al. (2000), is of the same order of magnitude as the optimal softening for a Dehnen sphere distribution (Dehnen, 1993). This distribution is similar to the power law distribution that we use, at least for what concerns the high central density peak, which is the key physical factor in the determination of the optimal softening. For an 80 000 particle distribution, $`ϵ_0`$ is about 15 times smaller than the mean inter-particle distance $`\mathrm{}`$ at the initial BH position $`R_00.87`$ (see Section 2.6). This value for $`ϵ`$ is small enough to avoid spurious effects in the force between a star and its neighbours, but is sufficient to inhibit very close encounters. The expression for $`\mathrm{}`$ can be obtained as:

$$\mathrm{}=n^{\frac{1}{3}}=\left(\frac{\rho }{m}\right)^{\frac{1}{3}}=\sqrt[3]{\frac{4\pi R^{3\alpha }}{NA\alpha }}$$
(7)
where $`n`$ is the star number density, and

$$\rho =\frac{1}{4\pi R^2}\frac{dM}{dR}=\frac{A\alpha }{4\pi }R^{\alpha 3}.$$

We used the expression in equation (4) for $`M`$, and the fact that the $`N`$ stars in the system have the same mass $`m=1/N`$.

One of the effects of softening is a damping in the BH infall at very small values of the galactocentric distance, more noticeable as $`N`$ increases. This can be explained with the fact that the inter-particle distance $`\mathrm{}`$ decreases as the BH approaches the Galactic centre (see equation 7). When $`\mathrm{}`$ becomes comparable to $`2ϵ`$, the role of softening in the force equation becomes dominant, since particles begin to ‘overlap’. With $`N=\mathrm{400\hspace{0.17em}000}`$, we get $`\mathrm{}=2ϵ_0`$ when $`R0.064`$, which is close to the value at which the damping arises, as Fig. 12 below clearly shows.

### 2.6 Initial condition

We generate the initial mass distribution according to the power law given by equation (4), with $`\alpha =1.2`$, which reproduces the mass distribution in the centre of the Galaxy, according to Mezger et al. (1999). The scale factor is $`A=4.2510^6\mathrm{M}_{}`$, corresponding to $`0.44`$ in the $`N`$-body standard units (Heggie & Mathieu, 1985), which are reported in Table 1. We use the standard units hereafter, unless other units are explicitly reported. The distributions that we generate are truncated at $`R`$ = 1.7 = 13.6 pc, with a total mass within this radius $`M_{tot}=1`$. The particles have equal mass $`m`$. Particles are assigned Maxwellian velocities, then the system is virialised to dynamical equilibrium. Then, before inserting the black hole (BH) particle, we let the system evolve for a few crossing times. The system reaches a stable configuration, whose mass profile is no more perfectly reproduced by Eq. 4. The best fit for $`A`$ and $`\alpha `$ on the mass profile of the stable configuration gives:

$`A`$ $`=`$ $`0.53`$ (8)
$`\alpha `$ $`=`$ $`0.9.`$
In fact, the mass profile having these coefficients diverges from the original one as the distance $`R`$ increases. On the other hand, in the region $`R<2`$, where we study the BH infall, the discrepancy between the two mass profiles is small. The relaxed profile values are within 10% of the initial profile values. Nevertheless, for consistency we will use the values in equation (8) for $`A`$ and $`\alpha `$ hereafter. This results in values of $`\mathrm{ln}\mathrm{\Lambda }10`$% smaller than the ones given by a mass profile with coefficients $`\alpha =1.2\mathrm{and}A=0.44`$.

The BH particle is placed at the half-mass radius $`R_{hm}0.87`$ with a circular orbit velocity, and its mass is $`M_{BH}=0.000528`$. The background particles number varies from 16 000 to 2 million. The low particle number simulations are performed with the PP code, the intermediate and high number simulations with the treecode and the PM code. Fig. 3 shows the particle range for each code. This allows us to span a large range in particle number, so that the influence of granularity in the BH motion towards the Galaxy centre can be studied.

In contrast to the other models, we choose physical units for the PM code simulations. The conversion factors from physical units to $`N`$-body units are shown in Table 1,

where l.u. denotes the unit length in $`N`$-body units, m.u. the unit mass, v.u. the unit velocity and t.u. the unit time.

The parameters of the PM calculations are chosen in the following way: the grid sizes are kept constant at

$`R_{\mathrm{system}}`$ $`=`$ $`140.0\mathrm{pc}`$
$`R_{\mathrm{out}}`$ $`=`$ $`9.6\mathrm{pc}`$ (9)
$`R_{\mathrm{core}}`$ $`=`$ $`2.4\mathrm{pc}`$
and are focussed on the center of mass of the “bulge” model, as sketched in Fig. 2. To change the resolution we alter the number of grid cells per dimension from $`32`$ up to $`128`$. With this choice the cell sizes listed in Table 2 are achieved.

To avoid a self-acceleration of the black hole we choose $`H_{0,0,0}=1.0`$ of the Greens-function, which worsens the energy conservation of the code (see Fellhauer et al. 2000) especially in simulations with high particle numbers per grid cell, but suppresses the forces of particles in the same cell on each other and on themselves.

To speed up the simulations, the time step in the PM code simulations should be as large as possible, but small enough to prevent spurious results. Therefore we started with a time step of 1,000 yr and reduced it to 200 and 50 yr. The results of the 200 yr and 50 yr time step do not differ from each other, therefore the global time step is chosen to be 200 yr. Conversely, the time step in the treecode and direct code simulations is variable and different for each particle. Time step values are in this case in the range 2–30 000 yr, with about 90% of them in the range 100–300 yr.

## 3 Results

We will now study the dependence of our results on the number of particles $`N`$ in Section 3.1, and compare the various $`N`$-body methods with identical initial realisations in Section 3.2. After having convinced ourselves that the various techniques produce consistent results, we continue by studying the effect of softening/cell size (Section 3.3) and black hole mass (Section 3.5) on the value of the Coulomb logarithm in the inner part of the Galaxy.

Our simulations aimed at several goals. 1) understanding the scaling of the system dynamics with respect to the number of particles $`N`$, and within this scaling, how results from different methods compare with each other. 2) How, at a fixed value of $`N`$, the softening parameter influences the dynamics, changing the value of $`\mathrm{ln}\mathrm{\Lambda }`$. The particle-mesh method does not make use of softening. The cell size in the PM code can be seen in this context as a softening length. In our framework, it is crucial to understand the relation between the PP code and treecode softening parameter and the PM code cell size. 3) We also study how the BH mass influences the infall time. We doubled and quadrupled the BH mass, and observed how this affects the value of $`\mathrm{ln}\mathrm{\Lambda }`$.

A resume of all the runs that we performed is reported in Table 3 for the PP code runs, Table 4 for the treecode runs, and finally Table 5 for the PM code runs. In all of our runs, the system remains in equilibrium during the whole BH infall, with no significant mass loss from stellar escapes, and a mass profile independent of time.

### 3.1 Dependence of $`\mathrm{𝐥𝐧}𝚲`$ on $`𝑵`$

In order to obtain a precise measure of $`\mathrm{ln}\mathrm{\Lambda }`$, ideally one would run a direct $`N`$-body simulation with $`N`$ of the order of the number of stars in the Galactic bulge, which amounts to $`10^8`$. Such high number makes a direct simulation unfeasible, and imposes the use of approximate methods instead. In order to evaluate the reliability of the approximate methods, we compared the PP code runs with the treecode runs. The PP code runs give a reliable picture of the system dynamics at low particle numbers, i.e. at high granularity. Using the treecode we can reach a much highernumber of particles, up to two million, which still is two orders of magnitude lower than the real system. A comparison of the results from the two methods allows us to estimate the validity of the treecode runs, up to 2 million particles. Then we can compare the treecode runs and the PM runs, in order to validate the results from the latter, which has the capability to simulate systems of about 100 million stars. In such a way we will eventually be able to study the infall of a BH into the Galactic centre in a simulation environment with a realistic value of $`N`$.

In Fig. 4 we show the evolution of the BH distance from the centre of mass of the system for three treecode simulations. $`N`$ varies from 80 000 to 400 000 and 2 million, with $`ϵ=ϵ_0=0.003735`$, corresponding to about 0.03 pc. In Fig. 5 we present a similar figure from PM code simulations. Here is $`N\{\mathrm{80\hspace{0.17em}000},\mathrm{400\hspace{0.17em}000},\mathrm{2\hspace{0.17em}000\hspace{0.17em}000}\}`$, with 32 cells per dimension, resulting in a cell size of about 0.69 pc.

Figs 4 and 5 show that increasing $`N`$ results in a much smoother motion of the BH in its infall towards the centre of the Galaxy. The BH infall rate is not much affected by a change in $`N`$. Accordingly, the value of $`\mathrm{ln}\mathrm{\Lambda }`$ for each of the two sets above is consistent, as values in Table 4 (first three rows) and Table 5 (rows with $`l=0.69`$) show.

In order to study further the extent of the influence of $`N`$ on the infall rate of the BH, and hence in $`\mathrm{ln}\mathrm{\Lambda }`$, we compare in Fig. 6 results from PM code simulations with increasing grid refinement, and extreme difference in $`N`$. To quantify the grid resolution, we use the cell length at intermediate refinement, which is the cell length pertaining to the physical region where the BH evolves for most of its infall. We measure this length in units of $`ϵ_0=0.003735`$, which makes the comparison with the softening parameter of the treecode easier. $`N`$ has no remarkable influence on the infall rate, except for the case where the cell size is $`l=0.15\mathrm{pc}5ϵ_0`$. In this case the simulation with $`N=\mathrm{80\hspace{0.17em}000}`$ ((data not reported in the figure), shows an incorrect BH infall, comparable to the case $`l=0.3210ϵ_0`$. This can be explained by the fact that in the low $`l`$, low $`N`$ case, the cells are so small, and the particles so few, that many cells in the PM grid are empty (see also the $`N_c`$ column in table 5, which gives the average number of particles per cell). When $`N_c1`$, the density field is incorrect, with many grid points having a null value, because the corresponding cell is empty. In this condition, the gravity field computed by the PM code becomes unreliable, affecting the numerical results, as in the simulation with $`N=\mathrm{80\hspace{0.17em}000}`$ and $`l5ϵ_0`$.

### 3.2 Comparison of the codes

In this section we compare the results obtained from the various codes, to check their consistency. The comparison of the PM results with the two other codes results is particularly critical, since the PM code computes forces using a different mathematical approach, i.e. a grid based force derivation vs a direct particle-particle computation for the PP code, or particle-multipole computation for the treecode. A consequence of this is a different parameter to tune the accuracy of the simulation, namely the cell size $`l`$ for the PM code, and the softening length $`ϵ`$ for the other two codes. We will study here how these two parameters influence the black hole infall.

In Fig. 7 we show the time evolution of the galactocentric BH distance $`R`$ simulated by the PP code, accompanied by a plot of the time evolution of $`\mathrm{\Delta }R/R_{PP}`$ for treecode and PM simulations, where $`\mathrm{\Delta }R=(RR_{PP})`$. The relative difference $`\mathrm{\Delta }R/R_{PP}`$ remains small for a large fraction of the infall, and the final discrepancy is mostly due to the small values of the quantities at that point, which are likely to amplify relative differences. As the following figures also show, the BH infall is predicted with very good consistency among the codes.

In Fig. 8 selected treecode runs with $`N=\mathrm{80\hspace{0.17em}000}`$ and increasing $`ϵ`$ are compared with the direct code runs having the same values of $`N`$ and $`ϵ`$. At the same time, the figure shows how the infall time increases (and implicitly how $`\mathrm{ln}\mathrm{\Lambda }`$ decreases), as $`ϵ`$ increases. Fig. 8 and Table 6 show that the results from the treecode, the PP code, and the PM code are in good agreement. The agreement of the results from the three methods, and the scaling of $`\mathrm{ln}\mathrm{\Lambda }`$ with $`ϵ`$, will be further studied quantitatively in Section 3.4.

In order to understand how the cell length $`l`$ of the PM code and the softening parameter $`ϵ`$ of the PP code and treecode relate with each other, we compare in Fig. 9 the results from the PM code and treecode simulations with $`\mathrm{80\hspace{0.17em}000}`$ particles. The BH infall as shown in Fig. 9 depends on the value of $`l`$ or $`ϵ`$. Remarkably, $`l`$ and $`ϵ`$ seem to play the same role not only qualitatively, but also quantitatively: in a PM run, a given value of $`l`$ induces an infall which is very similar to the infall, in a treecode run, with $`ϵ`$ assuming that same value. In Section 3.4 this relation will be studied further.

### 3.3 The effect of softening/grid

The influence of the softening parameter on the BH dynamics has been studied by running a number of simulations with the three codes. In Table 6 we report the value of $`\mathrm{ln}\mathrm{\Lambda }`$ obtained from our simulations. For the PP code and treecode simulations, we increase $`ϵ`$ from 0 to $`32ϵ_0=0.11950.96`$ pc. For the PM code, we increase $`l`$ from 2.5 $`ϵ_0`$ to 23 $`ϵ_0`$. In all cases is $`M_{BH}=0.52810^3\mathrm{62\hspace{0.17em}300}`$ M.

For the PP code and the treecode, we selected $`N=\mathrm{80\hspace{0.17em}000}`$ as a suitable value. The relaxation time is for this value of $`N`$ $`t_x0.1N/\mathrm{ln}NR/v_{typ}2000`$, about one order of magnitude larger than the typical BH infall time, so that the system is collisionless, and we can confidently use the treecode to simulate it. With this choice for $`N`$, the BH mass is $`M_{BH}/m42.3`$, (see Table 4). As a cross-check, we ran two PP runs with $`N=\mathrm{16\hspace{0.17em}000}`$ which, as expected, gave incorrect results (see Table 3). This is due to both a too small $`M_{BH}/m`$ ratio, and a too short relaxation time ($`t_x400`$ in this case). We did not increase $`ϵ`$ above $`0.12`$, since at this point $`ϵ`$ is already much bigger than $`b_{min}`$ (see Table 4), and the infall time is now close to $`t_x`$.

For the PM code simulations, we used $`N=2`$ million in order to have enough particles to fill all the cells, even for the simulations with a small $`l`$. As Table 5 shows, for $`l=0.076\mathrm{pc}2.5ϵ_0`$ the average number of particles per cell is already $`N_c=0.1`$. Since a PM simulation gives incorrect results for $`N_c1`$ (see also the discussion at the end of Section 3.1), we did not decrease $`l`$ below 2.5 $`ϵ_0`$.

The decrease of the value of $`\mathrm{ln}\mathrm{\Lambda }`$ as $`ϵ`$ or $`l`$ increases is clear from Table 6. In the next sub-section we focus on the relation between $`\mathrm{\Lambda }`$ and $`ϵ`$, and provide a fitting formula for $`\mathrm{ln}\mathrm{\Lambda }(ϵ)`$. We use hereafter $`ϵ`$ to refer either to the softening length of the PP and treecode, or the cell size of the PM code. As shown on Fig. 9 and discussed above, these two parameters play the same role even quantitatively in affecting $`\mathrm{ln}\mathrm{\Lambda }`$. In this respect, we refer to $`ϵ`$ as a generic accuracy parameter.

### 3.4 Determination of $`\mathrm{𝐥𝐧}𝚲`$

We will study here the relation between $`ϵ`$ and $`\mathrm{ln}\mathrm{\Lambda }`$. As just said above, in this context $`ϵ`$ will be used as the accuracy parameter, and it will refer to either the softening length used in the PP and treecode, or to the cell size in the PM code.

A mathematical expression for the relation between $`ϵ`$ and $`\mathrm{ln}\mathrm{\Lambda }`$ can be found by considering how softening affects two body scattering. The role of $`ϵ`$ is to prevent too close stellar encounters. In this respect, the effect of introducing a softening length is to increase the minimal impact parameter. Hence, we can define an effective impact parameter $`b_{eff}=b_{min}+ϵ`$, and we modify equation (3) to become:

$$\mathrm{ln}\mathrm{\Lambda }=\mathrm{ln}\frac{b_{max}}{b_{eff}}=\mathrm{ln}\frac{b_{max}}{b_{min}+ϵ}.$$
(10)
We will now fit this equation with the values reported in Table 6. In order to perform the fit, we change equation (10) in a more suitable form, as follows:

$$\mathrm{ln}\mathrm{\Lambda }=\mathrm{ln}b_{max}\mathrm{ln}(b_{min}+ϵ)=K\mathrm{ln}(a+ϵ).$$
(11)
We will refer to $`b_{max}`$ and $`b_{min}`$ as the theoretical values of the maximal and minimal impact parameters, as they can be obtained from equations (2) and (3), and $`K`$ and $`a`$ as the corresponding experimental values obtained with the fit.

The best fits for $`K`$ and $`a`$ with respect to simulation values are reported in Table 7 for all codes. Such fits have been performed with a fixed value for $`R_0`$, i.e. the $`R_0=R_{hm}`$. In fact, the not perfectly circular orbit of the BH results in an oscillatory behaviour for the BH galactocentric radius. In this case, having $`R_0`$ fixed could not be an appropriate choice for the fit. We checked whether having $`R_0`$ as a free parameter in the fit leads to different results in $`\mathrm{ln}\mathrm{\Lambda }`$. We obtained values for $`\mathrm{ln}\mathrm{\Lambda }`$ within the error bars in Fig. 10, and values for $`R_0`$ within $`R_{hm}\pm 0.05`$. We can conclude that, although the galactocentric BH radius does not decreases smoothly, but in an oscillatory fashion, having $`R_0`$ fixed to the actual initial BH radius in the simulations leads to correct fits for the value of $`\mathrm{ln}\mathrm{\Lambda }`$. With respect to the PM values, a further peculiarity is that when the BH enters the finest grid area, i.e. approximately at $`R=0.3`$, the value of $`l`$ decreases (see Section 2.3 and Fig. 2). This causes $`b_{eff}`$ to become smaller, increasing the value of $`\mathrm{ln}\mathrm{\Lambda }`$. In fact, a fit of the PM data limited to values of $`R>0.3`$ gives values of $`\mathrm{ln}\mathrm{\Lambda }`$ systematically higher by $`0.32\mathrm{\Delta }(\mathrm{ln}\mathrm{\Lambda })`$.

From the PP code value of $`K`$ in Table 7 we obtain for $`b_{max}`$ the experimental value $`b_{max}^E=e^K0.39`$. This value is quite smaller than what one would expect. Since $`b_{max}`$ has the meaning of the maximal impact parameter, a natural choice is to assign it a value of the order of the system size, which in our case would result in $`b_{max}=2`$. The maximal radius for dynamical friction in our system is then about one quarter of what it is customarily assumed. PM and treecode values are slightly higher, but still quite smaller than $`b_{max}=2`$. Also $`a`$ is smaller than the theoretical value $`b_{min}=G(M_{BH}+m)/v_{typ}^2=1.4110^3`$, by a factor 3. The $`a`$ value for all codes is perfectly consistent.

An explanation for the discrepancy between the values of $`b_{max}`$ and $`b_{max}^E`$ is that the BH, while moving to the Galactic centre, is off-centre with respect to the density peak (in fact the BH is spiralising towards it). With respect to the BH position, the density distribution is then asymmetric. This density peak has clearly a greater influence on the BH dynamics, contributing more than the other regions of the system to the dynamical friction on the BH. This leads to a value of $`b_{max}`$ affected by the galactocentric BH radius. This approach is studied in detail by Hashimoto et al. (2002), who propose the galactocentric radius as a value for $`b_{max}`$ in the context of the spiralisation of satellite galaxies.

In our simulations, the galactocentric radius varies from $`R0.9`$ at the beginning of a simulation, to $`R0`$ at the end of it. The value of $`b_{max}^E`$ that we find is within this range, and it can be interpreted as an order 0 estimate of a maximal impact parameter that depends on the galactocentric BH radius.

In order to explore this aspect further, we simulated the infall of the same BH, starting at the quarter mass radius $`R_{qm}0.43`$, for $`ϵ`$ ranging from 0 to $`16ϵ_0`$. What we expect is a smaller value of $`b_{max}^E`$, hence smaller values of $`\mathrm{ln}\mathrm{\Lambda }`$. All simulations are performed with the treecode, except for the $`ϵ=0`$ case, which is simulated with the PP code. Our results are in Fig. 11.

We can see there how smaller are the values of $`\mathrm{ln}\mathrm{\Lambda }`$ for the cases when the BH starts at the quarter mass radius. A fit on these data gives $`K1.1`$, which implies $`b_{max}^E0.33`$, which is smaller than the value of $`b_{max}^E`$ obtained for the BH starting from the half mass radius. Our findings support the argument of Hashimoto et al. (2002).

### 3.5 Varying black hole mass

We also studied the effect of a variable BH mass on the value of $`\mathrm{ln}\mathrm{\Lambda }`$. We simulated, using the treecode, the infall of a BH of mass two times and four times larger than the default mass $`M_0=0.52810^30.6210^5\mathrm{M}_{}`$. We studied this infall in both the 80 000 particles configuration, and the 400 000 particles configuration. In all cases, we used our standard value for $`ϵ`$, i.e. $`ϵ_0=0.003735`$. In Fig. 12 the distance $`r`$ of the BH from the centre of mass of the system is shown for all the cases mentioned above, together with the $`M_{BH}=M_0`$ cases. From equation (11) and Table 7, the appropriate value for $`\mathrm{ln}\mathrm{\Lambda }`$ in the above cases is:

$`\mathrm{ln}\mathrm{\Lambda }`$ $`=`$ $`K\mathrm{ln}(a+ϵ_0)\pm \mathrm{\Delta }=`$
$`0.64\mathrm{ln}(0.00088+0.003735)\pm 0.34.7\pm 0.3.`$
We also show in Fig. 12 the analytic curve $`r(t)`$, as given by equation (5), with $`\mathrm{ln}\mathrm{\Lambda }=4.7`$. An error bar gives, for each analytical curve, the spread corresponding to a variance $`\mathrm{\Delta }(\mathrm{ln}\mathrm{\Lambda })=0.3`$.

The results shown in Fig. 12 are consistent with the hypothesis that a variation in the BH mass has a little effect in the value of $`\mathrm{ln}\mathrm{\Lambda }`$. In fact, $`\mathrm{ln}\mathrm{\Lambda }`$ shows a logarithmic dependence on $`M_{BH}`$ through the parameter $`b_{min}`$, which depends linearly on $`M_{BH}`$ (see equations (2) and (3)). Assuming that also the experimental value $`a`$ depends linearly on $`M_{BH}`$, we obtain $`\mathrm{ln}\mathrm{\Lambda }4.6\pm 0.3`$, and $`\mathrm{ln}\mathrm{\Lambda }4.4\pm 0.3`$ respectively for the 2 $`M_{BH}`$ case and the 4 $`M_{BH}`$ case. This results in a small displacement towards the right of the corresponding analytic curves in Fig. 12, which does not affect the conclusions that can be drawn from the figure. The theoretical curve fits very well with the $`M=2M_0`$, $`N=\mathrm{400\hspace{0.17em}000}`$ case. The other simulation curves are within, or very close to, the error in $`r(t)`$ associated to the error in $`\mathrm{ln}\mathrm{\Lambda }`$. We can conclude that a variation in the mass of the infalling object has little influence in the value of $`\mathrm{ln}\mathrm{\Lambda }`$, which is important in view of extending this work to the case of the infall of a star cluster.

The fitting formula for $`\mathrm{ln}\mathrm{\Lambda }`$ vs $`ϵ`$ was obtained from simulations with $`M_{BH}=M_0`$. This formula predicts $`\mathrm{ln}\mathrm{\Lambda }`$ for the cases with $`M_{BH}>M_0`$ with a very good accuracy, showing that it can be applied in a more general context, in order to forecast the value of the Coulomb logarithm.

Fig. 12 also shows a damping in the BH infall at very small values of $`R`$, especially for the $`N=\mathrm{400\hspace{0.17em}000}`$ case. This effect, described in Section 2.5, is clearer in the $`N=\mathrm{400\hspace{0.17em}000}`$ case, since the particle density is higher in this case, compared to the $`N=\mathrm{80\hspace{0.17em}000}`$ case.

### 3.6 Comparison with related work

Milosavljević & Merritt (2001) study the dynamical evolution of two black holes, each one at the centre of a power law cusped galaxy core. They simulate the merging of the two galaxies, which results in the two black holes forming a hard binary at the centre of the merged galaxy. In Section 3 of their paper they discuss the value of $`\mathrm{ln}\mathrm{\Lambda }`$ in their simulations. They measure the decay rate of the two black holes, and compare this value with theoretical estimates. When they compare their experimental decay rate with an estimate for the case of the infall of an isolated black hole, they find a theoretical estimate about 6 times lower than the measured value, under the assumption that $`\mathrm{ln}\mathrm{\Lambda }1.6`$. If the value of $`\mathrm{ln}\mathrm{\Lambda }`$ is not theoretically pre-determined, and is instead obtained from the decay rate equation, the result is $`\mathrm{ln}\mathrm{\Lambda }10`$. Similarly, they compare the experimental value with an estimate for the case of two mutually spiralling spherical distributions of matter. In this case they assume $`\mathrm{ln}\mathrm{\Lambda }1.0`$, and obtain an estimate for the decay rate about a factor of 2 lower than the observed value. Determining $`\mathrm{ln}\mathrm{\Lambda }`$ from the measurement would give in this case $`\mathrm{ln}\mathrm{\Lambda }1.87`$. The values of $`\mathrm{ln}\mathrm{\Lambda }`$ that we find are between the two values above.

The value for $`\mathrm{ln}\mathrm{\Lambda }1`$ that they assume in their theoretical estimates, comes from a derivation that they present in appendix A of the same work. This derivation is based on results of Maoz (1993). Under the assumption that the stellar density obeys a power law centered on the BH position:

$$\rho (r)=\rho _0\left(\frac{r}{b_{min}}\right)^\alpha ,$$
(12)
they obtain $`\mathrm{\Lambda }1/\alpha 1`$, which actually implies $`b_{max}b_{min}`$, whereas it is customary to consider $`b_{max}b_{min}`$.

Their assumption in fact is valid only when the BH is close to the centre of the power law distribution. In their context this is true when: 1) the separation between the two BHs is much larger than the half mass radius of the two galaxies. In this case each BH is at the centre of its own galaxy, and at the same time its motion is not yet heavily perturbed by the other galaxy. 2) the BH binary has hardened, and occupies the centre of the merged galaxy.

During the transient phase, when the two BHs have not yet formed a binary, the density distribution that affects the motion of the BHs is double-cusped, with a BH in each of the two cusps. This is substantially different from the density distribution modelled by equation (12).

This qualitative argument would make the density distribution in equation (12) inapplicable during the transient phase, and could explain why Milosavljević & Merritt (2001) find a higher than expected value of $`\mathrm{ln}\mathrm{\Lambda }`$ in the transient. The analytical evaluation of $`\mathrm{ln}\mathrm{\Lambda }`$ according to the technique used by them is by no means trivial, when symmetry arguments cannot be straightforwardly applied. We will address this issue in future developments of the present work; the theory of linear response of Colpi & Pallavicini (1998) could be very useful in this context.

## 4 Applications to star clusters

Recent observations of the Galactic Centre have revealed a population of very young clusters with ages less than 10 Myr. The presence of such stars inside the inner pc of the Galaxy is puzzling, as the strong tidal field in the Galactic centre easily prevents star formation. The origin of these stars is therefore vividly debated (Gerhard 2001; McMillan & Portegies Zwart 2002; Kim & Morris 2002). Morris (1993) proposed that a star cluster at some distance from the Galactic centre could spiral-in due to dynamical friction (see also Gerhard (2001)). The efficiency of dynamical friction depends sensitively on the actual value of the Coulomb logarithm $`\mathrm{ln}\mathrm{\Lambda }`$.

### 4.1 Sinking of massive black holes in the Galactic centre

We performed $`N`$-body simulations for a large range of conditions. In Section 3.1 we varied the number of particles, in Section 3.3 we varied the size of the object, and in Section 3.5 we varied its mass. With direct $`N`$-body simulations we measured the actual value of the Coulomb logarithm $`\mathrm{ln}\mathrm{\Lambda }`$. We study the behaviour of $`\mathrm{ln}\mathrm{\Lambda }`$ for various types of $`N`$-body solvers and particle numbers. We also study the behaviour of $`\mathrm{ln}\mathrm{\Lambda }`$ as a function of the softening length $`ϵ`$. Only the direct $`N`$-body code can perform a true measurement of the Coulomb logarithm, because it is able to resolve even the smallest length scales and time scales. This, however, makes the direct code very slow and, even using the very fast GRAPE-6 special purpose device, we are able to perform simulations with only $`10^5`$ particles. This is a small number compared to the actual number of stars in the Galactic centre. With approximate methods (treecode and particle-mesh) we are able to increase the number of particles up to 2 million. The cost of this is a lower accuracy in calculating stellar motion below a typical length scale $`ϵ`$. We studied how this length scale influences $`\mathrm{ln}\mathrm{\Lambda }`$, by affecting the value of the minimal impact parameter.

### 4.2 Young dense clusters in the Galactic centre

The study of the dependence of $`\mathrm{ln}\mathrm{\Lambda }`$ on $`ϵ`$ described above is also of astronomical interest, because $`ϵ`$ can be interpreted as the typical length of a finite size infalling object. Based on this, our analysis of the dependence of $`\mathrm{ln}\mathrm{\Lambda }`$ on $`ϵ`$ can be seen as a first approach to the study of the infall of a star cluster of typical size $`ϵ`$ toward the Galactic centre. We found (see Fig. 10) that the value of $`\mathrm{ln}\mathrm{\Lambda }`$ decreases quite rapidly as $`ϵ`$ increases, with the logarithm argument $`\mathrm{\Lambda }1/ϵ`$. The typical size of the compact young clusters observed in the Galactic bulge is $`0.3`$ pc (Figer et al., 1999), which corresponds to $`ϵ10ϵ_0`$. With this value of $`ϵ`$, from equation (11) and Table 7, we obtain $`\mathrm{ln}\mathrm{\Lambda }2.9`$, about 60% less than the value for a point mass. The infall time is roughly doubled. For our choice of object mass, $`M`$ 62 300 M, and initial galactocentric radius, $`R_07`$ pc, we have an infall time that increases from $`6`$ Myr for the point mass, to $`12.5`$ Myr for an object of typical size $`10ϵ_00.3`$ pc.

We also studied the uncertainty associated with the maximal impact parameter $`b_{max}`$. We found that for an infall to the Galactic centre, the infalling object is mostly influenced by the density peak at the Galactic centre itself. A good choice for $`b_{max}`$ is then $`b_{max}\beta R_0`$, where $`R_0`$ is the initial galactocentric radius, and $`\beta 0.5`$

## 5 Conclusions

We simulated the evolution of a massive particle in a sea of lighter particles in a self gravitating system. The main goal of this work is to obtain an accurate value of the Coulomb logarithm ($`\mathrm{ln}\mathrm{\Lambda }`$). This helps us to understand the dynamics of the Galactic bulge and the rate at which intermediate mass black holes sink to the Galactic centre. We also study the effect of the finite size of the inspiraling object.

We ran both $`N`$-body particle-particle (PP) simulations, softened treecode simulations, and particle-mesh (PM) simulations. The comparative simulations are performed for 80 000 particles, and all result in the same value of $`\mathrm{ln}\mathrm{\Lambda }`$. For a point particle near the Galactic centre we find $`\mathrm{ln}\mathrm{\Lambda }=6.6\pm 0.6`$. In addition we measure the change in the Coulomb logarithm with respect to the softening parameter $`ϵ`$, which reveals $`\mathrm{\Lambda }1/ϵ`$. We argue that $`ϵ`$ can be interpreted as the typical length of a finite size object, such as a star cluster, so that $`\mathrm{ln}\mathrm{\Lambda }`$ as a function of $`ϵ`$ can be seen as a first approximation of the dependence of the Coulomb logarithm on the size of an infalling star cluster.

We also observed a value of the maximal impact parameter $`b_{max}`$ different from the customarily assumed value, which is proportional to the system size. We found that our results are more consistent to a value of $`b_{max}`$ linearly dependent on the BH galactocentric radius, which is in agreement with Hashimoto et al. (2002).

We performed simulations with up to 2 million particles using a treecode. The obtained value of $`\mathrm{ln}\mathrm{\Lambda }`$ does not depend on the number of particles. Apparently, 80 000 particles is already enough to eliminate any granularity for our choice of initial conditions. The results of the treecode, at the low $`N`$-limit, are in excellent agreement with the PP simulations, and we find the same scaling with respect to $`ϵ`$. Increasing the black hole mass reduces the time scale for spiral-in as was expected from theory (see McMillan & Portegies Zwart 2002).

Finally we compared the results of our PP and treecode simulations with a particle mesh (PM) method. We compared the methods for $`N`$ up to 2 million. The results of our PP, treecode, and PM calculations are in good agreement. The cell size in the PM model is directly comparable to the softening length $`ϵ`$ in the PP and tree methods.

## acknowledgements

We thank Douglas Heggie, Piet Hut, Jun Makino, and Steve McMillan for stimulating discussions. We acknowledge KNAW under grant 95-BTN-15, NWO under Spinoza grant 08-0, NOVA under grant V-37, DFG-grant FE564/1-1, EC through grant HPRI-1999-CT-00026. SPZ is a KNAW fellow.


text:
# Renormalization of Hamiltonian Field Theory; a non-perturbative and non-unitarity approach

## 1 Introduction

The power of Hamiltonian methods is well known from the study of non-relativistic many-body systems and from strongly-interacting few particle systems, even though a Lagrangian approach is usually chosen for relativistic theories. Hamiltonian methods for strongly-interacting systems are intrinsically non-perturbative and usually contain a Tamm-Dancoff type approximation, in the sense that, in practice, one has to limit the bound state as an expansion over states containing a small number of particles. This truncation of the Fock space gives rise to a new class of non-perturbative divergences, since the truncation does not allow us to take into account all diagrams for any given order in perturbation theory. Therefore renormalization issues have to be considered carefully. Two very different remedies for this issue are the use of light-front Tamm-Dancoff field theory (LFFT)  and the application of the coupled cluster method (CCM) . In the LFFT the quantization plane is chosen to coincide with the light front, therefore the divergences that plagued the original theories seem to disappear . Furthermore, not having to include interactions in boost operators allows a renormalizable truncation scheme . One of the most important difficulties in LFFT is the complicated structure of the renormalization process . In the standard form of CCM, on the other hand, the amplitudes obey a system of coupled non-linear equations which contain some ill-defined terms because of ultraviolet divergences. It has been shown  that the ill-defined amplitudes, which are also called critical topologies, can be systematically removed, by exploiting the linked-cluster property of the ground state. This can be done by introducing a mapping which transfers them into a finite representations without making any approximation such as a coupling expansion. Thus far this resummation method has been restricted to superrenormalizable theories due to its complexity.

A natural question in the renormalization within the Hamiltonian formalism arises, since one could also perform renormalization in a Lagrangian framework and finally construct the corresponding Hamiltonian by means of standard Legendre transformation. \[It should be noticed that this is not generally applicable if the Lagrangian contains higher-order time-derivatives.\] The Hamiltonian formalism, despite a certain lack of elegance, has the advantage that it is very economical, and one can use all the know-how of quantum many-body theory.

In the last decade extensive attempts have been made to give a prescription for renormalization within the Hamiltonian formalism . Commonly unitary transformations are used to decouple the high- and low-energy modes aiming at the partial diagonalization of the Hamiltonian. One of the elegant approach in this context has been the so-called similarity renormalization group (SRG) proposed by Glazek and Wilson  (and by Wegner  independently). The SRG resembles the original Wilsonian renormalization group formalism , since a transformation that explicitly runs the cutoff is developed. In an early paper, Wilson  exploited a similar transformation which was originally introduced by Bloch . There the Hamiltonian matrix is reduced by using a transformation which lowers the cutoff initially imposed on the individual states. Later, Wilson abandoned this formalism in favour of a Lagrangian one. The most important reason was that the Bloch transformation is ill-defined and produces fake divergences. These divergences emerge from denominators which contain a small energy difference between states retained and states removed by the transformation. The SRG  is designed to be free of such small energy denominators and to eliminate interactions which change the energies of the unperturbed states by a large amount. However, there are several issues in this approach: it is hard to incorporate loop expansions within the method, the SRG can not systematically remove interactions which change the number of particles (i.e, when the Hamiltonian is not diagonal in particle number space), and most importantly, the computations are complex and there is no an efficient non-perturbative calculable scheme.

In this paper we introduce a method for obtaining the low energy effective operators in the framework of a CCM approach. The transformation constructed avoids the small denominators that plague old-fashioned perturbation theory. Neither perturbation theory nor unitarity of the transformation are essential for this method. The method is non-perturbative, since there is no expansion in the coupling constant; nonetheless, the CCM can be conceived as a topological expansion in the number of correlated excitations. We show that introducing a double similarity transformation using linked-cluster amplitudes will simplify the partial diagonalization underlying renormalization in Hamiltonian approaches. However, a price must be paid: due to the truncation the similarity transformations are not unitary, and accordingly the hermiticity of the resultant effective Hamiltonian is not manifest. This is related to the fact that we have a biorthogonal representation of the given many-body problem. There is a long tradition of such approaches. The first we are aware of are Dyson-type bosonization schemes . \[Here one chooses to map the generators of a Lie algebra, such that the raising generators have a particularly simple representation.\] The space of states is mapped onto a larger space where the physically realizable states are obtained by constrained dynamics. This is closely related to CCM formalism, where the extended phase space is a complex manifold, the physical subspace constraint function has been shown to be of second class and the physical shell itself was found to be a Kähler manifold . The second is the Suzuki-Lee method in the nuclear many-body (NMT) problem , which reduces the full many-body problem to a problem in a small configuration space and introduces a related effective interaction. The effective interaction is naturally understood as the result of certain transformations which decouple a model space from its complement. As is well know in the theory of effective interactions, unitarity of the transformation used for decoupling or diagonalization is not necessary. Actually, the advantage of a non-unitarity approach is that it can give a very simple description for both diagonalization and the ground state. This has been discussed by many authors  and, although it might lead to a non-hermitian effective Hamiltonian, it has been shown that hermiticity can be recovered . Nevertheless non-hermiticity is negligible if the model space and its complement are not strongly correlated . Therefore defining a good model space can in principle control the accuracy of CCM.

To solve the relativistic bound state problem one needs to systematically and simultaneously decouple 1) the high-energy from low-energy modes and 2) the many- from the few-particle states. We emphasize in this paper that CCM can in principle be an adequate method to attack both these requirements. Our hope is to fully utilize Wilsonian Exact renormalization group  within the CCM formalism. Here the high energy modes will be integrated out leading to a modified low-energy Hamiltonian in an effective many-body space. Notice that our formulation does not depend on the form of dynamics and can be used for any quantization scheme, e.g., equal time or light-cone.

The organization of this paper is as follows. In section II, we discuss our approach and it’s foundation. Finally we conclude and present an outlook in section III.

## 2 Formalism

The discussion in this section is partially based on the work of Suzuki and Okamoto . Let us consider a system described by a Hamiltonian $`H(\mathrm{\Lambda })`$ which has, at the very beginning, a large cut-off $`\mathrm{\Lambda }`$. We assume that the renormalized Hamiltonian $`H^{\text{eff}}(\mathrm{\Lambda })`$ up to scale $`\mathrm{\Lambda }`$ can be written as the sum of the canonical Hamiltonian and a “counterterm” $`H_C(\mathrm{\Lambda })`$,

$$H^{\text{eff}}(\mathrm{\Lambda })=H(\mathrm{\Lambda })+H_C(\mathrm{\Lambda }).$$
(1)
Our aim is to construct the renormalized Hamiltonian by obtaining this counterterm. Now imagine that we restrict the Hamiltonian to a lower energy scale $`(\mu )`$, where we want to find an effective Hamiltonian $`H^{\text{eff}}(\mu )`$ which has the same energy spectrum as the original Hamiltonian in the smaller space. Formally, we wish to transform the Hamiltonian to a new basis, where the medium-energy modes $`\mu <k<\mathrm{\Lambda }`$, decouple from the low-energy ones, while the low-energy spectrum remains unchanged. We define two subspaces, the intermediate-energy space $`Q`$ containing modes with $`\mu <k<\mathrm{\Lambda }`$ and a low-energy space $`P`$ with $`\mu k`$. Our renormalization approach is based on decoupling of the complement space $`Q`$ from the model space $`P`$. Thereby the decoupling transformation generates a new effective interaction $`\delta H(\mu ,\mathrm{\Lambda })`$ containing the physics between the scale $`\mathrm{\Lambda }`$ and $`\mu `$. One can then determine the counterterm by requiring coupling coherence , namely that the transformed Hamiltonian has the same form as the original one but with $`\mathrm{\Lambda }`$ every where replaced by $`\mu `$. \[This is in contrast to the popular Effective Field Theory approach, where one includes all permissible couplings of a given order and fixes them by requiring observable computed be both cutoff-independent and Lorentz covariant.\] The operators $`P`$ and $`Q`$ which project a state onto the model space and its complement, satisfy $`P^2=P`$, $`Q^2=Q`$, $`PQ=0`$ and $`P+Q=1`$. We introduce an isometry operator $`G`$ which maps states in the $`P`$\- onto the $`Q`$\- space,

$$|q=G|p(|qQ,|pP).$$
(2)
The operator $`G`$ is the basic ingredient in a family of “integrating-out operators”, which passes information about the correlations of the high energy modes to the low-energy space. The operator $`G`$ obeys $`G=QGP`$, $`GQ=0`$, $`PG=0`$ and $`G^n=0`$ for $`n2`$. The rather surprising direction for $`G`$ to act in is due to the definition Eq. (5) below (cf. the relation between the active and passive view of rotations). In order to give a general form of the effective low-energy Hamiltonian, we define another operator $`X(n,\mu ,\mathrm{\Lambda })`$,

$$X(n,\mu ,\mathrm{\Lambda })=(1+G)(1+G^{}G+GG^{})^n.$$
(3)
($`n`$ is a real number.) The inverse of $`X(n,\mu ,\mathrm{\Lambda })`$ can be obtained explicitly,

$$X^1(n,\mu ,\mathrm{\Lambda })=(1+G^{}G+GG^{})^n(1G).$$
(4)
The special case $`n=0`$ is equivalent to the transformation introduced in Ref.  to relate the hermitian and non-hermitian effective operators in the energy-independent Suzuki-Lee approach. We now consider the transformation of $`H(\mathrm{\Lambda })`$ defined as

$$\overline{H}(n,\mu ,\mathrm{\Lambda })=X^1(n,\mu ,\mathrm{\Lambda })H(\mathrm{\Lambda })X(n,\mu ,\mathrm{\Lambda }),$$
(5)
where we have

$$H(\mathrm{\Lambda })\overline{H}(n,\mu ,\mathrm{\Lambda })H(\mu )+\delta H(\mu ,\mathrm{\Lambda }).$$
(6)
One can prove that if $`\overline{H}(n,\mu ,\mathrm{\Lambda })`$ satisfies the desirable decoupling property,

$$Q\overline{H}(n,\mu ,\mathrm{\Lambda })P=0,$$
(7)
or more explicitly, by substituting the definition of $`X(n,\mu ,\mathrm{\Lambda })`$ and $`X^1(n,\mu ,\mathrm{\Lambda })`$ from Eqs. (3-4),

$$QH(\mathrm{\Lambda })P+QH(\mathrm{\Lambda })QGGPH(\mathrm{\Lambda })PGPH(\mathrm{\Lambda })QG=0,$$
(8)
that $`H^{\text{eff}}(\mu )(n,\mu )=P\overline{H}(n,\mu ,\mathrm{\Lambda })P`$ is an effective Hamiltonian for the low energy degrees of freedom. In other words, it should have the same low-energy eigenvalues as the original Hamiltonian. The proof is as follows:

Consider an eigenvalue equation in the $`P`$ space with $`\{|\varphi (k)P\}`$,

$`P\overline{H}(n,\mu ,\mathrm{\Lambda })P|\varphi (k)=E_kPX^1(n,\mu ,\mathrm{\Lambda })X(n,\mu ,\mathrm{\Lambda })P|\varphi (k).`$ (9)
By multiplying both sides by $`X(n,\mu ,\mathrm{\Lambda })`$ and making use of the decoupling property Eq. (7), we obtain

$$H(\mathrm{\Lambda })X(n,\mu ,\mathrm{\Lambda })P|\varphi (k)=E_kX(n,\mu ,\mathrm{\Lambda })P|\varphi (k).$$
(10)
This equation means that $`E_k`$ in Eq. (9) agrees with one of the eigenvalue of $`H(\mathrm{\Lambda })`$ and $`X(n,\mu ,\mathrm{\Lambda })P|\varphi (k)`$ is the corresponding eigenstate. Now we demand that

$$H^{\text{eff}}(\mu )H(\mu )+H_C(\mu ).$$
(11)
This requirement uniquely determines the counterterm $`H_C`$.

In the same way we can also obtain the $`Q`$-space effective Hamiltonian, from the definition of $`\overline{H}(n,\mu ,\mathrm{\Lambda })`$. It can be seen that if $`G`$ satisfies the requirement in Eq. (8), then we have additional decoupling condition

$$P\overline{H}(n,\mu ,\mathrm{\Lambda })Q=0.$$
(12)
We will argue later that *both* of the decoupling conditions Eqs. (7) and (12) are necessary in order to have a sector-independent renormalization scheme. The word “sector” here means the given truncated Fock space. Let us now clarify the meaning of this concept. To maintain the generality of the previous discussion, we use here the well known Bloch-Feshbach formalism . The Bloch-Feshbach method exploits projection operators in the Hilbert space in order to determine effective operators in some restricted model space. This technique seems to be more universal than Wilson’s renormalization formulated in a Lagrangian framework. This is due to the fact that in the Bloch-Feshbach formalism, other irrelevant degrees of freedom (such as high angular momentum, spin degrees of freedom, number of particles, etc.) can be systematically eliminated in the same fashion.

Assume that the full space Schrödinger equation is $`H(\mathrm{\Lambda })|\psi =E|\psi `$ and for simplicity $`|\psi `$ has been normalized to one. We explicitly construct the energy dependent effective Hamiltonian in this formalism,

$`H^{\text{eff}}`$ $`=`$ $`P\overline{H}P+P\overline{H}Q{\displaystyle \frac{1}{EQ\overline{H}Q}}Q\overline{H}P,`$ (13)
where $`\overline{H}`$ can be a similarity transformed Hamiltonian. This equation resembles the one for Brueckner’s reaction matrix (or “G”-matrix) equation in nuclear many-body theory (NMT). In the same way for arbitrary operator $`O`$ (after a potential similarity transformation), we construct the effective operator

$`O^{\text{eff}}`$ $`=`$ $`P\overline{O}P+P\overline{H}Q{\displaystyle \frac{1}{EQ\overline{H}Q}}Q\overline{O}P+P\overline{O}Q{\displaystyle \frac{1}{EQ\overline{H}Q}}Q\overline{H}P`$ (14)
$`+`$ $`P\overline{H}Q{\displaystyle \frac{1}{EQ\overline{H}Q}}Q\overline{O}Q{\displaystyle \frac{1}{EQ\overline{H}Q}}Q\overline{H}P.`$
The $`E`$-dependence in Eqs. (13) and (14) emerges from the fact that the effective interaction in the reduced space is not assumed to be decoupled from the excluded space. However, by using the decoupling conditions introduced in Eqs. (7) and (12), we observe that energy dependence can be removed, and the effective operators become

$`H^{\text{eff}}`$ $`=`$ $`P\overline{H}P=(n,\mu ),`$
$`O^{\text{eff}}`$ $`=`$ $`P\overline{O}P=𝒪(n,\mu ).`$ (15)
The decoupling property makes the operators in one sector independent of the other sector. The effects of the excluded sector is taken into account by imposing the decoupling conditions. This is closely related to the folded diagram method in NMT for removing energy-dependence . (It is well-known in NMT that $`E`$-dependence in the $`G`$-matrix emerges from non-folded diagrams which can be systematically eliminated using the effective interaction approach). The above argument was given without assuming an explicit form for $`X`$ and thus the decoupling conditions are more fundamental than the prescription used to derive these conditions. We will show later that one can choose a transformation $`X`$, together with the model space and its complement, which avoids the occurrence of “small energy denominators”. We now show that Lorentz covariance in a given sector does not hinge on a special form of similarity transformation. We assume ten Poincaré generators $`L_i`$ satisfying

$$[L_i,L_j]=a_{ij}^kL_k,$$
(16)
where the $`a_{ij}^k`$ are the known structure coefficients. One can show that if the operators $`L_i`$ satisfy the decoupling conditions $`Q\overline{L}_iP=0`$ and $`P\overline{L}_iQ=0`$ then it follows that

$$[L_i^{\text{eff}},L_j^{\text{eff}}]=a_{ij}^kL_k^{\text{eff}}.$$
(17)
This leads to a relativistic description even after simultaneously integrating out the high-frequency modes and reducing the number of particles. Indeed we conjecture that requiring the decoupling conditions makes the effective Hamiltonian free of Lorentz-noninvariant operators for a given truncated sector regardless of the regularization scheme. However, one may still be faced with an effective Hamiltonian which violates gauge invariance (for e.g., when sharp cutoff is employed).

Note that the solution to Eq. (8) is independent of the number $`n`$. One can make use of Eq. (8) and its complex conjugate to show that for any real number $`n`$, the following relation for the effective low-energy Hamiltonian

$$(n,\mu )=^{}(n1,\mu ).$$
(18)
The case $`n=1/2`$ is special since the effective Hamiltonian is hermitian,

$$(1/2,\mu )=(P+G^{}G)^{1/2}H(\mathrm{\Lambda })(P+G)(P+G^{}G)^{1/2}.$$
(19)
Hermiticity can be verified from the relation

$$e^TP=(1+G)(P+G^{}G)^{1/2},$$
(20)
where,

$$T=\mathrm{arctan}(GG^{})=\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{2n+1}(G(G^{}G)^n\text{h.c}.).$$
(21)
Since the operator $`T`$ is anti-hermitian, $`e^T`$ is a unitary operator. From the above expression Eq. (19) can be written in the explicitly hermitian form

$$(1/2,\mu )=Pe^TH(\mathrm{\Lambda })e^TP.$$
(22)
As was already emphasized, renormalization based on unitary transformations is more complicated and non-economical. Thus we will explore a non-unitary approach. An interesting non-hermitian effective low-energy Hamiltonian can be obtained for $`n=0`$,

$$(0,\mu ,\mathrm{\Lambda })=PH(\mathrm{\Lambda })(P+QG).$$
(23)
This form resembles the Bloch and Horowitz type of effective Hamiltonian as used in NMT , and was used by Wilson in quantum field theory . In the context of the CCM, this form leads to the folded diagram expansion well known in many-body theory . It is of interest that various effective low-energy Hamiltonians can be constructed according to Eq. (15) by the use of the mapping operator $`G`$ which all obey the decoupling property Eq. (7) and Eq. (12). Neither perturbation theory nor hermiticity is essential for this large class of effective Hamiltonians.

The CCM approach is, of course, just one of the ways of describing the relevant spectrum by means of non-unitary transformations. According to our prescription the model space is $`P:\{|L|0,b_h,L\mu \}`$, where $`|0,b_h`$ is a bare high energy vacuum (the ground state of high-momentum of free-Hamiltonian) which is annihilated by all the high frequency annihilation operators $`\{C_I\}`$ (for a given quantization scheme, e.g., equal time or light-cone) , the set of indices $`\{I\}`$ therefore defines a subsystem, or cluster, within the full system of a given configuration. The actual choice depends upon the particular system under consideration. In general, the operators $`\{C_I\}`$ will be products (or sums of products) of suitable single-particle operators. We assume that the annihilation and its corresponding creation $`\{C_I^{}\}`$ subalgebras and the state $`|0,b_h`$ are cyclic, so that the linear combination of state $`\{C_I^{}|0,b_h\}`$ and $`\{{}_{h}{}^{}b,0|C_I\}`$ span the decoupled Hilbert space of the high-momentum modes, $`\{|H\}`$, where $`\mu <H<\mathrm{\Lambda }`$. It is also convenient, but not necessary, to impose the orthogonality condition, $`0|C_IC_J^{}|0=\delta (I,J)`$, where $`\delta (I,J)`$ is a Kronecker delta symbol. The complement space is $`Q:\{|L\left(|H|0,b_h\right)\}`$. Our main goal is to decouple the $`P`$-space from the $`Q`$-space. This gives sense to the partial diagonalization of the high-energy part of the Hamiltonian. The states in full Hilbert space are constructed by adding correlated clusters of high-energy modes onto the $`P`$-space, or equivalently integrating out the high-energy modes from the Hamiltonian,

$`|f=X(\mu ,\mathrm{\Lambda })|p=e^{\widehat{S}}e^{\widehat{S}^{}}|0,b_h{\displaystyle |L}=e^{\widehat{S}}|0,b_h{\displaystyle |L},`$ (24)
$`\stackrel{~}{f}|=L|{\displaystyle }{}_{h}{}^{}b,0|X^1(\mu ,\mathrm{\Lambda })=L|{\displaystyle }{}_{h}{}^{}0|e^{\widehat{S}^{}}e^{\widehat{S}},`$ (25)
where the operators $`X(\mu ,\mathrm{\Lambda })`$ and $`X^1(\mu ,\mathrm{\Lambda })`$ have been expanded in terms of independent coupled cluster excitations $`I`$,

$`\widehat{S}={\displaystyle \underset{m=0}{}}\widehat{S}_m\left({\displaystyle \frac{\mu }{\mathrm{\Lambda }}}\right)^m,\widehat{S}_m=\underset{I}{{\displaystyle ^{}}}\widehat{s}_IC_I^{},`$
$`\widehat{S}^{}={\displaystyle \underset{m=0}{}}\widehat{S}_m^{}\left({\displaystyle \frac{\mu }{\mathrm{\Lambda }}}\right)^m,\widehat{S}_m^{}=\underset{I}{{\displaystyle ^{}}}\widehat{s}_I^{}C_I.`$ (26)
Here the primed sum means that at least one fast particle is created or destroyed $`(I0)`$, and momentum conservation in $`PQ`$ is included in $`\widehat{s}_I`$ and $`\widehat{s}_I^{}`$. $`\widehat{S}_m(\widehat{S}_m^{})`$ are not generally commutable in the low-energy Fock space, whereas they are by construction commutable in the high-energy Fock space. It is immediately clear that states in the interacting Hilbert space are normalized, $`\stackrel{~}{f}|f={}_{h}{}^{}b,0|0,b_{h}^{}=1`$. We have two types of parameters in this procedure: One is the coupling constant of the theory ($`\lambda `$), and the other is the ratio of cutoffs ($`\mu /\mathrm{\Lambda }`$). The explicit power counting makes the degree of divergence of each order smaller than the previous one. According to our logic, Eq. (23) can be written as

$$(\mu )={}_{h}{}^{}b,0|X^1(\mu ,\mathrm{\Lambda })H(\mathrm{\Lambda })X(\mu ,\mathrm{\Lambda })|0,b_{h}^{},$$
(27)
with $`X(\mu ,\mathrm{\Lambda })`$ and $`X^1(\mu ,\mathrm{\Lambda })`$ defined in Eq. (24) and Eq. (25). We require that effective Hamiltonian $`(\mu )`$ obtained in this way remains form invariant or coherent . This requirement satisfies on an infinitely long renormalization group trajectory and thus does constitutes a renormalized Hamiltonian. Thereby one can readily identify the counter terms produced from expansion of Eq. (27).

It is well-know in many-body applications that the exponential Ansatz Eqs. (24) and (25) guarantees automatically proper size-exclusivity and conformity with the Goldstone linked-cluster theorem to all level of truncation. This parametrization does not manifestly preserve hermitian conjugacy. However, it is compatible with the Hellmann-Feynman theorem (HFT), in other words, demanding hermiticity will violate this theorem at any level of truncation. On the other hand, with this parametrization the phase space $`\{\widehat{s}_I,\widehat{s}_I^{}\}`$ for a given $`m`$ is a symplectic differentiable manifold. Thereby all the geometrical properties of the configuration space can be precisely defined. (In a more ordinary language, the canonical equations of motion with respect to phase space define a set of trajectories, which fill the whole dynamically allowed region of the phase space.) There is a deep connection between these three properties and we can not give up one without losing something else as well . The individual amplitudes for a given $`m`$, $`\{\widehat{s}_I^m,\widehat{s}^{}{}_{I}{}^{m}\}\{\widehat{s}_I,\widehat{s}_I^{}\}_m`$, can generally be functionals of the low- and high-energy field operators and have to be fixed by the dynamics of quantum system. This is a complicated set of requirements. However, we require less than that. Suppose that after a similar transformation of Hamiltonian, $`\overline{H}`$, we obtain an effective Hamiltonian of the form

$$\overline{H}=H(\text{low})+H_{\text{free}}(\text{high})+C_I^{}V_{IJ}C_J,$$
(28)
where $`V_{IJ}`$ is an arbitrary operator in the low frequency space. The $`I`$ and $`J`$ indices should be chosen such that the last term in Eq. (28) contains at least one creation- operator and one annihilation-operator of high frequency. By using Rayleigh-Schrödinger perturbation theory, it can be shown that the free high-energy vacuum state of $`H_{\text{free}}(\text{high})`$ is annihilated by Eq. (28) and remains without correction at any order of perturbation theory. Having said that, we will now consider how to find the individual amplitudes $`\{\widehat{s}_I,\widehat{s}_I^{}\}_m`$ that transfer the Hamiltonian into the form Eq. (28). We split the Hamiltonian in five parts:

$$H=H_1+H_2^{\text{free}}(\text{high})+V_C(C_I^{})+V_A(C_I)+V_B,$$
(29)
where $`H_1`$ contains only the low frequency modes with $`k\mu `$, $`H_2`$ is the free Hamiltonian for all modes with $`\mu <k<\mathrm{\Lambda }`$, $`V_C`$ contains low frequency operators and products of the high frequency creation operators $`C_I^{}`$ and $`V_A`$ is the hermitian conjugate of $`V_C`$. The remaining terms are contained in $`V_B`$, these terms contain at least one annihilation and creation operators of the high energy modes. Our goal is to eliminate $`V_C`$ and $`V_A`$ since $`V_B`$ annihilates the vacuum. The ket-state coefficients $`\{\widehat{s}_I\}_m`$ are worked out via the ket-state Schrödinger equation $`H(\mathrm{\Lambda })|f=E|f`$ written in the form

$$0|C_Ie^{\widehat{S}}He^{\widehat{S}}|0=0,I0.$$
(30)
The bra-state coefficients $`\{\widehat{s}_I,\widehat{s}_I^{}\}_m`$ are obtained by making use of the Schrödinger equation defined for the bra-state, $`\stackrel{~}{f}|H(\mathrm{\Lambda })=\stackrel{~}{f}|E`$. First we project both sides on $`C_I^{}|0`$, then we eliminate $`E`$ by making use of the ket-state equation projection with the state $`0|e^{\widehat{S}^{}}C_I^{}`$ to yield the equations

$$0|e^{\widehat{S}^{}}e^{\widehat{S}}[H,C_I^{}]e^{\widehat{S}}e^{\widehat{S}^{}}|0=0,I0.$$
(31)
Alternatively one can in a unified way apply $`e^{\widehat{S}}e^{\widehat{S}^{}}C_I^{}|0`$ on the Schrödinger equation for the bra-state and obtain

$$0|e^{\widehat{S}^{}}e^{\widehat{S}}He^{\widehat{S}}e^{\widehat{S}^{}}C_I^{}|0=0,I0.$$
(32)
Equation (30) and Eqs. (31) or (32) provide two sets of formally exact, microscopic, operatorial coupled non-linear equations for the ket and bra. One can solve the coupled equations in Eq. (30) to work out $`\{\widehat{s}_I\}_m`$ and then use them as an input in Eqs. (31) or (32).

It is important to notice that Eqs. (30) and (31) can be also derived by requiring that the effective low-energy Hamiltonian defined in Eq. (27), be stationary (i.e. $`\delta (\mu )=0`$) with respect to all variations in each of the independent functional $`\{\widehat{s}_I,\widehat{s}_I^{}\}_m`$. One can easily verify that the requirements $`\delta (\mu )/\delta \widehat{s}_I=0`$ and $`\delta (\mu )/\delta \widehat{s}_I^{}=0`$ yield Eqs. (30) and (31). The combination of Eqs. (30) and (31) does not manifestly satisfy the decoupling property as set out in Eqs. (7) and (12). On the other hand Eqs. (30) and (32) satisfy these conditions. Equations (30) and (32) imply that all interactions including creation and annihilation of fast particles (“$`I`$”) are eliminated from the transformed Hamiltonian $`(\mu )`$ in Eq. (27). In other words, these are decoupling conditions leading to the elimination of $`V_C`$ and $`V_A`$ from Eq. (29), which is, in essence, a block-diagonalization. Therefore it makes sense for our purpose to use Eqs. (30) and (32) for obtaining the unknown coefficients, losing some of the elegance of the CCM elsewhere.

So far everything has been introduced rigorously without invoking any approximation. In practice one needs to truncate both sets of coefficients $`\{\widehat{s}_I,\widehat{s}_I^{}\}_m`$ at a given order of $`m`$. A consistent truncation scheme is the so-called SUB($`n,m`$) scheme, where the $`n`$-body partition of the operator $`\{\widehat{S},\widehat{S}^{}\}`$ is truncated so that one sets the higher partition with $`I>n`$ to zero at a given accuracy $`m`$. Notice that, Eqs. (31) and (32) provide two equivalent sets of equations in the exact form, however after the truncation they can in principle be different. Eqs. (30) and (32) are compatible with the decoupling property at any level of the truncation, whereas Eqs. (30) combined with (31) are fully consistent with HFT at any level of truncation. Thus the low-energy effective form of an arbitrary operator can be computed according to Eq. (15) in the same truncation scheme used for the renormalization of the Hamiltonian. In particular, we will show that only in the lowest order ($`m=0`$), equations (31) and (32) are equivalent, independent of the physical system and the truncation scheme.

Although our method is non-perturbative, perturbation theory can be recovered from it. In this way, its simple structure for loop expansion will be obvious and we will observe that at lower order hermiticity is preserved. Now we illustrate how this is realizable in our approach. Assume that $`V_C`$ and $`V_A`$ are of order $`\lambda `$, we will diagonalize the Hamiltonian, at leading order in $`\lambda `$ up to the desired accuracy in $`\mu /\mathrm{\Lambda }`$. We use the commutator-expansion

$$e^SHe^S=H+[H,S]+\frac{1}{2!}[[H,S],S]+\mathrm{}.$$
(33)
Eq. (30) can be organized perturbatively in order of $`m`$, aiming at elimination of the high momenta degree of freedom up to the first order in the coupling constant, thus yields

$`m=0:0|C_I(V_C+[H_2,\widehat{S}_0])|0=0,`$
$`m=1:0|C_I([H_1,\widehat{S}_0]+[H_2,\widehat{S}_1]+[V_A,\widehat{S}_1]+[V_C,\widehat{S}_1])|0=0,`$
$`\mathrm{}`$
$`m=n:0|C_I([H_1,\widehat{S}_{n1}]+[H_2,\widehat{S}_n]+[V_A,\widehat{S}_n]+[V_C,\widehat{S}_n])|0=0,`$ (34)
where $`I0`$. Notice that $`\widehat{S}_0`$ is chosen to cancel $`V_C`$ in the effective Hamiltonian, hence it is at least of order of $`\lambda `$, consequently it generates a new term $`[H_1,\widehat{S}_0]`$ which is of higher order in $`\mu /\mathrm{\Lambda }`$ and can be canceled out on the next orders by $`\widehat{S}_1`$. The logic for obtaining the equations above is based on the fact that $`\widehat{S}_n`$ should be smaller than $`\widehat{S}_{n1}`$ (for sake of convergence) and that the equations should be consistent with each other. Since $`H_2,V_A,V_C\mathrm{\Lambda }`$ and $`H_1\mu `$, from Eq. (34) we have the desired relation $`\widehat{S}_n\frac{\mu }{\mathrm{\Lambda }}\widehat{S}_{n1}`$. The same procedure can be applied for Eq. (32) which leads to the introduction of a new series of equations in order of $`m`$,

$`m=0:0|(V_A[H_2,\widehat{S}_0^{}])C_I^{}|0=0,`$
$`m=1:0|([H_1,\widehat{S}_0^{}]+[H_2,\widehat{S}_1^{}]+[V_C,\widehat{S}_1^{}]+[V_A,\widehat{S}_1^{}][V_A,\widehat{S}_1])C_I^{}|0=0,`$
$`\mathrm{}`$
$`m=n:0|([H_1,\widehat{S}_{n1}^{}]+[H_2,\widehat{S}_n^{}]+[V_C,\widehat{S}_n^{}]+[V_A,\widehat{S}_n^{}][V_A,\widehat{S}_n])C_I^{}|0=0.`$ (35)
Alternatively, we can use Eq. (31) to yield the equations

$`m=0:0|\left([V_A,C_I^{}][[H_2,C_I^{}],\widehat{S}_0^{}]\right)|0=0,`$
$`m=1:0|\left([[V_A,C_I^{}],\widehat{S}_1][[H_2,C_I^{}],\widehat{S}_1^{}][[V_A,C_I^{}],\widehat{S}_1^{}]\right)|0=0,`$
$`\mathrm{}`$
$`m=n:0|([V_A,C_I^{}],\widehat{S}_n][[H_2,C_I^{}],\widehat{S}_n^{}][[V_A,C_I^{}],\widehat{S}_n^{}])|0=0.`$ (36)
It is obvious that at order $`m=0`$, Eqs. (35) and (36) are the same and $`\widehat{S}_0^{}=\widehat{S}_0^{}`$, which indicates that the similarity transformation at this level remains unitary. It should be noted that diagonalization at first order in the coupling constant introduces a low-energy effective Hamiltonian in Eq. (27) which is valid up to the order $`\lambda ^3`$. In the same way, diagonalization at second order in $`\lambda `$ modifies the Hamiltonian at order $`\lambda ^4`$ and leads generally to a non-unitarity transformation. In this way one can proceed to diagonalize the Hamiltonian at a given order in $`\lambda `$ with desired accuracy in $`\mu /\mathrm{\Lambda }`$ . Finally, the renormalization process is completed by introducing the correct $`Z(\mathrm{\Lambda })`$ factors which redefine the divergences emerging from Eq. (27).

## 3 Conclusion and outlook

In this paper we have outlined a strategy to derive effective renormalized operators in the Hamiltonian formulation of a field theory. An example of the application of the method can be found in Ref. . The effective low-frequency operator is obtained by the condition that it should exhibit decoupling between the low- and high-frequency degrees of freedom. We showed that the similarity transformation approach to renormalization can be systematically classified. The non-hermitian formulation gives a very simple description of decoupling, leading to a partial diagonalization of the high-energy part. In the recent paper , we showed that non-unitarity representation inherent in our formulation is in favour of economic computation and does not produce any non-hermiticity in the relevant terms. The techniques proposed are known from the coupled cluster many-body theory and invoke neither perturbation nor unitarity transformation. We showed that our formalism can be solved perturbatively. In this way, it was revealed that diagonalization at first order in coupling constant defines a correct low-energy effective Hamiltonian which is valid up to the order $`\lambda ^3`$. One can show that the non-hermiticity of the effective Hamiltonian is controllable and might appear in higher order which is beyond our approximation or in irrelevant terms which can be ignored in renormalization group sense.

One of the key features which has not yet been exploited is the non-perturbative aspect of the method; it may well be able to obtain effective degrees of freedom that are very different from the ones occurs at the high-energy scale. This is a promising avenue for future work.

## 4 Acknowledgment

One of the authors (AHR) acknowledges support from British Government ORS award and UMIST grant. The work of NRW is supported by the UK engineering and physical sciences research council under grant GR/N15672.