text:
# Ongoing Mass Transfer in the Interacting Galaxy Pair NGC 1409/10 1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract No. NAS5-26555.

## 1 Introduction

Galaxy interactions are clearly linked to bursts of star formation, and in more restricted ways, to triggering of nuclear activity. What physical mechanisms mediate these connections remains unclear, since there are multiple well-motivated and plausible processes. Among these, mass transfer between the galaxies in an interacting system has been long discussed, but evidence for its occurrence has remained largely circumstantial. Such transfer could play a particular role in dumping gas into the inner regions of a galaxy, since the angular-momentum barrier which restricts the inward transport of disk gas may be greatly reduced if the gas comes from another galaxy with the appropriate encounter geometry.

The conditions for mass transfer have been examined analytically as well as numerically. Sotnikova (1990) considered the fate of dense ISM clouds crossing between galaxies and heated by both conduction and radiation in a hot IGM, finding that the clouds could survive long enough to reach the second galaxy, and that a stream of clouds from galaxies in low-eccentricity orbits could reach a quasi-steady state configuration. Her numerical simulations (Sotnikova 1988) indicated that this situation is reached much faster when the spiral spin and orbit are parallel than the antiparallel case, in roughly an orbital timescale of the binary system.

Wallin & Stuart (1992) presented an extensive study of mass transfer from restricted $`n`$-body simulations, including over 1000 encounter geometries to map the systematics of mass transfer. The largest amounts finally captured by the second galaxy are found for encounters with the orbital plane close to the donor’s disk plane, with this favoritism enhanced by the strong dependence on orientation of periapsis for non-planar encounters. The intuitive increases in mass transfer for closer periapsis and larger mass ratio also hold. These simulations do not give much detail on how the companion captures the mass, but should be robust as to the totals bound in its potential well. These results fit with patterns in the 86 simulations described by Howard et al. (1993), where the particles were tracked but there was no modelled companion structure to interact in detail with the incoming particles.

The observational signatures of mass transfer can be hidden in the complex kinematics and star-formation behavior of gas-rich interactions and mergers. Thus, it has long been recognized that such effects are easiest to see in pairs of mixed morphology, since large amounts of cool gas in the early-type E/S0 members of these systems very likely originated in the gas-rich companion (as set forth by Domingue et al. 2003). This has been particularly discussed in the context of emission lines, far-infrared emission, and blue (post-starburst) optical colors for these galaxies, as possible evidence of past gas transfer (de Mello et al. 1995, 1997). This would also account for the Holmberg effect, a correlation between colors of pair members, even for early-type members (Holmberg 1958, Demin et al. 1984, Madore 1986, Reduzzi & Rampazzo 1995). Marziani et al. (2003) present kinematic evidence of ongoing or recent dumping of gas in the pair Arp 194, and stress that mass transfer itself is necessary but not sufficient for cross-fuelling of nuclear activity or starbursts.

An unusually clear example of current mass transfer appeared in the HST WFPC2 snapshot survey of AGN and starburst galaxies by Malkan, Gorjian, & Tam (1998). Their data for NGC 1409 showed a dust lane crossing between NGC 1409 and its close companion (and Seyfert galaxy) NGC 1410, also known as III Zw 55 (although there has been some ambuguity on the literature as to which member of the pair is which). This feature extends even beyond the region shown in Fig. 1 of Malkan et al. (where the galaxy identified there as the Seyfert galaxy NGC 1410 is in fact its earlier-type and more quiescent companion NGC 1409). Such cases highlight the value of dust absorption as a tracer of the interstellar medium, also exploited by several groups analyzing HST imagery of the nuclei of early-type galaxies (Martel. et al. 2000, Tran et al. 2001, Martini et al. 2003).

The extensive tidal envelope of stars around NGC 1409/10 provides backlighting to trace this feature and estimate the total mass and transfer rate involved. To allow such a study, two-color STIS imagery was obtained with higher signal-to-noise in each band. The results are described here, along with mapping in H$`\alpha `$ and \[N II\] which constrains the encounter geometry and timing through the velocity field.

We see the NGC 1409/10 pair deep into a strong encounter, as shown by the presence of a substantial stellar envelope extending well beyond the two galaxy cores, to a radius at least 15 kpc. The nuclei are projected at a separation of only 7 kpc, suggesting that there may have been significant tidal loss of stars from both galaxies. The spiral NGC 1410 in particular is rather small in linear size, with the nearly face-on spiral pattern spanning a diameter of 5.3 kpc (for $`\mathrm{H}_0=72`$ km s-1 Mpc -1 following the WMAP analysis of Spergel et al. 2003). This galaxy hosts a Seyfert nucleus and active disk star formation. In contrast, NGC 1409 is a barred early-type system, best described as type SB0. Because of the small angular extent of these structures, these types are not obvious from typical ground-based images (such as shown by Keel 1996).

## 2 Observations

### 2.1 STIS imaging

The STIS CCD imaging mode was used to observe NGC 1409/10 in two passbands, the “white-light” 50CCD and red F28X50LP bands. For study of extinction properties, the improved surface-brightness sensitivity compared to WFPC2, and fine angular sampling of this instrument, outweighed the nonstandard and overlapping passbands available with its filter complement. In this range of surface brightness, for single-orbit exposures split in two for cosmic-ray rejection, a STIS F28X50LP image has signal-to-noise ratio per pixel 2.5 times greater than would be obtained with the WFPC2 PC CCD and F606W filter, a combined result of improved sensitivity and broader passbands. Even using the WF chips and a pixel scale twice as coarse gives the advantage per pixel to STIS by about 1.8 times, plus its finer sampling of the PSF (slightly better than the Nyquist value across most of the red passband). This improvement does come with accepting a smaller field of view than the entire WFPC2 outline, and of variations in the PSF across the field. From the cycle 13 STIS instrument handbook (Kim Quijano et al. 2003), the fraction of encircled flux within a one-pixel radius varies by an extreme of 16% around the CCD field of view, with 5% of this variation occurring within the region analyzed here. This effect translates into a change in the pixel-to-pixel contrast of the dust features, one which is small compared to the photon noise in this application. PSF changes of similar amplitude also occur with WFPC2; the encircled energy within a one-pixel radius varies by 5% rms across its field.

The very broad STIS passbands also introduce the possibility of confusion between starlight and line emission, since multple strong lines from photoionized material are included. This confusion is not an important issue for these data unless there is a close connection between the ionized gas and dust even far from plausible ionizng sources. Beyond the Seyfert nucleus of NGC 1410, the strongest line in equivalent width is H$`\alpha `$, reaching values of 30 Å in the inner disk of NGC 1410. Folding a typical set of H II region lines ratios through the STIS response functions indicates that this level of line emission would contaminate local measures of the starlight at the 3% level in F28X50LP and less for the 50CCD mode. Larger contamination would require that there be substantial emisison-line structure within the aperture size of the WIYN spectra used for this scaling, and that such structure not be reflected in the starlight. Since there are no potential ionizing clusters seen near the extended dust feature, and the equivalent widths in the extended emission occurring in this region are much smaller than in the star-foring disk of NGC 1410, emission-line contamination is unlikely to be an issue in these results.

The white-light image (in the 50CCD mode) is shown in Fig. 1, to illustrate the geometry of the NGC 1409/10 pair, the extensive diffuse envelope of scattered stars, and the third companion spiral seen to the north of the strongly interacting pair. The stellar envelope provides smooth background light to trace the dust between and beyond the two pair members.

The observations, sequence root o5ev01, include full-orbit (2767 and 2815-second) sets of exposures, two in each filter for cosmic-ray rejection. The filter change introduces a small scale change between the two images, evaluated from the locations of 11 unresolved or compact objects. The angular size of the red F28X50LP pixels is larger than the white-light 50CCD pixels by a factor $`1+7.9\times 10^4`$, amounting to 0.8 pixel across the frame. This was corrected for differential analysis by transforming the red image to the coordinate system of the white-light data.

To analyze the extinction properties of the dust feature, a blue-light image was constructed as a scaled difference between the white- and red-light data sets. This allowed the use of the very sensitive STIS CCD for multicolor measurements, since the exact shape of the blue passband is less important than that this shape be known. The preflight estimates of the relative throughputs of the two filters in their overlap region had to be substantially revised, since using the initial throughput values led to unphysical negative flux in the difference. The F28X50LP filter was designed to replicate the white-light 50CCD configuration with only a grey shift longward of a sharp cutoff to the blue. To the extent that this was realized in fabrication, the count rate at each pixel may be expressed in terms of the white-light rate $`C_W`$ and that in the red F28X50LP filter $`C_R`$ to give a differential blue rate $`C_B`$ according to

$$C_W=C_B+C_R/f$$
$`(1)`$
where $`f`$ is the mean transmission ratio of red and white-light passbands longward of the cutoff wavelength. The in-flight calibration, now incorporated into the STSDAS $`crrefer`$ files, gives $`f=0.873`$ averaged across the 6000–9500 Å range. This fits with astrophysical constraints from the range of color and direction of color gradient in the early-type galaxy, which show that $`f`$ must lie in the range 0.83–0.93. Figure 2 shows the effective passbands of the two filters and their scaled difference.

An image in the blue passband was produced following equation (1). To avoid confusion with standard photometric passbands, the original and derived images will be referred to simply as white-light, red, and blue. To show the overall structure of the system, the central portion of the STIS field is shown in figure 3, comparing the blue and red images.

Some analysis techniques can make more appropriate use of the observed white and red passbands, since they are statistically independent and free of the correlated noise that can result from measuring subtle features in a difference image. Similarly, some points will be made more clearly using the disjoint passbands of the red and blue images. To allow approximate scaling of the extinction results from these broad passbands, the effective extinction was evaluated using a model old stellar population (appropriate for NGC 1409 and probably for most of the stellar envelope as well), in photon units. These very broad bands have the potential for nonlinear extinction behavior (that is, significant changes in effective wavelength with reddening). In practice, these numerical tests indicate that such shifts are still small for extinctions $`A_V<1`$ and slopes near that of a normal Galactic extinction law ($`A_V/E_{BV}=3.1`$). The blue passband behaves very much like $`V`$; its extinction ranges from $`1.010.97A_V`$ over the range $`A_V=04.0`$, varying in a nearly linear fashion. The red passband has extinction $`0.580.54A_V`$ over the same range, again changing almost linearly with $`A_V`$. The much broader white-light band shows more pronounced color behavior with increasing reddening, with extinction from $`0.680.59A_V`$ for $`A_V=04.0`$ in a roughly quadratic way, as the redder part of the passband becomes more important at large extinctions. However, the extinctions which matter in this work are small enough that these are not major complications.

### 2.2 Spectroscopic mapping

To help understand the kinematics and interaction history of NGC 1409/10, the velocity field in H$`\alpha `$ and \[N II\] emission was measured using the 3.5m WIYN telescope and Dense-Pak fiber array. As described by Barden, Sawyer, & Honeycutt (1998), the array includes a $`7\times 13`$ configuration of fibers in a roughly hexagonal packing covering a $`35\times 45`$” region, set for these observations with the long axis at PA 0. Four outlying fibers allow sky subtraction far from the galaxy centers. During 6/7 December 2000, two 30-minute exposures were obtained in each of four positions, with 60-minute exposures at each position obtained on the following night when the transparency was more stable. A final 60-minute exposure was obtained with an offset of 15” E and 16” N, to cover the fainter third galaxy just east to NGC 1410. This galaxy, centered at (2000) $`\alpha =`$03h 41m 11.60s, $`\delta =01^{}`$ 17’ 40$`\stackrel{}{\mathrm{.}}`$4 from the STIS coordinates, ended up not being well centered in any of the fibers. The fiber location nearest its core does show features at about the $`3\sigma `$ level that match H$`\alpha `$ and \[N II\] emission at $`cz=7808\pm 28`$ km s-1 (heliocentric), which is at least plausible given the values of 7710 and 7596 for NGC 1409 and 1410 respectively, and the evidence that the envelope of tidally stripped stars from the interaction encompasses this galaxy as well as the two brighter ones.

