text:
# Building a CCD Spectrograph for Educational or Amateur Astronomy

## 1 Introduction

The widespread availability of affordable, high-quality commercial CCD cameras for astronomy now makes it possible for educators and amateur astronomers to conduct sophisticated observing programs with small telescopes and limited resources. In recent years the potential for building inexpensive CCD spectrographs around these cameras has been recognized by several companies and individuals (e.g. Glumac & Sivo, 1999). Observations of the spectra of stars, planets, and nebulae open up an exciting new realm of inquiry for beginning astronomers, and the design and construction of the spectrograph itself can provide students with a hands-on education in the principles of optics and engineering. Here we describe one approach to designing and building a simple CCD spectrograph, based on IntroSpec, an educational instrument we recently built for the Harvard Astronomy Department.

We had four design goals for IntroSpec: (1) sufficient spectral range and resolution for a variety of projects in optical astronomy, (2) high throughput to enable observations of nearby nebulae and variable stars, (3) simple controls and a durable design suitable for student operation, and (4) a small telescope-borne weight. We decided upon a spectral range of $``$4500–6800 Å at $``$6 Å resolution. This wavelength range includes many important spectral features: H, Fe, Mg and Na absorption lines (prominent in stellar atmospheres), molecular absorption bands from CH4 and NH3 (visible in planetary atmospheres), as well as H, N, and O emission lines (characteristic of star-forming nebulae). Spectral coverage extending $``$1000 Å futher to the blue would be useful, but most inexpensive CCD chips have poor UV response. Our chosen resolution of 6 Å resolution allows detection of most of the commonly observed features in this spectral range and is adequate to resolve the closely spaced emission lines of H$`\alpha `$ and \[NII\].

We were able to meet our design goals on a budget of slightly under $5500, of which $``$$4000 went into buying a mid-range Meade CCD camera and paying for the time and materials of a skilled machinist (CBH).111These numbers do not include an expensive but desirable accessory: a Hg-Ne lamp for wavelength calibration. We purchased one from Oriel for an additional $450. Below we describe how IntroSpec’s design was tailored to the specific challenges posed by our budget and performance requirements. §2 reviews the general optical design constraints on any similar spectrograph, while §3 examines the specific choices we made for IntroSpec, with special attention to the instrument’s fiber-fed design.

## 2 Optical Design

We begin with a brief review of the layout of a generic spectrograph, then consider how to optimize its optical design.

### 2.1 Layout of a Generic Spectrograph

Figure 1 shows the light path in a simple CCD spectrograph. The aperture accepts light from a small portion of the telescope focal plane that may contain the image of a star or a portion of an extended nebula. This light is collimated with a lens, producing a parallel beam. The beam travels to a diffraction grating, which disperses the light into a spectrum described by the grating equation (e.g. Schroeder, 2000):

$$m\lambda =d(\mathrm{sin}\theta _1+\mathrm{sin}\theta _2),$$
(1)
where $`m`$ is the diffraction order, $`\lambda `$ is the wavelength of the light, $`d`$ is the grating line spacing, and $`\theta _1`$ and $`\theta _2`$ are the incoming and outgoing angles of the light with respect to the normal to the grating surface (see Schroeder for sign conventions).

A reflection grating such as that shown in Figure 1, folds the light path as well as dispersing the light. In this case the angle between the incident and diffracted beam typically must be $`>`$30$`\mathrm{°}`$ to allow the camera optics to stay clear of the collimated beam heading toward the grating. The grating normal generally points closer to the camera axis than to the collimator axis (Schroeder, 2000), leading to anamorphic magnification: a geometrical expansion of the dispersed collimated beam along the dispersion axis described by $`AM=\mathrm{cos}\theta _2/\mathrm{cos}\theta _1`$ (see Figure 1). Anamorphic magnification increases the spectral resolution by the same factor that the collimated beam is expanded, but requires a larger camera acceptance aperture to avoid vignetting. Camera-collimator angles exceeding 45$`\mathrm{°}`$ are rarely used for this reason.

The light dispersed by the grating is a superposition of collimated beams of different colors, all traveling at slightly different angles. The camera lens focuses each color onto a location on the CCD chip determined by its incoming angle, thereby forming a spectrum. This focusing process is exactly analogous to the way a telescope focuses incoming parallel light beams from stars at different angles on the sky onto different positions on the telescope focal plane: in the spectrograph, each point along the spectrum is an image of the light from the spectrograph entrance aperture at a different wavelength. A higher dispersion grating increases the angular spread between collimated beams of different colors, spreading the corresponding images further apart. For a detector of fixed size, higher dispersion yields higher spectral resolution at the expense of reduced spectral coverage.

### 2.2 Setting the Relative Scales of the System

Independent of the grating choice or the absolute scale of the spectrograph, we wish to specify the relative focal lengths of the spectrograph optics so that the image of the spectrograph entrance aperture is correctly sampled by the the CCD pixels. The width of this image in the dispersion direction defines the spectroscopic resolution element. Oversampling the image wastes CCD pixels and decreases spectral coverage, while undersampling it degrades spectral resolution. In the Nyquist sampling limit, $``$2 pixels per wavelength allow one to recover the peak and valley of an oscillating waveform. Typically, finer sampling up to 3–4 pixels per resolution element is advantageous in professional instruments with large CCD chips. However, for an economical instrument, CCD pixels for spectral coverage are at a premium, and one may prefer sampling near the Nyquist limit. With Nyquist sampling, IntroSpec’s design target of 6 Å resolution combined with 2300 Å spectral coverage requires a CCD chip with $`2\times 2300/6=767`$ pixels. We chose this combination of spectral resolution and spectral coverage knowing that the largest CCD within our budget had 768 pixels in the long dimension.

Once the desired sampling is chosen, the ratio of the focal lengths of the collimator and camera lenses, $`f_{\mathrm{coll}}`$ and $`f_{\mathrm{cam}}`$, is determined by the diameter ratio of the spectrograph entrance aperture and its image. For Nyquist sampling, the image of the spectrograph entrance aperture must be reduced to $``$2 CCD pixels, or $``$20 $`\mu `$m for a typical amateur CCD camera. The spectrograph entrance aperture has generally been sized to the expected diameter of a stellar image at the telescope’s focus, which is determined by atmospheric conditions and the telescope optics. This choice maximizes the constrast between the the object and the background sky light. For the 40 cm $`f`$/10 Knowles Telescope at Harvard, the local five-arcsecond seeing disk makes a 100 $`\mu `$m image at the telescope’s focus, implying that $`f_{\mathrm{coll}}/f_{\mathrm{cam}}`$ must equal $``$5 to reduce the image to $``$20 $`\mu `$m. For a 20 cm $`f`$/10 telescope, similar seeing would make a 50 $`\mu `$m point source image, requiring an $`f_{\mathrm{coll}}/f_{\mathrm{cam}}`$ ratio of $``$2.5. However, a larger entrance aperture (implying higher $`f_{\mathrm{coll}}/f_{\mathrm{cam}}`$ to achieve the same resolution) might be desirable for a 20 cm telescope to allow for image shifts due to imperfect guiding.

The acceptance angles of the collimator and camera lenses can also be estimated approximately independently of the absolute scale of the spectrograph. The collimator must be fast enough to accept an $`f`$/10 light cone from the telescope from every point within the spectrograph entrance aperture. This aperture typically takes the form of a long slit to enable spatial sampling perpendicular to the dispersion direction. The camera must be faster than the collimator: if $`f_{\mathrm{coll}}/f_{\mathrm{cam}}5`$, each monochromatic beam intercepted by the camera has a focal ratio of $``$$`f/2`$ and is elongated in the dispersion direction by anamorphic magnification (see §2.1). These beams must be accepted over a field of view that is extended in both the spatial direction (due to the slit length) and the dispersion direction (due to the angular divergence of the different wavelengths in the spectrum, §3.3).

### 2.3 Setting the Absolute Scale of the System

We have some flexibility in choosing the overall dimensions of the optics, limited by financial constraints and the desired spectroscopic performance. We must balance two considerations. First, larger optics and their mounts are generally more expensive, although the economies of mass-produced commercial optics can bend this rule. Second, to maintain a given spectral range and resolution on a fixed CCD format, one must keep the product of the grating dispersion (ruling density) and the focal length of the camera constant. The grating cost goes up nonlinearly with the ruling density, so the grating cost and availability must be balanced against the cost of the optics, particularly the camera lens.

The widespread availability of commercial 35 mm format single lens reflex (SLR) photographic lenses suggests an attractive solution. These lenses are inexpensive, offer excellent image quality, have sophisticated anti-reflection coatings, and are designed to work over a film format of 24 $`\times `$ 36 mm. CCD chips in typical amateur cameras fit easily within this format. Many photographic lenses can produce images of $``$20 $`\mu `$m diameter, or about 2 pixels for a CCD camera with $``$10 $`\mu `$m pixels. Excellent used lenses can be obtained for less than $100, particularly if an older manual focus lens is selected.

