text:
# A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy

## 1 Introduction

At a distance of $`4.2\pm 0.8`$ Mpc (Freeman et al.,, 1977), the Circinus galaxy hosts one of the nearest Seyfert II nuclei. The bolometric luminosity is $`4\times 10^{43}`$ ergs s-1, $`10\%`$ of which is contributed by star formation in the inner $`100`$ pc (Moorwood et al.,, 1996). The central engine radiates $`2\times 10^{42}`$ erg s-1 (2-10 keV) and is obscured by a n$`{}_{\mathrm{H}}{}^{}4\times 10^{24}`$ cm-2 column at energies lower than $`10`$ keV (Matt et al.,, 1999). Reflection by cool material dominates the X-ray spectrum, which includes a prominent ($`2`$ keV equivalent width) Fe K$`\alpha `$ fluorescence line (Matt et al.,, 1996). The optical spectrum includes coronal lines from photoionized gas (Oliva et al.,, 1994) and broad polarized H$`\alpha `$ emission, which is $`3300`$ km s-1 wide and betrays the existence of an obscured broad line region (Oliva et al.,, 1998). The mass of the central engine has been uncertain. At $`2.2\mu `$m wavelength, the central non-stellar source is $`<1.5`$ pc in radius, and its dynamical mass is $`<4\times 10^6`$ M(Maiolino et al.,, 1998). Although massive, the central engine constitutes a small fraction of the nuclear mass, which including stars and gas amounts to $`3.2\pm 0.8\times 10^8\mathrm{M}_{}`$ inside a radius of 140 pc (Curran et al.,, 1998).

The active galactic nucleus (AGN) of the Circinus galaxy drives a visible nuclear outflow on kpc scales, the position angle of which suggests an orientation for the central engine. Elmouttie et al., (1998) observe bipolar radio lobes, which lie at position angles (PA) of $`115`$ and $`295^{}\pm 5^{}`$. These are aligned with the minor axes of the galactic H I disk (PA$`=300^{}\pm 5^{}`$) (Freeman et al.,, 1977) and nuclear 12CO ring (PA$`=304^{}\pm 4^{}`$), which lies between radii of 140 and 600 pc and is inclined by $`78^{}\pm 1^{}`$ (Curran et al.,, 1998). Circinus also exhibits a kpc-scale, one-sided optical ionization cone (Marconi et al.,, 1994) with an opening angle of $`100^{}`$ and a mean PA$`295^{}`$, approximately along the outflow axis. The ionization cone comprises knots and radial filaments whose kinematics are indicative of outward motion, possibly in the form of bullet-like ejecta, at speeds on the order of 100 to 200 km s-1 (Veilleux & Bland-Hawthorn, 1997; Elmouttie et al.,, 1998; Wilson et al.,, 2000).

There is good evidence that extragalactic H2O masers ($`\lambda 1.3`$ cm) in NGC 4258, NGC 1068, and other galaxies (e.g., Miyoshi et al.,, 1995; Greenhill & Gwinn, 1997) trace edge-on, warped accretion disks bound by central engines $`10^6`$ M. Spectral-line interferometer data with milliarcsecond resolution straightforwardly provide accurate measurements of rotation curves, central engine masses, warp shapes, and disk orientations. The maser spectra in the systems that are the best examples exhibit complexes of “high-velocity” emission symmetrically bracketing “low-velocity” emission that is close to the systemic velocity of the host galaxy. In Circinus, H2O maser emission has been known for over 20 years (Gardner & Whiteoak, 1982), but only the relatively recent discoveries of spectral features blueward of the systemic velocity (Nakai et al.,, 1995; Braatz, Wilson, & Henkel, 1996; Greenhill et al.,, 1997) have led to speculation that the masers in Circinus also trace an accretion disk.

We present the first maps of the Circinus H2O maser source and report the detection of an accretion disk, bearing out the earlier speculation. We have also found that a fraction of the maser emission originates outside the disk, in what appears to be a wide-angle outflow within $`1`$ pc of the central engine and aligned with the ionization cone. The discovery of this qualitatively new phenomenon introduces a fresh element to the study of AGN, because the locations and velocities of dense non-disk material in such close proximity to a central engine can be mapped for the first time. As a result, it is possible to test in new ways proposed mechanisms for the driving and collimation of outflows in AGN (e.g., Bottorff et al.,, 1997, and references therein).

In Sections 2, 3, and 4, we summarize our Very Long Baseline Interferometer (VLBI) observations, post-correlation analysis, and proposed model of a warped disk and outflow. In Section 5 we quantify the model and evaluate ramifications for collimation and driving mechanisms, and discuss the relation of the warped disk to the outflow and the ionization cone.

## 2 Observations

We observed the $`6_{16}5_{23}`$ transition of H2O ($`\nu _{\mathrm{rest}}=22235.08`$ MHz) toward Circinus three times in 1997 and 1998 (Table 1) with the Australia Telescope Long Baseline Array (LBA)111The Australia Telescope Long Baseline Array is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO., which provided baselines of $`200`$ to 1000 km. Table 2 summarizes the characteristics of the individual antennas. Typical tracks were 18 h at Mopra and Hobart, 12 h at Tidbinbilla, and 14 h at Parkes. Approximately every 1.5 h we observed a continuum source for $`0.5`$ h (PKS 0537-448, PKS 1424-418, PKS 1144-379, or PKS 1921-293) to enable calibration of the interferometer and antenna pointing at Parkes and Tidbinbilla.

We recorded two 16 MHz bandpasses ($`215\mathrm{km}\mathrm{s}^1`$) at each station with an S2 VLBI terminal (Table 1). The S2 correlator at the Australia Telescope (AT) facility in Marsfield, NSW provided 1024 channels in each bandpass for cross-power and total-power spectra, following Fourier inversion of the correlation functions. The channel spacing was $`0.21`$ km s-1, which provided at least three channels across the half-power full-width of the narrowest spectral components (e.g., see Nakai et al.,, 1995; Greenhill et al.,, 1997).

## 3 Calibration and Imaging

We calibrated the amplitudes, delays, and phases of the data with standard VLBI techniques, using the AIPS package of the NRAO222The National Radio Astronomy Observatory is operated by Associated Universities, Inc., under cooperative agreement with the National Science Foundation. and the module ATLOD that is maintained separately by the AT.

We computed a time-series of amplitude calibration factors for each station by fitting total-power spectra for the line triplet near 565 km s-1 (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy) to template total-power spectra, constructed from data recorded at Mopra in 1997 July and Tidbinbilla in 1998 June. (We assume the radio astronomical definition of Doppler shift and heliocentric reference frame.) This calibration corrected for variations in antenna gain and atmospheric opacity. The template scans were calibrated with nominal system equivalent flux densities (SEFD) of 900 and 150 Jy, respectively, which provided a flux density scale accurate to $`30\%`$, given the elevations of the template scans and overall good weather conditions. The weak line strength and irregular spectral baselines at several stations in 1997 June precluded template fitting. Instead, amplitude calibrations were adjusted to maintain constant measured amplitude for the lines near $`565`$ km s-1 in the cross-power spectra, which is appropriate for relatively point-like emission. A total-power spectrum from Mopra, calibrated with the nominal SEFD, provided the absolute amplitude calibration.

We note that the quality of these amplitude calibrations was limited by the effects of interstellar scintillation, which in one instance has caused a doubling in line strength on time scales as short at 10 minutes (Greenhill et al.,, 1997). However, for the 1998 June and 1997 July epochs the effects of scintillation appear to have been relatively small. The total-power spectra constructed for the template fitting show variations in line strength of $`30\%`$ over time scales on the order of 30 minutes. Because we could not construct reliable templates for the 1997 June data, that epoch could have been more significantly affected by scintillation. This would reduce the dynamic range of images but not affect the centroid positions of individual emission features.

We used the calibrator scans to estimate delay, fringe rate, and fringe phase residuals induced by errors in the a priori station clocks and positions, and by electronics and the troposphere. We interpolated the delays and rates and applied them to scans of the maser. For the 1998 data, the calibration was accurate to $`2`$ ns and $`2`$ mHz in delay and rate, respectively. For the 1997 data, the accuracy for Mopra was worse, $`5`$ ns and $`5`$ mHz, because the maser clock was inadequately shielded against the (pointing-dependent) magnetic signature of the telescope structure. We also used observations of the calibrators PKS 1424-418 and PKS 1921-293 to estimate the time averaged (complex) bandpass response of each station. Overall, residual errors in calibration tended chiefly to scatter power in the synthesis images and reduce the dynamic range. Induced errors in the positions of emission features were ultimately much less than the size interferometer beam and did not affect the modeling.

The a priori position of the maser was another potential source of delay calibration error. Gardner & Whiteoak (1982) found that the emission lay within 10 to 15′′ of the optical center for the nucleus (Freeman et al.,, 1977), which leaves a very large uncertainty. We analyzed the time variation of fringe rate for the 565 km s-1 line triplet observed in 1997 July to improve the astrometric position. A $`\chi `$-squared minimization for a point-source model yielded a new position estimate, $`\alpha _{2000}=14^h13^m09\text{.}^s95\pm 0{}_{}{}^{s}.02`$, $`\delta _{2000}=65^{}20^{}21\text{.}^{\prime \prime }2\pm 0\text{.}^{\prime \prime }1`$, where the quoted uncertainties include random and systematic components. This position is offset from the a priori position by $`+1\text{.}^{\prime \prime }5`$ and $`+0\text{.}^{\prime \prime }1`$ in right ascension and declination, respectively. We corrected the 1997 July data and used the improved position estimate in correlation of the 1997 June and 1998 June epochs, after which the residual fringe rates were $`2`$ mHz (peak-to-peak) on even the longest baselines. This is consistent with the remaining $`0\text{.}^{\prime \prime }1`$ uncertainty.

To remove the effects of atmospheric path-length fluctuations at each epoch, we self-calibrated the emission within a few km s-1 of 565 km s-1 and applied the results to each spectral channel. We adjusted the weighting of data in the (u,v)-plane and the frequency-averaging of visibility data for each epoch, so that we could achieve a mean angular resolution of $`2.5`$ milliarcseconds (mas) and good signal-to-noise ratios ($`S/N`$). (We convolved the data with a 1.3 km s-1 wide boxcar, sampling every 0.65 km s-1 for 1997 June, a 0.87 km s-1 boxcar with sampling every 0.44 km s-1 for 1997 July, and a 0.44 km s-1 boxcar with sampling every 0.44 km s-1 for 1998 June.) The noise in the images ($`1\sigma `$) was 0.025 - 0.045 Jy, depending on the epoch and channel; the 1998 June observations were the most sensitive and best calibrated, with the broadest bandwidth (Table 1). To test the completeness of our census of emission detected in 1998 June, we also constructed images with natural weighting and achieved a sensitivity of 15 mJy over 0.44 km s-1. In these images, we detected two new emission clumps at the red and blue extrema of the measured spectrum. In the end, we detected emission between 179 and 699 km s-1(Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy).