The positions centered on NGC 1409/10 were dithered in a parallelogram pattern, with each leg offset by about 2”, to fill the gaps between 3” fiber apertures. The instrumental resolution with these fibers was 1.6 Å , well sampled by the 0.68-Å pixel scale. The spectral range observed was 6000-7400 Å . Velocity maps were produced from the H$`\alpha `$ and \[N II\] $`\lambda 6583`$ measures; for each observation, emission was detected in 21–28 fibers, for a total of 198 velocity measures from the higher-quality 60-minute exposures. The velocity maps were constructed on a 1” grid, with overlapping aperture data averaged at each pixel and numerical 0.1” subpixels used to track aperture outlines until the final averaging. Registration to the direct images used reconstruction of continuum images from the spectral measurements, giving positions of the galaxy nuclei in the DensePak coordinate system. The velocity field is shown in Fig. 4.

The nucleus of M32 was observed, in four exposures of 10 minutes each in a single fiber, to provide a reference old-population spectrum. These data are used in deriving a cross-correlation velocity scale and in subtracting a typical bulge contribution from the nucleus of NGC 1409.

## 3 The NGC 1409/10 Encounter: Geometry, Timing and Velocity Scale

The WIYN spectroscopic results give a velocity field (at least for ionized gas) showing the relative velocities and senses of rotation for NGC 1409 and 1410. The northeastern side of NGC 1410 is approaching, with the line of nodes in approximately position angle $`140^{}`$. For NGC 1409, the disk rotation is less clear, since much of the gas may be associated with the dust lane or other tidal debris. The velocity field in Fig. 4 shows the gradient across NGC 1409 approximately along the minor axis. The absorption-line results from these spectra are limited, but do show a gradient across the nucleus in the sense that the northerwestern side is receding. For the sense of orbital motion, the H$`\alpha `$ results and the central absorption-line velocity for NGC 1409 agree is showing NGC 1409 to be receding relative to the center of mass, with NGC 1410 approaching. The faint tails of starlight extending beyond the STIS field (see the inset to Fig. 1) indicate the rotation senses of the disks in the plane of the sky, since tidal tails act as a kind of material arm and wind up with the sense of disk rotation. The long and prominent one-sided tail wrapping north of NGC 1409 shows that the disk must share this sense of rotation in projection (clockwise), fitting with a polar passage, and putting the eastern side of the disk in front. Together, the morphology and kinematics suggest that the orbital plane is viewed from its eastern side, and not too nearly edge-on. The northern tidal tail also indicates that there was a previous close passage, through the disk plane of NGC 1409.

If the spiral pattern in NGC 1410 is trailing, the western side is in the foreground. Only the tidal tail provides compelling evidence as to which orientation about the plane of the sky the disk of NGC 1409 has, since its internal structure shows parts of a ring rather than spiral features. NGC 1409 probably lies behind NGC 1410 from our perspective, since the dust feature crosses in front of NGC 1409, something that would require very special geometry if it were to be seen far from the plane of the sky.

The observed velocity difference between nuclei is 181 km s-1. There are indirect arguments suggesting that we see this system within about 30 of the orbital plane, so that this is a modest underestimate of the relative velocity. NGC 1409 is the receding member. The emission-line velocity field does follow this trend; since dust and gas are generally associated, the dust is likely to show this same behavior, which fits with its apparently becoming bound to NGC 1409.

The dominant grand-design pattern in NGC 1410 is easiest to understand if this disk is experiencing a near-planar direct (prograde) encounter, which favors the orbital plane being viewed from its eastern side. This would be a geometry most favorable to loss of mass from its disk during the encounter. In contrast, NGC 1409 is undergoing a roughly polar encounter, which would favor the acquisition of material into a polar ring.

Fig. 5 shows a sketch of the inferred geometry of the disks and relative orbit, incorporating the constraints from velocity field, central velocities, and tidal structure. To put a timescale to the encounter and mass-transfer event, the characteristic orbital time (that is, assuming circular orbits) for the two galaxies can be estimated from $`T=2\pi D/v`$ for a deprojected linear separation $`D`$ and space-velocity difference $`v`$. Including the orbital inclination to the line of sight $`i`$ and phase angle $`\varphi `$ along the orbit, measured from the plane of the sky, gives $`T=2D_{}\mathrm{cos}i\mathrm{cos}\varphi /\mathrm{\Delta }v(1\mathrm{cos}^2i\mathrm{sin}^2\varphi )^{1/2}`$ in terms of the projected linear separation $`D_{}`$ and observed velocity difference $`\mathrm{\Delta }v`$ (e.g. Karachentsev 1987). The observed $`\mathrm{\Delta }v=181`$ km s-1 suggests that the correction for projection is not large, since for even more luminous pairs, the projected velocity differences are greater than this only 20% of the time (from Karachentsev 1987). Taking the typical value of 250 km s-1 as the maximum for galaxies at this combined luminosity (following Karachentsev, in his table 5), the projection effects suggest $`\mathrm{cos}i\mathrm{cos}\varphi >0.72`$, so that $`T=1.72.7\times 10^8`$ yr. The small projected separation of this pair, and the extensive halo of diffuse starlight, indicate that the members are well into the spiralling stage preceding a final merger, so that this circular-orbit timescale is only a rough scaling value.

## 4 An intergalactic dust bridge

The NGC 1409/10 system shows a striking lane of dust, crossing between the galaxies. The connection to the spiral pattern of NGC 1410 and polar wrapping around NGC 1409 make it clear that the direction of ISM transfer is from the spiral NGC 1410 to the early-type disk of NGC 1409.

The extent and structure of the dust bridge in the NGC 1409/10 system may be most simply shown in median-windowed images and color maps. Fig 6 (top) shows a median-windowed image, obtained upon dividing the original white-light image by a version median-filtered over a 2” box. From this display, the dust lane clearly crosses continuously between the galaxies, probably connecting to a feature that crosses the southern spiral arm of NGC 1410 and curving so as to suggest that it passes behind most of the starlight in NGC 1409. This linkage might then reappear as the near-polar dust structure seen closer to the nucleus of NGC 1409. This image also highlights non-axisymmetric structure in the smooth disk of NGC 1409, in the form of a narrow stellar bar and connected arcs. The dust lane between the galaxies has only moderate observed attenuation, with a maximum light loss in this area of 13%. It is not completely smooth, with multiple filaments and offsets particularly apparent as it curves north of NGC 1409. The dust lane passes in front of the northern part of NGC 1409, especially well seen where the dust crosses the bar feature. It may pass behind the southern arm of NGC 1410. The median-windowed images in each passband were the starting point for generating maps of residual intensity in the dust lanes, where the median-windowed images were interpolated roughly perpendicular to the dust lanes to make a background model, retaining the ratio of observed to model images.

Many of the same features appear in a color map (Fig. 6), obtained as the simple ratio of blue to red images. The signal-to-noise ratio in the dust bridge is reduced, but structures near the center of NGC 1409 are traced much more clearly since changes in the background illumination do not affect the color. In front of NGC 1409 there are two distinct dust features, one with a braided appearance wrapping across the bar, and a narrower lane crossing almost exactly in front of the nucleus.

To the extent that the dust properties and the distribution of starlight, in particular color gradients, are known, the combination of color and intensity information can be used to reconstruct the three-dimensional location and intrinsic optical depth of dust features. I neglect scattering into the observed beam at this stage, since the optical depth in the NGC 1409/10 features is modest, so that observational error dominates over this neglect. Let a pixel have residual intensity $`R`$ with respect to the unabsorbed starlight, as determined either by simple interpolation or modelling of the entire image structure. When a single, spatially resolved dust feature is responsible for the flux deficit, this implies a relation between its optical depth in the observed passband and its depth within the starlight distribution. If $`I(r)`$ is the luminosity density of starlight at location $`r`$ along the line of sight, the dust is located at the position $`r_0`$ given implicitly by

$$R=\frac{_0^{r_0}I(r)𝑑r+e^\tau _{r_0}^{\mathrm{}}I(r)𝑑r}{_0^{\mathrm{}}I(r)𝑑r}$$
$`(2)`$
This already rules out certain locations as unphysical, formalizing the immediate constraints that a dust feature must have optical depth at least great enough to give the observed flux deficit, and that it cannot lie any deeper within the starlight distribution than an opaque absorber which gives the appropriate foreground:total ratio. From equation (2), rearrangement shows that

$$\tau =\mathrm{ln}[\frac{_{r_0}^{\mathrm{}}I(r)𝑑r}{R_0^{infty}I(r)𝑑r+_0^{r_0}I(r)𝑑r}].$$
$`(3)`$
Data at multiple wavelengths can refine this relation further; if the ratio of optical depths at the two wavelengths is known for the dust, a further constraint is that the ratio of derived $`\tau `$ must match the reddening curve assumed for the dust, which specifies a corresponding spatial location $`r_0`$ as well.

A similar approach was used by Shaya et al. (1996) in reconstructing the locations of dust patches in NGC 1316. The most robust results, such as front versus backside locations and derived optical depth, do not depend on the model for the stellar luminosity distribution, since the integral ratio gives the relative depth of the absorbing structure into the starlight. This application may be illustrated for the special case of no color gradient and small $`\tau `$, so that color terms due to finite passbands can be neglected. Then the residual intensity $`R_\lambda `$ at each wavelength constrains the fraction of starlight in front of the dust at its optical depth $`\tau _\lambda `$ through

$$R_\lambda =X+e^{k\tau _\lambda }(1X)$$
$`(4)`$
where the extinction curve enters through $`k`$, the ratio of $`\tau _\lambda `$ to the value at a fiducial wavelength such as the $`V`$ band. An example using two bands is shown in Fig. 7, where error bounds for 1% precision in residual intensity $`R_\lambda `$ are shown restricting the possible values of $`X`$ and $`\tau _B`$ from data in the $`B`$ and $`R`$ passbands.

Use of these expressions should be done at high enough spatial resolution to avoid the averaging effects of dust structures within the resolution element, which artificially flattens the extinction curve and introduces nontrivial weighting of unresolved structures by their transmission. WFPC2 imaging of overlapping galaxy systems shows that linear resolution of a few tens of parsecs substantially reduces these effects (Keel & White 2001a,b). The STIS pixel size corresponds to 26 pc (at 105 Mpc from $`\mathrm{H}_0=72`$ km s-1 Mpc-1), satisfying this condition.

For disk dust in luminous spirals, there is evidence that the slope of the optical extinction curve is close to the typical Galactic value. Among the four backlit spirals studied using HST imagery by Keel & White (2001a,b), the only one showing a significant departure from a Galactic extinction curve is the late-type object silhouetted in front of NGC 1275, where the unusual environment may have altered the grain population. Dust in the arms of NGC 2207 gives similar results from HST images, albeit with significant errors from the less favorable geometry (Elmegreen et al. 2001). In the nearby spiral M31, star-by-star spectroscopy indicates that the optical behavior is much like the Milky Way (rather than, for example, the SMC slope associated with its lower metallicity), as found by Bianchi et al. (1996). These observations motivate the use of a galactic extinction law where we don’t have evidence to the contrary. However, the reddening behavior in the extended features of NGC 1409/10 may not fit this expectation.