To minimize vignetting, our design for IntroSpec calls for a camera lens faster than $`f`$/2 (§2.2), so an $`f`$/1.4 lens is a good choice. Commercial 35 mm SLR lenses as fast as $`f`$/1.4 are typically available with focal lengths between 35 mm and 85 mm, with 50 mm being by far the most common. These focal lengths pair conveniently with readily available diffraction gratings to meet our goals for spectral coverage and resolution. Table 1 describes three combinations of gratings and lenses with focal lengths between 35 and 85 mm. We assume a camera-collimator angle of 36, a central wavelength of 5650 Å, and a 6.9 mm CCD chip with 768 pixels (as found in IntroSpec). All of these lens-grating combinations approximately satisfy our design targets of 6 Å resolution and 2300 Å spectral coverage. Given the ease of obtaining both 50 mm focal length lenses and 600 line/mm gratings, we selected this configuration for IntroSpec. The spectrograph would have been somewhat smaller had we used a more expensive 35 mm focal length camera lens, but as Table 1 demonstrates, vignetting due to anamorphic magnification would be slightly greater.

## 3 IntroSpec

### 3.1 Optical Fibers Ease Design Constraints

IntroSpec’s design follows that of the generic CCD spectrograph in Figure 1, with one key difference: the spectrograph does not mount directly to the telescope. This change reflects the optical design constraints outlined in §2, which led us to prefer moderately large optical elements (a 50 mm focal length camera lens and a 250 mm focal length collimator), implying an instrument

whose longest dimension exceeds 0.5 meter. A commercial 40 cm telescope can support neither the size nor the weight of such an instrument. Instead of sacrificing optical performance, we opted for a bench-mounted, fiber-fed design. Light from the telescope enters a set of optical fibers at the telescope’s focus, and from there the fibers carry the light to IntroSpec, which sits on a table. Because IntroSpec’s weight is removed from the telescope, it can be packaged in a heavy metal box with an easy-to-remove lid, which protects the instrument while keeping it accessible to students. IntroSpec’s wavelength calibration is stable because the mechanical flexure encountered in a telescope-mounted spectrograph is eliminated.

### 3.2 General Overview of the Instrument

Figure 2 and Table 2 provide an overview of IntroSpec. Optical fibers are coupled to the telescope via a guider/telescope attachment, further discussed in §3.6. The fiber train consists of six fibers, each of which accepts light from a different position on the telescope focal plane; these are bundled in several layers of protective coverings. (§3.4 discusses the technical aspects of working with fibers in more detail.) Upon entering the spectrograph enclosure, the fiber train loops around a large black spool for strain relief, then the fibers are channeled into the fiber output, from which the bare fiber ends emerge in a vertical stack (Figure 3). The base of the fiber output assembly is adjustable to permit focusing the collimator lens. Each optical fiber is an independent spectrograph entrance aperture; by stacking these apertures perpendicular to the dispersion axis, we obtain six parallel spectra on the CCD chip (Figure 4).

IntroSpec uses two fibers of each of three diameters: 100, 200, and 300 $`\mu `$m. The seeing disk at the Knowles telescope is well matched to a 100 $`\mu `$m fiber, and this fiber diameter provides the highest spectral resolution, $``$6 Å FWHM. Larger fibers give the user the choice to trade off spectral resolution in exchange for gathering additional light from an extended source such as a nebula. Larger fibers are also more forgiving of guiding errors when spectral resolution is not critical. The different fiber diameters allow students to learn about the instrument and to make the connection between fiber size and spectral resolution. For point sources, one fiber can be used to acquire the target spectrum and the second fiber of the same diameter can be used to acquire a sky spectrum for background subtraction.

From the collimator on, IntroSpec looks much like the generic spectrograph of Figure 1. The collimator lens is a 50 mm diameter achromat with $`f_{\mathrm{coll}}=250`$ mm, yielding a collimated beam diameter of 25 mm with an f/10 input light cone. We selected a 588 line mm-1 grating blazed222The ruled diffraction gratings typically used in astronomical instruments are blazed to maximize the amount of light in the desired diffraction order. A blazed reflection grating has a triangular groove profile, with the reflective surface of the groove inclined to the grating surface by an angle known as the blaze angle. The blaze angle is chosen such that simple geometrical reflection from the groove face sends light not into the central maximum, but instead into the chosen diffraction order, coincident with the desired peak-efficiency wavelength. In the case of transmission gratings, the blaze angle is chosen such that light geometrically refracted from the groove face is coincident with the desired peak-efficiency wavelength in the selected diffraction order. We chose a reflection grating for IntroSpec because reflection gratings are somewhat less expensive than transmission gratings and the folding of the light path is convenient. for peak efficiency at 5610 Å. Richardson Grating Laboratory offers this grating in sizes up to $``$50 mm square. The grating face is precisely aligned to the rotation axis of a rotation stage, permitting micrometer adjustment of the grating angle to vary the spectral range. With a camera-collimator angle of 36$`\mathrm{°}`$, IntroSpec’s design range of 4500–6800 Å corresponds to an 8$`\mathrm{°}`$ angle between the grating normal and the camera axis. At this rotation angle, the anamorphic magnification factor is 1.1, so the monochromatic outgoing beams are oval in shape with dimensions $``$25 $`\times `$ 28 mm.

These beams travel $``$90 mm from the grating to the camera lens, where the dispersed colors slightly overfill the 38-mm maximum aperture of the camera lens, a 50 mm focal length $`f/1.4`$ Nikon lens (see the discussion of vignetting in §3.3). The lens focuses the spectrum onto a Meade Pictor 416 XT CCD camera, containing a $`768\times 512`$ KAF-0400 chip with 9 $`\mu `$m square pixels. Attaching the Nikon lens to the Pictor 416 XT required machining a special connector made from a Nikon bayonet mount and a Meade CCD camera mount. Finally, the cord shown exiting the spectrograph enclosure to the front right in Figure 2 connects the CCD camera to the Meade Pictor control box (not shown) for data and power transfer.

### 3.3 Image Quality & Vignetting

We have modeled IntroSpec’s optical performance using the commercially available optical design and analysis code ZEMAX. The optical prescription for the collimator achromat is included in the lens catalog supplied with ZEMAX. The precise optical prescription for the Nikon lens is not available, so for this lens we have substituted the scaled prescription for a 100 mm focal length, $`f/1.4`$ Double Gauss lens from a compilation of lens designs by Smith (1992). The Nikon lens used in IntroSpec is also of the Double Gauss design and should perform similarly.

Our raytraces for the model IntroSpec give images with RMS image diameters between 12 and 30 microns across the target 4500 to 6800 Å wavelength range when the camera lens is stopped down to $`f/1.8`$. The performance is considerably worse when the lens is opened to its full $`f/1.4`$ aperture. Our laboratory tests show that the actual Nikon lens produces CCD images with FWHM diameters of $``$20–25 $`\mu `$m when the grating is illuminated with collimated light, in good agreement with the raytraces.

Our raytraces also allow us to estimate the vignetting in the camera lens as a function of wavelength. As shown in Figure 5, when the camera lens is stopped down to $`f/1.8`$, the on-axis (5650 Å) beam is unvignetted, and $``$68% of the extreme (4500 Å or 6800 Å) beams are transmitted.

### 3.4 Throughput & Optical Fiber Losses

We estimate IntroSpec’s total throughput at about 12% at 5650 Å. Throughput is lower at the ends of the spectral range due to vignetting in the camera lens, the decrease of the grating efficiency off the blaze peak, and the wavelength-dependence of the quantum efficiency of the CCD camera.

Table 3 lists the individual throughput of each of Introspec’s optical elements, not including losses due to guiding errors or the finite apertures of the fibers. With high-quality components, the lens and grating efficiencies are hard to improve. However, current versions of the Meade Pictor CCD cameras use CCDs that have quantum efficiencies of $``$50% at 5650 Å, noticeably higher than the 35% quantum efficiency of the older model in IntroSpec.

Achieving the optical fiber throughput shown in Table 3 requires some attention to minimizing fiber losses. IntroSpec uses multimode, step-index fibers with a pure fused silica core and doped fused silica cladding, because this type of fiber provides the best overall throughput for the wavelengths of interest in astronomical spectroscopy.

Fiber losses come from three main sources: end reflection losses, internal transmission losses, and focal ratio degradation. End reflection losses are simply the Fresnel reflection losses encountered at every material interface; for fused silica, the refractive index is $``$1.46 at the center of the IntroSpec passband, resulting in a loss of $``$3.5% at each fiber end. These losses are hard to avoid as it is usually not practical to apply antireflection coatings to the tiny fiber ends. Internal transmission losses can be estimated using data from the fiber manufacturer and vary according to the OH content of the fused silica and the material preparation. For the 10 meter long fibers used in IntroSpec, internal transmission losses are negligible in the 4500 to 6800 Å range, although transmission losses become more important below 4000 Å. A discussion of the internal transmission of modern fused silica fibers can be found in Lu et al. (1998).