The synthesized beamwidths were $`3.7\times 2.0`$ mas at PA $`=77^{}`$ (1997 June), $`3.7\times 1.4`$ mas at PA $`=80^{}`$ (1997 July), $`3.6\times 1.9`$ at PA$`=84^{}`$ (1998 June 27), and $`4.0\times 1.4`$ mas at PA=$`79^{}`$ (1998 June 28). We fit a 2-D Gaussian model brightness distribution to each emission “maser spot” stronger than $`5\sigma `$ in each deconvolved image. Because the image pixels are partially correlated, we adopted a position uncertainty for each fitted Gaussian that was (1) the formal fitting error or (2) the statistical error dictated by the beamwidth divided by $`2\times S/N`$ (Reid & Moran, 1988; Fomalont, 1999), whichever was larger. To this uncertainty, we added in quadrature a term accounting for the effect of possible error in the absolute position of the 565 km s-1 emission ($`0\text{.}^{\prime \prime }1`$). This term contributed 35$`\mu `$as per 100 km s-1 of velocity offset from $`565\mathrm{km}\mathrm{s}^1`$. For 1998 June and 1997 July, the spreads in spot positions for individual clumps of emission are consistent with the computed total position uncertainties. The larger spreads observed in 1997 June are suggestive of an unmodeled systematic component of position uncertainty, probably stemming from the lower quality calibration and limited coverage of the (u,v)-plane for that epoch.

## 4 Results and Interpretation

At each epoch, the sky distribution of maser emission appears to comprise two populations (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy) that we will argue sample a roughly edge-on warped accretion disk as well as an outflow of material from the vicinity of the central engine. There is a thin, densely sampled, S-shaped locus that comprises redshifted emission to the west, and blueshifted emission to the east. The arms of the S exhibit velocity gradients along their lengths with increasingly large Doppler shifts (relative to the systemic velocity) toward the center. Most of the remainder of the maser emission comprises knots that lie outside the S and away from the putative disk. However, even when relatively close to the S, the line-of-sight velocity of the non-disk maser emission is offset by tens of km s-1. We propose that this non-disk emission originates in a wide-angle nuclear outflow.

The velocities of the most red and blueshifted maser spots, 179 km s-1 and 699 km s-1 (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy), symmetrically bracket the nominal systemic velocity of Circinus, $`439\pm 2`$ km s-1(Freeman et al.,, 1977; Curran et al.,, 1998). The overall spectrum and sky distribution changed from epoch to epoch, with emission features coming and going, probably because of natural variation in the intensities of individual emission features, as well as interstellar scintillation. To obtain the most complete assessment of the extent of dense molecular structure that underlies the maser emission in Circinus, we superposed the maps for two well separated epochs, 1997 July and 1998 June (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). To register the two maps, we aligned the positions of the common emission peak at 565.2 km s-1. A comparison of relative positions for other features that are present in both maps suggests the registration uncertainty is $`0.15`$ mas.

## 5 Discussion

### 5.1 The warped disk

We interpret the distribution of maser emission in Circinus in the context of the warped accretion disk in NGC 4258, which has been well modeled (e.g., Herrnstein, Greenhill, & Moran, 1996). In general, maser emission in warped disks is beamed toward the observer from regions where the gradient in line-of-sight velocity along the line of sight is close to zero and the disk is approximately tangent to the line of sight. This corresponds to two loci, the “midline,” where the orbital velocity is parallel to the line of sight, and an arc along the near side of the disk, where the orbital velocity is nearly perpendicular to the line of sight. The former is responsible for high-velocity emission that traces the rotation curve of the disk, and the latter is associated with low-velocity emission. The detailed spatial and velocity distribution of observable high- and low-velocity emission is also affected by the geometry of the warp, which (1) determines the solid angle into which emission is beamed, (2) shadows parts of the disk from X-ray emission by the central engine, which otherwise heats gas and probably supports maser emission (Neufeld, Maloney, & Conger, 1994; Neufeld & Maloney, 1995; Collison & Watson, 1995; Wallin & Watson, 1997), and (3) determines how the disk might overlie and thereby amplify background nonthermal continuum emission, such as from a jet (e.g., Herrnstein et al.,, 1997).

A warped disk model for Circinus (see Table 3) receives its strongest support from the following: (1) the elongated, antisymmetric angular distribution of maser spots that forms the shallow S in Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy, (2) the symmetric bracketing of the nominal systemic velocity by the most highly red and blueshifted emission, (3) the antisymmetric distribution of Doppler shifts along the S with velocity declining by roughly $`b^{0.5}`$, where $`b`$ is impact parameter (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy), and (4) the orientation of the inner disk (as seen in projection), more or less perpendicular to the axis of the known radio lobes (Elmouttie et al.,, 1998) and ionization cone (Marconi et al.,, 1994; Veilleux & Bland-Hawthorn, 1997; Wilson et al.,, 2000).

To establish which maser spots arise in the disk, we fit a second order polynomial to the sky position of the midline, which is identifiable because that is where the greatest Doppler shifts are observed (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy), and we calculated the angular offset between each maser spot and the midline. We classified emission that is offset by $`<1`$ mas (about one half beamwidth, north-south, in 1998 June) as disk emission (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). About 43% of the integrated emission we detected in 1998 June arises in the disk, and on average, maser spots in and out of the disk are comparably bright. The non-disk emission alone occupies the center $`260`$ km s-1 of the spectrum, while the disk emission lies on either side in velocity (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). There is some overlap in velocity, the extent of which depends on the chosen cutoff to the angular offset from the midline (i.e., here, 1 mas). This overlap may be a consequence of resolved disk thickness, azimuthal structure and a not quite edge-on orientation, or interaction between the disk and the outflow whereby the fastest spots (in projection) within the outflow lie closest to the disk.

The position of the dynamical center of the disk, its inner (observed) radius, and peak rotation speed may be estimated from the innermost red and blueshifted masers, for which we assume a common edge-on orbit. The peak observed rotation speed, $`V_{max}`$, is 260 km s-1 in an orbit at a position angle of $`29^{}`$ with radius, $`r_{\mathrm{in}}`$, of $`5.3\pm 0.1`$ mas ($`0.11\pm 0.02`$ pc, including the uncertainty in distance). The inferred disk center lies at 17.6 and 13.1 mas east and north of the origin in the maps in Figures A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy and A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy. Taking into account only gravity and assuming circular motion, we infer an enclosed mass, $`M(r<r_{in})`$, of $`1.7\pm 0.3\times 10^6`$ M or mean mass density of $`3.2\pm 0.9\times 10^8`$ M pc-3. We note that because the disk does not display identifiable low-velocity emission on the near side, it is difficult to estimate inclination directly (cf. NGC 4258; Miyoshi et al., 1995). Formally, the mass and density estimates depend on the inclination and are lower limits. However, we argue that the disk must be close to edge-on because the masers cluster tightly around a sharply defined disk midline.

The orientation of the innermost orbit on the disk is close to perpendicular to the flow axes of the known bipolar radio lobes (Elmouttie et al.,, 1998) and wide-angle CO outflow cone (Curran et al.,, 1999), the rotation axes of the H I galactic disk (Freeman et al.,, 1977) and the circumnuclear CO disk of radius 140 to 600 pc (Curran et al.,, 1998), and the principle axis of the most prominent \[O III\] filament in the ionization cone, for which the position angle is $`50^{}`$ (Veilleux & Bland-Hawthorn, 1997). However, the assumed $`90^{}`$ inclination of the disk at the inner radius is somewhat offset from the $`65^{}\pm 2^{}`$ inclination of the galactic disk (Freeman et al.,, 1977), the $`78^{}\pm 1^{}`$ inclination of the CO disk (Curran et al.,, 1998), and the $`65^{}`$ inclination inferred by Wilson et al., (2000) for a 40 pc-radius ring of H II regions centered on the central engine. Some of these alignments are sensible, such as the position angles of the inner accretion disk and the radio lobe axis. However, others may include an element of coincidence, such as the alignment between the circumnuclear CO disk and the accretion disk, for which agreement with the inner orbital axis is better than agreement with the outer orbital axis (Table 3). Overall, we suggest that the orientation of the accretion disk should be only weakly coupled via gravity to the surrounding large-scale dynamical structures, because the central engine contributes only $`0.5\%`$ of the $`3.2\times 10^8`$ M inside 140 pc and has a sphere of influence less than a few parsecs in radius (Curran et al.,, 1998).

The inferred peak rotation velocity as a function of impact parameter, $`b`$, measured from the estimated dynamical center, declines approximately as $`b^{0.5}`$ (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy), which is suggestive of a relatively low-mass disk, as is observed in NGC 4258 (Miyoshi et al.,, 1995). However, a disk that is both massive and warped in inclination angle could have a similarly steep rotation curve. Specifically, if a disk is edge-on at radius $`r_{}`$ but inclined by $`\delta i`$ at a radius $`r+\delta r`$, then it will display the same apparent rotation speed as a flat, edge-on, Keplerian disk if the fractional disk mass $`M_{\mathrm{disk}}/M`$ is $`\mathrm{tan}^2(\delta i)`$, where $`M`$ is the encircled mass at radius $`r_{}`$. In order to place a limit on the disk mass for Circinus, we assume that the disk is warped both in position and inclination angles, and that the warps are of the same order of magnitude, $`27^{}`$ between 0.11 and 0.40 pc (see Table 3). Hence, our data are consistent with a fractional mass of 26%, or $`4\times 10^5`$ M. In the absence of an actual measurement of the inclination warp, we adopt this mass estimate as an upper limit.