The multicolor comparison of reddening and extinction can, in principle, refine the extinction and mass estimates of the dust bridge in NGC 1409/10. It can also suggest whether the this bridge is spatially continuous with the other polar features around NGC 1409, if the fading of the dust extinction is accompanied by a greying trend which suggests that it is also wrapping behind most of the starlight. In practice, the reddening curve is steeper than expected for the STIS effective passbands and a Galactic extinction law, limiting what we can learn about its spatial location from the photometric results alone. This can be shown using the region where the dust passes in front of the bar of NGC 1409, and is thus plainly in the foreground of most of the starlight (with the difference in $`X`$ between here and the middle of the dust feature roughly a factor 9 from relative starlight intensities). For a given depth of the dust feature into the starlight and modest values of $`\tau `$, changing $`\tau `$ for different places in the dust feature will give a roughly linear locus in the color-residual intensity plane (Fig. 8). The slope of points in this part of the dust lane is inconsistent with “normal” reddening curves ($`A_V/E_{BV}=3.1`$), implying more reddening per unit extinction (or that the STIS passbands have not been correctly evaluated). The slope expected for a Galactic relation is the ratio of $`k`$ (as in equation 4) for the two bands, about 0.93 in this case. A linear relation was fitted to the data from this foreground region (region 1), incorporating the typical photometric errors from pixel scatter in both coordinates and using the procedure from Press et al. (1992) as implemented in the GSFC IDL library, giving a slope $`0.800\pm 0.019`$. The dust lane between galaxies and backlit by the diffuse envelope light (region 2) gives a slope $`0.835\pm 0.012`$. If the slope from region 1 represents the actual dust properties for these bands, the slope is equal to the ratio of $`k`$ values in the limit of small $`\tau `$. For finite $`X`$, the locus of values of $`R_1`$ and $`R_2`$ for changing $`\tau _1`$, as would be found in a single structured dust feature, follows a family of curves whose mean slope in a given range of $`R_1`$ steepens with $`X`$. This happens as nonlinear structure in the $`R\tau `$ relation occurs at higher values of $`R`$ for dust seen deeper into the starlight. In this case, the best fit occurs for $`X=0.4\pm 0.1`$, a slight improvement over the intuitive conclusion that dust is most visible for $`X<0.5`$.

This information then allows a correction to the apparent extinction, giving typical values of $`A_V`$ and hence column density. A map of $`R`$ at either wavelength becomes a map of $`\tau `$ according to

$$e^\tau =\frac{RX}{1X}$$
$`(5)`$
(and ideally the same for each wavelength observed to within errors of measurement and elimination of structure in the starlight). Integration over the resulting map of $`\tau `$ can give characteristic dust masses and column densities. In this case, I follow Domingue et al. (1999) in using a mix of graphite and silicate grain properties, which yields a column density of $`1.11\times 10^{10}\tau `$ in g cm-2. The mean map for the main dust lane as converted to $`\tau _B`$ was formed by averaging results from the white and red images. Region 1 is unaffected by this scaling, being in front of nearly all the starlight. The largest extinctions, scaled to $`\tau _B`$, reach 0.55 in region 2 and 0.35 in region 1 (where the extinction is visually more dramatic because of the brighter background intensity). This conversion gives a dust mass of $`7\times 10^5`$ solar masses in region 2, with about $`2\times 10^5`$ in region 1.

Additional dust features are of interest to the question of whether the dust wraps behind NGC 1409 and reappears as the polar lanes in one continuous structure. For region 3, where the lane appears faintly in the median-windowed image and might be turning behind much of the starlight, the $`R_WR_R`$ relation is still steeper than for region 2. How much steeper is poorly determined, since the fit has a slope formally greater than unity. This would then be farther back than the reference curve $`X=0.7`$ in Fig. 8. That would give this part of the lane similar optical depth (and mass) to region 2 and support the continuity of dust around NGC 1409. The structure in front of the central region of NGC 1409 is in front of the bulk of the starlight ($`X<<1`$) by inspection. It has $`\tau _B0.4`$ at its deepest points, much like regions 1 and 2. These results are consistent with the suspicion that there is a single continuous dust pane crossing from the disk of NGC 1410 and wrapping in front of the northwestern edge of NGC 14109, turning behind it, and reappearing in front on the western side near the nucleus to form a polar structure. The total dust mass involved is about $`2\times 10^6`$ solar masses, including the unseen portion which would lie behind the disk of NGC 1409. For a typical gas-to-dust ratio of 160 by mass (e.g. Sodroski et al. 1994), the accompanying gas phase material would total $`3\times 10^8`$ solar masses. The time spent crossing along this path would be of the same order as the orbital timescale, implying a mean mass-transfer rate of $`1.11.4`$ solar masses per year.

## 5 Star formation

The members of this pair differ as much in star formation rate as in morphology and color. NGC 1410 shows substantial recent star formation, including a population of luminous blue star clusters (Figs. 1 and 6). Of these, 10 are isolated enough for simple aperture photometry. Roughly transforming to the $`V`$ band gives them observed magnitudes $`22.924.4`$, or absolute magnitudes -12.2 – -10.7. These are luminous but by no means unprecedented, falling within the range of the brightest clusters seen in more local interacting systems (Whitmore et al. 1999, Keel & Borne 2003) in which only a handful of clusters appear brighter than $`M_V=12`$. Clusters are short-lived at this luminosity (a few $`10^7`$ years), indicating that their formation is ongoing at the timescale of the interaction.

Similarly, the line emission from the disk of NGC 1410 (away from the Seyfert nucleus) indicates brisk star formation. Equivalent widths in H$`\alpha `$ from 10–30 Å are found throughout the region within about 6” of the core (and excluding the active nucleus itself). This corresponds to a flux of $`7.6\times 10^{14}`$ erg cm-2 s-1, or an H$`\alpha `$ luminosity of $`9\times 10^{40}`$ erg s-1. Using the formulation by Kennicutt (1983) for a Salpeter initial mass function, this corresponds to a star-formation rate of 0.8 solar masses per year in the spiral disk of NGC 1410.

These points are in stark contrast to the lack of either associated line emission or stellar clusters in NGC 1409 (since even the brightest old globular clusters would be just at the threshold of these observations). The line emission showin in Fig. 4 is kinematically decoupled from the rotating disk of NGC 1409, and is therefore suspect as an indicator of any UV stellar flux. The broadband color, smooth image, and characteristic old stellar features in the spectrum of NGC 1409 limit its overall star-formation rate to only a few per cent of what we seen in NGC 1410. Whatever the fate of the gas reaching NGC 1409, it is not fuelling star formation at the epoch we observe. A star-formation rate of less than 0.05 solar mass per year would stand out in these data.

## 6 Nuclear activity

NGC 1410 is a well-known type 2 Seyfert galaxy, and as a member of such a strongly-interacting pair, has been included in many samples of AGN designed to test for links between nuclear activity and interaction. It is therefore ironic that in this clear instance of mass transfer, it is the Seyfert galaxy which acts as the donor, with a much more quiescent galaxy receiving the material. Nuclear activity in NGC 1409 is much weaker if present. The WIYN spectrum of its nucleus is compared to NGC 1410, with and without correction for the starlight continuum, in Fig. 9. The continuum was subtracted using M32 as a template, rebinned to match redshifts and broadened by a Gaussian of FWHM=11.4 Å (520 km s-1), the best fit for the combined velocity dispersion and central rotation gradient in the NGC 1410 spectrum. NGC 1409 shows typical LINER emission with \[N II\] $`\lambda 6583`$/H$`\alpha `$ near unity. Such emission, while often associated with other signs of genuinbe nuclear activity at low luminosity, is common enough not to be a particular indication of anything special happening in the core of NGC1409.

## 7 Summary

HST STIS images in two broad passbands have been used to trace a dust lane marking mass transfer in the interacting galaxy pair NGC 1409/10. The combined color and extinction behavior support the impression that there is a single feature crossing the 7-kpc gap between NGC 1410 and NGC 1409, then wrapping behind NGC 1409 and reappearing on the opposite side to become a small lane crossing its disk over the pole. The depth into the starlight derived from the two-color data and the residual intensity yield estimates of the dust mass, totalling about $`2\times 10^6`$ solar masses and likely accompanied by $`3\times 10^8`$ solar masses of gas. At the characteristic orbital velocities in this system, this suggests a mean rate of mass transfer slightly above one solar mass per year onto NGC 1409.

This mass-transfer calculation may be too simple, since the feature may be more exactly described as an isochrone than a pipeline. It is striking that only in a small region has material been launched away from NGC 1410, perhaps indicating that a special location in the system was needed to start such transfer. This fits with the fact that the dust lane appears to cross the spiral pattern of NGC 1410.

In light of the mechanisms discussed to accont for star formation and nuclear activity in interacting galaxies, it is ironic that the Seyfert galaxy NGC 1410 is the donor, rather than the recipient, of mass transfer. The recipeint NGC 1409 shows no signs of any fate of the gas reaching its disk. Limits to its rate of star formation stand at a few per cent of the inflow rate, and its nuclear activity is limited to a modest LINER. Either the flow has yet to actually reach the dense inner regions, or some additional trigger must set in after a critical density or mass of infalling gas is achieved. If flows such as these drive starbursts in interacting systems, such a cycle is needed to drive star-formation rates an order of magnitude greater than the estimated flow rate in NGC 1409/10.

This work was supported by NASA thought STScI grant STScI GO-8147.01-97A. Charles Proffitt helped in resolving the STIS filter throughput issues.


format: text

text:
# The sextic oscillator as a 𝛾-independent potential

## I Introduction

In the last few years there has been considerable interest in looking for analytic solutions of the Bohr Hamiltonian which describes the collective motion in nuclei in terms of shape variables ($`\beta `$, $`\gamma `$) bm . This was initiated with the introduction of the interacting boson model (IBM) ibm . This model present four different dynamical symmetries each one associated to a well defined nuclear shape. The model Hamiltonian provides a natural way of going from a phase to another one by changing systematically few parameters and, consequently allows to study shape phase transitions. Thus, the IBM phase diagram has been analyzed from different points of view Die ; Feng ; Cas ; Jo1 ; Ar3 . Specially important are the critical points since in these situations structural changes occur rapidly and is difficult design appropriate theoretical models. Recently, Iachello has proposed a new dynamical symmetry for describing the critical point at the transition from spherical to deformed $`\gamma `$-unstable shapes fi00 . This symmetry has been called E(5) and is expected to occur in nuclei in which the $`V(\beta ,\gamma )`$ potential depends only on the $`\beta `$ variable, and it has a relatively flat shape in the $`\beta `$ variable. Around this critical point of the phase transition the potential appearing in the Bohr Hamiltonian can be approximated by an infinite square well in the $`\beta `$ variable. In this case the Bohr Hamiltonian can be solved exactly in terms of Bessel functions, and various quantitative predictions can be obtained for various spectroscopic properties (ratios of the excitation energies and $`B(\mathrm{E2})`$ transitions), on the basis of which one can search for candidates for the E(5) symmetry among nuclei. Lately, other two dynamical symmetries X(5) fi01  and Y(5) fi03  to describe the critical point at the transition from spherical to axially deformed shapes and from axially deformed to triaxial shapes, respectively, have been proposed.

The introduction of the E(5) symmetry renewed the interest in studying exactly solvable $`V(\beta ,\gamma )`$ potentials. Most efforts have concentrated on $`\gamma `$-independent potentials, for which the most well-known example is the harmonic oscillator in the five-dimensional space wj56  and its extension which contains a term proportional to $`\beta ^2`$ elliott . More recently two more exactly solvable $`\gamma `$-unstable potentials have been discussed: the Coulomb and the Kratzer potentials fortunato . Here again the latter potential is an extension of the former in the sense that it contains a term proportional to $`\beta ^2`$, which formally changes the $`\tau `$ variable, the analogue of the $`l`$ angular momentum of radial potentials in three spatial dimensions. This formal changing of $`\tau `$ introduces a minimum of the potentials at $`\beta >0`$ in both cases elliott ; fortunato . The bound solutions of the harmonic oscillator and the Coulomb potentials (and their extensions) are given in terms of generalized Laguerre polynomials as70 .