Focal ratio degradation is the major controllable source of fiber losses. If a fiber is manufactured to high standards and not pinched in use, a light cone incident on the fiber axis exits the fiber with approximately the same focal ratio with which it entered. (The fibers we have chosen accept a numerical aperture of 0.22, corresponding to a focal ratio of $``$$`f/2.2`$, much faster than the required $`f/10`$.) However, stresses that locally deform the fiber will scatter light into a faster cone, resulting in focal ratio degradation. This degradation in turn leads to reduced system throughput if the optics are not enlarged to capture the faster cone of light. Our tests have shown that even if fibers are carefully mounted, protected from stress, and properly aligned, focal ratio degradation will scatter about 25% of the light incident at $`f/10`$ outside an $`f/10`$ light cone. Accounting for reflection losses at both ends of the fiber, the total fiber throughput under ideal conditions is then about 70% at $`f/10`$. Focal ratio degradation is much less important for input beams significantly faster than $`f/10`$. We recommend avoiding fibers with core diameters smaller than 100 $`\mu `$m because they are more difficult to handle and are more prone to focal ratio degradation.

### 3.5 Mounting Optical Fibers

The most common source of excess focal ratio degradation is external stress due to inappropriate fiber mounting techniques. For example, if the fibers are bonded into an oversized ferrule with a circumferential epoxy bond, the epoxy shrinkage upon curing can introduce serious stress. A much better technique is to bond the fibers onto precisely machined V-grooves (long grooves with a V-shaped cross-sectional profile) with small beads of silica-filled epoxy. Mixing epoxy and silica half and half by weight reduces the epoxy shrinkage and the viscous mixture does not completely surround the fiber (Figure 3).

To machine the V-grooves, we use a 45$`\mathrm{°}`$ double angle cutter in an R-8 or 3/4” straight shank stub milling machine arbor. Cutters with noticeable radii on the cutting teeth should be avoided, because a large radius will disturb the fiber alignment by allowing the fiber to contact the groove bottom. After the grooves are cut, we lightly bead blast them with glass beads; this process removes burrs that may stress or scratch the fibers and also creates a good surface for bonding.

Stress along the length of the fiber can be minimized by enclosing each fiber in a loosely fitting sleeve of teflon tubing (in addition to any plastic buffer the fiber may already be coated with); for convenience and extra protection, we also bundle these tubes within a larger fabric sheath. Several stages of strain relief at both ends of the fiber train protect the fibers from accidental tugs and maintain a large bend radius.

Accurate fiber alignment and polishing is also critical, both to avoid focal ratio degradation losses at the entrance fiber face and misdirection of the exiting light at the exit fiber face. Both the entrance and exit ends of the fibers must be polished flat to a good optical finish, and the polished ends must be perpendicular to the fiber axes (which are aligned perpendicular to the focal plane using the precisely machined V-grooves). We used a fiber polishing machine built for an earlier project; a similar machine could be constructed for $``$$300 in parts. The machine consists of a flat glass disk attached to an aluminum turntable rotating on a 72 rpm motor. A rigid arm holds a brass cylinder perpendicular to the glass disk. The brass cylinder has V-grooves cut along its length, and the fibers are temporarily attached to the V-grooves with Duco cement for polishing. A large drop of Duco cement covers the end of each fiber, which protrudes by $``$2 mm (most of this will be polished away). We attach disks of high quality self-adhesive abrasive paper in successively finer grades to the turntable, and the brass cylinder containing the fibers is swept across the abrasive paper. The sequence begins with 30 $`\mu `$m grit paper and ends with 0.3 $`\mu `$m grit paper. A thin solution of dishwashing soap in water provides adequate lubrication; clean water is used to rinse off grit residue before moving to the next smaller grit (the motor is isolated from water spills). Finally, we inspect the fiber ends under a microscope to verify that no scratches remain. The Duco cement can be removed with acetone to release the fiber from the polishing cylinder.

### 3.6 The Manual Guider

Acquiring and tracking astronomical targets requires some technique for viewing the focal plane. In a professional instrument this function is usually performed by a sensitive, high frame rate electronic camera. In IntroSpec we have provided optics that allow the observer to visually monitor and correct the telescope tracking while viewing a 12 arcminute field of view imaged around the spectrograph entrance aperture.

Figure 6 illustrates the primary components of the guider. A flat metal mirror intersects the focal plane of the telescope at a 45$`\mathrm{°}`$ angle. The mirror is machined in two pieces and the fibers protrude through the mirror at the intersection of the focal plane and the mirror’s surface. All the light that does not enter the fibers strikes the mirror and is reflected by 90$`\mathrm{°}`$ into the eyepiece extension. Inside the extension, two matched achromatic lenses relay the image of the telescope focal plane to a 26 mm Plossl eyepiece. The observer guides the telescope by centering the object of interest on one of the fiber ends so that most of the light disappears down the fiber.

The guider also serves as the mechanical interface between the telescope and the fibers, incorporating mechanisms for eyepiece rotation, strain relief, and fiber alignment. To minimize weight, nearly all of the parts were made from aluminum, including the mirror. We finished the mirror in two steps, first removing machining marks with a fine grit paste on a turntable and then polishing with white rouge compound and a hard felt wheel mounted on a lathe. The two halves of the mirror were registered with pins, allowing us to disassemble the mirror to mount the fibers in internal V-grooves without disturbing the polished surface.

## 4 Conclusion

IntroSpec saw first light on March 17, 1999 and is now in regular use in a Harvard astronomy course for non-majors. Students perform projects such as comparing the spectra of the two stars

in Albireo. Even under a bright urban sky, we have successfully used IntroSpec to observe objects ranging from stars and planets to emission nebulae and binary accretion systems. Figure 7 shows a few spectra taken with IntroSpec along with their exposure times. We have also used IntroSpec to demonstrate to students how a spectrograph works; the instrument’s easy-to-remove top cover allows students to peer directly inside. In another article (Kannappan & Fabricant, 2000) we provide a brief introduction to amateur spectroscopy with low-budget spectrographs like IntroSpec, including basics of data analysis and interpretation as well as a variety of spectroscopy project ideas.

IntroSpec was designed and built at the Harvard-Smithsonian Center for Astrophysics, with financial support from Harvard College and the Harvard Astronomy Department arranged by Jonathan Grindlay and Ramesh Narayan. We thank Professor Grindlay for his participation in planning and commissioning the instrument. The Knowles Telescope where IntroSpec is used was generously donated by C. Harry and Janet Knowles. We are grateful to John Geary, Steve Amado, Joe Zajac, Douglas Mar, John Roll, Warren Brown, and Anita Bowles for advice, resources, and technical assistance during IntroSpec’s construction phase. We also thank Tom Narita, John Raymond, Dimitar Sasselov, and the students, friends, and family who helped with commissioning runs.


format: text

text:
# Persistence length of a polyelectrolyte in salty water: a Monte-Carlo study

## I Introduction

Despite numerous theoretical studies of polyelectrolyte (PE), due to the long range nature of the Coulomb interaction, the description of their conformation is still not as satisfactory as that of neutral polymers. One of the longest standing problem is related to the electrostatic effect on the rigidity of a PE. In a water solution with monovalent ions, within the Debye-Hückel linear screening theory, the electrostatic interaction between PE charged monomers has the form:

$$V(r)=\frac{e^2}{Dr}\mathrm{exp}\left(\frac{r}{r_s}\right),$$
(1)
where $`r`$ is the distance between monomers, $`D`$ is the dielectric constant of water, $`e`$ is the elementary charge, and $`r_s`$ is the Debye-Hückel screening length, which is related to the ionic strength $`I`$ of the solution by $`r_s^2=4\pi l_BI`$. ($`l_B=e^2/Dk_BT`$ is the Bjerrum length, $`T`$ is the temperature of the solution).

The rigidity of a polymer is usually characterized by one parameter, the so called persistence length $`l_p`$. For a polyelectrolyte chain, besides the intrinsic persistence length $`l_0`$ which results from the specific chemical structure of the monomers and bonds between them, the total persistence length also includes an “electrostatic” contribution $`l_e`$ which results from the screened Coulomb interactions between monomers:

$$l_p=l_0+l_e.$$
(2)
Because the interaction (1) is exponentially screened at distances larger than $`r_s`$, early works concerning the structure of the PE assumed that $`l_e`$ is of the order of $`r_s`$. However, this simple assumption was challenged by the pioneering works of OdijkOdijk  and Skolnick and FixmanSF  (OSF), who showed that Debye-Hückel interaction can induce a rod-like conformation at length scales much larger than $`r_s`$. Their calculation gives

$$l_e=l_{OSF}=\frac{\eta _0^2}{4Dk_BT}r_s^2,$$
(3)
where $`\eta _0`$ is the linear charge density of the PE. Because $`l_er_s^2`$, it can be much larger than $`r_s`$ at weak screening (large $`r_s`$).

Although the idea that electrostatic interaction enhances the stiffness of a PE is qualitatively accepted and confirmed in many experiments, the quadratic dependence of $`l_e`$ on the screening length $`r_s`$ is still the subject of many discussions. In the work of OSF, the bond angle deflection was assumed to be small everywhere along the chain, what is valid for large $`l_0`$. They suggested that if $`l_0`$ is not small but $`l_e`$ is large enough (week screening), their assumption is still valid. Ref. BJ, , however, has questioned this assumption especially when $`l_0`$ is so small that the bond angle deflection is large before electrostatics comes into play and rigidifies the chain.