The relatively high particle densities inferred from the presence of H2O maser emission — $`10^8`$ to $`10^{10}`$ cm-3 for a H2O fractional abundance of $`10^5`$ to $`10^4`$ (Elitzur, 1992) — are also consistent with a relatively massive disk. For a uniform particle density and a planar disk, the disk mass enclosed inside radius $`r`$ is $`1.4\times 10^6n_{10}\gamma r_1^3`$ M, where $`r_1`$ is the radius in units of 0.1 pc, $`n_{10}`$ is density in units of $`10^{10}`$ cm-3, and $`\gamma `$ is the height of the disk as a fraction of its radius. Hydrostatic equilibrium dictates a lower limit on $`\gamma `$ that is on the order of $`c_s/v_{\mathrm{rot}}`$, where $`c_s`$ is the sound speed ($`1.6`$ km s-1 for gas at 400 K). Our upper limit on disk mass inside 0.40 pc ($`4\times 10^5`$ M), is consistent with a $`4\times 10^9`$ cm-3 upper limit on density, for $`v_{\mathrm{rot}}140`$ km s-1 and hydrostatic equilibrium, which agrees with expectations.

Estimates of the Toomre Q-parameter for gas densities sufficiently high to support maser action ($`>10^8`$ cm-3) indicate that self-gravity is important and probably causes clumping of disk material. We note that the distributions of disk masers on the sky (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy) and in the position-velocity diagram (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy) exhibit substantial substructure, which may be a signature of clumping. This clumping could limit the maser gain paths in the disk and thus be responsible for the comparable strengths of the disk and non-disk emission. The Toomre $`Q`$-parameter is $`c_s\mathrm{\Omega }/\pi G\mathrm{\Sigma }`$, where $`\mathrm{\Omega }`$ is the angular rotation speed, and $`\mathrm{\Sigma }`$ is the mass surface density. Instabilities arise for $`Q<1`$. By assuming the disk thickness is $`\eta `$ times the hydrostatic limit, we obtain $`Q0.6M_6/\eta n_{10}r_1^3`$, where $`M_6`$ is the central engine mass in units of $`10^6`$ M. For $`n=4\times 10^9`$ cm-3, $`r=0.40`$ pc, and hydrostatic equilibrium ($`\eta =1`$), $`Q0.04`$. The Q-parameter could be even smaller if the apparent thickness of the disk (which is on the order of a few tenths of a mas) is real and thereby much greater than the hydrostatic limit (which is about tens times smaller). However, it is difficult to distinguish between measurable disk thickness and nonaxisymmetric structure (e.g., Maoz & McKee, 1998), because the disk is highly inclined and the angular resolution of our observations is an order of magnitude larger than the putative thickness.

### 5.2 The clumpy wide-angle outflow

The association of water maser emission with a bipolar wide-angle wind in an AGN is qualitatively new and distinct from the reported connection to jet activity in some galaxies \[i.e., NGC 1068 (Gallimore et al.,, 1996), NGC 1052 (Claussen et al.,, 1998), NGC 3079 (Trotter et al.,, 1998), and Mrk 348 (Peck et al.,, 2001)\]. Less than 1 pc from a central engine, outflowing material is often assumed to be photoionized, but away from the flow axis, gas columns may be large enough ($`10^{22\pm 1}`$ cm-2) to maintain a molecular component with internal pressures ($`p/k`$) that are $`10^{11}`$ to $`10^{13}`$ K cm-3 (Emmering, Blandford, & Shlosman, 1992; Königl & Kartje, 1994; Kartje, Königl, & Elitzur, 1999). Maser emission can be readily stimulated in this gas by X-ray heating (Neufeld, Maloney, & Conger, 1994; Neufeld & Maloney, 1995) or shock heating (Elitzur, Hollenbach, & McKee, 1989; Kaufman & Neufeld, 1996), if turbulence and gradients in line-of-sight velocity permit.

In Circinus, blueshifted “outflow masers” lie south and east of the central engine, and redshifted counterparts lie chiefly to the west (Figures A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy, A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy, & A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). The line-of-sight velocities of these masers are offset from the systemic velocity by up to $`\pm 160`$ km s-1 (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy), though we note a few maser spots that are somewhat more redshifted and that lie almost due west of the dynamical center in 1997 July. To start, we infer that this emission originates outside the disk because (1) it lies up to $`1`$ pc from the disk center, while the disk midline can be traced only between radii of 0.11 to $`0.40`$ pc, (2) many of the non-disk masers lie close to the disk rotation axis, as seen in projection, (3) the line-of-sight velocities of the non-disk masers that instead lie close to the disk differ by tens of km s-1 from the velocities of adjacent disk material, and (4) the non-disk masers display no systematic gradient in line-of-sight velocity as a function of position angle relative to the disk center (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy), as would be expected were they to lie in a disk sufficiently warped to include them. Having excluded a disk origin, we infer that the masers in question are associated with outflow because it is unlikely that material close to the axis of an accretion disk, and also within 1 pc of a supermassive blackhole would be infalling. (We note that it is also unlikely that the non-disk masers are associated with a relaxed central cluster of mass losing stars, principally because the distribution of masers is not at all centered on the inferred position of the central engine.)

We adopt a bipolar, cone-like model geometry to describe the outflow (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). Emission with different Doppler shifts is well mixed on the sky (e.g., Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). This is suggestive of relatively short gain paths peppering the outflow volume, where emitting regions that move along very different trajectories lie along proximate lines of sight (i.e., the masers probably do not just lie along the limbs of the outflow). Within this volume, the masers presumably do lie toward the front because of the opacity of the photoionized gas that probably lies between the maser-emitting regions (see below).

Although a compelling general case may be made for outflow, we cannot explain the observed angular distributions and Doppler shifts of the outflow masers in detail. The chief hurdles are (1) the prevalence of blueshifted emission to the south and east of the disk center and redshifted emission on the reverse side, and (2) the difference in angular distributions on opposite sides of the accretion disk, where some areas are devoid of visible maser emission while their counterparts (mirrored on the other side of the disk) are not (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). The blue-red asymmetry in the outflow could indicate that the flow is inclined to the line of sight, despite the putative edge-on orientation of the accretion disk. However, it is difficult to test this hypothesis because the redshifted outflow is poorly sampled on the sky and almost entirely defined by maser emission at low latitudes (i.e., in close proximity to the disk). Gas dynamic or magnetic coupling between outflowing and rotating disk material could be responsible for the observed Doppler shifts at low latitudes, and in position-velocity space, for some portion of the structure we observe near (the representation of) the disk midline (see Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). Regarding the second hurdle, the difference in maser distribution above and below the disk could be a consequence of differential extinction. However, if the axis of the outflow is close to the plane of the sky, then this effect would be small. Alternatively, the angular distributions could depend chiefly on where velocity coherent maser gain paths happen to exist. In general, the longest paths are probably associated with limb-brightening and motions that are closer to the plane of the sky than to the line of sight (e.g., Elitzur, Hollenbach, & McKee, 1992). This observation may answer the question, why do we not see highly blueshifted outflow material close to the line of sight toward the central engine. However, it is not clear whether it can explain the absence of a broad distribution on the sky of redshifted outflow masers.

Because water masers require high gas densities, their presence broadly distributed in the outflow argues that it must be inhomogenous or “clumpy.” Were the outflow smooth, the associated mass loss would be prohibitively large, $`6\times 10^3r_1^2n_{10}v_2`$ M yr-1, where $`v_2`$ is velocity in units of 100 km s-1. For a 100 km s-1 flow that is $`10^9`$ cm-3 at a radius of 0.1 pc, the flow would consume a mass equivalent to the central engine and accretion disk in only a few $`\times 10^3`$ years. However, the appearance of the flow, as traced by maser emission, is at best sparse and irregular. If this is representative and not merely a consequence of anisotropic beaming of maser radiation, then the mass loss would be orders of magnitude smaller. For instance, if the rate of loss balances the rate of accretion, then the mean density of the outflow at 0.1 pc radius is $`10^4`$ cm-3, where we express the accretion rate as $`1.7\times 10^4L_{43}ϵ^1`$, and where $`L_{43}`$ is the bolometric luminosity in units of $`10^{43}`$ ergs s-1 and $`ϵ`$ is the accretion efficiency of the central engine, probably $`0.1`$. The apparent mean clump-interclump density contrast is on the order of $`10^5`$, and the shielding column in the interclump medium is on the order of $`10^{22}`$ cm-2 for a path length comparable to the disk radius, as required (Kartje, Königl, & Elitzur, 1999). (We note that the flow could instead contain pockets of molecular gas with a range of sizes and densities that span orders of magnitude, only some of which would be conducive to observable maser action. This would affect the estimates of density contrast and shielding column, but it would not reduce the imperative that the flow be inhomogeneous.)

### 5.3 Collimation by the warped disk

The angular distribution of the outflow masers is consistent with the hypothesis that excitation depends on a direct line of sight to the central engine. In the shadow of the disk (i.e., in regions where this path is blocked) no maser emission was detected (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy). We suggest that this shadowing is a physical blockage of both radiation and mechanical energy through which the warped accretion disk collimates the nuclear outflow. Channeling on such small scales is consistent with the inference by Wilson et al., (2000) based on HST imaging, that the Circinus ionization cone is collimated within 2 pc of the central engine. In the proposed model, on each side of the maser disk, one radial edge of the outflow is determined by the position angle of the disk at the inner radius, and the other edge is fixed by the concave surface of the disk warp along which the outflow brushes, up to the outer radius of the disk.

We estimate the outer radius to be $`0.40`$ pc based on the close packing of disk masers along the midline up to this radius, the angular distribution of outflow masers at larger radii, and the extent of the ionization cone. First, the midline of the warped disk is traced in the maps (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy) and the position-velocity diagram (Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy) by a relatively continuous locus of maser spots up to a radius of $`0.40`$ pc. Beyond 0.40 pc, an extrapolation of the midline does pass close to some other masers, but they do not lie on the Keplerian rotation curve. Second, the only way in which the southwesternmost maser clump near position (-16, -8) mas in Figure A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy can achieve a direct line of sight to the central engine is if the disk is truncated at $`0.40`$ pc. Third, the southern edge of the ionization cone, observed in \[O III\] light (Marconi et al.,, 1994; Veilleux & Bland-Hawthorn, 1997; Wilson et al.,, 2000), lies at a position angle of $`120^{}`$, which is close to the position angle of the disk midline at a radius of $`0.40`$ pc. Were the warped disk to extend substantially beyond this proposed outer radius, the disk would block a larger solid angle of radiation and create a different southern limit for the ionization cone. We speculate that an outer radius of $`0.40`$ pc may also be associated with the presence of a maser spot cluster surrounding the outer tip of the redshifted midline, forming a backward “C” visible near the origin in Figures A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy and A Warped Accretion Disk and Wide Angle Outflow in the Inner Parsec of the Circinus Galaxy. These masers move with sub-Keplerian velocities and could correspond to material lost from the outer disk where it is roughened by interaction with outflow and radiation.