Similarly to the three-dimensional case, these examples, together with the infinite square well, practically exhaust those exactly solvable potentials that contain a $`\beta ^2`$ term. Here we propose another potential which has this property, and although it is not exactly solvable in the classical sense, it has a number of features that make it an ideal potential to be used in the Bohr Hamiltonian. This is the sextic oscillator, which belongs to the class of quasi-exactly solvable (QES) potentials qes . These potentials have the property that their solutions can be obtained in closed form for a number of energy eigenvalues, i.e. for the lowest few values of the $`n`$ principal quantum number, which are also the lowest in energy. This is clearly sufficient for potentials appearing in the Bohr Hamiltonian: it is rarely necessary to consider more than a few levels with the same angular momentum, and these can be obtained from the lowest few solutions of the sextic oscillator. Furthermore, the sextic oscillator has a more flexible shape than other solvable potentials, as depending on its parameters, it can have a minimum at $`\beta =0`$ or at $`\beta =\beta _{\mathrm{min}}>0`$, and in addition, it can also have a local maximum at $`\beta _{\mathrm{max}}<\beta _{\mathrm{min}}`$.

The paper is structured as follows. In Sect. 2 the Bohr Hamiltonian is revised together with its solutions for $`\gamma `$-independent potentials. In Sect. 3 the lowest energy solutions of the Bohr Hamiltonian for the sextic oscillator potential are worked out. Section 4 is devoted to show the simple use of the sextic oscillator as a flexible $`\gamma `$-independent potential. Finally, in Sect. 5 we present preliminary applications and discuss future extensions.

## II The Bohr Hamiltonian for $`\gamma `$-independent potentials

Let us first consider the Bohr Hamiltonian describing the collective motion of a deformed nucleus in the five-dimensional space determined by the $`\theta _i`$ Euler angles ($`i=1,2,3`$) and the intrinsic $`\beta `$ and $`\gamma `$ variables is bm

$$H=\frac{\mathrm{}^2}{2B}\left(\frac{1}{\beta ^4}\frac{}{\beta }\beta ^4\frac{}{\beta }+\frac{1}{\beta ^2\mathrm{sin}3\gamma }\frac{}{\gamma }\mathrm{sin}3\gamma \frac{}{\gamma }\frac{1}{4\beta ^2}\underset{k}{}\frac{Q_k^2}{\mathrm{sin}^2(\gamma \frac{2}{3}\pi k)}\right)+V(\beta ,\gamma ).$$
(1)
In what follows we assume that the potential in (1) depends only on $`\beta `$, i.e. $`V(\beta ,\gamma )=U(\beta )`$. For these $`\gamma `$-independent potentials the wavefunctions can be separated into two parts

$$\mathrm{\Psi }(\beta ,\gamma ,\theta _i)=f(\beta )\mathrm{\Phi }(\gamma ,\theta _i),$$
(2)
which satisfy the following differential equations:

$`\left({\displaystyle \frac{1}{\mathrm{sin}3\gamma }}{\displaystyle \frac{}{\gamma }}\mathrm{sin}3\gamma {\displaystyle \frac{}{\gamma }}+{\displaystyle \frac{1}{4}}{\displaystyle \underset{k}{}}{\displaystyle \frac{Q_k^2}{\mathrm{sin}^2(\gamma \frac{2}{3}\pi k)}}\right)\mathrm{\Phi }(\gamma ,\theta _i)=\mathrm{\Lambda }\mathrm{\Phi }(\gamma ,\theta _i),`$
$`\mathrm{\Lambda }=\tau (\tau +3),\tau =0,1,2,\mathrm{}`$ (3)
$$\left(\frac{1}{\beta ^4}\frac{}{\beta }\beta ^4\frac{}{\beta }+\frac{\mathrm{\Lambda }}{\beta ^2}+u(\beta )\right)f(\beta )=ϵf(\beta ).$$
(4)
Here we have introduced $`ϵ=\frac{2B}{\mathrm{}^2}E`$ and $`u(\beta )=\frac{2B}{\mathrm{}^2}U(\beta )`$. Note that the $`\tau `$ values determine the allowed angular momenta $`J`$ too bes . By setting $`\varphi (\beta )=\beta ^2f(\beta )`$ we obtain an equation which has the form of a radial Schrödinger equation

$$\frac{\mathrm{d}^2\varphi }{\mathrm{d}\beta ^2}+\left(\frac{(\tau +1)(\tau +2)}{\beta ^2}+u(\beta )\right)\varphi =ϵ\varphi .$$
(5)
Note that this is different from the Eq. (6) in fi00  in that it contains no linear derivative term due to the different definition of $`\varphi (\beta )`$. This choice also implies that the factor corresponding to the $`\beta `$ volume element in the integration of functions of the type $`f(\beta )`$ in (4) has been transferred to the solutions of the type $`\varphi (\beta )`$ in (5). Thus, in the integration of these no factor arising from the volume element has to be included. The complete solution of the problem implies the solution of Eq. (3) too, this was solved in Ref. bes .

## III The sextic oscillator

The sextic oscillator with a centrifugal barrier is defined qes  as

$$H=\frac{\mathrm{d}^2}{\mathrm{d}x^2}+\frac{(2s1/2)(2s3/2)}{x^2}+\left(b^24a(s+\frac{1}{2}+M)\right)x^2+2abx^4+a^2x^6,$$
(6)
where $`x[0,\mathrm{})`$ and $`M`$ is a non-negative integer. This potential is quasi-exactly solvable, which means that for any non-negative integer value of $`M`$, $`M+1`$ of its solutions can be obtained in an algebraic way. The (unnormalized) solutions are written in the form

$$\varphi _n(x)=P_n(x^2)(x^2)^{s\frac{1}{4}}\mathrm{exp}\left(\frac{a}{4}x^4\frac{b}{2}x^2\right)n=0,1,2,\mathrm{}$$
(7)
where $`P_n`$ is a polynomial of order $`n`$. Obviously, normalizability requires $`a0`$, while $`a=0`$ reduces the problem to the exactly solvable harmonic oscillator.

The simplest solutions are obtained for $`M=0`$ and $`M=1`$ qes . For $`M=0`$ only one nodeless (i.e. ground-state) solution appears at $`E_0^{(M=0)}=4bs`$, with the corresponding wavefunction being

$$\varphi _0^{(M=0)}(x)(x^2)^{s\frac{1}{4}}\mathrm{exp}\left(\frac{a}{4}x^4\frac{b}{2}x^2\right).$$
(8)
For $`M=1`$ two solutions appear, one nodeless, and another with one node for $`x>0`$. These correspond to the ground-state and the first excited state, respectively, at energies $`E_0^{(M=1)}=4bs+\lambda _{}(s)`$ and $`E_1^{(M=1)}=4bs+\lambda _+(s)`$, where

$$\lambda _\pm (s)=2b\pm 2(b^2+8as)^{1/2}$$
(9)
are the roots of the equation $`\lambda ^24b\lambda 32as=0`$. The corresponding wavefunctions are

$$\varphi _n^{(M=1)}(x)\left(1\frac{\lambda }{8s}x^2\right)(x^2)^{s\frac{1}{4}}\mathrm{exp}\left(\frac{a}{4}x^4\frac{b}{2}x^2\right),$$
(10)
and the $`\lambda =\lambda _{}(s)`$ and $`\lambda =\lambda _+(s)`$ choice has to be made for $`n=0`$ and $`n=1`$, respectively qes . (Note that $`a0`$ and $`s0`$ imply $`\lambda _{}(s)0`$, so the polynomial part of (10) is nodeless.) It has to be mentioned that the solutions for $`M=0`$ and $`M=1`$ belong to different sextic potentials if $`s`$ is the same, as the coefficient of the quadratic term is different then. We shall see, however, that with appropriate combinations of $`s`$ and $`M`$ it is possible to solve sextic potentials that differ only in the strength of the centrifugal term.

The normalization of the wavefunctions can also be given in closed form. For this one has to evaluate integrals of the type

$`I^{(A)}`$ $``$ $`{\displaystyle _0^{\mathrm{}}}x^A\mathrm{exp}\left({\displaystyle \frac{a}{2}}x^4bx^2\right)`$ (11)
$`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Gamma }\left({\displaystyle \frac{A+1}{2}}\right)a^{\frac{A+1}{4}}\mathrm{exp}\left({\displaystyle \frac{b^2}{4a}}\right)D_{\frac{A+1}{2}}\left({\displaystyle \frac{b}{a^{1/2}}}\right)`$
$`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Gamma }\left({\displaystyle \frac{A+1}{2}}\right)(2a)^{\frac{A+1}{4}}U({\displaystyle \frac{A+1}{4}},{\displaystyle \frac{1}{2}};{\displaystyle \frac{b^2}{2a}}),`$ (12)
where $`D_p(z)`$ is the parabolic cylinder function and $`U(\alpha ,\beta ;z)`$ is one of the forms of the confluent hypergeometric function as70 .

Larger values of $`M`$ can also be considered (e.g. for $`M=2`$ three different solutions are obtained for the three roots of a cubic algebraic equation for $`\lambda `$), but $`M=0`$ and $`M=1`$ are sufficient for our purposes in this paper. A complete study including more solutions, explicit closed forms for the normalization factors and applications to actual nuclei is underway.

## IV Application as a $`\gamma `$-independent potential

In order to cast (6) in a form similar to (5) we have to write $`x=\beta `$ and $`s=\frac{\tau }{2}+\frac{5}{4}`$ (remember that $`\tau 0`$). In order to keep the quadratic term at a constant value (once the $`a,b`$ parameters are fixed) we also have to prescribe

$$s+M+\frac{1}{2}=\frac{1}{2}\left(\tau +2M+\frac{7}{2}\right)c=const.$$
(13)
With this the sextic oscillator Hamiltonian can be brought to the form of (4) with $`u(\beta )`$ being

$$u^\pi (\beta )=(b^24ac^\pi )\beta ^2+2ab\beta ^4+a^2\beta ^6+u_0^\pi ,$$
(14)
where the index $`\pi =\pm `$ is included to distinguish the potential for even/odd $`\tau `$’s which is slightly different as explained below. In Eq. (14) $`c^\pi `$ are the constants obtained in Eq. (13) for even/odd values of $`\tau `$ and we have introduced a constant $`u_0^\pi `$ for convenience as it will be discussed below.

Equation (13) implies that increasing/decreasing $`M`$ in one unit has to come with decreasing/increasing $`\tau `$ in two units. Thus, once the values of the $`(a,b)`$ parameters are fixed, the sequence of $`(M,\tau )`$ values $`(K,0),(K1,2),(K2,4),\mathrm{}`$ correspond to solutions of Eq. (14) with $`c^+=\frac{7}{4}+K`$. In the same way, the sequence of $`(M,\tau )`$ values $`(K,1),(K1,3),(K2,5),\mathrm{}`$ correspond to solutions of Eq. (14) with $`c^{}=\frac{9}{4}+K`$. Consequently, the potential for $`\tau `$-even and $`\tau `$-odd states is slightly different due to the fact that the coupling coefficient $`b^24ac^\pm `$ of the quadratic term is different in the two cases due to the choices for $`c^+`$ and $`c^{}`$ that are necessary to separate the $`(\tau +1)(\tau +2)\beta ^2`$ term in a uniform way. This situation can be handled using different strategies. One possibility is setting $`u_0^+=u_0^{}=0`$ and considering $`b^2>10a`$, which minimizes the deviation of the quadratic terms compared to the quartic and sextic terms. Another possibility is introducing a relative energy shift between the $`\tau `$-even and $`\tau `$-odd potentials by setting $`u_0^+`$ and $`u_0^{}`$ such that the potential minima are at the same energy. We shall discuss this possibility after analyzing qualitatively the spectrum and the potential shapes.