A significant progress was made by Khokhlov and Khachaturian (KK) who proposed a generalized OSF theoryKK  for the case of flexible polyelectrolyte (small $`l_0`$). It is known that in the absence of screening ($`r_s\mathrm{}`$), the structure of a polyelectrolyte can be conveniently described by introducing the concept of electrostatic blobs. A blob is a chain subunit within which the electrostatic interaction is only a weak perturbation. The blob size $`\xi `$ is related to the number of Kuhn segments $`g`$ within one blob as $`\xi =l_0g^{1/2}`$. The condition of weak Coulomb interaction suggest that the electrostatic self energy of a blob, $`(\eta _0gl_0)^2/D\xi `$ is of the order of $`k_BT`$. This leads to $`\xi (Dk_BTl_0^2/\eta _0^2)^{1/3}`$. At length scale greater than $`\xi `$, Coulomb interaction plays important role and the string of blobs assumes a rod-like conformation, with the end-to-end distance proportional to the number of blobs.

Using this blob picture, KK proposed that OSF theory is still applicable for a flexible PE provided one deals with the chain of blobs instead of the original chain of monomers. This means, in Eq. (3), one replaces the bare linear charge density $`\eta _0`$ by that of the blob chain $`\eta =\eta _0gl_0/\xi `$. The intrinsic persistence length $`l_0`$ should also be replaced by $`\xi `$. As a result, the total persistence length of the flexible PE reads:

$$l_{p,KK}=\xi +\frac{\eta ^2}{4Dk_BT}r_s^2.$$
(4)
Thus, in KK theory, despite the flexibility of the PE, its electrostatic persistence length remains quadratic in $`r_s`$. Small $`l_0`$ only renormalizes the linear charge density from $`\eta _0`$ to $`\eta `$.

Note that $`r_s`$ is implicitly assumed to be larger than the blob size $`\xi `$ in KK theory (weak screening). For strong screening $`r_s<\xi `$, there are no electrostatic rigidity and the chain behaves as flexible chain with the Debye-Hückel short range interaction playing the role of an additional excluded volume interaction.

A number of variational calculations have also been proposed to describe more quantitatively the structure of flexible chain. These calculations, although based on different ansatz, have the same basic idea of describing the flexible charged chain by some model of noninteracting semiflexible chain and variationally optimizing the persistence length of the noninteracting system. Surprisingly, while some of these calculations support the OSF-KK dependence $`l_er_s^2`$ such as Refs. Witten, ; Netz, ; Thirumalai2, , other calculations found that $`l_e`$ scales linearly with $`r_s`$ insteadBJ ; Dawson ; Thirumalai . However, because variational calculation results depend strongly on the variational model Hamiltonian, none of these results can be considered conclusive.

Computer simulationsBB ; Reeds ; Seidel ; Ullner2 ; Ullner3 ; Kremer  also have been used to determine the dependence $`l_e`$ on $`r_s`$ and to verify OSF or variational theories. Some of these papers claim to support the linear dependence of $`l_p`$ on $`r_s`$. The simulation of Ref. Kremer,  concludes that the dependence of $`l_p`$ on $`r_s`$ is sublinear. Thus, the problem of the dependence $`l_e(r_s)`$, despite being very clearly stated, still remains unsolved for a flexible PE. More details about the present status of this problem can be found in Ref. Ullner4, .

In this paper, we again use computer simulations to study the dependence of $`l_e`$ on $`r_s`$. The longest polyelectrolyte simulated in our paper contains 4096 charged monomers, four times more than those studied in previous simulations. This allows for better studying of size effect on the simulation result. Furthermore, we use a more refined analysis of the simulation result, which takes into account local fluctuations in the chain at short distance scale. Our results show that OSF formula quantitatively describes the structure of a polyelectrolyte.

The paper is organized as follows. The procedure of Monte-Carlo simulation of a polyelectrolyte using the primitive freely jointed beads is described in the next section. The data for the end-to-end distance $`R_{ee}`$ is given. In Sec. III, we analyze this data using the scaling argument to show that it is consistent with OSF theory. In Sec. IV, we analyze the data for the case of large $`r_s`$, where excluded volume effect is not important, in order to extract $`l_e`$ and again show that it obeys OSF theory in this limit. In Sec. V, we use the bond angle correlation function to calculate $`l_e`$ and to confirm the result of Sec. IV. The good agreement between $`l_e`$ calculated using different methods further suggests that OSF theory is correct in describing a polyelectrolyte structure. We conclude in Sec. VI.

Several days after the submission of our paper to the Los Alamos preprint archiveNguyen , another paperEveraers  with Monte-Carlo simulations for PE molecules in the same range of lengths appears in the same archive. Results of this paper are in good agreement with our Sec. III.

## II Monte-Carlo simulation

The polyelectrolyte is modeled as a chain of $`N`$ freely jointed hard spherical beads each with charge $`e`$. The bond length of the PE is fixed and equal to $`l_B`$, where $`l_B=e^2/Dk_BT`$ is the Bjerrum length which is about 7Å at room temperature in water solution. Thus the bare linear charge density of our polyelectrolyte is $`\eta _0=e/l_B`$. Because we are concerned about the electrostatic persistence length only, the bead radius is set to zero so that all excluded volume of monomers is provided by the screened Coulomb interaction between them only. For convenience, the middle bead is fixed in space.

To relax the PE configuration globally, the pivot algorithmpivot  is used. In this algorithm, in an attempted move, a part of the chain from a randomly chosen monomer to one end of the chain is rotated by a random angle about a random axis. This algorithm is known to be very efficient. A new independent sample can be produced in a computer time of the order of $`N`$, or in other words, uncorrelated samples are obtained every few Monte Carlo (MC) steps (one MC step is defined as the number of elementary moves such that, on average, every particle attempts to move once). To relax the PE configuration locally, the flip algorithm is used. In this algorithm, a randomly chosen monomer is rotated by a random angle about the axis connecting its two neighbor (if it is one of the end monomers, its new position is chosen randomly on the surface of a sphere with radius $`l_B`$ centered at its neighbor.) In a simulation, the number of pivot moves is about 30% of the total number of moves. The usual Metropolis algorithm is used to accept or reject the move. About $`1÷2\times 10^4`$ MC steps are run for each set of parameters ($`N`$, $`r_s`$), of which 512 initial MC steps are discarded and the rest is used for statistical average (due to time constrain, for $`N=4096`$, only 2000 MC steps are used). Two different initial configurations, a Gaussian coil and a straight rod, were used to ensure that final states are indistinguishable and the systems reaches equilibrium.

The simulation result for the end-to-end distance $`R_{ee}`$ of a polyelectrolyte for different $`N`$ is plotted in Fig. 1 as a function of the screening radius $`r_s`$ of the solution. At very small $`r_s`$, Coulomb interactions between monomers are strongly screened and the chain behaves as a neutral Gaussian chain with $`R_{ee}=l_B\sqrt{N1}`$. At very large $`r_sN`$, Coulomb interactions between the monomers are not screened and $`R_{ee}`$ is saturated and equal to that of an unscreened PE with the same number of monomers (see the arrows in Fig. 1).

Three different methods are used to verify the validity of OSF theory for flexible PE: i) study of the scaling dependence of $`R_{ee}`$ on $`r_s`$ in whole range of $`r_s`$, ii) extraction of $`l_e`$ in the large $`r_s`$ limit and iii) analysis of the bond correlation function. In the next three sections, we discuss these methods in details together with their limitations. Comparison with previous simulations is also made to explain their results which so far have not supported either of the theories.

## III Scaling dependence of $`R_{ee}`$ on $`r_s`$

Let us first describe theoretically how the chain size should behave as a function of the screening radius $`r_s`$ when $`r_s`$ increases from 0 to $`\mathrm{}`$.

When $`r_sl_B`$, the Coulomb interaction is strongly screened. Because there are no other interaction present in our chain model of freely jointed beads, the chain statistic is Gaussian. Its end-to-end distance $`R_{ee}`$ is proportional to the the square root of the number of bonds and independent on $`r_s`$:

$$R_{ee}^2=l_B^2N.$$
(5)
When $`r_sl_B`$, the chain persistence length is dominated by the Coulomb contribution $`l_pl_e`$. If $`N`$ is very large such that the chain contour length $`Nl_B`$ is much larger than $`l_e`$ then the chain behaves as a linear chain with $`Nl_B/l_e`$ segments of length $`l_e`$ each and thickness $`r_s`$. The excluded volume between segments is $`vl_e^2r_s`$, and the end-to-end distanceKK :