The putative outer radius of the accretion disk constrains the longevity of activity in Circinus at the observed luminosity, $`4\times 10^{43}`$ ergs s-1. The wide-angle outflow restricts the solid angle (as seen from the central engine) over which material can fall in and replenish the accretion disk. The outflow blocks 80% to 90% of the surroundings. If the supply of new material is thereby reduced, then for the accretion rate ($`1.7\times 10^4L_{43}ϵ^1`$) that corresponds to an efficiency of 10%, the accretion disk will be exhausted in $`6`$ Myr for the previously estimated maximum disk mass of $`4\times 10^5`$ M.

### 5.4 The Origin of the Outflow

Two possible origins for the molecular outflow are: (1) radial ejection of ionized broad line clouds that form dust and molecules as they cool, and (2) uplift and ejection of warm molecular disk material. The association of masers and cooled broad line cloud material is largely circumstantial. It rests on the observation that the densities and temperatures required for maser emission are consistent with adiabatic expansion of material with intial densities and temperatures on the order of $`10^{11}`$ cm-3 and $`10^4`$ K (e.g., Kaspi & Netzer, 1999; Korista & Goad, 2001, and references therein), assuming an adiabatic constant of $`\frac{5}{3}`$. However, it is not clear whether maser action is possible in freely expanding clouds, where velocity gradients and turbulent motions may be significant, or in the shocks that may arise within or ahead of the clouds. The observed line-of-sight outflow velocities, which are smaller than the Keplerian rotation velocity and thereby the escape velocity at all radii, are a secondary concern. Because the motions of the observable masers may lie somewhat preferentially toward the plane of the sky (see Section 5.2), the too small velocities may be a projection effect.

In order to describe self-consistently the geometry and dynamics of broad line regions, Emmering, Blandford, & Shlosman (1992) modeled a hydromagnetic wind created by the uplift and ionization of clumps from a cool molecular accretion disk. The magnetic field, anchored in the disk, confined the ionized clouds and imparted rotation to the outflow. Over time, the clouds moved outwards with increasingly radial trajectories as the azimuthal field component or cloud coupling to the field weakened (see, e.g., Bottorff et al.,, 1997). Königl & Kartje (1994) expanded this model by showing that the hydromagnetic disk wind could be clumpy, dusty, and molecular over a broad range of latitudes. They suggested that such a wind might be responsible in general for the obscuration observed toward type-2 AGN, rather than oft posited geometrically thick tori. Kartje, Königl, & Elitzur (1999) applied the model to the stimulation of maser emission in the uplifted clumps once they rise from the disk to be irradiated by the central engine. (We note that in this picture, the disk is not directly the source of maser radiation.)

The hydromagnetic wind model is somewhat more attractive than the simple ejection model for broad line clouds principally because in the former the clumps are magnetically confined. They do not need to expand by orders of magnitude and remain coherent to achieve temperatures and densities conducive to maser action. However, we note that in the limit where outflowing material moves largely radially well after uplift, the models are similar. The simple ejection model as stated finesses the origin of the broad line clouds, but if the clouds originate in an accretion disk, then the models differ chiefly in whether the disk material remains neutral after it is removed from the disk.

However, the hydromagnetic wind model, as applied by Kartje, Königl, & Elitzur (1999) makes three predictions that are not easily reconciled with the observations. First, the masers should lie away from the rotation axis of the disk. In contrast, we observe outflow masers that are apparently quite close to the rotation axis. If this material is removed from the observed accretion disk (between radii of $`0.1`$ and $`0.4`$ pc), then it would have to travel almost vertically as it is uplifted. Second, masers in a hydromagnetic wind should exhibit a sub-Keplerian rotation curve (i.e., in a position-velocity diagram the masers would at least in part lie above the disk rotation curve). In fact, the signature of rotation is weak if present at all, and for any given radius, the Doppler shifts of outflow material are smaller than for the disk. Third, the wind should support maser emission close to the loci where dust sublimation occurs, resulting in a more or less flattened X-pattern on the sky for relatively edge-on orientations. However, the actual distribution of masers is dichotomous, featuring a sharply defined, thin, S-shaped disk and a broad outflow component.

We propose a modified wind model. First, the outflow begins well inside the observed inner radius of the disk and thereby can reach high latitudes at radii of 0.1 to 0.4 pc more readily. Very high velocity maser emission, offset up to $`400`$ km s-1 from the systemic velocity, has recently been discovered in Circinus (Greenhill et al.,, 2003). If the emission originates in the accretion disk, it provides evidence for a resevoir of molecular material $`0.03`$ pc from the central engine. Second, the uplifted material moves by and large radially, perhaps because it decouples quickly from the magnetic field threading the accretion disk. Third, the outflow is sparsely populated with clumps of material, so that the warped accretion disk at larger radii is not significantly shadowed and can emit maser radiation, as is observed.

## 6 Conclusions and Comment

We find that the H2O maser emission in the Circinus Galaxy traces a warped, edge-on accretion disk at radii between $`0.11\pm 0.02`$ pc and $`0.40`$ pc. The peak observed disk rotation velocity is $`260`$ km s-1, and the encircled mass is $`1.7\pm 0.3\times 10^6\mathrm{M}_{}`$. The implied ratio of bolometric luminosity to Eddington luminosity is $`0.2`$. The rotation curve is nearly Keplerian and the disk mass is probably $`<24\%`$ of the central mass. Given the (low) rotation speed of the disk, this is substantial enough to cause clumping due to self-gravity, which may explain observed substructure in the disk. The inner radius of the maser disk is comparable to the well resolved maser disk in NGC 4258 (0.14 pc), even though the X-ray luminosity of Circinus is $`10`$ times larger (Makishima et al.,, 1994; Reynolds et al.,, 2000). If the recently reported very high velocity emission in Circinus originates in the disk (Greenhill et al.,, 2003), then the paradox is more striking still. However, it is possible the inner radii of maser emission are not dictated by dust sublimation. The discrepancy may be due to the geometries of the warped disks. For example, in NGC 4258 maser emission at particularly small radii may not be directed toward us.

A second population of masers in the Circinus AGN traces a wide-angle outflow up to $`1`$ pc from the central engine with line-of-sight velocities of up to $`\pm 160`$ km s-1 with respect to the systemic velocity. These masers are the first direct evidence (i.e., provided by imaging) of dusty, high-density, molecular material in a nuclear outflow, at such small scales. The mechanism that drives the wind is uncertain but may be hydromagnetic. The wind is observed in the regions around the disk that are not shadowed by the warp, from which we suggest that the disk channels the outflow and occults radiation from the central engine. The position angles of the edges of the outflow correspond well to those of the outflow and ionization cone observed on scales of hundreds of parsecs (Veilleux & Bland-Hawthorn, 1997; Curran et al.,, 1999). Here too, Circinus may provide the first direct evidence for collimation of an AGN outflow and illumination pattern by a warped accretion disk. In the classical paradigm for AGN, a geometrically thick tori has been cast in this role, even though it has been difficult to explain how such structures can be supported vertically (e.g., Krolik & Begelman, 1988). For Circinus and perhaps other type-2 AGN, sub-parsec, warped accretion disks are a simpler explanation of a breadth of observations.

It is an open question whether thick tori are required at all to explain the X-ray obscuring columns observed in type-2 AGN. Warped disks of the type found in NGC 4258 and Circinus could easily supply the needed column densities, but they are relatively flat. These two disks subtend on the order of 10 to 20% of the sky as seen from their respective central engines and do not readily explain the high proportion of obscured AGN out of the total (e.g., Miller & Goodrich, 1990). As suggested by Königl & Kartje (1994), dusty hydromagnetic winds may be a general source of obscuration. Alternatively, structures on the small spatial scales of broad line regions (Risaliti et al.,, 2002) or the $`>100`$ pc scales of warped CO gas disks (Schinnerer et al.,, 2000) and other dusty galactic structures (Malkan et al.,, 1998) may be important.

The Circinus H2O maser is the third one for which detailed VLBI study has detected disk structure and a declining rotation curve. As has been demonstrated for the NGC 4258 and NGC 1068 masers, the Circinus maser makes possible unusually detailed study of the structure and dynamics of a region $`<1`$ pc from what is most probably a supermassive black hole. We anticipate that higher angular resolution VLBI observations covering the recently recognized, full $`900`$ km s-1 breadth of the maser will improve our understanding of the structure of the $`0.1`$ pc radius warped accretion disk and the mechanism that drives the wide-angle outflow apparently emanating from still closer to the black hole.

We thank W. Wilson and R. Ferris for their support of the S2 VLBI correlator, and H. May and S. Amy for network and software support. We acknowledge the helpful comments of G. Ball, M. Elitzur, A. Konigl, and an anonymous referee. We are grateful for the dedication of the staff at the LBA antennas. We acknowledge the use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Figure Captions


format: text

text:
# Fluctuating Filaments Under Tension - From Flexible Chains to Rigid Rods

## I Introduction

Recent experimental studies of the mechanical properties of biological materials, ranging from individual DNA moleculesbustamante  to single actin filaments and their networksalbert ; fred ; sackmann , led to revival of interest in the worm-like chain model of semiflexible biopolymers. In particular, the force-extension curves predicted by this model in the long chain limit (i.e., when the length of the chain is much larger than the persistence length) were shown to agree quantitatively with experimental results on dsDNA at large elongationsms1 . However, many interesting questions concerning worm-like chains under tension are still open and there are delicate issues concerning even the linear response of the worm-like chain to stretching forcefrey2 . Also, while the equilibrium distribution function of the extension was calculated for the worm-like chainfrey1 ; sinha ; dhar  and fluctuating ribbonjcp-kr  models, only approximate (using the Gaussian approximation) calculations of the distribution function for stretched chains were reported to datewinkler . Such distributions have been measured in experiments on stretched DNA moleculesthesis  and provide a much more stringent test of the theory than the force-extension curves. Finally, even though all mechanical experiments involve biopolymers attached to surfaces, none of the present theories incorporates wall effects on force-extension curves and on the corresponding distributions. While such effects are negligible for nearly fully stretched chains, they are expected to be important in the weak force regime.