Let us analyze the spectrum obtained in a simple case. Taking $`M=1`$ we obtain, as explained in the preceding section, two solutions for $`\tau =0`$ (one with no nodes, $`n=0`$, and the other one with one node, $`n=1`$). In the notation introduced in Ref. fi00  the label $`\xi `$ is our $`n+1`$. Thus the two solutions of Eq. (10) with $`s=5/4`$ are $`\varphi _n^{(M=1)}(\beta )`$ with $`n=0`$ and $`1`$ and correspond to $`\varphi _{\xi ,\tau }=\varphi _{1,0}`$ and $`\varphi _{\xi ,\tau }=\varphi _{2,0}`$ respectively in Ref. fi00  notation. Inspecting the energy eigenvalues, the corresponding energies are $`E_1(M=1,\tau =0)=E_{1,0}=7b2(b^2+10a)^{1/2}+u_0^+`$ and $`E_2(M=1,\tau =0)=E_{2,0}=7b+2(b^2+10a)^{1/2}+u_0^+`$, respectively. The same potential is obtained by taking $`M=0`$ and $`\tau =2`$. In this case we have a single solution, Eq. (8) with $`s=9/4`$, $`\varphi _0^{(M=0)}(\beta )`$ with no nodes and with energy $`E_1(M=0,\tau =2)=9b+u_0^+`$. This corresponds to $`\varphi _{\xi ,\tau }=\varphi _{1,2}`$ in the notation of Ref. fi00 . A similar analysis can be performed for the solutions with odd $`\tau `$-values. For $`M=1`$ there are two solutions with $`\tau =1`$, Eq. (10) with $`s=7/4`$, which correspond to $`\varphi _{\xi ,\tau }=\varphi _{1,1}`$ and $`\varphi _{\xi ,\tau }=\varphi _{2,1}`$ for $`n=0`$ and $`n=1`$ respectively. The corresponding energy eigenvalues are $`E_{1,1}=9b2(b^2+14a)^{1/2}+u_0^{}`$ and $`E_{2,1}=9b+2(b^2+14a)^{1/2}+u_0^{}`$. Again the same potential is obtained for $`M=0`$ and $`\tau =3`$. In this case there is a single solution with no nodes, Eq. (8) with $`s=11/4`$, which corresponds to $`\varphi _{1,3}`$ and has an energy $`E_{1,3}=11b+u_0^{}`$. In Fig. 1 a schematic spectrum is shown with indication of the relevant quantum numbers. In Fig. 2 the corresponding wavefunctions with the notation $`\varphi _{\xi ,\tau }`$ are presented.

Now we analyze the different potential shapes that can be produced by different election of parameters in Eq. (14). From (14) we find that the shape of the potential $`u^\pi (\beta )`$ depends on the sign of $`b^24ac^\pi `$ and $`b`$, which set the coefficients of the quadratic and quartic terms. (The coefficient of the leading sextic term is always positive.) When $`b^2>4ac^\pi `$ and $`b>0`$ hold (i.e. for $`b>2(ac^\pi )^{1/2}`$), the potential has a minimum at $`\beta =0`$ and it increases monotonously with $`\beta `$. When $`b^2<4ac^\pi `$, irrespective of the sign of $`b`$ (i.e. for $`2(ac^\pi )^{1/2}<b<2(ac^\pi )^{1/2}`$), a minimum appears for $`\beta >0`$, while for $`b^2>4ac^\pi `$ and $`b<0`$ (i.e. for $`b<2(ac^\pi )^{1/2}`$), first a maximum appears and then a minimum as $`\beta `$ increases. In all three cases the exact location of the extremal point(s) can be obtained from the real and positive solutions of

$$(\beta _0^\pi )^2=\frac{1}{3a}[2b\pm (b^2+12ac^\pi )^{1/2}].$$
(15)
Due to the relatively small difference in $`c^+`$ and $`c^{}`$, the $`\tau `$-even and $`\tau `$-odd potentials have the same types of extrema at about the same $`\beta `$, except for some peculiar combinations of $`a`$ and $`b`$. Assuming that there are no complications of this kind, we can now return to the question of renormalizing the minima of the $`\tau `$-even and $`\tau `$-odd potentials. For $`b>2(ac^\pi )^{1/2}`$, $`\pi =+,`$ the minima of the two potentials will be $`u_0^+`$ and $`u_0^{}`$ at $`\beta =0`$, so they coincide if $`u_0^+=u_0^{}`$ holds. For $`b<2(ac^\pi )^{1/2}`$ we can equate the minima of $`u^+(\beta )`$ and $`u^{}(\beta )`$ if we set $`u_0^+=0`$ and

$$u_0^{}=(b^211a)(\beta _0^+)^2(b^213a)(\beta _0^{})^2+2ab[(\beta _0^+)^4(\beta _0^{})^4]+a^2[(\beta _0^+)^6(\beta _0^{})^6],$$
(16)
where the $`\beta _0^\pi `$ are obtained from (15) with the choice of the “+” sign. With this the two potentials have their minima at the same energy, but they take on different values at the origin. Illustrations of the possible potential shapes are displayed in Fig. 3. Obviously, $`u_0^{}`$ in (15) also has to be added to the energies of the $`\tau `$-odd potential. Fig. 4 show the relative position of the energy levels $`E_{\xi ,\tau }`$ as either $`a`$ or $`b`$ is varied and the other parameter is kept at a fixed value.

The electric quadrupole transition rates can also be determined analytically by calculating the matrix elements of the transition operator wj56 ; fi00

$$T^{(\mathrm{E2})}=t\alpha _{2\mu }=t\beta \left[D_{\mu ,0}^{(2)}\mathrm{cos}\gamma +2^{1/2}(D_{\mu ,2}^{(2)}+D_{\mu ,2}^{(2)})\mathrm{sin}\gamma \right].$$
(17)
The radial integrals that appear in the $`\beta `$ variable in the matrix elements of $`T^{(\mathrm{E2})}`$ can again be determined using (11). In order to obtain the total matrix elements, one has to calculate also the components depending on $`\gamma `$ and the Euler angles $`\theta _i`$. This can be done following the techniques described in Ref. bes . These parts introduce certain selection rules not only for the angular momenta, but also for $`\tau `$.

## V Discussion

In order to compare the main characteristics of the sextic oscillator as a $`\gamma `$-unstable potential with those of other potentials of this kind, we present calculations for a particular value of the parameters, $`a=40000`$ and $`b=200`$. These numbers were chosen such that the resulting energy spectrum approximates that of the 134Ba nucleus, the first candidate for E(5) symmetry casten . We stress that our aim is not to reproduce the experimental data, rather to get a qualitative picture about the general performance of the model. The potentials $`u^\pm (\beta )`$ are displayed in the middle panel of Fig. 3, while the energy eigenvalues are shown in Fig. 5, together with the corresponding experimental energy levels. Fig. 5 also shows the calculated and the experimental $`B(\mathrm{E2})`$ values for transitions between the energy levels. Note that electric quadrupole transitions which change $`\tau `$ with more than one unit are zero if we use the transition operator (17), but finite $`B(\mathrm{E2})`$ strengths can be obtained if we apply terms of the next order (see e.g. frank ).

In Table 1 we summarize the ratio of the most important energy eigenvalues and those of the most characteristic $`B(\mathrm{E2})`$ transition rates obtained from the sextic oscillator with parameters $`a=40000`$, $`b=200`$, the infinite square well potential fi00  and the numerically solved $`\beta ^4`$ potential beta4  together with the corresponding experimental values for 134Ba, whenever available. It is seen that the energy ratios corresponding to the E(5) symmetry systematically fall between the values of the $`\beta ^4`$ potential and the sextic oscillator. The situation is less obvious for the ratio of the $`B(\mathrm{E2})`$ values: here the sextic oscillator and the infinite square well seem to yield similar ratios, while the numbers obtained from the $`\beta ^4`$ potential are systematically higher. This might be due to the fact that the sextic oscillator potential goes to infinity steeper than the $`\beta ^4`$ potential, so the asymptotic behaviour of its wavefunctions can be closer to that of the wavefunctions of the infinite square well. Comparing the results with the experimental data for 134Ba we can conclude that, at least in this case, the sextic oscillator allows a better approximation than the other essentially parameter-free potentials. We expect that this conclusion will be general due to the flexible nature of the sextic potential whose shape is governed by two parameters. In fact this potential can be used not just at the critical point but it can be useful to model the full shape phase transition from spherical to deformed $`\gamma `$unstable nuclei by changing the parameters $`a`$ and $`b`$.

Before closing, we mention some aspects of the sextic oscillator that might give further help in the analysis of nuclei near critical points. First we note that with $`M=2`$ in (6) the analysis can be extended to further states, such as $`\varphi _{1,4}`$, $`\varphi _{1,5}`$, $`\varphi _{2,2}`$, $`\varphi _{2,3}`$, $`\varphi _{3,0}`$ and $`\varphi _{3,1}`$. Second, the potential shape which contains both a local maximum and a minimum at $`\beta >0`$ might be useful in the description of nuclei with the so-called X(5) symmetry, which is thought to occur in the shape phase transition between the spherical and the axially deformed domain fi01 . Third, there are further quasi-exactly solvable potentials both with confining and non-confining nature qes , which can also be considered in the Bohr Hamiltonian.

###### Acknowledgements.

This work was supported by the OTKA grant No. T37502 (Hungary) and by the Spanish MCyT under project No. BFM2002-03315


text:
# Nonperturbative renormalisation of composite operators with overlap quarksPresented by J. B. Zhang

## 1 INTRODUCTION

Lattice QCD is a unique tool with which we can compute physical observables non-perturbatively from first principles. Renormalisation of lattice operators is an essential ingredient needed to deduce physical results from numerical simulations. In this paper we study the renormalisation properties of composite bilinear operators with the overlap quark action.

In principle, the renormalisation of a quark bilinear can be computed by lattice perturbation theory. However, lattice perturbation theory converges slowly and the higher-order corrections may not be small, thus introducing a large uncertainty in the calculation of the renormalised matrix elements in some continuum scheme. To overcome these difficulties, Martinelli et al.  have proposed a promising non-perturbative renormalisation procedure. The procedure allows a full non-perturbative computation of the matrix elements of composite operators in the Regularisation Independent (RI) scheme . The matching between the RI scheme and $`\overline{MS}`$, which is intrinsically perturbative, is computed using only continuum perturbation theory and has been carried out to higher orders, which is well behaved. This method has been shown to be quite successful in reproducing results obtained by other methods, such as using chiral Ward Identities (WI) . The method has also been successfully applied to determine renormalisation coefficients for various operators using the Wilson  and staggered actions , domain-wall fermions , as well as the quark mass renormalisation constant for overlap fermions . The purpose of the current work is to study the application of this non-perturbative renormalisation procedure to the renormalisation of the quark field and the flavour non-singlet fermion bilinear operators in the case of overlap fermions.

We shall use Neuberger’s overlap fermions  which have lattice chiral symmetry at finite cut-off. As a result, many chiral-symmetry relations are still valid  and the quark propagator  preserves the same structure as in the continuum. The use of the overlap action entails many theoretical advantages , such as no additive mass renormalisation, no $`𝒪(a)`$ error, and no mixing among operators of different chirality. The latter is very helpful for computing weak matrix elements.

## 2 NON-PERTURBATIVE RENORMALIZATION METHOD

In this section we review the nonperturbative renormalisation method of Ref. , which we use to compute the renormalisation constants of quark bilinears in this paper. The method imposes renormalisation conditions non-perturbatively, directly on quark and gluon Green’s functions, in a fixed gauge, for example the Landau gauge.

Let $`S(x,0)`$ denotes the quark propagator on a gauge-fixed configuration from a source $`0`$ to all space-time points $`x`$. The momentum space propagator is defined as the discrete Fourier transform over the sink positions

$$S(p,0)=\underset{x}{}\mathrm{exp}\left(ipx\right)S(x,0).$$
(1)
We use periodic boundary conditions in spatial directions and an anti-periodic boundary condition in time direction. The dimensionless lattice momenta are