$$R_{ee}^2=l_e^2\left(\frac{v}{l_e^3}\right)^{2/5}\left(\frac{Nl_B}{l_e}\right)^{6/5}\{\begin{array}{cc}r_s^{4/5}\hfill & \text{if }l_er_s\hfill \\ r_s^{6/5}\hfill & \text{if }l_er_s^2\hfill \end{array}.$$
(6)
At larger $`r_s`$ where $`l_e`$ becomes comparable to the PE contour length, the excluded volume effect is not important. In this case, the chain statistics is again Gaussian and

$$R_{ee}^2l_e^2\frac{Nl_B}{l_e}\{\begin{array}{cc}r_s\hfill & \text{if }l_er_s\hfill \\ r_s^2\hfill & \text{if }l_er_s^2\hfill \end{array}.$$
(7)
Finally, at even larger $`r_s`$ when $`l_p`$ is greater than $`Nl_B`$, the chain becomes a straight rod with length independent on $`r_s`$:

$$R_{ee}^2l_B^2N^2.$$
(8)
If $`l_er_s^2`$, the transition from the scaling range of Eq. (6) to Eq. (7) happens at $`r_sl_BN^{1/4}`$, while the transition from the scaling range of Eq. (7) to Eq. (8) happens at $`r_sl_BN^{1/2}`$. On the other hand, if $`l_er_s`$, both transitions from the scaling range of Eq. (6) to Eq. (7) and from the scaling range of Eq. (7) to Eq. (8) happen at $`r_sl_BN`$. This means, there is no scaling range of Eq. (7) in this theory.

Thus, one can distinguish between the OSF result, $`l_pr_s^2`$, and the variational result, $`l_pr_s`$ by plotting the exponent $`\alpha =\mathrm{ln}[R_{ee}^2]/\mathrm{ln}r_s`$ as a function of $`\mathrm{ln}r_s`$. The schematic figure of this plot is shown in Fig. 2. OSF theory gives plateaus at $`\alpha =6/5`$ and 2, and when $`r_s>l_BN^{1/2}`$, $`\alpha `$ drops back to 0. Variational theories, on the other hand, would suggest one large plateau at $`\alpha =4/5`$ up to $`r_sl_BN`$.

The simulation results for $`\alpha `$ are shown in Fig. 3 for different $`N`$. One can see that as $`N`$ increases, the agreement with OSF theory becomes more visible. Note that the plateaus in Fig. 2 are scaling ranges, and relatively sharp transitions between plateaus are valid only for $`N\mathrm{}`$. For a finite $`N`$, the plateaus may be too narrow to be observed and can be masked in the transition regions. This explains why one cannot see the plateau at $`\alpha =2`$ in our results. Nevertheless, the tendencies of $`\alpha `$ to develop a plateau at $`\alpha =6/5`$, then to grow higher toward $`\alpha =2`$ at larger $`r_s`$ and finally to collapse to zero when approaching relatively small $`r_s=l_B\sqrt{N}`$ are clearly seen for large $`N`$. Thus, generally speaking, the curves agree with OSF theory much better than with variational theories (where $`\alpha `$ is supposed to be about 4/5 and to decrease to zero only when $`r_sl_BN`$, i.e. at much larger $`r_s`$ than what observed).

Fig. 3 also shows one reason why similar simulations done by other groups do not support OSF theory. All of these simulations are limited to 512 charges. As one can see from Fig. 3, the curves for $`N512`$ do not permit to discriminate between the two theories as clearly as the case $`N=2048`$ or 4096. Only when $`N`$ becomes very large can scaling ranges with $`\alpha >1`$ show up and one observes better agreement with OSF result.

## IV Large $`r_s`$ limit

In this section, we attempt to extract directly from the simulation data the persistence length in order to compare with OSF theory. To do this, one notices that even a chain with excluded volume interaction behaves as a Gaussian chain when its contour length is very short such that it contains only a few Kuhn segments. In this case, one can use the Bresler-Frenkel formulaLandau  to describe the relationship between the end-to-end distance $`R_{ee}`$ and the chain persistence length $`l_p`$:

$$R_{ee}^2=2Ll_p2l_p^2[1\mathrm{exp}(L/l_p)],$$
(9)
where $`L`$ is the contour length of the chain.

For our polyelectrolyte, this formula can be used for large $`r_s`$ when the persistence length is of the order of $`R_{ee}`$ or larger. However, one cannot use the bare contour length $`L_0=(N1)l_B`$ in the Eq. (9) because the chain where OSF theory is supposed to be applicable is not the bare chain but an effective chain which takes into account local fluctuations. The contour length $`L`$ of this effective chain is

$$L=Ne/\eta $$
(10)
where $`\eta `$ is the renormalized linear charge density of the PE.

In KK theory, the effective chain is the chain of electrostatic blobs, and the normalized charge density is $`\eta =\eta _0gl_0/\xi `$. However, the standard blob picture can only be used to describe flexible weakly charged chains where the fraction of charged monomers is small so that the number of monomers, $`g`$, within one blob is large and Gaussian statistics can be used to relate its size and molecular weight. Because, for a given number of charged monomers, Monte-Carlo simulation for weakly charged polyelectrolyte is extremely time consuming, all monomers of our simulated polyelectrolyte are charged. In this case, the neighbor-neighbor monomers interaction equals $`k_BT`$. This makes $`g1`$ and the standard picture of Gaussian blobs does not apply. Thus, in order to treat our data, we assume that both $`l_p`$ and $`\eta `$ are unknown quantities.

To proceed further, one needs an equation relating $`\eta `$ and $`l_p`$, and in order to verify OSF theory, we could use their formula

$$l_p=\eta ^2r_s^2/4Dk_BT,$$
(11)
for this purpose. Thus, we could substitute Eq. (10) and (11) into Eq. (9), and solve for $`\eta `$ using $`R_{ee}^2`$ obtained from simulation. If OSF theory is valid, the obtained values of $`\eta `$ should be a very slow changing function of $`r_s`$. In addition, in the limit $`N\mathrm{}`$, they should also be independent on $`N`$.

The OSF equation (3), however, was derived for the case $`r_sL`$ while in our simulation, the ratio $`r_s/L`$ is not always small. Therefore, instead of Eq. (11), we use the more general Odijk’s finite size formulaOdijk

$$l_p=\frac{\eta ^2r_s^2}{12Dk_BT}\left[3\frac{8r_s}{L}+\left(5+\frac{L}{r_s}+\frac{8r_s}{L}\right)e^{L/r_s}\right].$$
(12)
for the persistence length $`l_p`$. When $`Lr_s`$, the term in the square brackets is equal to 3 and the standard OSF result is recovered. On the other hand, when $`r_sL`$, the persistence length $`l_p`$ saturates at $`\eta ^2L^2/72Dk_BT`$.

Below, we treat our Monte-Carlo simulation data with the help of Eq. (9) using Eq. (10) and (12) for $`L`$ and $`l_p`$. The results for $`\eta `$ are plotted in Fig. 4 for different PE sizes $`N`$.

As one can see, at large $`r_s`$, $`\eta `$ changes very slowly with $`r_s`$, and as $`N`$ increases, tends to saturate at an $`N`$ independent value.

It should be noted that the lines $`\eta (r_s)`$ in Fig. 4 unphysically start to drop below certain values of $`r_s`$. This is because at smaller $`r_s`$, the electrostatics-induced excluded volume interactions between monomers become so strong that the right hand side of Eq. (9) (which is derived for a Gaussian worm like chain) strongly underestimates $`R_{ee}`$.

Even though the picture of Gaussian blobs does not work for our chain, $`\eta `$ can still be calculated analytically in the limit $`Lr_s`$ ($`N\mathrm{}`$). Indeed, let us assume that the effective chain is straight at length scale smaller than $`r_s`$ (which is a reasonable assumption because all the analytical theories so far suggested that the PE persistence length scales as $`r_s`$ or $`r_s^2`$). Thus, the self energy of the chain can be written as $`E=L\eta ^2\mathrm{ln}(r_s/l_B)/D`$. At length scale smaller than $`r_s`$, the polyelectrolyte behaves as a neutral chain under an uniform tension

$$F=E/L=\eta ^2\mathrm{ln}(r_s/l_B)/D.$$
(13)
The average angle a bond vector makes with respect to the axis of the chain, therefore, is:

$`\mathrm{cos}\theta `$ $`=`$ $`{\displaystyle \frac{_0^\pi \mathrm{exp}(Fl_B\mathrm{cos}\theta /k_BT)\mathrm{cos}\theta \mathrm{sin}\theta d\theta }{_0^\pi \mathrm{exp}(Fl_B\mathrm{cos}\theta /k_BT)\mathrm{sin}\theta d\theta }}`$ (14)
$`=`$ $`\mathrm{coth}{\displaystyle \frac{Fl_B}{k_BT}}{\displaystyle \frac{k_BT}{Fl_B}}.`$
The charge density $`\eta `$ can be calculated as

$$\eta =\eta _0/\mathrm{cos}\theta .$$
(15)
At weak screening $`r_sl_B`$, it can be estimated analytically:

$$\eta \eta _0\left[1+\frac{1}{\mathrm{ln}(r_s/l_B)}+\mathrm{}\right],$$
(16)
where the expansion terms of the order of $`1/\mathrm{ln}^2(r_s/l_B)`$ and higher were neglected.

The more accurate numerical solution of Eq. (13), (14) and (15) for $`\eta `$ is plotted in Fig. 4 by the thick solid line. One can see that the values $`\eta (r_s)`$ calculated experimentally using OSF theory with growing $`N`$ converge well to the theoretical curve for $`N=\mathrm{}`$. Remarkably, the theoretical estimate for $`\eta `$ does not use any fitting parameters. This, once again, strongly suggests the OSF theory is valid for flexible PE as well.

The Bresler-Frenkel formula, Eq. (9), is also used to extract the persistence length in Ref. Kremer,  where the authors concluded that the dependence of $`l_e`$ on $`r_s`$ is sublinear. The authors, however, used in Eq. (9) the bare contour length $`L`$, or in other words $`\eta =\eta _0`$, for the calculation of $`l_e`$. As one can see from Fig. 4, this leads to 20-30% overestimation of the contour length of the effective chain where OSF theory is supposed to apply. To show that this overestimation is crucial, let us treat our data similarly to Ref. Kremer,  using $`\eta =\eta _0`$. We plot the resulting dependence of $`l_e/r_s`$ on $`r_s`$ (similarly to Fig. 4 of Ref. Kremer, ) and compare it with our own results using corrected $`\eta `$.

The case $`N=1024`$ is shown in Fig. 5. Obviously, the two results are different qualitatively. While the upper curve follows Eq. (12) with slightly decreasing $`\eta `$, the lower curve shows sublinear growth of $`l_e`$ with $`r_s`$ ($`l_e/r_s`$ is a decreasing function of $`r_s`$). This sublinear dependence observed in Ref. Kremer,  is clearly a manifestation of their overestimation of the PE length $`L`$ which should be used in Eq. (9).

Note that the true $`l_e/r_s`$ curve should also eventually decrease to zero because $`l_e`$ saturates to the constant value $`\eta ^2L^2/72Dk_BT`$ when $`r_sL`$ \[See Eq. (12)\]. But, according to Eq. (12), this decay starts only at very large $`r_s`$ where $`r_s/L0.25`$. The deviation from $`l_er_s^2`$ at large $`r_s`$ seen in Fig. 5 is due to both the violation of the inequality $`r_sL`$ and to the slight decrease of $`\eta `$ with $`r_s`$.

## V Bond angle correlation function

Another standard procedure used in literature is to calculate the persistence length of a polyelectrolyte as the typical decay length of the bond angle correlation function (BACF) along the contour of the chain, assuming the later is exponential.

$$f(|s^{}s|)=\mathrm{cos}[\mathrm{}(𝐛_s,𝐛_s^{})]\mathrm{exp}\left(\frac{|s^{}s|}{l_p}\right).$$
(17)
Here $`𝐛_s`$ and $`𝐛_s^{}`$ are the bond numbered $`s`$ and $`s^{}`$ respectively and $`\mathrm{}(𝐛_s,𝐛_s^{})`$ is the angle between them. The symbol $`\mathrm{}`$ denotes the averaging over different chain conformations. To improve averaging, the pair $`s`$ and $`s^{}`$ are also allowed to move along the chain keeping $`|s^{}s|`$ constant.

We argue in this section that this method of determining persistence length actually has a very limited range of applicability. At either small or large $`r_s`$, the results of persistence length obtained from BACF are not reliable. In the range where this method is supposed to be applicable, we show that the obtained $`l_p`$ are close to those obtained in Sec. IV above.

For small $`r_s`$, excluded volume plays important role and, strictly speaking, it is not clear whether BACF is exponential, and if yes, how one should eliminate excluded volume effect and extract $`l_p`$ from the decay length. According to Ref. Ullner3, , the decay is not exponential in this regime.

The procedure of determining the persistence length using BACF becomes unreliable at large $`r_s`$ as well. To elaborate this point, in Fig. 6a, we plot the logarithm of the bond angle correlation function $`f(x)`$ along the PE contour length for a $`N=512`$ and $`r_s=50l_B`$, typical values of $`N`$ and $`r_s`$ where the excluded volume due to Coulomb interactions is small. There are three regions in this plot. In region A at very small distance along the PE contour length, monomers are within one electrostatic blobs from each other and the effects of Coulomb interaction are small. The bond angle correlation in this region decays over one bond length $`l_B`$. At larger distance along the PE contour length, the region B, the decay is exponential and a constant decay length seems well-defined. Finally, at distance comparable to the chain’s contour length, one again observes a fast drop of the BACF (region C). This end effect is due to the fact that the stress at the end of the chain goes to zero and the end bonds become uncorrelated. The persistence length of interest can be defined as the decay length in region B.

Problem arises, however, at large enough $`r_s`$ when the region C (the end effect) becomes so large that region B is not well defined. In this case the obtained decay length underestimates the correct persistence length. As one can see from Fig. 6, region C can be quite large. It occupies 40% of the available range of $`x`$, even though the screening length is only 10% of the contour length in this case.

There is an even more strict condition on how large $`r_s`$ is when the method of BACF loses its reliability. If $`l_e`$ is larger than $`L`$, the decrease of $`\mathrm{ln}f(x)`$ in region B is less than unity. When this happens, an exponential decay is ambiguous.

Because of all these limitations, in this section we use BACF to calculate $`l_e`$ only in the very limited range of $`r_s`$ where excluded volume is not important and $`l_e`$ is not much larger than $`L`$ (the decrease in region B is greater than 0.1). The obtained $`l_{BACF}`$, which is measured along the chain contour, is compared to $`\eta l_e/\eta _0`$ obtained using the Bresler-Frankel formula in the previous subsection. (The factor $`\eta /\eta _0`$ is needed because $`l_{BACF}`$ is measured along the real PE contour while $`l_e`$ is measured along the renormalized PE contour.) The results are shown in the Table 1. The two persistence lengths are within 20-25% of each other. This reasonably good agreement between two different methods shows that our calculations are consistent. It further strengthens the conclusion of two previous sections that OSF theory is correct in describing flexible polyelectrolytes.

## VI Conclusion

In this paper, we use extensive Monte Carlo simulation to study the dependence of the electrostatic persistence length of a polyelectrolyte on the screening radius of the solution. Not only did we simulate a much longer polyelectrolyte than those studied in previous simulations in order to show the scaling ranges, we also used a refined analyses which take into account local fluctuations to calculate the persistence length. These improvements result in a good support for OSF theory. They also help to explain why previous simulations failed to support OSF theory.

In order to describe our numerical data we used a modified OSF theory in the framework of ideas of KK. Linear charge density $`\eta `$ was corrected to allow for short range fluctuations. In our case this is a relatively small correction to $`\eta _0`$ because we deal with a strongly charged PE. When one crosses over to sufficiently weakly charged PE linear charge density becomes strongly renormalized and matches KK expressions. We confirmed that corrections of $`\eta `$ do not affect $`r_s^2`$ dependence of persistence length which was predicted by OSF for $`l_0r_sL`$. In other words, we confirm KK idea that at large $`r_s`$ all effects of flexibility of PE are limited to a renormalization of $`\eta `$. At $`r_s`$ comparable to contour length $`L`$ we found a good agreement of the numerical data with OSF formula modified for this case \[Eq. (12)\], which is derived in Ref. Odijk, . Again all effects of local flexibility are isolated in the small correction to the linear charge density $`\eta `$.

###### Acknowledgements.

The authors are grateful to A. Yu. Grosberg M. Rubinstein, M. Ullner and R. Netz for useful discussions and comments. This work is supported by NSF No. DMR-9985785. T.T.N. is also supported by the Doctoral Dissertation Fellowship of the University of Minnesota.


text:
# Simple analytic evaluation of the integrals in the Gaussian model of Bose-Einstein correlations

## Abstract

A generating function for products of Hermite polynomials is used to significantly simplify the evaluation of the integrals $`G_n(p,p^{})`$ occurring in the Gaussian model of multiple particle production. These integrals are crucial for studies of multiparticle effects in Bose-Einstein correlations.

PACS numbers 25.75.Gz, 13.65.+i
Bose-Einstein correlations, Gaussian model.

1. Bose-Einstein correlations in multiple particle production processes are a field of intense study. For recent reviews see -. An important subfield is the study of multiparticle effects on multiplicity and momenta distributions. Here exactly soluble models are a convenient starting point. A soluble model, which is particularly popular, is the Gaussian model  \- . In this model all the integrations are Gaussian and consequently can be performed analytically. The generic integral is

$$G_n(p,p^{})=\rho (p,k_1)𝑑k_1\rho (k_2,k_3)\mathrm{}𝑑k_{n1}\rho (k_{n1},p^{}),$$
(1)
where

$$\rho (p,p^{})=\frac{1}{\sqrt{2\pi \mathrm{\Delta }^2}}\mathrm{exp}\left[\frac{K^2}{2\mathrm{\Delta }^2}\frac{1}{2}R^2q^2\right]$$
(2)
with

$$K=\frac{p+p^{}}{2};q=pp^{}.$$
(3)
The constants $`R^2`$ and $`\mathrm{\Delta }^2`$ are real and positive. For simplicity we start with the one-dimensional case. In the Gaussian model the transition to three dimensions is trivial. It is seen that the integrals $`G_n`$ are Gaussian and consequently can be performed analytically. One needs, however, a closed formula valid for any $`n`$. This is more difficult to obtain. Zimányi and Csörgö  following some suggestions from Pratt  constructed a system of recurrence relations for three series of parameters and after a rather lengthy calculation got a formula for the integrals, which can be written in a fairly simple form . Wiedemann  applied successfully the MATHEMATICA program directly to the integrals $`G_n`$. In the present note we point out that a very short derivation of the formula for the integrals $`G_n`$ can be given, if a mathematical identity for Hermite polynomials is made use of. This work is an extension of the method used in .

2. An identity, valid for $`|z|<\frac{1}{2}`$, is

$$\mathrm{\Phi }(x,y,z)\underset{n=0}{\overset{\mathrm{}}{}}\frac{z^n}{n!}H_n(x)H_n(y)=\frac{1}{\sqrt{14z^2}}\mathrm{exp}\left[4z\frac{z(x^2+y^2)xy}{14z^2}\right].$$
(4)
A convenient starting point to check this formula is the definition of the Hermite polynomials in the form 222In this edition of the tables the power of two is mistakenly written as $`\frac{n}{2}`$ instead of $`n`$. In order to check the formula it is enough to expand the binomial into powers of $`z`$, integrate term by term and compare the coefficients with the formulae given in .

$$H_n(z)=\frac{2^n}{\sqrt{\pi }}_{\mathrm{}}^+\mathrm{}(z+it)^ne^{t^2}𝑑t.$$
(5)
After substituting it into (4) the summation over $`n`$ gives an exponential and it is enough to perform a double Gaussian integral to get the right-hand side. The function $`\mathrm{\Phi }(x,y,z)`$ is defined by this formula.

3. For any orthonormal set of functions $`\psi _n(p)`$ and any set of nonnegative numbers $`\lambda _n`$ satisfying the normalization condition

$$\underset{n}{}\lambda _n=1,$$
(6)
the function

$$\rho (p,p^{})=\underset{n}{}\psi _n(p)\lambda _n\psi _n^{}(p^{})$$
(7)
can be interpreted as a well-defined single-particle density matrix in the momentum representation. The functions $`\psi _n`$ are its eigenfunctions and the numbers $`\lambda _n`$ are the corresponding eigenvalues. Let us choose in particular for the functions $`\psi _n`$ the energy eigenfunctions in the momentum representation of a one-dimensional harmonic oscillator and for the numbers $`\lambda _n`$ a geometric progression :

$$\psi _n(p)=\sqrt{\frac{\alpha }{\sqrt{\pi }2^nn!}}\mathrm{exp}\left[\frac{\alpha ^2p^2}{2}\right]H_n(\alpha p);\lambda _n=(1z)z^n.$$
(8)
It is assumed here that $`|z|<1`$, so that the eigenvalues $`\lambda _n`$ are correctly normalized. The parameter $`\alpha `$ is positive. Its interpretation in terms of the mechanical parameters of the oscillator is of no importance here. This choice corresponds to a single particle density matrix which using (4) can be written in the form

$$\rho (p,p^{})=\frac{\alpha (1z)}{\sqrt{\pi }}e^{\frac{1}{2}\alpha ^2(p^2+p^2)}\mathrm{\Phi }(\alpha p,\alpha p^{},\frac{z}{2}),$$
(9)
or substituting the definition of the function $`\mathrm{\Phi }`$ and replacing the parameters $`p,p^{}`$ by the parameters $`K,q`$ according to (3)

$$\rho (K,q)=\frac{\alpha }{\sqrt{\pi }}\sqrt{\frac{1z}{1+z}}\mathrm{exp}\left[\alpha ^2K^2\frac{1z}{1+z}\frac{\alpha ^2q^2}{4}\frac{1+z}{1z}\right].$$
(10)
This density matrix differs only by notation from the density matrix (2). The correspondence between the parameters is

$$\alpha =\sqrt{\frac{R}{\mathrm{\Delta }}};z=\frac{v1}{v+1}.$$
(11)
Following  we have introduced the notation

$$v=2\mathrm{\Delta }R1,$$
(12)
where the inequality follows from the uncertainty principle.

4. In order to calculate the function $`G_n(p,p^{})`$ let us note first that it is the $`p,p^{}`$ matrix element of the $`n`$-th power of the matrix $`\rho `$. Therefore, it has the same eigenfunctions and its eigenvalues are $`n`$-th powers of the corresponding eigenvalues of matrix $`\rho `$. Thus

$$G_n(p,p^{})=\underset{k=0}{\overset{\mathrm{}}{}}\psi _k(p)\psi _k(p^{})\lambda _k^n=\frac{\alpha (1z)^n}{\sqrt{\pi }}e^{\frac{1}{2}\alpha ^2(p^2+p^2)}\mathrm{\Phi }(\alpha p,\alpha p^{},\frac{z^n}{2}).$$
(13)
Substituting the definition of function $`\mathrm{\Phi }`$ and changing variables from $`p,p^{}`$ to $`K,q`$ we find

$$G_n(K,q)=\frac{\alpha (1z)^n}{\sqrt{\pi (1z^{2n})}}\mathrm{exp}\left[\alpha ^2K^2\frac{1z^n}{1+z^n}\frac{1}{4}\alpha ^2q^2\frac{1+z^n}{1z^n}\right].$$
(14)
It is convenient to rewrite this formula in the form

$$G_n(K,q)=\frac{c_n}{\sqrt{2\pi \mathrm{\Delta }^2}}\mathrm{exp}\left[\frac{K^2}{2\mathrm{\Delta }_n^2}\frac{1}{2}R_n^2q^2\right].$$
(15)
By comparison with the previous form $`\mathrm{\Delta }_n=a_n\mathrm{\Delta }^2;R_n^2=a_nR^2`$ with

$$a_n=\frac{1}{v}\frac{1+z^n}{1z^n}=\frac{1}{v}\frac{(v+1)^n+(v1)^n}{(v+1)^n(v1)^n}$$
(16)
and

$$c_n=\frac{2^n\sqrt{v}}{\sqrt{(v+1)^{2n}(v1)^{2n}}}.$$
(17)
5. Let us add some comments.

* As mentioned, the transition from the one-dimensional case to $`d`$-dimensions is very simple. The density matrix (2) gets replaced by a product

$$\rho (\stackrel{}{p},\stackrel{}{p^{}})=\underset{i=1}{\overset{d}{}}\frac{1}{\sqrt{2\pi \mathrm{\Delta }_i^2}}\mathrm{exp}\left[\frac{K_i^2}{2\mathrm{\Delta }_i^2}\frac{1}{2}R_i^2q^2\right].$$
(18)
Consequently, the integrals $`G_m`$ from (1) factorize and one obtains for each $`G_m`$ a product of expressions (15), in general with different values of $`\mathrm{\Delta }^2`$ and $`v`$ for each factor.
* For the harmonic oscillator the energy eigenfunctions in the momentum representation and in the coordinate representation have the same form. Besides the replacement of $`p`$ by $`x`$, the interpretation of the parameter $`\alpha `$ changes, but this is irrelevant for the present analysis. Thus, the results given in the present paper can be used, after some changes in notation, also for the study of Bose-Einstein correlations in ordinary space.
* Since the energy eigenvalues of the harmonic oscillator are equidistant, the corresponding Boltzmann factors form a geometrical progression. Thus, formula (9) as well as the corresponding formula in the coordinate representation, can be used to calculate the canonical density matrix for the harmonic oscillator. According to Landau and Lifszyc , this calculation for the diagonal elements of the density matrix in the coordinate representation has been first performed by F. Bloch in year 1932.

The author thanks Andrzej Bialas for collaboration and stimulating discussions.


text:
# TeV gamma rays from the blazar H 1426+428 and the diffuse extragalactic background radiation

## 1 Introduction

Many nonthermal extragalactic objects representing different classes of AGNs are considered as potential sources of TeV photons. First of all this concerns the BL Lac population of blazars in general; two nearby representatives of this class, Mkn 421 and Mkn 501 with redshifts of $`z`$=0.030 and $`z`$=0.034, respectively, are firmly established as TeV $`\gamma `$-ray emitters. Of special interest are the so-called “extreme synchrotron blazars”, BL Lac objects with flat spectra of synchrotron emission and high X-ray to radio flux ratios (e.g. Costamante et al. 2001). The observations of two such objects, Mkn 501 (Pian et al. 1998) and 1ES 2344+514 (Giommi et al. 2000), by BeppoSAX showed that in a flaring state the synchrotron peak of “extreme blazars” can reach 100 keV. Remarkably, both these objects are also reported as TeV blazars (see e.g. Catanese & Weekes 1999). Most probably, this is not a mere coincidence; since the high synchrotron peak is an indicator of acceleration of electrons to ultrarelativistic energies, the “extreme blazars” are obvious candidates for TeV emission. Consequently, some of these objects were intensively monitored with imaging atmospheric Cherenkov telescopes (IACT) over the last several years. Possible detections of TeV signals have been claimed by different groups for 1ES 2344+514 (Catanese et al. 1998), 1ES 1959+650 (Nishiyama et al. 1999), and PKS 2155-304 (Chadwick et al. 1999). Most recently, statistically significant TeV signals have been reported also for H 1426+428, by the VERITAS (Horan & the VERITAS coll. 2001) and HEGRA (Götting et al. 2001; Aharonian 2001) collaborations. The relatively large redshift ($`z`$=0.129) of H 1426+428, as well as the fact that it is the third 100 keV blazar ever found (Costamante et al. 2001), make this object an extremely important target for future multiwavelength studies of the nonthermal processes of particle acceleration and radiation in AGN jets, as well as for probing the DEBRA in the near infrared region.

In this letter we report the HEGRA observations of H 1426+428 performed in 1999 and 2000, and briefly discuss some apparent astrophysical implications suggested by the obtained results.

## 2 Observations and results of analyses

For the following analysis, 14.5 hrs of data from 1999 and 29.9 hrs from 2000 were available, taken with the HEGRA IACT system; see Table 1 for the dates of observations. The median zenith angle was 17°, resulting in a $`\gamma `$-shower peak detection rate at an energy of $`700\text{GeV}`$ for Crab-like spectra (Konopelko et al. 1999).

Shower reconstruction and cuts were described e.g. in Aharonian et al. (1999a). The methods applied in this analysis differ slightly, depending on the tasks which were performed. For the search for new sources we applied tight cuts and stereo algorithm #3 (Hofmann et al. 1999) to achieve high sensitivity. For the spectral analysis of the detected source, somewhat relaxed cuts were used which provide a larger $`\gamma `$-acceptance; these cuts are less susceptible to systematic errors, and the spectral reconstruction was tested in great detail using Monte Carlo studies and $`\gamma `$-ray data from known TeV sources (Aharonian et al. 1999a, b). The background was estimated using a set of control regions in the field of view (FOV); the setup provides seven times more background data than on source data under identical acceptance conditions (Aharonian et al. 2001). The optimum angular cut (for the signal search) was derived on the basis of the angular resolution – using $`\gamma `$-ray events from the Crab nebula – and the background level. Table 2 summarizes cut parameters, resulting event numbers, and significances for the signal search and spectral analysis.

Figure 1 shows the event distribution, obtained from the signal search, as a function of the squared angular distance to the source position. The excess significance, calculated according to Li & Ma (1983), amounts to $`5.8\sigma `$ ($`2.4\sigma `$ in 1999, $`5.3\sigma `$ in 2000). While there is an indication that the source was in a higher state in May 2000, the total data set is also statistically consistent with constant signal accumulation.

The event distribution in the FOV was used to reconstruct the source position as seen in TeV $`\gamma `$-rays; for details of the procedure see Pühlhofer et al. (1997) and Aharonian et al. (2001). Figure 2 shows the sky area around H 1426+428. The positions of X-ray selected QSOs and galaxies were obtained from the NASA/IPAC Extragalactic Database (NED). The median angular resolution of the HEGRA system is approx. 0$`\stackrel{}{.}`$1; the given statistics allowed the reconstruction of the center of TeV emission with a statistical error of 50″(1 $`\sigma `$). A confusion of the TeV source with other known X-ray sources is excluded. We note that the very prominent X-ray source GB 1428+4217 is located at an angular distance of 41$`\stackrel{}{.}`$2 to H 1426+428, far outside the area shown in Fig. 2.

The energy reconstruction, which has a single event resolution of $`\mathrm{\Delta }E/E=20\%`$, and the spectral evaluation were described in detail in Aharonian et al. (1999a, b). The energy spectrum was derived from the raw, background subtracted photon count spectrum using a so called effective area, which is a response function depending upon the reconstructed energy and zenith angle; the effective area was adjusted regularly according to the varying detector conditions and different system setups (3/4/5 telescopes included in the system).

Figure 3 shows the reconstructed differential energy spectrum of H 1426+428. To first order, the spectrum can be described by a pure power law $`\mathrm{d}\mathrm{\Phi }/\mathrm{d}E=(2.0\pm 1.3_{\mathrm{stat}}\pm 0.1_{\mathrm{syst}})\times 10^{12}(E/\mathrm{TeV})^{(2.6\pm 0.6_{\mathrm{stat}}\pm 0.1_{\mathrm{syst}})}\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1\mathrm{TeV}^1`$ (the systematic errors do not include the 15% error on the energy scale). The spectral slope is similar to the one measured for the Crab nebula (Aharonian et al. 2000). The measured integral flux (derived using loose cuts) above 1 TeV is $`(1.7\pm 0.5_{\mathrm{stat}}\pm 0.1_{\mathrm{syst}})\times 10^{12}\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1`$ which corresponds to 10% of the Crab nebula flux. A comparison of the total photon count rates of H 1426+428 and Crab nebula after tight cuts yields 6%; within statistical errors, both results are compatible.

## 3 Discussion

For any reasonable DEBRA model (see Fig. 4 a), the TeV radiation from H 1426+428 is expected to arrive with a drastically modified spectrum. Reasons are the strong energy dependence of the mean free path of $`\gamma `$-rays in the intergalactic medium, $`\mathrm{\Lambda }(E)`$, and the fact that the optical depth $`\tau =czH_0^1\mathrm{\Lambda }^1(E)`$ significantly exceeds unity at all $`\gamma `$-ray energies above 300 GeV for H 1426+428 with a redshift of $`z`$=0.129. The conclusion about strong absorption at all energies covered by HEGRA does not depend on the large systematic uncertainties in the reported DEBRA fluxes. This is demonstrated in Fig. 4 b where we show the intergalactic absorption factors calculated for two reference DEBRA models. Model 1 is close to the recent calculations of Primack et al. (2001) based on the so-called Kennicut initial mass function (IMF), and fits quite well the DEBRA fluxes reported recently by Wright & Johnson (2001) at 1.25, 2.2 and 3.5 $`\mu \mathrm{m}`$ wavelengths. Model 2 is designed to match the NIR fluxes reported by Matsumoto (2000) and Cambresy et al. (2001) who claim very high fluxes at wavelengths shorter than $`2\mu \mathrm{m}`$. Despite significant quantitative differences, the two DEBRA models result in similar absorption features in the TeV spectrum: both models predict a very strong steepening of the spectrum at energies below 1-2 TeV, but an almost energy independent absorption at energies around several TeV (see Fig. 4 b). This effect is common for all realistic model curves describing the DEBRA in the near IR (e.g. Primack et al. 2001), since they reflect the characteristic shape of the starlight spectrum which in the wavelength band between 1 and several microns behaves as $`\nu I(\nu )\lambda ^\beta `$ with $`\beta 1`$. The corresponding number density of background photons $`n(ϵ)ϵ^1`$ results in a nearly constant optical depth for $`\gamma `$-rays at energies around several TeV (Aharonian et al. 1999b).

The X-ray spectrum of H 1426+428 which extends up to 100 keV favours an intrinsic pure power law spectrum in the relevant energy range, with an index of 1.9 (a detailed discussion is beyond the scope of this paper and will be given elsewhere). The observed spectrum – after absorption – is therefore not expected to obey a power law. The results presented in Fig. 3 show that the HEGRA spectral points are indeed better described by an intrinsic power law spectrum modified by intergalactic absorption. However, a power law fit with an adjustable spectral index results in a reduced $`\chi ^2/\mathrm{dof}=4.4/3=1.47`$, while absorbed spectra – with a fixed source spectrum $`E^{1.9}`$ – yield $`\chi ^2/\mathrm{dof}=3.7/4=0.9`$ and $`\chi ^2/\mathrm{dof}=1.0/4=0.25`$, for the DEBRA models 1 and 2 (see Fig. 4), respectively. Hence, more accurate spectral measurements, as well as detailed theoretical studies (concerning both the formation of the $`\gamma `$-ray spectra in the source and their subsequent deformation during the passage through the extragalactic photon fields), are required.

Nevertheless, the sheer fact of detection of TeV $`\gamma `$-rays from H 1426+428 should lead to a revision of the current conceptual view of TeV blazars, according to which the synchrotron (X-ray) peak in the SED dominates over the inverse Compton (TeV) peak (see e.g. Fossati et al. 1998). Indeed, although the detected energy flux of $`\gamma `$-rays is only $`4\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, corrected for intergalactic absorption this flux may well exceed $`10^{10}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. For comparison, the X-ray flux measured by BeppoSAX (Costamante et al. 2001) and ASCA (T. Takahashi, private communication) is well below $`10^{10}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. Since the corrected TeV luminosity seems to exceed the level anticipated from the current models of TeV blazars by far, this result may have a crucial impact on the further development of models for TeV blazars.

H 1426+428 is the most distant source detected in TeV $`\gamma `$-rays so far. The HEGRA collaboration intends to continue observing this source extensively in 2002 in the hope to decisively improve the statistics of the measurement.

###### Acknowledgements.

The support of the German ministry for research and technology (BMBF) and of the Spanish Research Council (CICYT) is gratefully acknowledged. GR acknowledges receipt of a Humboldt fellowship. We thank the Instituto de Astrofísica de Canarias for the use of the site and for supplying excellent working conditions at La Palma.