In previous workKats ; jcp-kr , we have discussed how to generate via Monte-Carlo equilibrium ensembles of fluctuating ribbons with spontaneous curvature and/or twist. This method was based on using the intrinsic geometry of the filament, in which representation the energy is purely quadratic. This enables a direct generation of independent samples distributed according to the Gibbs equilibrium ensemble. When external forces act on the filaments, the energy is no longer a quadratic function of the curvatures, and the problem is now not only nonlinear but also nonlocal. In this paper we discuss a biasing method which succeeds in efficiently generating Monte-Carlo samples for the simplest case of an external force, namely that of the stretched filament. This problem was studied by Marko and Siggia, (MS), ms1 ; MS94  and others odijk ; ha  who calculated the force-extension curve for long worm-like chains. We will reproduce these results as well as extend the analysis to the case of nearly stiff filaments with length comparable to or shorter than the persistence length, and in addition calculate the distribution of extensions for various forces. Lastly, we will show how we can also easily incorporate steric constraints; we will demonstrate this in the context of a restriction of the filament to a half-space by the presence of a wall.

In sections 2 and 3 we review the general theory of fluctuating ribbons (we define a ribbon as a filament with an asymmetric cross–section) and the Frenet algorithm for computer simulations of such ribbons, respectively. In section 4 we describe the biased Monte-Carlo simulations of filaments under tension and present the resulting force-extension curves and distribution functions. Wall effects are considered in section 5. In section 6 we discuss our results and outline directions for future research. In Appendix A we expand the Euler angles in terms of the principal curvatures and then use this expansion in Appendix B in order to calculate the partition function of a stretched filament. In part 1 of this appendix we derive an analytic expression for the partition function and obtain a force-extension relation for arbitrary force in the limit of large persistence length (expansion about the rigid rod limit). In part 2 of the appendix we derive asymptotic expressions for the above quantities in the limit of large force (for arbitrary rigidity). Finally, in Appendix 3 we present the calculation of the renormalization of the elastic moduli of the various Fourier modes of the curvature.

## II General theory

The linear elasticity of ribbons is described by the elastic energyHelixpre

$$E_{el}=\frac{1}{2}\underset{k=1}{\overset{3}{}}b_k_0^L𝑑s\left(\omega _k\omega _{0k}\right)^2,$$
(1)
where $`s`$ is the distance measured along the contour of the ribbon, the coefficients $`b_1`$ and $`b_2`$ are the moduli associated with bending along the two principal axes of inertia of the cross section (they differ if the cross section is not circular), and $`b_3`$ is the twist modulus. In this paper we treat $`\left\{b_i\right\}`$ as given material parameters of the ribbon. The functions $`\left\{\omega _k(s)\right\}`$ and $`\left\{\omega _{0k}(s)\right\}`$ are the generalized curvatures in the deformed and the stress-free states of the ribbon, respectively. These curvatures completely determine the three dimensional conformation of the ribbon through the generalized Frenet equationsHelixpre

$$\frac{d𝐭_i(s)}{ds}=\underset{j,k}{}\epsilon _{ijk}\omega _j(s)𝐭_k(s).$$
(2)
Here $`𝐭_3`$ is the unit tangent to the centerline and the unit vectors $`𝐭_1`$ and $`𝐭_2`$ are oriented along the principal axes of the cross section ($`\epsilon _{ijk}`$ is the antisymmetric unit tensor). Note that since these equations describe a pure rotation of the triad $`\left\{𝐭_i(s)\right\}`$ as one moves along the contour of the ribbon, the constraint $`\left|𝐭_3\right|=\left|d𝐫(s)/ds\right|=1`$ is automatically satisfied in this intrinsic coordinate description.

The physical frame of the ribbon $`\left\{𝐭_i(s)\right\}`$ can be related to the original Frenet description of space curves in terms of a unit tangent (which coincides with $`𝐭_3`$), normal $`𝐧`$ and binormal $`𝐛`$, provided that in addition to the Frenet triad $`\{𝐛(s),𝐧(s),𝐭_3(s)\}`$ associated with the centerline of the ribbon, one introduces an angle $`\alpha (s)`$ that describes the rotation of two of the principal axes of the cross section $`𝐭_1`$, $`𝐭_2`$ with respect to the binormal $`𝐛`$ and normal $`𝐭`$ (see Fig. 1 in the second of Refs. \[Helixpre,  \] ),

$$𝐭_1=𝐛\mathrm{cos}\alpha +𝐧\mathrm{sin}\alpha ,𝐭_2=𝐛\mathrm{sin}\alpha +𝐧\mathrm{cos}\alpha .$$
(3)
The three dimensional configuration of the centerline is determined by the Frenet equations

$$\frac{d𝐛}{ds}=\tau 𝐧,\frac{d𝐧}{ds}=\kappa 𝐭_3\tau 𝐛,\frac{d𝐭_3}{ds}=\kappa 𝐧$$
(4)
where the curvature $`\kappa (s)`$ and torsion $`\tau (s)`$ are treated as known functions of position along the contour. The three generalized curvatures $`\left\{\omega _k(s)\right\},`$ can be expressed in terms of these parameters and the angle $`\alpha (s),`$

$$\omega _1=\kappa \mathrm{cos}\alpha ,\omega _2=\kappa \mathrm{sin}\alpha ,\omega _3=\tau +d\alpha /ds.$$
(5)
## III The Frenet algorithm

Since the energy is a quadratic form in the deviations $`\delta \omega _k`$ of the generalized curvatures from their values in the stress free state, Eq. (1) is valid as long as the characteristic length scale of the deformation is much larger than the diameter of the ribbonLove . When a ribbon undergoes fluctuations in the presence of a thermal bath at temperature $`T,`$ this energy determines the statistical weight of the configuration $`\left\{\omega _k\right\}`$. The statistical average of any functional of the configuration $`A(\left\{\omega _k\right\})`$ is defined as the functional integralHelixpre

$$A\left(\left\{\omega _i\right\}\right)=\frac{D\left\{\delta \omega _i\right\}A\left(\left\{\omega _i\right\}\right)e^{E_{el}\left\{\delta \omega _i\right\}/kT}}{D\left\{\delta \omega _i\right\}e^{E_{el}\left\{\delta \omega _i\right\}/kT}}.$$
(6)
where $`k`$ is the Boltzmann constant. From the form of the energy, Eq. 1, it follows that the distribution of $`\delta \omega _k`$ is Gaussian, with zero mean ($`\delta \omega _i(s)=0`$) and second moments given by

$$\delta \omega _i(s)\delta \omega _j(s^{})=\frac{kT}{b_i}\delta _{ij}\delta (ss^{}).$$
(7)
In order to study the distribution functions of the various fluctuating quantities (e.g., the extension $`R`$), we developed what we entitled the Frenet algorithm, which is based on the combination of statistical methods and differential geometry and is described in the following. Examination of Eq. (7) shows that the decoupling of the “noises” $`\left\{\delta \omega _k\right\}`$ in the intrinsic coordinate representation permits an efficient numerical generation of independent samples drawn from the exact canonical distribution. The Gaussian distribution of the $`\delta \omega _k`$ means that each $`\delta \omega _k(s)`$ can be directly generated as one of a string of independent random numbers drawn from a distribution symmetric about the origin with width $`\sqrt{kT/(b_k\mathrm{\Delta }s)}`$ where $`\mathrm{\Delta }s`$ is the discretization step length. The remaining task is to construct the curve using the Frenet equations with $`\omega _k(s)=\omega _{0k}(s)+\delta \omega _k(s)`$. The Frenet equations are best integrated by stepping the basic triad $`t_k`$ forward in $`s`$ through a suitable small rotation. In this way, the orthonormality of the triad is guaranteed to be preserved up to machine accuracy. To construct this rotation matrix, we begin with Eqs. (2) and, defining the three vectors $`v^x=(t_1^x,t_2^x,t_3^x)`$, and so forth, we can write this equation as

$$\frac{\mathrm{\Delta }v^i}{\mathrm{\Delta }s}=Av^i,$$
(8)
where $`A`$ is an antisymmetric matrix with elements $`A_{ij}=_k\epsilon _{ijk}\omega _k`$. We now discretize Eq. (8) as

$$v^i(s+\mathrm{\Delta }s)=Ov^i(s),$$
(9)
where the orthogonal matrix $`O`$ is

$$O=\left(1+\frac{\mathrm{\Delta }s}{2}A\right)\left(1\frac{\mathrm{\Delta }s}{2}A\right)^1.$$
(10)
The choice of discretization step length ($`L/\mathrm{\Delta }s`$ is the number of points) depends on the amplitude of the random noise, $`\sqrt{kT/(b_k\mathrm{\Delta }s)}.`$ We found that $`400`$ points gave sufficient accuracy (i.e., was limited by the statistical noise determined by the number of independent samples) for the calculation of the distribution function of the extension, when the smallest of the moduli satisfied $`kTL/b_k^{\mathrm{min}}1.`$

This algorithm is optimal for studying the fluctuations of a free ribbon. However, if the ribbon is acted upon by external forces and/or constraints, the algorithm must be modified. This is described in the next section.

## IV Biasing the Monte-Carlo

Consider a filament with a circular cross–section ($`b_1=b_2`$) which is characterized by a bending persistence length

$$a=\frac{b_1}{kT}.$$
(11)
The simplest example of a filament under tension is that of a constant force acting on one end and directed along the $`x`$-axis, with the second end constrained to lie at the origin but otherwise free to rotate about it. The energy of the filament thus picks up an additional contribution $`Fx`$, where $`x`$ is the projection of the end to end vector on the direction of the force. The problem is that $`x`$ is a complicated nonlinear function of all the $`\delta \omega `$, and so the external force destroys the Gaussian structure of the energy on which the generation of the $`\delta \omega `$ depends. One can simply incorporate the extra contribution to the energy by treating it as a weighting factor, or bias, and suitably normalizing all averages. Thus, for example, one can calculate the mean value of $`x`$ via

$$\overline{x}=\frac{xe^{fx}}{e^{fx}}$$
(12)
where the sum is over configurations generated via the Frenet algorithm described above and $`fF/kT`$ is the force in units of inverse length. For small forces, this works well. However, as the force increases, the probability of generating highly extended configurations which are heavily favored by the bias, is so small that the sampling becomes ineffective. One therefore needs a way to bias the generation of configurations towards those that are extended, so they can be sampled appropriately.