$$p_i=\frac{2\pi }{N_s}(n_i\frac{N_S}{2}),p_t=\frac{2\pi }{N_t}(n_t\frac{1}{2}\frac{N_t}{2})$$
(2)
for an $`N_s^3\times N_t`$ lattice, where $`n_i`$ ($`n_t`$) may in principle lie in the range $`0N_s(N_t)`$. In practice, however, only a subset of this range is used. In this paper the momentum range was restricted to those momenta for which $`n_i=0,1,2`$ and $`n_t=1,2,3,4`$.

We also define the square of the absolute momentum,

$$\left(pa\right)^2=a^2\underset{\mu =0}{\overset{3}{}}p_\mu p_\mu .$$
(3)
### 2.1  Three Point Function

Consider the two-fermion operators

$$O_\mathrm{\Gamma }(x)=\overline{\psi }(x)\mathrm{\Gamma }\psi ,$$
(4)
where $`\mathrm{\Gamma }`$ is the Dirac gamma matrix

$$\mathrm{\Gamma }\{1,\gamma _\mu ,\gamma _5,\gamma _\mu \gamma _5,\sigma _{\mu \nu }\},$$
(5)
and the corresponding notation will be {S, V, P, A, T} respectively. The three point function with the operator insertion at position $`0`$ and the propagator from $`y`$ to $`0`$ and then from $`0`$ to $`x`$ is given by

$`G_O(x,0,y)`$ $`=`$ $`\psi (x)O_\mathrm{\Gamma }(0)\overline{\psi }(y)`$ (6)
$`=`$ $`S(x,0)\mathrm{\Gamma }S(0,y),`$
where $`S(0,y)`$ is the quark propagator from $`y`$ to $`0`$. It is the inverse of the Dirac operator.

The Fourier transform of the three-point function is given by

$`G_O(pa,qa)`$ $``$ $`{\displaystyle d^4xd^4ye^{i(pxqy)}G_O(x,0,y)}`$ (7)
$`=`$ $`S(p,0)\mathrm{\Gamma }\left(\gamma _5S^{}(q,0)\gamma _5\right).`$
From this, one can define the amputated three-point function

$$\mathrm{\Lambda }_O(pa,qa)=S(pa)^1G_O(pa,qa)S(qa)^1,$$
(8)
where

$$S(pa)=S(p,0),$$
(9)
which is translational invariant and a $`12\times 12`$ matrix.

Finally, we define a projected vertex function

$$\mathrm{\Gamma }_O(pa)=\frac{1}{\mathrm{Tr}(\widehat{P}_O^2)}\mathrm{Tr}\left(\mathrm{\Lambda }_O(pa,pa)\widehat{P}_O\right),$$
(10)
where $`\widehat{P}_O=\mathrm{\Gamma }`$ is the corresponding projection operator.

### 2.2 RI/MOM Renormalisation Condition

The renormalised operator $`O(\mu )`$ at scale $`\mu `$ is related to the bare operator

$$O(\mu )=Z_O(\mu a,g(a))O(a),$$
(11)
and the renormalisation condition is imposed on the vertex function $`\mathrm{\Gamma }_O(pa)`$ at a scale $`p^2=\mu ^2`$,

$$\mathrm{\Gamma }_{O,\mathrm{ren}}(pa)|_{p^2=\mu ^2}=\frac{Z_O}{Z_\psi }\mathrm{\Gamma }_O(pa)|_{p^2=\mu ^2}=1,$$
(12)
to make it agree with the tree-level value of unity. Here $`Z_\psi `$ is the field or wavefunction renormalisation

$$\psi _{\mathrm{ren}}=Z_\psi ^{1/2}\psi .$$
(13)
In order to obtain the renormalisation constant $`Z_O`$ for the operator $`O`$, one needs to know $`Z_\psi `$ first. It has been suggested that $`Z_\psi `$ be obtained from the renormalisation of vector or axial-vector currents. For example, if one uses the conserved vector current, then one expects $`Z_{V^C}=1`$. Therefore, from the renormalisation condition Eq. (12), one obtains,

$$Z_\psi =\frac{1}{48}\mathrm{Tr}\left(\mathrm{\Lambda }_{V_\mu ^C}(pa)\gamma _\mu \right)|_{p^2=\mu ^2}.$$
(14)
However, in this work, we will obtain $`Z_\psi `$ directly from the quark propagator. It can be defined from the Ward Identity (WI) as

$$Z_\psi =i\frac{1}{12}\mathrm{Tr}\left[\frac{S(pa)^1}{/p}\right]|_{p^2=\mu ^2}.$$
(15)
To avoid derivatives with respect to a discrete variable, we have used

$$Z_\psi ^{}=i\frac{1}{12}\frac{\mathrm{Tr}_{\mu =1,4}\gamma _\mu (p_\mu a)S(pa)^1}{4_{\mu =1,4}(p_\mu a)^2}|_{p^2=\mu ^2},$$
(16)
which, in the Landau Gauge, differs from $`Z_\psi `$ by a finite term of order $`\alpha _s^2`$ , and the difference is less than 1% at the typical scale $`\mu `$ 2 GeV. It has been pointed out in Ref.  that the systematic error due to the definition of lattice momentum $`p`$ is much larger than 1%. We will not distinguish $`Z_\psi ^{}`$ and $`Z_\psi `$ later on.

## 3 NUMERICAL RESULTS

We work on a $`16^3\times 32`$ lattice with lattice spacing, $`a`$=0.093 fm. The gauge configurations are created using a tadpole improved plaquette plus rectangle (Lüscher-Weisz ) gauge action through the pseudo-heat-bath algorithm. A total of 50 configurations are considered. The lattice parameters are summarised in Table 1. The lattice spacing is determined from the static quark potential with a string tension $`\sqrt{\sigma }=440`$ MeV .

The gauge field configurations are gauge fixed to the Landau gauge using a Conjugate Gradient Fourier Acceleration  algorithm with an accuracy of $`\theta \left|_\mu A_\mu (x)\right|^2<10^{12}`$. We use an improved gauge-fixing scheme  to minimise gauge-fixing discretisation errors.

The massive overlap Dirac operator  is defined so that at tree-level there is no mass or wavefunction renormalisation ,

$$D(m_0)=\rho +\frac{m_0}{2}+(\rho \frac{m_0}{2})\gamma _5ϵ(H).$$
(17)
Here $`ϵ(H)=H/\sqrt{H^2}`$ is the matrix sign function and $`H`$ is taken to be the hermitian Wilson-Dirac operator, i.e., $`H=\gamma _5D_\mathrm{W}`$. Here $`D_\mathrm{W}`$ is the usual Wilson fermion operator, but with a negative mass parameter $`\rho =1/2\kappa 4`$ in which $`\kappa _c<\kappa <0.25`$.

Our numerical calculation begins with an evaluation of the inverse of $`D(m_0)`$ for each gauge configuration in the ensemble. We use a 14th order Zolotarev approximation  to the sign matrix function $`ϵ(H_w)`$, and in the selected window of $`x[0.031,2.5]`$ of $`ϵ(x)`$, the approximation is better than $`3.3\times 10^{10}`$ . In the calculations, $`\kappa =0.19163`$ is used, which gives $`\rho a=1.391`$. We calculate 15 quark masses by using a shifted version of a Conjugate Gradient solver. The bare quark masses $`m_0a`$ are chosen to be $`0.0250`$, $`0.0278`$, $`0.0334`$, $`0.0389`$, $`0.0445`$, $`0.0501`$, $`0.0584`$, $`0.0668`$, $`0.0834`$, $`0.1001`$, $`0.1252`$, $`0.1669`$, $`0.2086`$, $`0.2503`$, $`0.2921`$. The $`m_0`$ values in physical units are $`53`$, $`59`$, $`71`$, $`83`$, $`94`$, $`106`$, $`124`$, $`142`$, $`177`$, $`212`$, $`266`$, $`354`$, $`442`$, $`531`$, and $`620`$ MeV respectively.

### 3.1 Axial and vector currents

Let us consider first the vector and axial vector currents. Since each obeys a chiral WI  their renormalisation constants are finite. In Fig. 1, we show $`\mathrm{\Gamma }_A`$ and $`\mathrm{\Gamma }_V`$, calculated in the RI scheme for different bare quark masses as a function of the lattice momentum $`(pa)^2`$. We find they are weakly dependent on the mass, and almost scale independent for $`(pa)^20.5`$.

For the overlap fermion, which satisfies the Ginsparg-Wilson  relation and preserves chiral symmetry on finite lattice. In Ref.  it has been proved that the renormalisation constants for naive vector and axial vector currents satisfy

$$Z_A=Z_V.$$
(18)
In Fig. 2 we show the quantities $`\mathrm{\Gamma }_A\mathrm{\Gamma }_V`$ and and $`\frac{1}{2}(\mathrm{\Gamma }_A+\mathrm{\Gamma }_V)`$, after extrapolating to the chiral limit. We can see from upper part of the figure that there is no signal of effects from chiral symmetry breaking, since $`\mathrm{\Gamma }_A\mathrm{\Gamma }_V`$ is tending to zero at moderate and high momenta. The effects of spontaneous chiral symmetry breaking are visible at low momenta and they are damped at higher momenta. At the momenta of interest, there also seems to be no significant splitting due to non-perturbative effects with the difference between $`\mathrm{\Gamma }_A`$ and $`\mathrm{\Gamma }_V`$ being less than $`1\%`$ at $`(ap)^2=0.7`$ and smaller for momenta above this. In the lower part of Fig. 2, we plot $`\frac{1}{2}(\mathrm{\Gamma }_A+\mathrm{\Gamma }_V)`$ against $`(ap)^2`$, and it can be used for the extraction of both $`Z_A/Z_\psi `$ and $`Z_V/Z_\psi `$ to increase the statistical accuracy.

### 3.2 Pseudoscalar and scalar densities

The pseudoscalar and scalar densities differ from the axial and vector currents in that their renormalisation is not scale independent. For overlap fermions with chiral symmetry the RI scheme preserves the well known $`\overline{MS}`$ relations

$$Z_S=Z_P\mathrm{and}Z_m=\frac{1}{Z_S},$$
(19)
and the quantities $`Z_S/Z_P`$, $`Z_SZ_m`$ and $`Z_mZ_P`$ are expected to be scale independent.

Fig. 3 and Fig. 4 show $`\mathrm{\Gamma }_S`$ and $`\mathrm{\Gamma }_P`$ versus $`(pa)^2`$ with different masses. Unlike in the case of the vector and axial vector currents, they are strongly dependent on the quark mass $`m_0`$ for light quarks. When we want to extract the chiral limit value, we cannot use a simple linear or quadratic fit. In Fig. 3 and Fig. 4, we can clearly see the pole effect in the small $`(pa)^2`$ region.

As discussed in detail in Ref. , due to the contribution from the zero-modes, we may fit the $`\mathrm{\Gamma }_S`$ to the form

$$\mathrm{\Gamma }_S=\frac{c_{1,S}}{(am_0)^2}+c_{2,S}+c_{3,S}(am_0)^2$$
(20)
at a fixed momentum with $`Z_\psi /Z_S`$ being given by $`c_{2,S}`$. The fitted value of $`c_{2,S}`$ is plotted in Fig. 5 as $`Z_\psi /Z_S`$. We find that the coefficient of the pole term, $`c_{1,S}`$ is very small, between $`10^3`$ to $`10^4`$, this means that the pole term only provides a large contribution at very small $`m_0`$, i.e., near the chiral limit.

For $`\mathrm{\Gamma }_P`$, we fit to the form

$$\mathrm{\Gamma }_{P,\mathrm{latt}}=\frac{c_{1,P}}{(am_0)^2}+\frac{c_{2,P}}{(am_0)}+c_{3,P}+c_{4,P}(am_0)^2,$$
(21)
with $`c_{3,P}`$ being equal to $`Z_\psi /Z_P`$. The quadratic mass pole is due to zero-mode effects in $`\overline{q}q`$. The fitted value of $`c_{3,P}`$ is plotted in Fig. 5 as $`Z_\psi /Z_P`$.

