text:
# The symbiotic star CH Cygni. III. A precessing radio jet

## 1 Introduction

Symbiotic stars occupy an extreme and relatively poorly understood region of the binary classification scheme. The name was coined by Merril (1941) to describe stars which appeared to have a combination spectrum: that of high excitation lines usually associated with a hot, ionised nebula superimposed on a cool continuum with prominent absorption features consistent with a late–type star. At present they are understood to be interacting binaries (with orbital periods of a few to tens of years) consisting of a cool giant losing material mostly via the stellar wind and a hot, luminous compact object which ionises a portion of the cool component wind (Seaquist et al., 1984). Such a state of affairs represents the so–called quiescent phase, which can be interrupted by periods of activity. The active phases start with an eruption of the hot star, an event indicated photometrically by an increase of the star’s brightness by 2-6 mag, and/or spectroscopically by high velocity broad emission features from the central star. Both radio and Hubble Space Telescope (HST) imaging can directly resolve the remnants of such dramatic events (Hack & Paresce 1993; Eyres et al. 1995; Kenny et al. 1996; Richards et al. 1999; Eyres et al. 2001b).

The symbiotic star CH Cygni displays particularly complex behaviour. Optical spectroscopic studies by Mikolajewski et al. (1987) showed that the orbits of the stars within the system are likely to be coplanar and eclipsing, with eclipses separated by around 5700 days. Later studies (Mikolajewski et al. 1990b and references therein) confirmed this period.

Further spectral studies of the system (Hinkle et al., 1993) suggested that CH Cyg is probably a triple-star system consisting of an inner 756-day period binary which is orbited by an unseen G-K dwarf on a 5300 day orbit. Skopal et al. (1996) uncovered 756 day interval eclipses which show that all three stars in the system are likely to be in coplanar orbits.

Each of the three outbursts seen since 1978 was accompanied by high velocity broad emission features consistent with mass outflows. During 1984-85 the material was ejected at $``$600-2500 km s-1 (Mikolajewski & Tomov, 1986), whilst the 1992–95 active phase was characterized by sporadic and in part bipolar outflow at $``$1000-1600 km s-1 (Leedjärv & Mikolajewski 1996; Skopal et al. 1996; Iijima 1996) and finally, during the recent, 1998–2000, outburst mass outflows at about 1000 km s-1 were observed (Eyres et al. 2001a, hereafter Paper II). The outflows may be correlated with a significant increase of the radio emission and the radio light curves during these periods fit well with the optical ones (Kenny et al., 1996). The 1984 mass ejection has been linked to the emergence of bipolar emission (Taylor et al., 1986) which has been attributed to high velocity ($``$1000 km s-1) jets. Non–thermal emission features have been discovered in these jets from radio observations (Crocker et al. 2001, hereafter Paper I), explained by shocked regions that arise when the high–speed ejecta interact with existing wind material.

This paper aims to extend our understanding of the jets in this interesting system by developing theories and models to explain the change in their morphology in the years following their emergence.

In this paper we present analysis of radio data from the VLA and MERLIN, along with HST observations.

## 2 Observations

CH Cygni has been observed on several occasions since 1985, in several different radio wavelengths and optical domains, further details of which are given in Paper I. All observations have been presented previously with the exception of the 2000 November 10? VLA data which is shown here for the first time.

Radio data reduction was carried out using the AIPS software package (Greisen, 1999). The flux density of the phase reference source was determined, and the data calibrated to correct phase and amplitude instabilities. These corrections were applied to the observations of CH Cygni and a map produced by deconvolving the instrumental response (the so-called “dirty beam”) from the map using the well-established CLEAN algorithm (Högbom, 1974).