This task is not simple, since the connection between the $`\delta \omega `$’s through which the filament is parameterized and the value of $`x`$ is very complicated. We can gain insight into this by transforming to Fourier space and looking at the expectation value of the moments of the $`n`$’th Fourier component of the curvature, $`\delta \stackrel{~}{\omega }(n)`$, using the above algorithm, Eq. 12, for small to intermediate forces where it is still effective. For an extensional force, we do not expect spontaneous curvature to develop, and indeed the first moments of the Fourier components all vanish. The second moments, however, do exhibit an interesting structure, as seen in Fig. 1, where the expectation value of the squared amplitude of the low $`n`$ Fourier modes are presented for the case $`a=3`$, $`L=2\pi `$. For $`F=0`$, the second moments are independent of $`n`$, as expected. When the filament is subjected to a tensile force we see that the small $`n`$ modes are suppressed, with the degree of suppression increasing with $`F`$ and decreasing with $`n`$. In the case of large $`a`$, where the fluctuations are in any case small, we observe that only the $`n=0`$ mode is significantly affected by the force (not shown). The suppression of the various curvature modes can be calculated analytically when fluctuations are small, namely either in the large $`FL`$ limit, or for any $`F`$ provided $`\chi \sqrt{FL^2/(4akT)}1`$ (see Appendix). The calculation confirms the above-mentioned trends. We learn from this that the high Fourier modes (or equivalently the small-scale fluctuations) of the filament are essentially unaffected by the force, whereas the long-wavelength fluctuations are significantly damped.

This suggests the possibility of biasing the Monte-Carlo by selectively dampening the small-$`n`$ modes. The idea is to generate the $`\delta \stackrel{~}{\omega }`$’s in Fourier space, with appropriate amplitudes for the various modes, and construct the real-space $`\delta \omega (s)`$ via a Fast Fourier Transform. Of course, this bias has to be corrected explicitly in the weighting in the ensemble average. We do this in a self-consistent fashion, choosing the amplitude of each Fourier mode to reproduce the measured square magnitude of the mode. In practice, it is sufficient to do this only for the first few modes, which are the ones most severely impacted by the force. Also, in order to reduce the noise in the measurement of the mode amplitudes, we only calculate the averages after every block of $`N/10`$ realizations and update accordingly, where $`N`$ is the total number of realizations we intend to generate for these parameters. The force dependence of the mode amplitudes is obtained by a stepwise increase of $`F`$, where for each value of the force we use as a starting point the amplitudes measured at the previous value.

There is one other degree of freedom of the problem that also must be treated appropriately, namely the initial orientation of the filament at its constrained end. This is because the force not only straightens the filament, but also orients it toward its direction. We have employed two different algorithms to deal with this problem. The simpler approach is to simply bias the choice of initial orientation, and correct for this in the weighting. The zero force distribution of the orientation angle, $`\theta `$, of the tangent vector with respect to the direction of the force is $`P(\theta )=sin\theta d\theta =d(\mathrm{cos}\theta )`$. In the presence of the force we modified this to $`P(cos\theta )=exp(\alpha \mathrm{cos}\theta )d(\mathrm{cos}\theta ),`$ where one can choose $`\alpha `$ to reproduce the measured mean of $`cos\theta `$.

A more sophisticated approach is to explicitly perform the integral over $`\theta `$. We generate a given configuration and calculate its end-to-end distance $`R`$ and angle $`\theta _{\text{e}nd}`$ that the tangent to the end of the filament makes with the direction of the end-to-end vector $`𝐑`$. In order to evaluate the contribution of this realization to the average elongation, Eq. 12, one can perform a rigid rotation of this configuration about the end point and calculate the contribution to the statistical sum, from each angle of rotation $`\theta `$. The resulting integrals over $`\theta `$ can be expressed in terms of $`R`$ and $`\theta _{\text{e}nd}`$ as

$$_0^\pi 𝑑\theta sin\theta \{Rcos(\theta +\theta _{\text{e}nd}),1\}e^{fR\mathrm{cos}(\theta +\theta _{\text{e}nd})}$$
(13)
where the curly braces indicate the integrand for the numerator and denominator of the ensemble sum, respectively. These integrals can easily be evaluated numerically, using standard techniques. This second method is more efficient, but more involved to implement.

Using these techniques, we have measured the force-extension curve $`\overline{x}`$ vs. $`f`$. We present the results for $`a=.75`$, $`L=2\pi `$ in the left box of Fig. 2 together with the interpolation formula of MS