As the evidence for the above fitting form, we consider the resulting values for $`Z_\psi /Z_S`$ and $`Z_\psi /Z_P`$. Fig. 5 shows a comparison between the extracted values of these two quantities. As chiral symmetry would predict for $`Z_\psi /Z_S`$ and $`Z_\psi /Z_P`$, the two quantities coincide at moderate and large momenta.This provides an excellent test of both the fitting method to extract the poles and the chiral properties of overlap fermions.

### 3.3 The Tensor current

In Fig. 6 we show $`\mathrm{\Gamma }_T`$ versus $`(pa)^2`$ with different masses. We can see that at moderate and large $`(pa)^2`$, $`\mathrm{\Gamma }_T`$ is not sensitive to the quark masses. The chiral limit value is obtained by a linear fit in quark mass as for the cases of vector and axial vector currents, and the results will be presented in the next subsection.

### 3.4 Running of the renormalisation constants

The renormalised operators are defined as

$$Z_OO_{\mathrm{bare}}=O_{\mathrm{ren}}.$$
(22)
Requiring the bare operator to be independent of the renormalisation scale $`\mu ^2`$ gives the RG equation,

$$\mu ^2\frac{d}{d\mu ^2}O_{\mathrm{ren}}=\frac{1}{Z_O}\mu ^2\frac{dZ_O}{d\mu ^2}O_{\mathrm{ren}}=\frac{\gamma _O}{2}O_{\mathrm{ren}}.$$
(23)
The solution of $`Z_O(\mu ^2)`$ can be written in the following form

$$Z_O\left(\mu ^2\right)=\frac{C_O\left(\mu ^2\right)}{C_O\left(\mu _{}^{}{}_{}{}^{2}\right)}Z_O\left(\mu _{}^{}{}_{}{}^{2}\right).$$
(24)
The coefficient $`C_O\left(\mu ^2\right)`$ can be found in Appendix B of Ref.  and references within.

In this work, the value of $`\alpha _s`$ was calculated at three loops using a lattice value of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ taken from Ref.  as $`\mathrm{\Lambda }_{\mathrm{QCD}}=238\pm 19\mathrm{M}\mathrm{e}\mathrm{V}`$ .

Both $`Z_A`$ and $`Z_V`$ should be scale independent, but this is not the case for $`Z_\psi `$. Fig. 7 shows both $`\mathrm{\Gamma }_A`$ and the scale invariant (SI) quantity (the data after removing the renormalisation group running) calculated as described above:

$$\mathrm{\Gamma }_A^{SI}\left((ap)^2\right)=\mathrm{\Gamma }_A\left((ap)^2\right)/C_\psi \left((ap)^2\right).$$
(25)
The quantity $`C_\psi `$ is normalised so that $`C_\psi (1)=1`$. As can be seen, in this case the renormalisation group running is coming from $`Z_\psi `$ alone. It is small, but it actually improves the scale independence of the data. The remain scale dependence of this data is very small and a plausible explanation for this is a small $`(ap)^2`$ error. Indeed, when a linear fit of the SI data versus $`(ap)^2`$ is performed, for $`0.8<(ap)^2<2.0`$, the gradient is $`0.01`$.

In the case of $`Z_\psi /Z_S`$, the two renormalisation constants are both running with $`(ap)^2`$. Fig. 8 shows both $`\mathrm{\Gamma }_S`$ and the scale invariant (SI) quantity after three loop running. The modification to the raw data, which is taken after the mass-pole has been subtracted gives a satisfactory result. The linear fit of the SI data versus $`(ap)^2`$ in the range of $`0.8<(ap)^2<2.3`$, gives a gradient of about 0.02.

The SI result of $`\mathrm{\Gamma }_T`$ which comes from the two loop running of $`Z_T`$  and three loop running of $`Z_\psi `$ is plotted in Fig. 9, The “bare” data is taken by a simple linear extrapolation to the chiral limit. The linear fit to the SI data gives gradient about 0.01.

With the interpretation that the remaining scale dependence is due to $`𝒪(ap)^2`$ effects, the correct way to extract the renormalisation coefficients is to first construct the SI quantity as described above, and then fit any remaining scale dependence  to the form

$$y=c_1+c_2(ap)^2,$$
(26)
for a range of momenta that is chosen to be “above” the region for which condensate effects are important. We choose the range $`0.8<(ap)^2<2.3`$. We then apply the renormalisation group running formula (an inverse operation of Eq. (25) ) to $`c_1`$ in Eq. (26) to get the ratio of the renormalisation constant which is now free of $`(ap)^2`$ error in RI scheme.

### 3.5  Extracting $`Z_\psi `$ from the propagator

Here we will use Eq. (16) to extract the wave function renormalisation constant $`Z_\psi `$. We will use two different definitions of lattice momenta, i.e. , the discrete lattice momentum defined in Eq. (2) and the kinematic lattice momentum $`q_\mu `$ which was introduced in Eq. (A9) of Ref. .

#### 3.5.1 Discrete lattice momentum $`p_\mu `$

We here use the discrete lattice momentum $`p_\mu `$ in Eq. (16) to calculate $`Z_\psi ^{}`$. By using Eq. (16), we can get the corresponding quantity of $`Z_\psi ^{}`$ (here we use the symbol $`A(p)`$ for $`Z_\psi ^{}`$). The resulting $`A(p)`$ is strongly dependent on the scale. We use three loop perturbative running in RI scheme and finally extract the value of $`Z_\psi `$ from $`A(p)`$ by fitting $`A(p)`$ to Eq. (26) for a range of momenta for which the condensate effects are unimportant. Here we chose the range $`(ap)^2(0.70,2.32)`$, which corresponds to $`p`$ (1.78, 3.23) GeV. $`c_1`$ can be identified as the SI value of $`Z_\psi `$. In Fig. 10 we show a linear fit for $`A(p)`$ versus $`(pa)^2`$ as in Eq. (26).

#### 3.5.2 Kinematic lattice momentum $`q_\mu `$

Another way to calculate $`Z_\psi `$ is to use the kinematic lattice momentum $`q_\mu `$ defined in Eq. (A9) of Ref. . The calculations are similar to those using $`p_\mu `$. The result is plotted in Fig. 11. The fitted value of $`Z_\psi `$ is 1.179(5), with the difference to the value calculated by $`p_\mu `$ being less than 4%. In both calculations, we attribute the remaining $`(pa)^2`$ dependence to the possible $`𝒪((pa)^2)`$ error and use the linear fit Eq.( 26) to remove it. When using the kinematic lattice momentum $`q_\mu `$ to calculate $`Z_\psi `$, the $`𝒪((pa)^2)`$ errors are much smaller than in the case of using the discrete lattice momentum $`p_\mu `$.

### 3.6 Matching to $`\overline{MS}`$ scheme

In order to confront experiment, one frequently prefers to quote the final results in the $`\overline{MS}`$ scheme at a certain scale. For light hadrons, the popular scale is 2 GeV.

The perturbative expansion of the ratio $`Z^{\overline{MS}}/Z^{RI}`$ to two-loop order is given by

$$R=\frac{Z^{\overline{MS}}}{Z^{RI}}=1+\frac{\alpha _s}{4\pi }(Z^{RI})_0^{(1)}+(\frac{\alpha _s}{4\pi })^2(Z^{RI})_0^{(2)}+\mathrm{}.$$
(27)
The numerical values of the matching coefficients, $`Z_0^{(1)}`$ and $`Z_0^{(2)}`$ in Eq. (27) for $`Z_\psi `$ and $`Z_S`$ can be found in Appendix C of Ref. , we have not found the numerical value for $`Z_T`$ in the literature. The final results for the renormalisation constants are listed in Table 2.

## 4 SUMMARY

In this work we performed the non-perturbative renormalisation of composite operators with overlap fermions. The main results are: in $`\overline{MS}`$ scheme at $`\mu `$ = 2 GeV, $`Z_\psi `$ = 1.153(5)(22), $`Z_S`$ = 0.970(4)(20), $`Z_A`$ = 1.070(5)(21), and $`Z_T`$ = 1.167(3)(23), where the first error is statistical and the second error comes from using the two different lattice momenta to get $`Z_\psi `$.

Support for this research from the Australian Research Council is gratefully acknowledged.


text:
# Quantum dots based on spin properties of semiconductor heterostructures

## Abstract

The possibility of a type of semiconductor quantum dots obtained by spatially modulating the spin-orbit coupling intensity in III-V heterostructures is discussed. Using the effective mass model we predict confined one-electron states having peculiar spin properties. Furthermore, from mean field calculations (local-spin-density and Hartree-Fock) we find that even two electrons could form a bound state in these dots.

Low-dimensional semiconductor structures constitute one of the most suitable scenarios to observe quantum phenomena in condensed matter physics. In this sense, the atomic-like properties of quantum dots and the conductance quantization in narrow wires are seminal examples widely studied . At present, the available fabrication techniques of such semiconductor structures rely on charge-based confinement at length scales where quantum effects are unavoidable. More precisely, confinement is obtained through the energy band discontinuities of different materials or by directly creating a potential landscape by means of metallic electrodes on a semiconductor heterostructure.

The aim of this work is to point out, using a theoretical model, the feasibility of a class of mesoscopic semiconductor structures where the spin plays a relevant role in the mechanism that determines the confining properties, without being a pure spin-dependent confinement . The mechanism also differs from conventional electrostatic confinement in that it does not rely on the direct coupling of the electron charge with a horizontal electric field.

In III-V non-magnetic semiconductor heterostructures there are two relevant interactions concerning the spin: the spin-orbit coupling and the hyperfine interaction between carriers and polarized nuclei. The latter mechanism has a characteristic energy scale usually lower than that of the spin-orbit coupling .

In semiconductors, spin-orbit coupling stems from the relativistic effect caused by the electric field due to the lack of inversion symmetry of certain alloys like III-V heterostructures. Depending on the particular origin of the electric field two contributions are distinguished; the electric field created by the bulk inversion asymmetry (BIA) of the material gives rise to the Dresselhaus term , while the asymmetry in the profile of the heterostructure, structural asymmetry (SIA), generates the so-called Bychkov-Rashba term . This latter mechanism constitutes the basis of the proposed electronic nanostructures.

As mentioned above, the intensity of the Bychkov-Rashba spin-orbit coupling depends on the effective electric field perpendicular to the plane of the quantum well. This implies that it can be, in principle, engineered by adjusting the specific fabrication settings of each heterostructure. Another possibility to control the Rashba interaction, that has been demonstrated experimentally , is by means of an external electric field perpendicular to the plane of the quantum well.

Exploiting the tunability of the Bychkov-Rashba mechanism we consider quantum dots whose main ingredient is a space-modulated spin-orbit intensity. This spatial modulation could be obtained with biased small electrodes on top of the semiconductor heterostructure. However, such an arrangement would not only modify the vertical electric field but it would also generate an in-plane gradient of electric potential, leading to hybrid quantum dots where both mechanisms are present: horizontal electric potential gradient and spin-orbit intensity gradient. A more refined method to obtain space modulated spin-orbit intensities requires a deeper insight in the origin of the Rashba intensity. In this respect, de Andrada et al  discussed the relevance of the barrier layers in the determination of the spin-orbit intensity. More recently, Grundler  has shown that for an InAs-based quantum well an applied vertical electric field of the order of $`10^3`$ V/cm can lead to a noticeable control of the spin-orbit intensity, whereas the built-in electric field of the sample is of the order of $`10^5`$ V/cm. This high sensitivity to a relatively weak electric field has been attributed to the induced variations in the wavefunction penetration into the barriers. The strong relevance of the barrier layers in the determination of the Rashba intensity supposes a great advantage that could help to create a quantum well having a space-modulated spin-orbit coupling through space-engineered wavefunction barrier penetration (using inhomegeneous barrier layer composition, thickness, …). This would lead to a fixed (not tunable) substantial step in the spin-orbit intensity capable of inducing the confinement.