Observations made on 1999 August 12 with the Hubble Space Telescope’s Wide Field and Planetary Camera 2 through the F502N filter, corresponding to the \[O III\] line at $`\lambda \lambda `$4956 and 5007Å) (Biretta et al., 2001), are used here as these, with a diffraction limted resolution of around 0.05 arcseconds, offer direct comparison with the spatial resolution of the radio data (see Table 1. Each observation was made up of two sub–exposures, allowing cosmic ray subtraction to be carried out (see Paper II and references therein).

## 3 Initial data analysis

A preliminary examination of radio maps of CH Cygni (taken with the VLA and MERLIN), as well as the HST optical images (using the F502N \[OIII\] filter) show that the position angle of the extended regions changes drastically over the 15 years between the earliest (1985) and most recent (2000) observations.

A more quantitative analysis requires a measurement of the position angle of the central regions of the radio nebula. A problem inherent in this method is that the geometry of the restoring beam (in the case of the radio observations) may dominate over any natural geometry of the region. This arises as a result of the $`(u,v)`$ coverage of the instrument, which depends upon the altitude of the source above the horizon.

To avoid these problems, the radio maps were recreated using the methods of Högbom (1974) but a circular restoring beam was imposed upon them, as laid out in Table 1.

Two dimensional Gaussians were fitted to the bright peak in the centre of each map, and the position angle (anticlockwise through east from north) of the semi–major axis of the resulting ellipse was taken as the position angle of the central material. The uncertainties are those produced by the AIPS task ’IMFIT’, and are derived from the ratio of the axes of the fitted ellipse. These measurements are shown in Table 2.

Since the details of the alignment of the CH Cygni system are not known, it was assumed that the orbital plane of the system is being viewed edge–on, as eclipses are seen within the system (Mikolajewski et al. 1987; Skopal et al. 1996). The extension might be expected to precess about an axis perpendicular to the orbital plane if this is associated with the collimating mechanism. Assuming that the rate of precession is constant, the position angle, $`\theta _{\mathrm{PA}}`$ should follow a simple sinusoid so that at time $`t`$ the position angle is given by

$$\theta _{\mathrm{PA}}=\theta _{\mathrm{mid}}+\theta _\mathrm{o}\mathrm{cos}\left(2\pi \times \frac{tt_0}{P}\right).$$
(1)
Such a curve was fitted, giving a precession of period $`P=8373\pm 800`$ days about a position angle on the sky of $`\theta _{\mathrm{mid}}=146.7\pm 3.4`$ degrees. The opening angle of the precession cone was $`\theta _\mathrm{o}=37.1\pm 4.7^{}`$.

This precession period is longer than the $`5294\pm 117`$ day period quoted for the length of the orbit of the outer giant of the system by Hinkle et al. (1993). The precession period we see could be a resonance of this, the 756 day orbital period, or both. However, the period we find should be treated with some caution, as our data set does not appear to cover a full period of precession.

## 4 Evidence of precession within the radio jets

There is a large deviation of some of the points from the fitted curve. A flaw inherent in the method of fitting gaussians is that it will naturally pick up the position angle of material at a certain distance from the core – that distance being determined by the resolution of the map. Hence what is being picked up is the position angle of material that was ejected some time before the observation was made.

If the material in the extension is being ejected ballistically from a precessing source then it would be expected to trace out a conical helix of a form seen in many precessing objects (Zaninetti & van Horn, 1988). Expansion has been detected in the jet since the first radio jets were seen by Taylor et al. (1986).

Precessing jets are frequently seen in astronomical objects, one of the best studied being SS433. Its extended nebula exhibits a characteristic “S” profile that has been attributed to a precessing jet (Hjellming & Johnston, 1981) with great success. Code derived from the work of R. E. Spencer (priv. comm.) was adapted for CH Cygni by removing the relativistic behaviour of the jets. This code gives the on–the–sky positions of discrete knots in the jet at a given time, as determined by a set of parameters. These patterns were compared to the distribution of the 40 brightest CLEAN components in each map at a given frequency and the deviation minimised through use of the Downhill Simplex Method (Nelder & Mead, 1965).

The process was carried out using the set of 5 GHz maps taken by the VLA in the highest resolution “A” configuration, as these provided the most extensive coverage of the system over time.

### 4.1 Simple ballistic model

The simplest model of radio jet precession is one in which material is ejected continuously, at a constant velocity, from the poles of a precessing object. The locus described by the points along the resulting curved path at a time $`t`$ is described (see Figure 2) in parametric form by

$$x=v_0(t\tau )\mathrm{cos}\theta _0$$
(2)
$$y=v_0(t\tau )\mathrm{sin}\theta _0\mathrm{sin}(\varphi )$$
(3)
$$z=v_0(t\tau )\mathrm{sin}\theta _0\mathrm{cos}(\varphi )$$
(4)
where $`x`$, $`y`$ and $`z`$ are the coordinates of a given knot of material, ejected from the central source (at the origin) at a time $`\tau `$ with velocity $`v_0`$. The opening angle of the precession cone is $`\theta _0`$, measured from the $`x`$ axis, and the azimuthal angle of the jet (relative to the $`y`$ axis) is $`\mathrm{\Phi }=\mathrm{\Omega }\tau `$ where $`\mathrm{\Omega }`$ is the precession rate.

Such a model has been used to fit curves to the high–velocity jets of the microquasar SS433 by Spencer (1984). For this model, the parameters were: $`P`$, the precession period of the jet ($`=\frac{2\pi }{\mathrm{\Omega }}`$); $`\tau _0`$, a fixed time corresponding to a point when the precession angle of the jets was along the line of sight (the $`y`$ axis in Figure 2); $`\theta _{\mathrm{mid}}`$ the position angle (on the sky) of the axis of precession; $`\theta _0`$, the opening angle of precession cone; $`v_0`$, the velocity of the material ejected in the jet and $`i`$, the inclination angle of the precession axis to the line of sight.

The size of the initial simplex was set up so that it had parameters centred around $`P=5294`$ days (Hinkle et al., 1993), $`T_0=`$JD2446275 (Mikolajewski et al., 1990a), $`\theta _{\mathrm{mid}}=140^{}`$, $`\theta _0=30^{}`$ (Section 3), $`V=1230`$ kms-1 (Paper I) and $`i=90^{}`$. The process was allowed to run until the standard errors of the evaluated fits became smaller than 0.01. The “best–fit” parameters from this run were then fed in as the starting values for a second run. This process was repeated so that four runs were carried out. The final set of parameters is shown in Table 3. Results very similarb to this were found if the starting period was set to 8373 days, as found in Section 3. The resulting ephemerides, resulting from material being ejected in discrete bullets every 50 days and following the model described above, are shown in Figures 3 and 4.

The model fits reasonably well to the geometry seen in the maps, especially in those resulting from observations closely following the emergence of the radio jet in 1986. In these images (JD2446509 and JD2446576) the jets are two–sided and well extended, making fitting more robust.

In later epochs, when the jet is more one–sided, the curves still fit the morphology seen in the radio maps well. The fact that such a fit is achieved so long after the initial outburst (and following at least two other outbursts) demonstrates that the behaviour of the bipolar jets is not far from a simple ballistic model of the kind applied successfully to SS433 and other sources.

### 4.2 Disrupted jet model

The greatest deviations of the simple model from the observations are seen in the 1991 map (which has a low peak brightness and high background noise level) and in the motion of the NW component in the 1987 and 1988 maps. If the component is the same object in all of the maps, then it is moving in a fashion entirely different to a simple, ballistic, bullet of material. By inspection, it appears to be decelerating, perhaps due to interaction with dense material.

The synchrotron theory put forward in Paper I required ejecta to be impacting upon previously ejected stellar wind material, in order to cause shocks and local condensations in the magnetic field. This might be expected to cause some degree of deceleration and disruption.

Direct evidence of deceleration is provided by optical spectroscopy (Skopal et al., 2001), which detects the velocity of material within a few hundred R of the central stars. Terminal velocities of up to 2500 km s-1 are seen. Radio maps, which observe the motion of material much further out (from a few AU up to 100s of AU), never detect on–the–sky velocities of more than 1500 km s-1, assuming a distance to CH Cygni of 268 pc (Viotti et al., 1998). This is evidence of a large degree of deceleration between these two scales.

The actual offsets of the apparently–slowed component from the central radio peak at the different observing epochs are shown in Table 4.

A component moving ballistically, as described in Section 4.1, would exhibit a constant position angle ($`\theta _{\mathrm{PA}}`$) and a radial distance from the core ($`R`$) that increased linearly with time. If there is simple deceleration, $`\theta _{\mathrm{PA}}`$ would remain constant while $`R`$ changed according to an environment–dependent model.

Compared with the $`1230`$ km s-1 velocity seen in the main portion of the jets, the north western emission is almost stationary on the sky. It is possible that the material has indeed almost stopped moving, perhaps due to interactions with the surrounding environment. The jet may have become entirely disrupted through shock–interactions and so would no longer obey the simple ballistic model.

For a jet that has a time–variable ejection direction, such as a precessing bipolar jet, shocks will form between the jet material and the surrounding medium not only in the main direction of jet propagation but also as a result of the transverse interactions between the two regions (Dyson, 1987).

The transverse velocity of the jet becomes an important factor at large distances from the point of ejection.

A colimated outflow will only survive until the transverse shock has crossed the diameter of the jet. To simplify things, the approximation of Raga et al. (1993) is used so that

$$t_\mathrm{d}\sqrt{\frac{r_\mathrm{j}}{\beta v_0\mathrm{\Omega }\mathrm{sin}\theta _0}}.$$
(5)
Most of these values can be estimated from work already carried out. The jet radius, $`r_\mathrm{j}`$, can be estimated by measuring the size of the anomalous region in the 1986 VLA U band map. The region was found to be approximately circular with a mean diameter of 0.18 arcsec. Assuming a distance of 268 pc this corresponds to a linear size of $`6.9\times 10^{12}`$ m. The initial jet ejection velocity was found in Paper I to be $`1.23\times 10^6`$ m s-1 which is in agreement with that found by fitting to the the purely ballistic ejection model. The rate of precession can be found from the fitted period of 6519 days to be $`\mathrm{\Omega }=1.1\times 10^8`$rad s-1. The opening angle, $`\theta _0`$ was estimated previously to be $`35^{}\pm 1^{}`$.

The main unknown is the relative density parameter, $`\beta `$. Following Raga et al. (1993), $`\beta `$ is assumed to be constant, so that the ambient medium is entirely homogeneous. This may not be the case, as a constant spherical outflow from a star is expected to produce an ambient density that varies as $`\rho _\mathrm{e}r^2`$, and the presence of a jet would be expected to cause disruptions in the local region. However, without more knowledge of this region, a homogeneous medium is assumed for simplicity. Raga & Noriega–Crespo (1993) derived a ratio of $`\rho _\mathrm{e}/\rho _\mathrm{j}0.05`$ from line ratios in the Herbig–Haro object HH 111 V. If this applies to CH Cygni, then $`\beta `$ is approximately equal to 0.2.

Applying these parameters to Equation 5 gives a value for the time taken for the disruptive shock to cross the jet as $`t_\mathrm{d}7\times 10^7`$ s. The distance the material has covered in this time (along the $`x`$ axis in Figure 2) will be, for the parameters discussed above, $`x_\mathrm{d}=311`$ AU. The actual distance of the observed anomalies in the radio jet in 1986–1988 is indeed around 300 AU.

To determine whether this disruption is occurring, direct imaging of the shocked region is necessary. This has been done in the case of the Herbig–Haro objects HH46/47 (Reipurth & Heathcote, 1991) and HH34 (Reipurth & Heathcote, 1992). H$`\alpha `$ maps were subtracted from \[S II\] maps (at wavelengths 6717Å and 6731Å) to produce images in which the leading bow shock and trailing jet shock could be seen separately. The ratio of the intensity of emission in the H$`\alpha `$ and \[S II\] was then used to derive the relative densities in the jet material and surrounding medium (Raga & Noriega–Crespo, 1993). Unfortunately, whilst H$`\alpha `$ emission maps of CH Cygni exist (Paper II), no maps in \[S II\] are available.

### 4.3 Models with reflective symmetry

If the geometry of the jets is a result of them being gravitationally attracted to the outer giant, they would be expected to exhibit a morphology in which the position angle of the material is a function of distance from the ejection source, and in which there is reflective symmetry about the central core (Zaninetti & van Horn, 1988). Such a model was applied successfully to 3C449 by Lupton & Gott (1982). Such a pattern does not fit the observed extensions in the CH Cygni system, implying that the precession is due to effects within the collimating mechanism, rather than post–ejection gravitational influences.

### 4.4 Ballistic model with time–dependent ejection velocity

CH Cygni is observed to undergo irregular outbursts with associated variable mass outflow rates. Such a state of affairs renders invalid the assumption that what is being observed is a simple, continuous jet with constant ejection velocity. A jet that results from truly irregular outflows is difficult to model analytically using the methods described within this work.

Recent work, as described in Paper II, has discovered what may be two types of mass outflow from the system. One is optically thin, at a velocity measured spectroscopically to be 1200 km s-1 and the other is optically thick and highly irregular. There is evidence that the increase in mass–transfer from the cool component to the compact component that occurs near to times of periastron leads to higher luminosity of the compact component and, as a consequence, greater radiatively driven mass outflow.

Assuming a distance of 268 pc (Viotti et al., 1998), the movement between 1986 and 1988 shown in Table 4 represents an on–the–sky velocity of only $`200`$ km s-1. Extrapolating this motion backwards gives the time of ejection of this material, assuming ballistic motion, as around December 1978. This time is near the beginning of the large outburst that lasted between 1977 and 1984. Minor outflows have been observed spectrally during the outburst (Hack et al. 1988; Anderson et al. 1980) at velocities of a few hundred km s-1. These may coincide with the ejection of this anomalous packet of material.

Recent HST observations (from October 1999) made by Corradi et al. (2001) show the existence of loops of material in the \[OII\] nebula at a distance of $`1.5`$ and $`5`$ arcsec. Echelle spectroscopy of these loops places their kinematical axis at around 100 km s-1. Allowing for projection effects (i.e. assuming they are tilted at 37 to the plane of the sky) their kinematical ages are approximately 15.2 years for the inner loop and 51 years for the outer one, giving the times of ejection of the loops as August 1984 and February 1949 respectively. The 1984 event fits in well with the emergence of the radio jets first observed by Taylor et al. (1986) and the model proposed in this section, but no records of a 1949 outburst could be found (Mikolajewski et al., 1990a).

## 5 Possible causes of precession

Precession of bipolar jets is a common phenomenon in astrophysics (Fukue & Yokoo, 1986). The precession is usually attributed to motion of the collimating mechanism of the jets.

A simple collimating structure is a circumstellar torus or disk, such as that postulated by Fukue (1983). Flickering in the optical brightness of CH Cygni has been observed on several occasions (Wallerstein 1968; Mikolajewski et al. 1990b) and has been explained by the orbital motion of a bright point on a disk surrounding the compact component.

There are several mechanisms by which such a disk could be induced to precess, and four of these scenarios are examined here. Only one of these is able to cause precession at the rate that is observed, given suitable estimates for the appropriate parameters.

### 5.1 Precession of a tilted accretion disk

The gravitational influence of the nearby cool component, or perhaps the possible third star in the system, could force precession in this disk, and hence the jets. Following the work of Merritt & Petterson (1980), the angular speed of a torus or disk, with radius $`r`$, undergoing forced precession as a result of the influence of a companion is

$$\mathrm{\Omega }=\frac{3}{4}\sqrt{\frac{G}{m_\mathrm{h}}}\frac{m_\mathrm{h}m_\mathrm{c}}{m_\mathrm{h}+m_\mathrm{c}}\frac{r^{\frac{3}{2}}}{a^3}\mathrm{cos}\alpha $$

where $`m_\mathrm{h}`$ is the mass of the hot component, $`m_\mathrm{c}`$ is the mass of the cool component, $`a`$ is their separation and $`\alpha `$ is the angle between the equatorial plane of the torus and the orbital plane of the stars. Estimates for the parameters, from the work of Skopal (1997), are $`m_\mathrm{h}=0.2`$M, $`m_\mathrm{c}=0.93`$M and $`a=1.5`$ AU. This relation assumes that the disk is centred on the hot component, and the companion star causing the forced precession is the cool giant on the 756 day orbit. It would also be sensible to assume that $`\alpha <45^{}`$ (Mikolajewski et al., 1990b) and $`0<r<a`$. To obtain a precession period of 8373 days requires a disk radius of around 0.7 AU. This is a very large accretion disk, and it would be hard to imagine how such a disk could form purely via wind accretion mechanisms, making this scenario unlikely.

If instead it is assumed that CH Cygni is a binary system on a 5294 day orbit, with $`a=7.9`$AU the disk radius is found to be around 21 AU, around twice the size of the proposed semimajor axis of the outer giant’s orbit.

### 5.2 Slaved disk

An attempt to explain the observed 164 day period precession of the bipolar jets in SS433 was made by van den Heuvel et al. (1980), by assuming that the compact component was surrounded by an accretion disk that precessed at a rate governed by the apsidal motion of the parent star.

Using Equations 12, 13 and 16 of that paper, taking reasonable stellar and orbital parameters of CH Cygni from Skopal (1997) and assuming that the coefficient of apsidal motion is 0.005 (Schwarzschild, 1958) suggests that the precession period of a slaved disk would be around 12000 years. This result is a strong function of the radius of the cool component, so an increase in this star’s size by a factor of 4 would reduce this precession period to about 15 years.

To evaluate the possibility of such a model, precise eclipse timings need to be made in the future.

### 5.3 Irradiated disk

An accretion disk that is irradiated by a central source can become warped as a result of a torque that arises due to radiation that is re–emitted from the disk (Livio & Pringle, 1997).

Such a warped disk might be expected to precess. In contrast to the gravitationally induced precession mentioned above, this is “self precession”, in that no companion star is required. Following the work of Livio & Pringle (1996), the period of this precession is of the same order as the time it takes for the warping to develop in the disk. Using sensible parameters for the mass of a possible accretion disk (around $`10^7`$ M), this timescale is of the order of $``$2000 years, making this an unsuitable model for precession.

### 5.4 Magnetic precession

The inner region of an accretion disk surrounding a magnetic star, such as a white dwarf (Mikolajewski et al., 1990b), is subject to a twisting torque and a warping torque. These arise as a result of the interaction between the surface current of the disk and the horizontal magnetic field that is a feature of the dipole field of the star (Lai, 1999).

The warping torque arises through the twisting of the field that threads the disk vertically. The precessing torque is a result of the azimuthal screening current that occurs if the disk is diamagnetic. In this case, the vertical field lines cannot penetrate the disk. An azimuthal current is induced in the disk. The radial component of the central star’s magnetic field interacts with this current to set up a vertical force acting on the disk. This force is uneven across the disk, so a net torque results. This torque acts to cause the disk (and hence the collimation mechanism of the jets) to precess about the magnetic axis of the star.

Lai (1999) gives the angular precession rate of such a disk in Equation 2.35 of that paper. Presently, most of the required quantities are unknown, so only an order of magnitude estimate is possible. The magnetic moment of the white dwarf was estimated to be $`2\times 10^{28}`$ T m3 using the method of Postnov (1985), and the surface density of the disk was assumed to be of the order 10 kg m-2 (Lin et al., 1988). The angle between the disk axis and stellar rotation axis was taken to be equal to the opening angle of the precession cone, 35, found previously. The magnetic axis of the system was assumed to lie along the disk axis.

The magnetically induced precession period is then directly related to the radial distance, $`r`$, of the warped torus from the white dwarf. Since this varies as $`r^7`$, small changes in $`r`$ greatly affect the period. A warped torus at around 2 R from the star will precess with a period of around 6000 days.

A magnetic collimation and ejection mechanism for the jets gives a possible explanation for their one–sidedness (Chagelishvili et al., 1996). Jets powered by a central magnetic rotator (Mikolajewski et al., 1990b) can become inherently asymmetric, if the solution to the magnetohydrodynamic equations for each hemisphere of the star are asymmetric. Such conditions arise when there is no change in the relative sign of the poloidal and toroidal fields across the equator of the star. The different magnetic states above and below the equator can then easily lead to a condition in which only one supports the collimation and ejection of material via a jet (Wang et al., 1992).

## 6 Conclusion

It is certain that the bipolar jets of CH Cygni exhibit precession, with a period controlled by the motion of the mechanism resposible for the collimation of the outflow. In addition, the variability of the mass–loss rate from the cool component causes the material to be ejected at a rate that is highly variable (between $`500`$ and $`2000`$km s-1), leading to complications in predicting the outflow geometry. The activity of the bipolar ejection seems to be tied to the state of the hot–component, with fast jets seen at or soon after times of optical outburst. This activity then decreases during times of quiescence.

Although the precession period found here is similar to the orbital period of the outer giant in the system, the current time resolution of the observations does not rule out shorter periods. A precession that resulted from the motion of the inner, symbiotic, pair and had a similar period to its 756 day orbit would need observations to be taken at least every year for it to be detected.

The simple ballistic model fits the geometry of the nebula extremely well. Knots seen in the larger–scale radio maps (the VLA in C, X and U band) that do not fit the simple ballistic model can be explained as regions of jet disruption caused by sideways motion of the ejected material brought about as a result of the precession. The predicted distance at which this disruption would occur is in close agreement with the observations, although it relies upon several assumptions about the nature of the ambient material. Confirmation of this model would require simultaneous optical images in both H$`\alpha `$ and \[S II\].

An alternative explanation is that the anomalous material was ejected at a much slower velocity and at an earlier time than the majority of the gas in the bipolar nebula. This would indicate the presence of a minor ejection event prior to the main ejection of the jets first seen by Taylor et al. (1986). The variable ejection velocity model discussed here and the Echelle spectroscopy of Corradi et al. (2001) suggest the emission of material with a velocity of no more than a few 100s of km-1 just before the main ejection of the high–speed material.

The cause of the precession is unknown but, given realistic estimates for the masses and separation of the stellar components, warping of the collimating accretion disk by a magnetic white dwarf is the only mechanism that can give precession periods similar to those found by model fitting. The other possible causes give periods that are several orders of magnitude too large. The magnetic explanation also allows for the one–sided nature of the emission seen in all VLA radio maps following the initial outburst in 1984.

## Acknowledgments

MC is supported by a grant from the Particle Physics and Astronomy Research Council (PPARC). The contribution of AS was supported by a grant of the Slovak Academy of Science, number 2/1157/01. The VLA is operated by the National Science Foundation operated under cooperative agreement with Associated Universities, Inc. MERLIN is a national facility operated by the University of Manchester on behalf of PPARC. Ultraviolet spectral data are based on INES data from the IUE satellite, operated by ESA.


format: text

text:
# Global Optimization on an Evolving Energy Landscape

## Abstract

Locating the global minimum of a complex potential energy surface is facilitated by considering a homotopy, namely a family of surfaces that interpolate continuously from an arbitrary initial potential to the system under consideration. Different strategies can be used to follow the evolving minima. It is possible to enhance the probability of locating the global minimum through a heuristic choice of interpolation schemes and parameters, and the continuously evolving potential landscape reduces the probability of trapping in local minima. In application to a model problem, finding the ground–state configuration and energy of rare–gas (Lennard–Jones) atomic clusters, we demonstrate the utility and efficacy of this method.

Introduction: Global optimization problems  can often be formulated in terms of finding the minimum (or maximum) of a multidimensional potential energy surface (PES). Such problems, which occur in a variety of areas, are of considerable practical and theoretical interest . The “energy landscape”  paradigm is particularly useful when the potential energy function is continuously varying with the physical configurations relevant to the problem . An example of such a situation is the protein–folding problem , namely determining the native configuration of complex molecules given their atomic composition. A simpler variant is the determination of the ground state configuration of atomic or molecular clusters .

In this Letter we propose a new homotopy method to study such problems by a controlled deformation of the potential energy surface. If $`V_f`$ is the potential energy hypersurface under consideration, we study the landscapes

$$V(\alpha )=(1\alpha )V_i+\alpha V_f,$$
(1)
with $`\alpha `$ a parameter. Given a choice of initial potential, $`V_i`$, this is a 1–parameter family of potential energy surfaces which smoothly evolves from $`V_i`$ into $`V_f`$ as $`\alpha `$ varies from $`01`$.

The minima of the landscapes continuously change with $`\alpha `$, and in order to track them, one of two strategies are possible. Varying the interpolation parameter $`\alpha `$ in a finite number of steps, a standard technique such as conjugate gradient (CG) minimization  can be employed at each $`\alpha `$. On the other hand, one can consider $`\alpha `$ as a time–dependent function such that the PES evolves according to

$$V(t)=(1h(t))V_i+h(t)V_f,$$
(2)
where $`h(t)`$ is suitably chosen with $`h(0)=0`$, and $`lim_{tT}h(t)1`$. Over a timescale $`T`$, therefore, the potential deforms from the initial to the desired potential energy surface, and the evolving minima can be tracked, for example, by following the damped dynamics in this potential via molecular dynamics (MD) simulation.

In the present work we follow both these strategies, and show how homotopic deformation facilitates location of the global minimum in a model problem. Similar (so–called “continuation”) homotopic methods have frequently been employed in related situations, as for example in finding roots of polynomial equations in several variables  or in the mean–field dynamics in attractor neural–networks .

Different global optimization methods frequently find optimal solutions by elimination, by seeking lower and lower minima. Trapping in local minima—and escape from these minima—is a major practical issue. A number of different strategies have been suggested in order to engineer escape from local minima. These include both techniques to allow for large excursions in the phase space by the use of temperature or similar auxiliary parameters (such as simulated annealing  and its variants ) as well as methods that deform the potential energy surface. The diffusion equation method  and the distance scaling method  fall in this latter class. Other methods utilize both strategies, as for example the stochastic tunneling method  where simulated annealing is performed on a surface where the barriers are exponentially reduced so as to facilitate escape from local minima, the landscape paving technique , or the basin hopping technique  which replaces the potential surface by a set of piecewise flat regions.

The present technique is in the class of optimization methods that exploit potential surface deformation to avoid trapping in local minima. The interpolation parameter $`\alpha `$, or the switching functions $`h(t)`$ smoothly convert one PES into another. The intermediate potentials are qualitatively not very different from the asymptotic potential in terms of the number of minima and maxima, although the relative depths and curvatures are quite different. As we discuss below, this feature contributes to efficiency of the present technique in locating minima. The lowest energy achieved when an ensemble of suitably compact initial configurations is evolved is taken as the ground state prediction of this method.

Application:  The problem of minimum energy configuration determination for $`N`$ particle atomic clusters is computationally hard, and the validity of a global solution cannot, typically, be verified. Existing data for global minima  are usually the “lowest minima as yet located” in all but the simplest cases. A variety of global optimization techniques have been applied to this problem  with differing degrees of success.

For the most extensively studied such systems, namely model rare–gas clusters, the potential energy surface (PES) is an additive pairwise Lennard-Jones interaction,

$$V_f=\underset{i<j}{}V(r_{ij})=\underset{i<j}{}4ϵ[(\frac{\sigma }{r_{ij}})^{12}(\frac{\sigma }{r_{ij}})^6]$$
(3)
where $`r_{ij}`$ is the distance between particles $`i`$ and $`j`$, and $`ϵ,\sigma `$ are the standard Lennard-Jones parameters. The potential energy landscape varies greatly with cluster size. Notable difficult optimization problems in this regard are, for example, 38, 75, or 98 atom clusters, where the potential energy surface has the so–called multiple funnel structure .

In the implementation of the MD approach we proceed as follows. $`V_i`$ is taken to be a pairwise sum of harmonic terms $`V(r_{ij})=(r_{ij}2^{1/6}\sigma )^2/2`$. We perform molecular dynamics simulations  of the $`N`$ particle system, with an additional damping term for each particle,

$$m\ddot{\stackrel{}{r}_i}+\gamma \dot{\stackrel{}{r}_i}+\frac{V(t)}{\stackrel{}{r}_i}=0,i=1,\mathrm{},N$$
(4)
where $`m`$ is the mass of the particle and $`\gamma `$ is the damping coefficient. The internal timescale of interparticle vibrations depends on the parameters $`m,\sigma `$ and $`ϵ`$. For a given switching function $`h(t)`$ (we have explored a variety of such functions listed in Table I) the adiabatic time scale is set by the parameter $`\zeta `$; the entire system dynamics thus has two external time scales $`\zeta ^1`$ and $`m\gamma ^1`$. In the limit $`\gamma \mathrm{}`$, our procedure reduces to a steepest descent minimization on the evolving potential. The dynamics of the system is followed until a stationary configuration is reached.

In order to quantitatively assess the efficiency of this procedure, we define the measure

$$𝒫_g=\frac{\text{Number of ground state configurations}}{\text{Total number of condensates}},$$
(5)
a condensate being a configuration such that all atoms are within a single cluster. This is clearly a function of $`\gamma `$ and $`\zeta `$. For the ground state energy, comparison is made to the existing benchmark calculations already available for Lennard–Jones clusters .

In the CG approach, $`V_i`$ is taken to be $`\beta _{j=1}^N(\stackrel{}{r}_j\stackrel{}{r}_j^0)^2`$, $`\stackrel{}{r}_j^0`$ being the (random) initial position for the $`j`$th atom. This choice of $`V_i`$ ensures that the initial configuration is the exact global minimum for the potential energy surface, Eq. (1) with $`\alpha =0`$; $`\beta `$ is a constant that tunes the curvature of the PES. The parameter $`\alpha `$ is then varied from 0 to 1 in $`N_s`$ discrete steps; the result of the CG minimization (we follow the standard method ) at each step is taken to be the starting configuration for the CG minimization at the next value of $`\alpha `$. In this latter approach, therefore, the attempt is to allow the global minimum itself to evolve homotopically.

Results: The present application is intended to be illustrative rather than exhaustive. We have systematically studied different cluster sizes up to $`N=40`$ and in all cases the calculated ground–state energy and configuration matches existing results exactly. This includes the difficult case of the 38-atom cluster which is an interesting and important test for any optimization method . The number of minima increases exponentially with cluster size; for LJ7 there are 4 minima, while for LJ55 the number exceeds $`10^{10}`$ . Detailed results, which clarify some aspects of the present technique are presented for the cases of $`N`$=19,22 (MD) and $`N`$=38 (CG).

In the MD version of the present technique, in the absence of switching, namely in the sudden limit $`V(t)=V_f`$, the system quickly settles into the nearest available minimum based on the level of damping introduced. By starting from an ensemble of initial conditions, a variety of different minima are reached but the probability of finding the true ground state is essentially zero for large clusters. With an adiabatic switch , the results are dramatically different. The continuous evolution of the potential energy landscape is a key factor in permitting escape from local minima. Only asymptotically does the system come to rest, but until then, there is always residual kinetic energy due to which the system avoids being trapped by small barriers. Shown in Fig. 1 is the typical variation of potential energy (in units of $`ϵ`$), which is nonmonotonic once the adiabatic switching is incorporated. Regardless of the actual form of the switching, more than 85% of all initially random configurations condense, except in the case where the switching is applied to the repulsive term of the potential. Representative data is given in Table I.

As emphasized, the adiabatic optimization proposed here is heuristic. The optimal choice for the parameters $`\gamma ,\zeta `$ for a given cluster size depend on a number of features such as the interaction potential parameters and the inherent time-scales. By scanning over reasonable values of the parameters, it is possible to determine regions in parameter space with a higher than average probability of reaching the ground state. It also appears that adiabaticity is crucial since the probability of reaching the ground state increases substantially with decreasing $`\zeta `$: $`𝒫_g`$ is shown versus $`\zeta `$ for the 19–atom case in Fig. 2.

In the CG method of following minima during homotopy, the probability of reaching the global minimum is enhanced through the modification of the PES curvature. Since $`V_i`$ adds a uniform positive curvature at the intermediate stages it effectively suppresses or eliminates many barrier. To perform some benchmarking of the advantage this gives, we present, in Table II, data pertaining to finding the global minimum for LJ38 comparing the present method and the basin-hopping technique. The three lowest minima are at energies -173.928, -173.252 and -173.134 respectively. In either method, all particles are initially placed randomly inside a sphere of radius $`(\frac{N}{4})^{1/3}`$. In $`V_i`$, the parameter $`\beta `$= 100. In our implementation of the basin-hopping algorithm, coordinate displacements are random in the interval \[-0.3,0.3\] and the temperature is taken to be 2. An overall confining potential of the form $`V_c=_i\mathrm{exp}(20(r_ia)),a=1+(\frac{N}{4})^{1/3}`$ was added to prevent dissociation, and a standard Polak–Ribiere algorithm was used for the conjugate–gradient minimization  with tolerance set between 10-5 and 10-7. The average computational effort required is a product of the number of trials needed in order to get to the ground state on average and the number of function and derivative calls per trial. In our implementation of the algorithms, we find that the reduction in computational effort in locating the global minimum through the homotopy method is about 40%. The relative efficiencies can, however, vary depending on the actual choice of the various adjustable parameters in the two techniques. In either the MD or the CG version, configurations that do not reach the global minimum still typically tend to find the lowest energy states, so that a by-product of this methodology is a considerably detailed map of the low excitation regime of the cluster. This feature, however, is not unique to the present method.

Summary: We have presented here a method for global optimization which relies on the guided evolution of the underlying landscape. The methodology for finding minima on this surface can vary, and in the examples presented here, we have used both the conjugate gradient technique as well as damped molecular dynamics. (Dynamics in the landscape has been incorporated in other techniques, for example in genetic algorithms .) As in other methods, apart from the global minimum, we also obtain a detailed picture of the excitation spectrum.

Within the context of cluster geometry determination itself, several issues need to be addressed. The adiabatic method can be shown to locate ground states even when there are bifurcations along the deformation pathway . Is it possible to design more efficient homotopic deformations? What is the role of $`V_i`$ in controlling the efficiency?

The application here, though in some ways a model problem, has all the complications that arise in more general optimization problems. The success of this simple technique is therefore encouraging.

Acknowledgment: This research is supported by a grant from the Department of Science and Technology, India.

Table I: Representative results using the MD version of the homotopy method for LJ19 and LJ22 as examples of magic and nonmagic clusters. The parameters are $`\zeta =\gamma `$ = 0.5, for 1000 trials starting from random initial conditions.

|  | $`𝒫_g`$ | |
| --- | --- | --- |
| $`h(t)`$ | LJ19 | LJ22 |
| 1 (No Switching) | 0.0 | 0.0 |
| $`1\mathrm{exp}(\zeta t)`$ | 0.0081 | 0.0083 |
| $`\mathrm{sin}\frac{\pi \zeta t}{2T}`$ | 0.0065 | 0.0181 |
| $`\frac{\zeta t}{T}`$ | 0.0044 | 0.0138 |
| $`(\mathrm{tanh}(\zeta t10)+1)/2`$ | 0.0124 | 0.0176 |
| $`1\mathrm{exp}(\zeta t)\mathrm{cos}^2(3\zeta t)`$ | 0.0111 | 0.0041 |

Table II: Comparative analysis of the homotopy method and basin–hopping. For each method, $`N_r`$ initial configurations are evolved to find the global minimum in 100 instances for the LJ38 cluster. $`N_j,j=0,1,2`$ are the number of times the lowest three minima are found in the two methods; the number of function and derivative calls needed (per initial condition) are also indicated to give an estimate of the computational effort involved.

| Optimization | $`N_r`$ | $`N_0`$ | $`N_1`$ | $`N_2`$ | Function | Derivative |
| --- | --- | --- | --- | --- | --- | --- |
| Method |  |  |  |  | Calls | Calls |
| Basin | 937674 | 100 | 239 | 941 | 3495 | 154 |
| Hopping |  |  |  |  |  |  |
| Homotopy | 195690 | 100 | 4 | 173 | 9260 | 475 |
| Method |  |  |  |  |  |  |


text:
# 1 Introduction

## 1 Introduction

Supersymmetric gauge theories with four conserved supersymmetries is conjectured to have a non-perturbative self-equivalence  and has been studied for many reasons. These theories have the most symmetry of gauge theories at arbitrary couplings, including conformal invariance. The well-known holographic description of IIB string theory via $`N=4`$ gauge theory , or vice versa, generates information in the latter at the non-perturbative corner and is a well-defined example of the gravity/gauge correspondence. The construction of the S-matrix in this work has manifest S-duality and unitarity at each order in the expansion.111The S-duality invariance within the AdS/CFT correspondence was explored to all orders in , and leads to the composite operator correlations in the gauge theory. It is possible to translate information from amplitudes to composite operator correlations.

Spontaneously broken $`N=4`$ supersymmetric gauge theories have a mass parameter (and a BPS self-dual set of states) and an expansion in this variable, $`\text{ }\text{ }\text{ }\text{ }/m^2`$, is readily available. S-duality of $`N=4`$ gauge theory requires that the S-matrix of gauge bosons and composite operators be invariant via,

$`\tau {\displaystyle \frac{a\tau +b}{c\tau +d}},`$ (1.1)
with the coupling,

$`\tau ={\displaystyle \frac{\theta }{2\pi }}+i{\displaystyle \frac{4\pi }{g^2}}.`$ (1.2)
This structure, together with the perturbative coupling dependence and instantons requires that the correlations be built from a classified set of automorphic functions, Eisenstien series and cusp form(s); for a review see . One interesting feature of S-duality is that knowledge of the perturbative scattering is enough to construct the instanton series; together, these two series are redundant. S-duality also might lead to exact solutions for correlators. This approach has manifest S-duality in the expansion, and is non-perturbative.

Following along the lines of the recursion formulae in scalar field theory, we now consider amplitudes and correlators in in the broken gauge theory. The tensor manipulations and integrals have already been done in . The primary difference is that $`N=4`$ gauge theory contains a conjectured S-duality, and this symmetry of the amplitudes helps control the coupling structure. $`N=8`$ was previously examined in the derivative expansion  and new insights regarding the finiteness was found; $`N=8`$ contains an $`SO(8)`$ $`N=4`$ subtheory.

## 2 Couplings and Modular Invariance

S-duality requires the quantum generating functional of the S-matrix (gauge bosons) and correlations to be invariant under $`SL(2,z)`$ transformations. The lagrangian is,

$`={\displaystyle d^4x\mathrm{Tr}\left(F^2+\psi \psi +\varphi \text{ }\text{ }\text{ }\varphi +[\varphi ^i,\varphi ^j]^2\right)},`$ (2.1)
with gauge group G, e.g. SU(N). The toroidal modular parameter in general is a matrix,

$`\tau _{ij}={\displaystyle \frac{\theta }{2\pi }}+{\displaystyle \frac{4\pi i}{g^2}}|_{ij},`$ (2.2)
and for simplicity we take the gauge group to be $`SU(2)`$. The complications from enlarging $`SU(2)U(1)`$ to $`SU(N)U(1)^{N1}`$ is in the modular functions and trace structures involved in the sewing.

We examine the gauge boson scattering through the quantum generating functional of the S-matrix. This functional is built from the operators,

$`\mathrm{Tr}F^n,\mathrm{Tr}F^n\mathrm{Tr}F^m,\mathrm{}`$ (2.3)
together with covariant derivatives. Label the operators by the dimension, $`𝒪_{(d_i)}^i`$; the dimensionality is labeled with the expectation value of a scalar field giving the correct dimensions,

$`\varphi ^{4d_i}𝒪_{(d_i)}^i.`$ (2.4)
Because of S-duality the terms in the quantum generating functional must be modular invariant functions. These have been classified , and separable into normalizable and non-normalizable ones, and convergent versus non-convergent. The former condition is not relevant because it is tantamount to integrating over the coupling constants over the inequivalent vacua. The nonconvergent modular forms could be zeta function regularized; however, this results in a modular anomaly (an example is produced below). Any automorphic function may be represented as an infinite sum of the independent ones for complex values of $`s`$,

$`E_s^{(q,q)}={\displaystyle \underset{p,q}{}}{\displaystyle \frac{\tau _2^s}{(p\tau +q)^{sq}(p\overline{\tau }+q)^{s+q}}}`$ (2.5)
They are divergent for $`\mathrm{Re}s<1`$. The scattering of gravitons and the bounds on the planar limit of $`N=4`$ strong coupling correlators indicate that only half-integral or integral values of $`s`$ enter into the description of the scattering. For non-vanishing $`q`$ these functions transform with weight $`q`$. They satisfy the $`SL(2,Z)`$ invariant Laplacian, and potentially leads to differential relations in the terms of the S-matrix expansion.

The set of functions we consider consist of the ring,

$`{\displaystyle E_{s_i}^{(q_i,q_i)}}`$ (2.6)
with $`s_i=s=n/2`$, $`n`$ an integer, and $`q_i=0`$ with $`|q_i|s_i`$. The amplitudes containing the fermions have non-zero $`q`$ ($`N_\psi =q/2`$).

In $`k`$-space, the amplitude for four gauge bosons is an expansion,

$`{\displaystyle \mathrm{Tr}^{2n}F^4h_n(\tau ,\overline{\tau })}`$ (2.7)
with the derivatives distributed among the four field strengths (the color structures are implied, and multi-trace structures are found by taking a product of lower trace ones). The general term is,

$`{\displaystyle \underset{i=1}{\overset{p^A}{}}}{\displaystyle \underset{j=1}{\overset{n_i^A}{}}}_{\mu _{\sigma ^{}(i,j)}}{\displaystyle \stackrel{m_i^A}{}}A_{\mu _{\sigma (i,j)}}{\displaystyle \underset{i=1}{\overset{p^\varphi }{}}}{\displaystyle \underset{j=1}{\overset{n_i^\varphi }{}}}_{\mu _{\sigma ^{}(i,j)}}`$ (2.8)
$`{\displaystyle \stackrel{m_i^\varphi }{}}\varphi _{a_{\sigma (i,j)}}{\displaystyle \underset{i=1}{\overset{p^\psi }{}}}{\displaystyle \underset{j=1}{\overset{n_i^\psi }{}}}_{\mu _{\sigma ^{}(i,j)}}{\displaystyle \stackrel{m_i^\psi }{}}\psi _{\alpha _{\sigma (i,j)}}`$ (2.9)
and may be grouped into gauge invariant operators, in the functional,

$`{\displaystyle \underset{i=1}{\overset{p^A}{}}}{\displaystyle \underset{j=1}{\overset{n_i^A}{}}}_{\mu _{\sigma ^{}(i,j)}}{\displaystyle \stackrel{m_i^A}{}}A_{\mu _{\sigma (i,j)}\nu _{\stackrel{~}{\sigma }(i,j)}}{\displaystyle \underset{i=1}{\overset{p^\varphi }{}}}\left({\displaystyle \underset{j=1}{\overset{n_i^\varphi }{}}}_{\mu _{\sigma ^{}(i,j)}}\right)`$ (2.10)
$`{\displaystyle \stackrel{m_i^\varphi }{}}\varphi _{a_{\sigma (i,j)}}{\displaystyle \underset{i=1}{\overset{p^\psi }{}}}\left({\displaystyle \underset{j=1}{\overset{n_i^\psi }{}}}_{\mu _{\sigma ^{}(i,j)}}\right){\displaystyle \stackrel{m_i^\psi }{}}\psi _{\alpha _{\sigma (i,j)}}`$ (2.11)
with the color indices implied. Unitarity is built into the sewing, and may also be seen via expanding the usual Feynman diagrams at small $`k^2/m^2`$. In order to generate the unitarity cuts of the massive modes, the gauge multiplet and the dyons, one has to resum the terms in the derivative expansion. Alternatively, from the expected cuts in the amplitudes, we may impose more conditions on coefficients beyond the recursion.

The prefactors $`h`$ are invariant, and model the coupling dependence,

$`h({\displaystyle \frac{a\tau +b}{c\tau +d}},{\displaystyle \frac{a\overline{\tau }+b}{c\overline{\tau }+d}})=h(\tau ,\overline{\tau }).`$ (2.12)
Our task is to find the ring of functions that enter into the $`h`$ functions and then to utilize the sewing to build the S-matrix. In this manner, unlike the $`\varphi ^n`$ theory, the coupling structure in $`\tau `$ is simplified and also allows for a determination of the complete set of instanton corrections to the amplitudes! One of the two series, perturbative and non-perturbative, is redundant. The perturbative derivative series is invariant under S-duality. In agreement with the scattering functions, we examine functions $`h_k`$ in the ring with $`s=m/2`$, an integer or half-integer (this will be elaborated further). The Eisenstein functions admit the representation,

$`E_s^{(q,q)}={\displaystyle \underset{p,q}{}}{\displaystyle \frac{\tau _2^s}{(p\tau +q)^{sq}(p\overline{\tau }+q)^{s+q}}}`$ (2.13)
and expansion, found via a Possion resummation,

$`E_s^{(q,q)}=2\zeta (s)\tau _2^s+\alpha _s\tau _2^{1s}+𝒪(e^{\tau _2}).`$ (2.14)
The power terms are associated with perturbation theory, and the exponential series the instanton corrections. The ring $`S_{n/2}^q`$ of modular invariant functions is the set,

$`\left\{{\displaystyle \underset{j=1}{\overset{p}{}}}E_{s_j}^{(q_j,q_j)}\right\}`$ (2.15)
with $`_js_jn`$ and $`q_j=q`$. These functions converge individually for $`s>1`$. The regularized versions could be utilized, for example, $`\overline{E}_1=\mathrm{ln}\tau _2\eta (\tau ,\overline{\tau })`$, but they possess a modular anomaly as a result. The basis was already commented on, and follows from examining compatibility with the perturbative coupling series.

A n-gauge boson scattering amplitude has the coupling series,

$`\left(g^2\right)^{n/21},\left(g^2\right)^{n/2},\mathrm{},\left(g^2\right)^{n/22+n_{max}/2}.`$ (2.16)
We have labeled the maximum coupling with the integer $`n_{max}`$, as the ring of functions related to this expansion always truncate. It is impossible to generate an infinite series in $`g^2`$ with $`sn/2`$. The series of terms may be regrouped into,

$`g^{n1}\left(g^2\right)^{n_{max}/2}[\left({\displaystyle \frac{1}{g^2}}\right)^{n_{max}/2},\mathrm{},\left({\displaystyle \frac{1}{g^2}}\right)^{n_{max}/2+1}]=g^{n1}\left(g^2\right)^{n_{max}/2}h_n(\tau ,\overline{\tau })`$ (2.17)
with $`k=n_{max}`$, with $`n_{max}=n+2`$; there are no inverse powers of $`g`$ in perturbation theory. The anomalous, from the modular point of view, front factor may be removed with a field redefinition,

$`Ag^2Axgx.`$ (2.18)
Then the Lagrangian is,

$`{\displaystyle d^4x\frac{1}{g^2}\mathrm{Tr}\left(A+\frac{1}{g^2}A^2\right)^2}`$ (2.19)
and the $`g^2`$ factor in front is absorbed by the implied metric $`\sqrt{g}`$ in the coordinate redefinition. The coordinate redefinition is a slight improvement on the statement of S-duality of $`N=4`$ super gauge theory, similar to Einstein frame in IIB superstring theory.

The amplitudes have the expanded form by varying the quantum generating functional,

$`A(k_1)\mathrm{}A(k_n)={\displaystyle \underset{q}{}}h_q^{(n)}(\tau ,\overline{\tau })f_q(k_1,\mathrm{},k_n)`$ (2.20)
with $`h_q`$ elements in the finite ring $`S_{n/2+1}^{(0)}`$. The form is similar for the fermions, and $`n+2`$ is changed in general to $`n_A+n_\varphi +n_\psi /2+2`$. Via the expansion of the modular functionis, the maximum loop is $`n_{max}1`$ for the $`n`$-point function.

The four-point function in massless $`N=4`$ has the property that the $`A^4`$ term receives corrections at tree-level and one-loop only. In the spontaneously broken theory the modular ring has $`n_{max}=6`$, $`n/2+14`$. This corresponds to a truncation of the $`A^4`$ after five loops. The relevant ring contains,

$`E_3,E_{3/2}^2,|E_{3/2}^{(1,1)}|^2,E_2`$ (2.21)
$`|E_{3/2}^{1/2,1/2}|^2,|E_{3/2}^{3/2,1/2}|^2,|E_{3/2}^{+3/2,3/2}|^2`$ (2.22)
where only the convergent ones are included. These separate functions leaves unknown parameters, which are determined by the derivative expansion, up to five loops if all functions are included. This is simpler than an infinite series, as the four-point with no derivatives has only six functions. The structure is similar at higher order in derivatives and number of gauge bosons, and the instanton corrections are determined by the perturbative terms.

## 3 Sewing and Amplitudes

The sewing relation, as presented in $`\varphi ^n`$ theory, may be used to generate all of the coefficients in the derivative expansion. The approach is the same as in the scalar theory, with the added complication of many tensors. However, all loop integrals may be performed and the amplitudes are generated via algebraic computations. The parameter $`L`$ counts $`L1`$ loops in these diagrams, and for simplicity we keep the parameter as $`L`$. We have to sum over all the vertices, with $`n_A`$, $`n_\varphi `$, and $`n_\psi `$ lines, together with multi-derivatives acting on them. The general vertex is listed in the previous section.

The recursion is illustrated in figure (1) of . The left hand side takes the form, with $`m_l+m_r=n`$-point (numbers of lines to the left and right of the quantum vertex),

$`{\displaystyle \underset{L}{}}\left[{\displaystyle \underset{i=1}{\overset{m_l}{}}}\left({\displaystyle \underset{j=1}{\overset{m_i^{}}{}}}_{\mu _{\sigma (i,j)}}A_{\mu _i}(k_i)\right)\right]\left[{\displaystyle \underset{i=m_l+1}{\overset{m_r}{}}}\left({\displaystyle \underset{j=1}{\overset{\stackrel{~}{m}_i^{}}{}}}_{\mu _{\stackrel{~}{\sigma }(i,j)}}A_{\nu _i}(k_i)\right)\right]`$ (3.1)
$`\left[{\displaystyle \underset{i=1}{\overset{L}{}}}\left({\displaystyle \underset{j=1}{\overset{m_i^{}}{}}}_{\mu _{\sigma ^{}(i,j)}}{\displaystyle \underset{j=1}{\overset{m_i^A}{}}}A_{\mu _i}\right)\right]\left[{\displaystyle \underset{i=1}{\overset{L}{}}}\left({\displaystyle \underset{j=1}{\overset{\stackrel{~}{m}_i^{}}{}}}_{\mu _{\stackrel{~}{\sigma }^{}(i,j)}}{\displaystyle \underset{j=1}{\overset{\stackrel{~}{m}_i^A}{}}}A_{\stackrel{~}{\mu }_i}\right)\right]\times t^{\mu _\sigma ,\mu _i,m_i^{}}t^{\stackrel{~}{\mu }_\sigma ,\stackrel{~}{\mu }_i,\stackrel{~}{m}_i^{}}`$ (3.2)
$`+\mathrm{permutations},`$ (3.3)
with $`m_i^A=\stackrel{~}{m}_i^A=L`$; the $`\tau `$ is within the tensor $`t`$, as well as the coefficients of the derivative structure. We have to sum over all possible combinations, as iterated from lower to higher orders in derivatives and numbers of external states. Internal fermions and scalars must be included in the sum.

The right-hand side of the equation contains the quantum vertices,

$`t^{\mu _\sigma ,\mu _i,m_i^{}}{\displaystyle \underset{i=1}{\overset{p^A}{}}}{\displaystyle \underset{j=1}{\overset{n_i^A}{}}}_{\mu _{\sigma ^{}(i,j)}}{\displaystyle \stackrel{m_i^A}{}}A_{\mu _{\sigma (i,j)}}`$ (3.4)
The integrals have to be performed in the former. Contractions of the internal fields generate the propagators, and due to the tensor structure the integrals are more complicated. However, all internal momenta or derivatives may be extracted from within the integral with the use of the identity,

$`^\mu \mathrm{}^L=L\left[(2d)+m^2_m^2\right]\times ({\displaystyle \frac{k^\mu }{k^2}})\mathrm{}^L.`$ (3.5)
The massive propagator is

$`\mathrm{}_{\mu \nu }=(x^2)^{d/2+1}K_{d/2}(mx)P_{\mu \nu },`$ (3.6)
in terms of the modified Bessel function and the tensor structure (gauge dependent) of the propagator. This identity is easily proved in x-space, and its application to the gauge field correlations is, upon iteration,

$`{\displaystyle e^{ixk}\underset{j=1}{\overset{n}{}}_{\mu _j}\mathrm{}^L}={\displaystyle \underset{\sigma ,\sigma ^{}}{}}{\displaystyle \eta _{\mu _\sigma \nu _\sigma }k_{\mu _\sigma ^{}}(\frac{1}{k^{2n_1+2n_2}})(2)^{n_21}\times \left[(2d)+m^2_m^2\right]^n\mathrm{}^L}.`$ (3.7)
We have to sum over all pairs $`(n_1,n_2)`$ such that $`n_1+n_2=n`$ and all partitions of $`n`$ to the $`\eta `$ and $`k`$’s. The gauge boson scattering is,

$`\left[{\displaystyle \underset{i=1}{\overset{L}{}}}\left({\displaystyle \underset{j=1}{\overset{m_i^{}}{}}}_{\mu _{\sigma ^{}(i,j)}}{\displaystyle \underset{j=1}{\overset{m_i^A}{}}}A_{\mu _i}\right)\right]\left[{\displaystyle \underset{i=1}{\overset{L}{}}}\left({\displaystyle \underset{j=1}{\overset{\stackrel{~}{m}_i^{}}{}}}_{\mu _{\stackrel{~}{\sigma }^{}(i,j)}}{\displaystyle \underset{j=1}{\overset{\stackrel{~}{m}_i^A}{}}}A_{\stackrel{~}{\mu }_i}\right)\right]`$ (3.8)
$`={\displaystyle \underset{\sigma ,\sigma ^{}}{}}{\displaystyle \eta _{\mu _\sigma \nu _\sigma }k_{\mu _\sigma ^{}}(\frac{1}{k^{2n_1+2n_2}})(2)^{n_21}\times \left[(2d)+m^2_m^2\right]^n\mathrm{}^L\times \underset{\rho (i),\stackrel{~}{\rho }(j)}{}\eta _{\mu _{\stackrel{~}{\rho }(i)}\nu _{\stackrel{~}{\rho }(j)}}}.`$ (3.9)
The tensor structure follows from the deriavtive diagram, and the last term comes from the Pfaffian,

$`{\displaystyle \underset{i}{\overset{L}{}}}A_{\mu _i}{\displaystyle \underset{j}{\overset{L}{}}}A_{\stackrel{~}{\mu }_j}.`$ (3.10)
The expression may seem complicated, but all loop integrals have been performed. The algebraic complexity is also less than typical loop diagram calculations, which at intermediate steps may contain millions of terms.

The integrals are evaluated to,

$`{\displaystyle d^dx}e^{ixk}\mathrm{}^L=(k^2)^{d/2L(d/21)}{\displaystyle \underset{n}{}}(k^2/m^2)^n\alpha _n^{(L)}`$ (3.11)
in dimensional regularization, and

$`(k^2)^{d/2L(d/21)}{\displaystyle \underset{m,n}{}}({\displaystyle \frac{k^2}{\mathrm{\Lambda }^2}})^m({\displaystyle \frac{k^2}{m^2}})^n\alpha _{m,n}^{(L)}`$ (3.12)
in a momentum cutoff scheme. $`N=4`$ is finite, and divergences at intermediate steps in the calculation should drop out in the final equations.

Equating the former with the quantum vertices generates a recursion in the couplings and tensor $`ttt`$, with the coupling structure found from the ring of Eisenstein functions up to number coefficients. The products of Eisenstein function iterate from $`n`$ and $`m`$-point to the higher derivative terms in $`t`$ (together with the enlarged basis in the latter). The iteration formulae allow for a computation of the scattering amplitude via polynomial equations, with no integrals. Similar results are available for the correlations of the composite operators.

## 4 Discussion

The scattering amplitudes of spontaneously broken $`N=4`$ supersymmetric gauge theory have been examined in the derivative expanson. The gauge theory may be iteratively constructed via sewing, and is equivalent to the sewing in the Feynman diagram expansion. Unitarity is obvious in this construction. Amplitudes may be constructed without any integrals; the integrals involved are of the free-field type and may all be integrated. Unitarity is obvious in this construction.

Scalar and gauge field theory has also been examined in ; S-duality has been implemented in the $`N=4`$ gauge theory and order by order the gauge boson expansion is invariant under $`SL(2,Z)`$ duality. S-duality in the amplitudes has the advantage that from the perturbative structure, all of the instanton corrections are computable and follow from the former. The instanton expansion via the terms, $`e^{2n\tau _2}`$, are connected to the perturbative terms and are redundant; the perturbative terms generate an instanton calculus via an $`SL(2,Z)`$ completion.

The derivative expansion, and the concise expressions which follow from it, should enable further progress several areas. The intermediate coupling regime is accessible as well as some possible information regarding the verification of duality.

## Acknowledgements

GC thanks the DOD, 444025-HL-25619, for support.