$$fa=\frac{\overline{x}}{L}+\frac{1}{4(1\overline{x}/L)^2}\frac{1}{4}$$
(14)
We see that for this value of $`a`$ ($`aL)`$ the agreement is very good. In the right half of Fig. 2 we present the results for $`a=48`$, again with $`L=2\pi `$. Here the MS formula does not work for small and intermediate values of the force, since we are not in the long-chain limit ($`La`$) for which it was derived. On the other hand, for this case analysis can also be done, as presented in the Appendix, and the analytic formula derived there is in perfect agreement with the simulations (see the right box of Fig. 2). Notice that the MS expression systematically overestimates the extension of the filament and approaches our results only in the infinite force limit. The physical basis of the deviation is the large orientational entropy which is not accounted for in the MS calculation.

Inspection of the right box of Fig. 2 suggests that both the simulation and our analytical approximation (Eq. 31) yield a finite initial slope of the extension versus force curve in the rigid rod limit ($`aL`$). Since the effective spring constant is proportional to the inverse of the initial slope of this curve this suggests that the linear response of the rigid rod to small stretching forces does not depend on the persistence length $`a`$ and is a function of its length $`L`$ only. The above expectation is confirmed by calculations reported in Appendix B where we show that the corresponding spring constant is $`3kT/L^2`$. Obviously, this behavior differs from the flexible chain limit ($`aL`$) for which the spring constant scales asPG  $`kT/(aL).`$ It also differs from the predictions of Ref. \[frey2,  \] where it was argued that the effective spring constant diverges in the rigid rod limit, as $`akT/L^3`$. Notice that such a divergence is not expected for a single filament under tension since there is a finite contribution to the susceptibility (and thus to the response to the external force) due to orientational fluctuations, even in the rigid rod limit. The difference between our result and that of Ref. \[frey2,  \] stems from the fact that the latter considers a network in which filaments are not free to rotate and the only angular averaging is over the randomly distributed orientations of different filaments in the network.

Not only can we calculate the average $`\overline{x}`$, but we can also measure the entire distribution $`P(x)`$. We show this in Fig. 3, for $`a=0.75`$ and $`6`$. In the former case, the persistence length is much shorter than length of the filament and the force-free distribution is nearly Gaussian. The peak of the distribution moves to higher values of the extension and its width decreases with increasing force. In the latter case (persistence length approximately equal to the length of the filament), the force-free distribution is nearly flat up to $`x\pm 5,`$ as expected for random sampling of orientations by a nearly rigid rod (due to thermal fluctuations, the effective length of this rod is somewhat shorter than its contour length, $`2\pi `$). The application of small force (in the range $`F/kT1`$) orients the filament and leads to dramatic change of the distribution. At larger forces, the distribution changes slowly with the force, in a manner similar to the case $`a=0.75`$. The transition from orientation-dominated to stretching-dominated behavior of nearly rigid rods is clearly observed in the right box of Fig. 2, where pronounced deviations from the Marko-Siggia formula are observed in the weak force regime (orientational effects are not taken into account in the MS theory).

## V Steric Constraints

It is fairly straightforward to incorporate steric constraints into our calculations. As an example, we consider the case where the filament is confined to the upper half space $`x>0`$. After generating each configuration according to the procedure outlined in the previous section, we simply reject any configuration which violates the constraint. Our method of integrating over the initial orientation angle $`\theta `$ can also be employed. We have simply to compute the limits of the integral over $`\theta `$ which can be done straightforwardly given the configuration. We present the results for the force-extension curve in Fig. 4 and for the distribution $`P(x)`$ in Fig. 5. We see that, as one might expect, the effect of the constraint is most pronounced at small $`F`$, where the straightening/aligning of the curve is smallest. For vanishing $`F`$, the average extension is positive as a result of the constraint. The constraint basically carves out a “hole” in the distribution at small $`x`$, resulting in a very asymmetric distribution for small force. For $`F/kT>1`$ the distribution function is practically unaffected by the presence of the wall and coincides with that shown in Fig. 3.

## VI Discussion

In this work we used biased Monte-Carlo simulations, based on an extension of our Frenet algorithm, to calculate the force-extension relations and the distribution functions of the extension for filaments of arbitrary persistence to contour length ratio. We showed that in the presence of external force, curvature fluctuations are no longer Gaussian and that their long wavelength Fourier modes are progressively suppressed by the force. This effect was taken into account in a self-consistent manner in the simulations by gradually increasing the force and adjusting the corresponding Fourier amplitudes of curvature fluctuations. Analytic results were obtained both for small fluctuations about the rigid rod limit, and in the limit of large force. In the former case, we calculated the force-induced renormalization of the length-scale-dependent bending moduli that govern the various Fourier modes of curvature. As expected, the largest effect is on the zeroth (constant curvature) mode; the corresponding modulus increases quadratically and linearly with the force in the weak and strong force limits, respectively. The underlying physical picture is reminiscent of the suppression of excluded volume by tensile force in the self-avoiding random walk model of polymers - while large scale fluctuations are strongly suppressed, small scale fluctuations (within the so called Pincus blobspincus ) are not affected by the force.

For flexible filaments with $`aL`$ our simulation reproduces the Marko-Siggia force-extension curves. However, our simulations and analytic results show that as the rigid rod limit $`aL`$ is approached, the MS expression increasingly overestimates the response of the filament at finite values of the force, and approaches the correct result only asymptotically, in the limit of infinite force. The discrepancy can be attributed to the neglect of orientational effects (the constant curvature mode is not taken into account) in the MS derivation which is strictly valid only in the $`a/L`$ $`0`$ limit. These orientational effects are however seen in the short-chain limit of the Gaussian model investigated recently by Winklerwinkler . The large force limit of the Gaussian model for short chains is however very different from the worm-like chain model studied herein. For long chains, the Gaussian model can reproduce either the weak or the strong force limit, by appropriate choice of the persistence length, but not both. This is an outcome of the length-scale dependence of the effective moduli discussed above and in Appendix C.

Examination of the distribution functions for the extension shows that increasing rigidity leads to a more rapid shift of the peak of the distribution and narrowing of its width, with external force. The largest deviations from Gaussian statistics of the distributions are observed for stiff filaments in the weak to intermediate force range ($`F/kT1`$), where fluctuations of initial orientation are not yet completely suppressed and where wall effects play an important role. While it may be difficult to approach the rigid rod limit in mechanical experiments on dsDNA that has a persistence length of $`500`$ Å, it can be readily accessed by experiments on actin filamentsalbert  for which the persistence length of order $`10`$ $`\mu `$m.

Finally, we would like to comment on possible extensions of the present work. Since we only considered here filaments with symmetric cross section and without spontaneous curvature, twist rigidity did not affect the distribution of extension and played no role in our considerations. However, since twist rigidity affects even the equilibrium distributions of ribbons with strongly asymmetric cross sectionsjcp-kr , it is expected that both twist rigidity and spontaneous twist will affect the deformation of such ribbons as well. Perhaps even more interesting is the effect of spontaneous curvature. Indeed, we have shown that a helical filament with sufficiently large twist to bending rigidity and radius to pitch ratios, undergoes a sequence of stretching instabilities under the action of tensile force, even in the absence of any thermal fluctuationsprl-kr . In the presence of thermal fluctuations or mechanical noise these instabilities are expected to lead to coexistence between different shapes of the deformed objectbenedek .

###### Acknowledgements.

DAK and YR acknowledge the support of the Israel Science Foundation. We also would like to thank J. Fineberg and G. Cohen for stimulating discussions.

## Appendix A Euler Angles

Since we are interested in three dimensions, we work in terms of the Euler angles $`\theta `$, $`\varphi `$ and $`\psi .`$ For convenience we will orient the force along the $`x`$-axis, i.e. along the direction $`\theta =\pi /2,`$ $`\varphi =0`$. The energy due to the force is then

$$E_F=Fx=F_0^L𝑑st_3^x=F_0^L𝑑s\mathrm{sin}(\theta (s))\mathrm{cos}(\varphi (s))$$
(15)
The Euler angles satisfy the differential equationsring :

$`{\displaystyle \frac{d\theta }{ds}}`$ $`=\omega _1\mathrm{sin}\psi +\omega _2\mathrm{cos}\psi `$
$`{\displaystyle \frac{d\varphi }{ds}}\mathrm{sin}\theta `$ $`=\omega _1\mathrm{cos}\psi +\omega _2\mathrm{sin}\psi `$
$`{\displaystyle \frac{d\psi }{ds}}\mathrm{sin}\theta `$ $`=(\omega _1\mathrm{cos}\psi \omega _2\mathrm{sin}\varphi )\mathrm{cos}\theta +\omega _3\mathrm{sin}\theta `$ (16)
We need to expand to quadratic order about the $`s=0`$ values $`\theta _0`$ and $`\varphi _0`$, and to linear order about $`\psi _0`$, (since the energy does not involve $`\psi `$ directly) giving

$`{\displaystyle \frac{d\theta }{ds}}`$ $`=\kappa \lambda (\psi \psi _0),`$
$`{\displaystyle \frac{d\varphi }{ds}}`$ $`={\displaystyle \frac{\lambda +\kappa (\psi \psi _0)}{\mathrm{sin}\theta _0}}{\displaystyle \frac{\mathrm{cos}\theta _0}{\mathrm{sin}^2\theta _0}}\lambda (\theta \theta _0),`$
$`{\displaystyle \frac{d\psi }{ds}}`$ $`={\displaystyle \frac{\mathrm{cos}\theta _0}{\mathrm{sin}\theta _0}}\lambda +\omega _3,`$ (17)
where $`\kappa `$ and $`\lambda `$ are the rotated curvatures

$`\kappa `$ $`\omega _1\mathrm{sin}\psi _0+\omega _2\mathrm{cos}\psi _0,`$
$`\lambda `$ $`\omega _1\mathrm{cos}\psi _0+\omega _2\mathrm{sin}\psi _0,`$ (18)
Thus, to second order in the curvatures

$`\theta (s)`$ $`=\theta _0+{\displaystyle _0^s}𝑑s^{}\kappa (s^{})+{\displaystyle \frac{\mathrm{cos}\theta _0}{2\mathrm{sin}\theta _0}}\left({\displaystyle _0^s}𝑑s^{}\lambda (s^{})\right)^2{\displaystyle _0^s}𝑑s^{}\lambda (s^{}){\displaystyle _0^s^{}}𝑑s^{\prime \prime }\omega _3(s^{\prime \prime }),`$
$`\varphi (s)`$ $`=\varphi _0+{\displaystyle \frac{1}{\mathrm{sin}\theta _0}}{\displaystyle _0^s}𝑑s^{}\lambda (s^{}){\displaystyle \frac{\mathrm{cos}\theta _0}{\mathrm{sin}^2\theta _0}}\left({\displaystyle _0^s}𝑑s^{}\kappa (s^{})\right)\left({\displaystyle _0^s}𝑑s^{}\lambda (s^{})\right)`$
$`+{\displaystyle \frac{1}{\mathrm{sin}\theta _0}}{\displaystyle _0^s}𝑑s^{}\kappa (s^{}){\displaystyle _0^s^{}}𝑑s^{\prime \prime }\omega _3(s^{\prime \prime }).`$ (19)
## Appendix B Force-Extension Curve

We will first compute the force-extension curve, via the partition function. We do this in two stages, depending on the value of the parameter $`\chi \sqrt{FL^2/4a}`$.

### B.1  Small $`\chi `$ ($`\chi 1,`$ arbitrary $`fL)`$

Here we write

$$E_F=FL\mathrm{sin}\theta _0\mathrm{cos}\varphi _0+F\delta $$
(20)
where $`\delta `$ is given by

$`\delta `$ $`={\displaystyle _0^L}ds[\mathrm{cos}\theta _0\mathrm{cos}\varphi _0{\displaystyle _0^s}\kappa \mathrm{sin}\varphi _0{\displaystyle _0^s}\lambda {\displaystyle \frac{1}{2}}\mathrm{sin}\theta _0\mathrm{cos}\varphi _0[\left({\displaystyle _0^s}\kappa \right)^2+\left({\displaystyle _0^s}\lambda \right)^2]`$
$`\mathrm{cos}\theta _0\mathrm{cos}\varphi _0{\displaystyle _0^s}ds^{}\lambda (s^{}){\displaystyle _0^s^{}}\omega _3\mathrm{sin}\varphi _0{\displaystyle _0^s}ds^{}\kappa (s^{}){\displaystyle _0^s^{}}\omega _3]`$ (21)
We expand $`e^{F\delta /kT}`$ to quadratic order, yielding

$`Z`$ $`={\displaystyle 𝑑\mathrm{\Omega }_0𝒟\kappa 𝒟\lambda \mathrm{exp}\left[\frac{a}{2}_0^L(\kappa ^2+\lambda ^2)+fL\mathrm{sin}\theta _0\mathrm{cos}\varphi _0\right]}`$
$`\times \left[1f\delta +{\displaystyle \frac{f^2}{2}}\delta ^2\right]`$ (22)
where $`𝑑\mathrm{\Omega }_0\frac{1}{4\pi }_0^\pi \mathrm{sin}\theta _0d\theta _0_0^{2\pi }𝑑\varphi _0`$. The functional integrals are all Gaussian and can be easily done by going to Fourier space. For example,

$$\kappa (s)=\kappa _0+\underset{n=1}{\overset{\mathrm{}}{}}\kappa _n^c\mathrm{cos}(\frac{2\pi ns}{L})+\underset{n=1}{\overset{\mathrm{}}{}}\kappa _n^s\mathrm{sin}(\frac{2\pi ns}{L})$$
(23)
so that

$$_0^L\left(_0^s𝑑s^{}\kappa \right)^2𝑑s=\frac{L^3}{3}\kappa _0^2+L^3\underset{n=1}{\overset{\mathrm{}}{}}\frac{(\kappa _n^c)^2+3(\kappa _n^s)^2}{2(2\pi n)^2}+mixedterms$$
(24)
where the mixed terms will vanish latter upon functional integration and are omitted in the following. Defining

$$𝒟_0\kappa 𝒟\kappa e^{\frac{a}{2}_0^L\kappa ^2𝑑s}=𝒟\kappa e^{\frac{aL}{2}\kappa _0^2\frac{aL}{4}_{n=1}^{\mathrm{}}\left((\kappa _n^c)^2+(\kappa _n^s)^2\right)}$$
(25)
and carrying out the functional integration by integrating over all the Fourier components

$$𝒟\kappa d\kappa _0\underset{i=1}{\overset{\mathrm{}}{}}d\kappa _i^c\underset{j=1}{\overset{\mathrm{}}{}}d\kappa _j^s,$$
(26)
we get

$$𝒟_0\kappa _0^L\left(_0^s𝑑s^{}\kappa \right)^2𝑑s=𝒩\frac{L^2}{2a}$$
(27)
where $`𝒩=𝒟_0\kappa `$ is the (normalization) integral over the exponential factor only. Similarly,

$$𝒟_0\kappa \left[_0^L\left(_0^s\kappa \right)𝑑s\right]^2=𝒟_0\kappa \left[\frac{\kappa _0L^2}{2}+\underset{n=1}{\overset{\mathrm{}}{}}\frac{\kappa _n^sL^2}{2\pi n}\right]^2=𝒩\frac{L^3}{3a}.$$

The functional integration with respect to $`\lambda `$ is carried out in the same manner. Putting this all together, we get

$$Z=𝒩^2𝑑\mathrm{\Omega }_0e^{fL\mathrm{sin}\theta _0\mathrm{cos}\varphi _0}\left[1\frac{fL^2}{2a}\mathrm{sin}\theta _0\mathrm{cos}\varphi _0+\frac{f^2L^3}{6a}\left(\mathrm{cos}^2\theta _0\mathrm{cos}^2\varphi _0+\mathrm{sin}^2\varphi _0\right)\right]$$
(28)
The integrals over solid angle are as follows:

$`{\displaystyle 𝑑\mathrm{\Omega }_0e^{fL\mathrm{sin}\theta _0\mathrm{cos}\varphi _0}}`$ $`={\displaystyle \frac{\mathrm{sinh}fL}{fL}}`$
$`{\displaystyle 𝑑\mathrm{\Omega }_0e^{fL\mathrm{sin}\theta _0\mathrm{cos}\varphi _0}\mathrm{sin}\theta _0\mathrm{cos}\varphi _0}`$ $`={\displaystyle \frac{\mathrm{cosh}fL}{fL}}{\displaystyle \frac{\mathrm{sinh}fL}{(fL)^2}}`$
$`{\displaystyle 𝑑\mathrm{\Omega }_0e^{fL\mathrm{sin}\theta _0\mathrm{cos}\varphi _0}\mathrm{cos}^2\theta _0\mathrm{cos}^2\varphi _0}`$ $`={\displaystyle \frac{1+\mathrm{cosh}fL}{(fL)^2}}{\displaystyle \frac{2\mathrm{sinh}fL}{(fL)^3}}`$
$`{\displaystyle 𝑑\mathrm{\Omega }_0e^{fL\mathrm{sin}\theta _0\mathrm{cos}\varphi _0}\mathrm{sin}^2\varphi _0}`$ $`={\displaystyle \frac{\mathrm{cosh}fL1}{(fL)^2}}`$ (29)
Collecting all terms we finally get

$$Z=𝒩^2\left\{\frac{\mathrm{sinh}fL}{fL}\frac{L}{6a}\left[\mathrm{cosh}fL\frac{\mathrm{sinh}fL}{fL}\right]\right\}$$
(30)
From this, we can differentiate to get the average projection of the end-to-end vector on the $`x`$-axis, $`\overline{x}=d(\mathrm{ln}Z)/df`$

$$\overline{x}=\frac{L}{\mathrm{tanh}fL}\frac{1}{f}\frac{2\chi ^2}{3f}\left[\frac{1}{\mathrm{tanh}fL}\frac{fL}{\mathrm{sinh}^2fL}\right]$$
(31)
At large $`fL`$ (but still small $`\chi `$), this goes as $`LL^2/(6a)1/f`$. For small $`fL`$ this projection goes linearly with $`f`$ as $`\overline{x}fL^2/3(1L/(3a))`$ and we conclude that in the limit $`aL`$ the effective “spring constant” is $`3kT/L^2`$.

### B.2 Large $`fL`$ (arbitrary $`\chi `$)

We start by expanding $`\theta _0`$, $`\varphi _0`$ about the values $`\theta _0^{}`$, $`\varphi _0^{}`$ which minimize $`E_F`$. To lowest order,

$`\theta _0^{}`$ $`={\displaystyle \frac{\pi }{2}}{\displaystyle \frac{1}{L}}{\displaystyle _0^L}𝑑s{\displaystyle _0^s}𝑑s^{}\kappa `$
$`\varphi _0^{}`$ $`={\displaystyle \frac{1}{L}}{\displaystyle _0^L}𝑑s{\displaystyle _0^s}𝑑s^{}\lambda `$ (32)
so that the expansion reads

$`E_F`$ $`=FL{\displaystyle \frac{F}{2L}}\left[\left({\displaystyle _0^L}𝑑s{\displaystyle _0^s}𝑑s^{}\kappa \right)^2+\left({\displaystyle _0^L}𝑑s{\displaystyle _0^s}𝑑s^{}\lambda \right)^2\right]+{\displaystyle \frac{FL}{2}}\left[\left(\theta _0\theta _0^{}\right)^2+\left(\varphi _0\varphi _0^{}\right)^2\right]`$
$`+{\displaystyle \frac{F}{2}}{\displaystyle _0^L}𝑑s\left[\left({\displaystyle _0^s}𝑑s^{}\kappa \right)^2+\left({\displaystyle _0^s}𝑑s^{}\lambda \right)^2\right]`$ (33)
Let us now consider the functional integral over $`\kappa `$, the integral over $`\lambda `$ being exactly the same. We have

$`Z_\kappa `$ $`={\displaystyle 𝒟\kappa e^{\frac{a}{2}_0^L𝑑s\kappa ^2\frac{f}{2}_0^L𝑑s\left(_0^s𝑑s^{}\kappa \right)^2+\frac{f}{2L}\left[_0^L𝑑s_0^s𝑑s^{}\kappa \right]^2}}`$
$`={\displaystyle }𝒟\kappa \mathrm{exp}\{{\displaystyle \frac{aL}{2}}\kappa _0^2{\displaystyle \frac{aL}{4}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}((\kappa _n^c)^2+(\kappa _n^s)^2)`$
$`{\displaystyle \frac{fL^3}{2}}[{\displaystyle \frac{\kappa _0^2}{12}}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}({\displaystyle \frac{(\kappa _n^c)^2}{2(2\pi n)^2}}+{\displaystyle \frac{(\kappa _n^s)^2}{2(2\pi n)^2}}{\displaystyle \frac{2\kappa _0\kappa _n^c}{(2\pi n)^2}})]\}`$ (34)
We can diagonalize this expression by shifting $`\kappa _n^c`$, defining

$$C_n\kappa _n^c\frac{b_n}{a_n}\kappa _0$$
(35)
with

$$a_n\frac{aL+b_n}{2};b_n\frac{fL^3}{(2\pi n)^2}$$
(36)
yielding an effective modulus of the zero mode

$$a_0=a+\frac{fL^2}{12}\underset{n=1}{\overset{\mathrm{}}{}}\frac{b_n^2}{a_n}=\frac{a\chi }{\mathrm{tanh}\chi }$$
(37)
We can now calculate $`Z_\kappa `$:

$`Z_\kappa /Z_{\kappa ,0}`$ $`=(a_0/a)^{1/2}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{aL}{2a_n}}`$
$`=(a_0/a)^{1/2}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left(1+{\displaystyle \frac{\chi ^2}{(\pi n)^2}}\right)^1`$
$`=(a_0/a)^{1/2}{\displaystyle \frac{\chi }{\mathrm{sinh}\chi }}`$
$`=\left[{\displaystyle \frac{2\chi }{\mathrm{sinh}2\chi }}\right]^{1/2}`$ (38)
where $`Z_{\kappa ,0}`$ is the $`\kappa `$ partition function at $`f=0`$. The entire partition function is then

$$Z=𝒩^2\frac{e^{fL}}{2fL}\frac{2\chi }{\mathrm{sinh}2\chi }$$
(39)
For small $`\chi `$, this yields

$$Z𝒩^2\frac{e^{fL}}{2fL}\left(1\frac{2\chi ^2}{3}\right),$$
(40)
which agrees precisely with the large $`f`$ expansion of Eq. (30).

The average extension in this second region ($`fL1`$) is given by

$$\overline{x}=L\frac{1}{2f}\frac{\chi }{f\mathrm{tanh}2\chi }$$
(41)
which for small $`\chi `$ gives $`\overline{x}L1/fL^2/(6a),`$ in accord with the large $`fL`$ limit of the first region ($`\chi 1`$). At large $`fL`$, it agrees with the Marko–Siggia result, Eq. (14). A uniform approximation for the force-extension relation is then

$$\overline{x}\frac{1}{2f}\frac{\chi }{f\mathrm{tanh}2\chi }+\frac{L}{\mathrm{tanh}fL}\frac{2\chi ^2}{3f}\left(\frac{1}{\mathrm{tanh}fL}\frac{fL}{\mathrm{sinh}^2fL}1\right).$$
(42)
This result is confirmed by simulations, as shown in data for $`a=48`$, Fig. 2. It even works well when the persistence length slightly shorter than the chain length, as can be seen in the data for $`a=6`$ in the right half of Fig. 4. Since our analysis is limited to the small fluctuation (stiff rod and/or large force) regime, we expect it to break down for very short persistence lengths in the small force limit. This is clearly observed in the case of $`a=0.75`$, shown in the left half of Fig. 2.

## Appendix C Moduli

We now turn to the computation of the effective moduli of the various modes. We first treat the zero modes, starting as above with the small $`\chi `$ limit. We proceed by calculating a modified partition function, $`\stackrel{~}{Z}`$ with a modified modulus of the zero modes, $`\stackrel{~}{a}`$, replacing the bare modulus $`a`$. Then

$$\frac{L}{2}\kappa _0^2+\lambda _0^2=\frac{d\mathrm{ln}\stackrel{~}{Z}}{d\stackrel{~}{a}}|_{\stackrel{~}{a}=a}$$
(43)
Tracking the changes induced in the partition function calculation above, we find

$$\stackrel{~}{Z}=\frac{a}{\stackrel{~}{a}}𝒩^2\left\{\frac{\mathrm{sinh}fL}{fL}fL^2\left[\left(\frac{1}{3\stackrel{~}{a}}+\frac{1}{6a}\right)\left(\frac{1}{4\stackrel{~}{a}}+\frac{1}{12a}\right)\right]\left[\frac{\mathrm{cosh}fL}{fL}\frac{\mathrm{sinh}fL}{(fL)^2}\right]\right\}$$
(44)
giving

$$\mathrm{ln}\stackrel{~}{Z}=\mathrm{ln}a\mathrm{ln}\stackrel{~}{a}+2\mathrm{ln}𝒩+\mathrm{ln}\frac{\mathrm{sinh}fL}{fL}fL^2\left(\frac{1}{12\stackrel{~}{a}}+\frac{1}{24a}\right)\left(\frac{\mathrm{cosh}fL}{\mathrm{sinh}fL}\frac{1}{fL}\right)$$
(45)
so that

$$\frac{d\mathrm{ln}\stackrel{~}{Z}}{d\stackrel{~}{a}}|_a=\frac{1}{a}+\frac{fL^2}{12a^2}\left[\frac{\mathrm{cosh}fL}{\mathrm{sinh}fL}\frac{1}{fL}\right]$$
(46)
This implies an effective zero mode modulus of

$$a_0=a+\frac{fL^2}{12a}\left(\frac{1}{\mathrm{tanh}fL}\frac{1}{fL}\right)$$
(47)
At small and large forces, this reads

$$a_0\{\begin{array}{cc}a+f^2L^3/(36a)& fL1\hfill \\ a+fL^2/(12a)& fL1\hfill \end{array}$$
(48)
At large $`f`$, the effective zero mode modulus was calculated above, Eq. (37). The small $`\chi `$ limit of this agrees with the large force limit of Eq. (48), as it must. A uniform approximation for $`a_0`$ reads

$$a_0=\frac{fL^2}{12a}\left(\frac{e^{fL}}{\mathrm{sinh}fL}\frac{1}{fL}\right)+a\frac{\chi }{\mathrm{tanh}\chi }$$
(49)
This result is confirmed by the simulations as seen in Fig. 6.

Now we turn to the finite wave-number modes. The interesting thing to note here is that the cos modes renormalize differently than the sin modes, due to the coupling of the zero mode to the cos modes. At large force, the effective modulus of the sin modes is the trivial

$$a_n^s=a+\frac{fL^2}{(2\pi n)^2}$$
(50)
so that the long wavelength modes are the most effected by the force. The effective modulus of the cos modes follows from

$$(\kappa _n^c)^2=C_n^2+\frac{b_n^2}{a_n^2}\kappa _0^2=\frac{1}{a_n}+\frac{b_n^2}{a_n^2}\frac{1}{a_{\text{eff},0}L}$$
(51)
with

$$a_n^c=\frac{2}{L(\kappa _n^c)^2}$$
(52)
This expression is fairly messy, so let us examine its various limits. For large enough mode number, the effective modulus is $`a+fL^2/(2\pi n)^2`$ as for the sin modes. For intermediate forces, all modes are large enough for this to be the case. For strong forces, however, the effective modulus is $`L\sqrt{fa}/4`$ up to some mode number, at which point it crosses over to the previous result.