We model structures confining electrons from the conduction band of the semiconductor and we use the effective mass approximation. We also suppose that the motion perpendicular to the well’s plane is restricted to the first transversal subband, so that the electrons effectively move in a two-dimensional region. These simplifications might be quite drastic in some cases and, therefore, the scope of this work is limited to a discussion of the feasibility of the effects rather than to a detailed quantitative description. In this model we include the standard Bychkov-Rashba term

$`_R={\displaystyle \frac{\lambda _R(r)}{\mathrm{}}}\left(p_y\sigma _xp_x\sigma _y\right),`$ (1)
where $`p_x`$ and $`p_y`$ represent the components of the in-plane electron’s momentum and the $`\sigma `$’s are the corresponding Pauli matrices.

In Eq. (1) the inhomogeneity of the spin-orbit intensity is included through the spatial dependence of the coupling constant $`\lambda _R(r)`$, which is assumed of radial type ($`r`$ and $`\theta `$ are used to label the standard polar coordinates in the plane). To characterize the quantum dot we consider the following step-like variation of the spin-orbit intensity

$$\lambda _R(r)=\lambda _e+(\lambda _i\lambda _e)\frac{1}{1+e^{\frac{rR_0}{\sigma }}},$$
(2)
where $`\lambda _e`$ and $`\lambda _i`$ represent the external and internal constant values of spin-orbit intensity, respectively. The quantum dot radius is given by $`R_0`$ and $`\sigma `$ is a small diffusivity introduced to avoid discontinuities in the numerical solutions.

The Bychkov-Rashba interaction breaks the rotational symmetry of the Hamiltonian. Nevertheless, a symmetry combining spatial and spin degrees of freedom is maintained when the Rashba intensity is a purely radial function. Namely, the Hamiltonian is invariant under a combined rotation in real and spin spaces by the same angle, i.e., it commutes with the generalized angular momentum $`J_z=L_z+S_z`$. Note that, at variance with conventional circular quantum dots (based on electric-potential confinement ), the Hamiltonian no longer commutes with $`L_z`$ and $`S_z`$ separately. As a consequence of the $`J_z`$-symmetry the single-particle solutions can be expressed, in a general form, as two-component spinors of type

$$𝚿_{nj}(r,\theta )=\left(\begin{array}{c}\varphi _{nj}(r)e^{i(j1/2)\theta }\\ \varphi _{nj}(r)e^{i(j+1/2)\theta }\end{array}\right),$$
(3)
where $`n`$ and $`j`$ are quantum numbers characterizing the states, with values $`n=1,2,\mathrm{}`$, $`j=\pm 1/2,\pm 3/2,\mathrm{}`$. Since the angular parts are analytical, we can reduce the Schrödinger equation to one-dimensional radial equations for the two spinor components,

$`{\displaystyle \frac{\mathrm{}^2}{2m^{}}}\left[{\displaystyle \frac{\varphi _{}}{r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\varphi _{}}{r}}{\displaystyle \frac{(j1/2)^2}{r^2}}\varphi _{}\right]\lambda _R(r)\left[{\displaystyle \frac{\varphi _{}}{r}}+{\displaystyle \frac{j+1/2}{r}}\varphi _{}\right]=\epsilon _{nj}\varphi _{}`$ (4)
$`{\displaystyle \frac{\mathrm{}^2}{2m^{}}}\left[{\displaystyle \frac{\varphi _{}}{r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\varphi _{}}{r}}{\displaystyle \frac{(j+1/2)^2}{r^2}}\varphi _{}\right]+\lambda _R(r)\left[{\displaystyle \frac{\varphi _{}}{r}}{\displaystyle \frac{j1/2}{r}}\varphi _{}\right]=\epsilon _{nj}\varphi _{}`$ (5)
where $`\varphi _{}`$ and $`\varphi _{}`$ denote $`\varphi _{nj}(r)`$ and $`\varphi _{nj}(r)`$, respectively.

The above model has been analyzed using two alternative numerical procedures: a) solving the coupled one-dimensional equations for the radial components of the eigenspinors, Eqs. (5); and b) solving the full two-dimensional problem in Cartesian coordinates, without imposing any symmetry restriction, in a uniform grid. Before discussing the numerical results it has to be noted that systems with the above spin-orbit interaction are still invariant under time reversal. Consequently, a two-fold degeneracy, known as Kramers degeneracy, must hold. From the radial equations it can easily be seen that the substitutions: $`\varphi _{}\varphi _{}`$, $`\varphi _{}\varphi _{}`$ and $`jj`$ lead to a new state (Kramers conjugate) with opposite generalized angular momentum and the same energy eigenvalue ($`\epsilon _{nj}=\epsilon _{nj}`$). The expectation values $`S_z`$ and $`L_z`$ also change sign with the Kramers transformation. It is worth stressing that in the usual quantum dots (with an electric potential confinement of circular symmetry) the degeneracies are 4 and 2; nonzero (zero) angular momentum states being four-fold (two-fold) degenerate. Therefore, the reduction in degeneracy from 4 to at most 2 is a genuine effect of the spin-orbit interaction.

Since the spin orientation in Eq. (3) depends on the position $`𝐫`$, the eigenstates will show characteristic spin textures. The distribution of vertical and parallel spin read

$`\sigma _z(r)`$ $`=`$ $`|\varphi _{nj}(r)|^2|\varphi _{nj}(r)|^2`$ (6)
$`\stackrel{}{\sigma }_{}(r,\theta )`$ $`=`$ $`2\varphi _{nj}(r)\varphi _{nj}(r)\left[\mathrm{cos}(\theta )\widehat{x}+\mathrm{sin}(\theta )\widehat{y}\right].`$ (7)
Note that while the vertical component is a pure radial function, the in-plane (horizontal) spin points in the radial direction. The fact that the horizontal spin follows the position vector can be easily explained from the combined rotational invariance in spin and position space: consider two points in the dot having the same radius and different azimuthal angle $`\theta `$; due to the symmetry of the Hamiltonian with $`J_z`$, the local in-plane spin orientation in the second point must be that of the first point rotated the angle needed to go from the first point to the second in the real space. This property of the eigenstates is of particular interest for determining the role of the geometrical phase (Berry’s phase) in mesoscopic structures .

From the numerical simulations we observe that a necessary requirement for confined states is that the spin-orbit intensity within the quantum dot must be larger than that of the outer region, i.e., $`\lambda _i>\lambda _e`$. The reason for this behaviour can be understood from bulk considerations: the bulk energy bands

$`\epsilon _{k\pm }^{\mathrm{bulk}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m^{}}}\left(k\pm {\displaystyle \frac{m^{}}{\mathrm{}^2}}\lambda _e\right)^2{\displaystyle \frac{m^{}}{2\mathrm{}^2}}\lambda _e^2`$ (8)
display a constant offset $`m^{}\lambda _e^2/2\mathrm{}^2`$ that depends on the spin-orbit intensity (for a schematic representation of the bulk bands see Fig. 1). Therefore a bulk uniform sample with a big value of spin-orbit intensity has an energy origin for the eigenstates lower than that of a sample having a smaller intensity. In a finite system this energy mismatch can lead to bound states in the region of higher intensity. Nevertheless, it has to be noted that for each spin orientation the energy mismatch is only confining for a certain range of momenta. This means that the confining mechanism is spin-selective through the motion of the particles. In the case of circular symmetry, it implies that electrons rotating in one sense are confined if they have a given spin orientation, while electrons rotating in the opposite sense should have the opposite spin orientation to be confined.

We stress that the SO-dependent offset in Eq. (8) is the real cause of confinement since this term acts as an electrostatic potential when the SO intensity is varying in space. Actually, a different SO Hamiltonian in which this offset were not included would not lead to confined states.

Another question of relevance refers to the magnitude of the step in the spin-orbit intensity needed to obtain confined states. For a typical size of dot radius of the order of 10 atomic effective units simulation shows that an intensity in the range $`\lambda _i1.52.0\times 10^9`$ eVcm is required to obtain a single-particle bound state when the external intensity is zero ($`\lambda _e=0`$). It is also important to specify the external value $`\lambda _e`$, since the minimum step required to achieve confined states strongly depends on this quantity. For instance, when the external intensity takes a value of the order of $`\lambda _e2\times 10^9`$ eVcm the step needed is reduced and the required internal value is in the range $`\lambda _i\lambda _e0.40.7\times 10^9`$ eVcm.

Quantum dots satisfying the required conditions stated above show a level structure analogous to that of conventional quantum dots. For instance, a dot characterized by spin-orbit parameters in the range of the biggest Rashba intensities experimentally measured , $`\lambda _i4\times 10^9`$ eVcm and $`\lambda _e2\times 10^9`$ eVcm, shows a total of 14 single-particle bound states distributed in two subbands (see Fig. 1). The first one containing 10 states characterized by the quantum numbers $`n=1`$, $`j=9/2,7/2,\mathrm{},7/2,9/2`$ and the second one having the remaining 4 states $`n=2`$, $`j=3/2,1/2,1/2,3/2`$. Figure 2 shows the spatial representation of the lowest energy state of this dot corresponding to $`n=1,j=1/2`$ and its conjugated state $`n=1,j=1/2`$. The upper panels show the radial profile of the spinorial components for both states, the only differences between the two cases are the exchange of the spin index and the sign change in one of the components, in agreement with symmetry properties derived above. Middle panel shows the spatial distribution of the density which has radial symmetry and is common for both conjugated states. Finally, the lower panel displays the radial distribution of the $`S_z`$ spin density for the two states, showing that conjugated states posses opposite spin characters.

Up to now, we have only presented results concerning the single-particle properties of the dots. There is also the possibility of the spin-orbit-based confinement to be robust against electron-electron interaction and, thus, able to support multielectron bound states. To treat the interacting problem we impose no symmetry restriction in space and we solve numerically the full two-dimensional problem including the interaction within density functional theory, using a generalization of the local-spin-density approximation (LSDA) for non-collinear spin densities (see Ref.  for a previous application of this formalism). Within this scheme, we find that even two-electron states can be confined to the dot. Figure 3 shows the ground state density for an interacting two-electron dot characterized by a radius of 15 units and the same spin-orbit parameters of Fig. 2. The density displays a circular ring-like shape although the spin-orbit term breaks the rotational invariance. As an additional check of this bound two-electron state, we have also solved the Hartree-Fock equations, finding a similar total energy and density distributions also confined to the dot region. Extending the calculations to treat three or more interacting electrons with the same spin-orbit term does not lead to confined states.

In summary, we have proposed a mechanism to create electronic spatial confinement based on the spin properties of III-V semiconductor structures and, more specifically, on the Bychkov-Rashba spin-orbit interaction and its strong dependence on the barrier layer composition. The eigenstates confined to these dots show common characteristics with conventional semiconductor quantum dots such as a discrete level structure of the single-particle states and the possibility of confined two-electron states. The peculiar origin of the confinement reflects in the spin properties of the single-particle eigenstates. Not any material is suitable to create this type of structures since, as stated above, Rashba spin-orbit intensities of the order $`10^9`$ eVcm and a range of tunability of the same magnitude are required. In this sense, InGaAs-based wells  could be enough to produce the confining effect at mesoscopic scale. However, since the bulk-intensity value $`\lambda _e`$ has a strong relevance in the determination of the minimum conditions for confinement, InAs-based wells reveal as the best candidates to fabricate these dots. This is because of their large built-in Rashba intensities, as large as $`2\times 10^9`$ eVcm, and their tunability range of the same order.

This work was supported by Grant No. BFM2002-03241 from DGI (Spain).