text:
# The Influence of Triaxial Halos on Collisionless Galactic Disks

## 1 Introduction

Galaxies are not spherical objects. Although it is difficult to measure the intrinsic shapes of astronomical bodies from their projected light distributions, inferences can be make from their observed ellipticities. Many studies of the distribution of observed axis ratios of galaxies and galaxy clusters (e.g. Ryden, 1992, 1996; Tremblay & Merritt 1996) have shown that they are well represented by a population of intrinsically triaxial objects. Similar statistical analysis performed by Alam & Ryden (2002) with a much larger sample of galaxies from the Sloan Digital Sky Survey Early Data Release has strengthened the claim that the observed distributions of ellipticities cannot be fit by projections of oblate or prolate spheroids.

This type of analysis is only available in the cases where the mass distribution creates light. The shapes of dark matter halos, therefore, must be inferred from even more indirect means. Proxy measures of the shape of our own halo offer widely varying results. The tidal stream of the Sagittarius dwarf galaxy implies a nearly spherical halo (Ibata et al. 2001), in contrast to the results of star counts, which give a minor to major axis ratio of $``$ 0.6 (Siegal et al. 2002). The dark matter halos created in numerical studies of structure formation are non-axisymmetric systems. $`N`$-body simulations of dissipationless collapse in the cold dark matter (CDM) scenario by Dubinski & Carlberg (1991) showed that the resulting virialized halos deviated substantially from spherical symmetry and from axisymmetry, with average axis ratios $`c/a=0.5`$ and $`b/a=0.7`$, where $`a`$, $`b`$ and $`c`$ are the major, intermediate, and minor axes respectively. Simulations of larger cosmological volumes, which more adequately reproduce the hierarchical nature of structure formation in CDM, have also shown that dark matter halos are triaxial (Barnes & Efstathiou 1987; Warren et al. 1992). Analysis of the currently favored $`\mathrm{\Lambda }`$CDM cosmology has shown similar results (Jing & Suto 2002).

The simulations cited above consist only of a collisionless dark matter component and include no gas dynamics. The introduction of a gas and stellar disk in the equatorial plane of the halo would certainly have effects on the internal structure of the halo. Fully self-consistent cosmological simulations of dissipational structure formation have problems reproducing more than general galactic attributes (e.g. Murali et al. 2002), and the resulting halo shape has not been the focus of these efforts. However, Dubinski (1994) showed that the adiabatic growth of a disk-like potential in a triaxial halo did not influence $`c/a`$, while the intermediate to major axis ratio was only slightly increased to $`b/a0.7`$.

Models and simulations of disk galaxies, however, usually employ spherical halos, which do not reflect these results from CDM simulations. The purpose of most current investigations of galaxy-scale simulations has been modeling transformations of galaxies through major mergers (e.g. Barnes 1992; Dubinski, Mihos, & Hernquist 1999), minor encounters (e.g. Struck 1997; Thakar & Ryden 1998) and gas-dynamical effects (e.g. Mihos & Hernquist 1996; Springel 2000). The triaxiality of the halo does not immediately effect the results of these simulations and the implications of using triaxial halos have not been explored.

Early studies of the stability of disk galaxies (e.g. Ostriker & Peebles 1973; Hohl 1976; Efstathiou et al. 1982) demonstrated that a cold, rotationally supported disk is dramatically unstable to bar formation. This fact was used as justification for the existence of heavy, spherical, dark matter halos enveloping disk galaxies. This dynamically hot component of the galaxy inhibits bar formation. It has also been shown that a dense core gives a galaxy an inner Lindblad resonance, which does not allow swing amplification of waves through the center of the galaxy and thus acts as a bar suppressant (Sellwood 1989)111Efstathiou et al. (1982) did perform two simulations with dense centers but found that it did not make their disks stable. Sellwood proposed that the low resolution of their simulations was responsible for the results of their “bulge” simulations.. In all cases the stabilizing mass component has been assumed to have spherical symmetry. If these galaxies were enveloped by triaxial halos, they would be more susceptible to bar instabilities. A systematic exploration of the stability of disks with non-axisymmetric halos has not been conducted. Curir & Mazzei (1999) presented fully self-consistent gas-dynamical simulations of disk galaxies with triaxial halos, but only used two halo models and did not employ bulges. The resolution in their simulations was severely limited ($`N_{disk}=3,000`$), far less than the number of particles used by Efstathiou et al. (1982) which Sellwood (1989) found to be insufficient.

The notion that the potential field of a galaxy may contain non-axisymmetric features is not a new one; it is, in fact, readily apparent from the bars seen in many galactic disks. The effect of an oval distortion, like a bar, on a gas disk has been proposed as a driving mechanism for spiral density waves (Lin 1970). This idea has been explored numerically as well (Sanders & Huntley 1976; Sanders 1977; Huntley et al. 1978). These early simulations generally were limited to two-dimensional, massless disks driven by an analytic distortion to the symmetric potential, but were successful in creating significant spiral response in the disk. However, bar forcing does not reproduce tightly bound spiral patterns. And in comparison to Sanders & Huntley (1976), the response of the outer disk is weak when a more realistic bar potential, one that drops off faster with radius, is used in the simulation (Sanders & Tubbs 1980). Observationally spiral patterns are as common in barred galaxies as galaxies without bars. Also, there is evidence that bars and spiral arms have different pattern speeds, showing that one feature may not drive the other (Sellwood & Sparke 1988). A triaxial halo would not be subject to these concerns since the potential of the dark matter would dominate the outer regions of the disk, especially low mass disks which are stable to bar formation.

The bulge-to-disk ratio of a spiral galaxy, as well as the tightness of the spiral pattern, are two integral criteria of the Hubble classification scheme. On the average, there exists a smooth increase in pitch angle with later Hubble types (Kennicutt 1981). Such a correlation, however, is little more than a consistency check between classification parameters, and there exists significant scatter in the measured pitch angle for a given Hubble type.

Even with this scatter, the correlation between arm pattern and galaxy type has prompted several theoretical studies on the origin of the this relation (see Kennicutt & Hodge 1982, and references therein). Kennicutt (1981) found that the correlation between pitch angle and maximum rotational velocity of the disk was as good as the theoretical model predictions, suggesting that the mass distribution and its resulting rotation curve is an important determinant in producing the shape of spiral arms. A more massive or more concentrated bulge will induce higher rotation velocities and loosely follow the correlation with Hubble type as well.

Our simulations include a collisionless stellar disk and a stiff, rotating halo. Most simulations include a bulge, which is also stiff. Our use of dissipationless simulations with halos which are not fully self-consistent allows us to cover a significant range of parameter space for all three galaxy components; the disk, bulge, and halo.

The structure of this paper is as follows: §2 presents our initial conditions as well as the $`N`$-body techniques used in this study. §3 shows the analysis and results of our fiducial simulation. §4 describes our suite of different galaxy models. In §5 we investigate the stability of our models with different halo and bulge masses, as well as different halo rotation rates. §6 presents our results for the spiral morphology of our simulations. We present results for models which follow the bulge-to-disk ratios of the Hubble sequence and models which vary widely in rotation velocity. We compare our results to observations of spiral galaxies.

## 2 Numerical Methods

### 2.1 Initial Conditions

The positions and velocities of the disk particles are initialized in the method outlined by Hernquist (1993). The density distribution of the disk is described by

$$\rho (R,z)=\frac{M_d}{4\pi h^2z_0}\mathrm{exp}\left(\frac{R}{h}\right)\text{sech}^2\left(\frac{z}{z_0}\right)$$
(1)
where $`R`$ and $`z`$ are cylindrical coordinates, $`h`$ and $`z_0`$ are the scale length and scale height of the disk, and $`M_d`$ is the total mass of the disk. The circular velocity and velocity dispersions are obtained from

$$v_c=\sqrt{|a_R|R}$$
(2)
$$\sigma _R^2=\sigma _{R,0}^2\mathrm{exp}(R/h)$$
(3)
$$\sigma _z^2=\pi G\mathrm{\Sigma }(R)z_0$$
(4)
$$\sigma _\varphi ^2=\sigma _R^2\frac{\kappa ^2}{4\mathrm{\Omega }^2}$$
(5)
where $`a_R`$ is radial acceleration, $`\mathrm{\Sigma }`$ is the surface density of the disk, $`\kappa `$ is the epicyclic frequency and $`\mathrm{\Omega }`$ is the angular rotation rate. The radial acceleration in equation (2) is the total acceleration contributed by the self-gravity of the disk, the dark matter halo, and a spherical bulge. The normalization of the radial velocity dispersion, $`\sigma _{R,0}^2`$, is obtained by requiring $`\sigma _R`$ at a specified radius $`R_0`$ to be a multiple of the critical dispersion at that radius:

$$\sigma _{R=R_0}=Q\frac{3.36G\mathrm{\Sigma }(R_0)}{\kappa (R_0)}$$
(6)
where $`Q`$ is Toomre’s parameter for gravitational stability in a rotationally supported disk. The disk models in this paper are defined to have $`Q=1.5`$ at $`R=2.43h`$, which in the Milky Way galaxy would be approximately the solar radius. The rotational velocity of each particle is then given by

$$\overline{v_\varphi }^2v_c^2=\sigma _R^2\left(1\frac{\kappa ^2}{4\mathrm{\Omega }^2}2\frac{R}{h}\right).$$
(7)
The different halo models used in this paper can be represented by Dehnen’s (1993) $`\gamma `$-models;

$$\rho (r)=M_h\frac{3\gamma }{4\pi }\frac{r_0}{r^\gamma (r+r_0)^{4\gamma }}$$
(8)
where $`\gamma `$ is the interior logarithmic slope of the density profile and $`r_0`$ is the scale length of the halo. A model with $`\gamma =1`$ corresponds to a model proposed by Hernquist (1990) for bulges and elliptical galaxies. A spherical Hernquist model is used here to represent the bulges of our galaxies. Dubinski & Carlberg (1991) showed that the Hernquist model is a good fit for CDM halos, and is close to the general density profile for dark halos found by Navarro, Frenk, & White (1997), which also asymptotes to $`r^1`$ at small radii. The $`\gamma `$-models can be generalized to triaxial systems by replacing the radial coordinate $`r`$ with $`m`$, where

$$m^2=\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}$$
(9)
where $`a`$, $`b`$, and $`c`$ are once again the principal axes of the system.

In practice, the disk-particle positions are sampled from the density function in equation (1), the accelerations on each particle from the different galaxy components are computed, and then the velocity dispersions and the rotational velocity of each particle are calculated. The peculiar velocities of each particle are sampled from Gaussian probability distributions with dispersions given by equations (3)-(5). Because $`\kappa `$ and $`\mathrm{\Omega }`$ in equations (5)-(7) are assumed to be azimuthally symmetric, the velocities for the disk particles are computed assuming that the halo is spherical; therefore when the disk models are placed inside a triaxial halo they are initially out of equilibrium.

### 2.2 N-body Techniques

To compute the self-gravity of the disk particles, we have employed GADGET (Springel, Yoshida, & White 2001), a tree code based on the hierarchical algorithm of Barnes & Hut (1986). GADGET employs a spline-softened force calculation that is expanded to quadrupole order for particle-cell interactions. The timesteps are continuously varying and individual for each particle. The spatial and time resolutions are set by the smoothing length, $`ϵ`$, and the timestep accuracy parameter, $`\alpha `$. For simulations with $`N=10^5`$, which are the majority of the simulations in this paper, these parameters were set to $`ϵ=0.057`$ and $`\alpha =0.04`$.

The gravitational contribution of the halo is computed by the self-consistent field (SCF) method of Hernquist & Ostriker (1992). This method is based on the expansion of the galaxy potential in an orthonormal set of basis functions. The accuracy of the potential calculated is driven by the number of terms in the expansion. The numbers $`n`$ and $`l`$ determine the number of radial and angular basis functions used. The integration time, however, is a linear function of the number of terms as well. We found that the best balance of accuracy and speed was achieved with $`n`$=16 and $`l`$=10, which were the values used in all the simulations. The zeroth-order basis function was chosen to be the Hernquist profile, and the coefficients of the expansion were computed by Monte Carlo sampling of the halo density distribution, $`\rho (m)`$. Since the coefficients only have to be calculated once for each halo model, accurate calculations can be made by initially sampling a large number of particles from each halo profile, $`10^7`$ particles. The SCF method has been used frequently to model triaxial systems (e.g. Merritt & Quinlan 1998; Holley-Bockelmann et al. 2001). When comparing the force on a set of test particles calculated from the SCF coefficients with the force calculated from direct summation over the particles used for the coefficients, the rms error was $`6\times 10^3`$. This error level was maintained in models where $`\gamma `$ was not equal to unity.

The contribution to the potential of the galaxy from the bulge can be computed analytically. Therefore, the total acceleration of each particle is given by

$$𝐚_i=𝐚_{tree}\frac{M_b}{(r+r_b)^2}\widehat{𝐫}\underset{n,l,m}{}A_{nlm}\mathrm{\Phi }_{nlm}.$$
(10)
The first term is the acceleration returned by GADGET from walking the tree of disk particles. The middle term is the acceleration from the bulge, with the subscript $`b`$ denoting bulge properties. The third term in equation (10) is the acceleration computed from the SCF expansion of the halo potential. The particle timesteps are inversely proportional to the particle acceleration, $`\mathrm{\Delta }t_i=\alpha /|𝐚_i|`$.

## 3 Simulating a Milky Way Galaxy

### 3.1 Parameters of the Fiducial Simulation

Table 1 lists the physical parameters of the fiducial simulation. Scaling our $`N`$-body units to physical values reasonable for the Milky Way galaxy as in Hernquist (1993), we get $`h=3.5`$ kpc and $`M_d=5.6\times 10^{10}M_{}`$, making the unit of time 13.1 Myr. With these scalings, the value of the rotational speed at the solar radius, $`R_{}=2.43h`$, is nearly 220 km s-1. The fiducial halo rotation rate, $`\mathrm{\Omega }_h`$ =0.06, is equivalent to a pattern speed of 4.5 km s-1 kpc-1 in these unit. (For comparison, the simulations of Bekki & Freeman (2002), in which a massless gas disk is embedded in a rotating triaxial halo, use a pattern speed of 3.8 km s-1 kpc-1.) Figure 1 shows the contributions of each galaxy component to the total circular velocity. Figure 1 also shows the values of the $`\mathrm{\Omega }`$, $`\mathrm{\Omega }+\kappa /2`$, and $`\mathrm{\Omega }\kappa /2`$ resonances as a function of radius. Although the simulations of Dubinski (1994) suggest that the equatorial plane of the halo can be as elliptical as $`b/a=q_b=0.7`$, we chose a more moderate value of $`q_b=0.85`$ for our fiducial simulation. The minor axis ratio, $`c/a=q_c`$, was set to 0.7. Figure 2 plots the fractional difference in the potential along the $`x`$ and $`y`$ axes of halos with different values of $`q_b`$. The dotted line represents the asymmetry of the potential from the entire fiducial model. At five scale lengths, the difference in the halo potential is $`3\%`$, while the difference in the total system potential is $`2\%`$.

To test our numerical method, we ran a simulation with a spherical halo but still computed the halo’s contribution to the potential through the SCF method. Outside of the shape of the halo, which was set to $`q_b=q_c=1`$, all the parameters listed in Table 1 were the same. The simulation was run to a time of 500, which corresponds to a physical time of 6.55 Gyr. Figure 3 shows the rotation curve, the surface density, and the vertical particle density for the initial and final states of the simulation. These curves show remarkably little deviation from their initial values, proving that the integration of the SCF into GADGET works well for this type of simulation.

### 3.2 Analyzing the Simulations

The fiducial model was integrated to a time of 6.55 Gyr as well. The influence of the triaxial halo of the density structure and rotation curve can be seen in Figure 4. As in the simulation with the spherically symmetric halo, there is little change in the azimuthally averaged structure of the disk. The maximum of the rotation curve is slightly lower than the initial velocity, caused by a slight heating of the disk as rotational energy is transformed into thermal. It is interesting to note that the vertical structure of the disk does not change even after this lengthy integration time. The disk is not vertically heated by the constant perturbation of the halo potential.

Structure in the disk is generated by the influence of the halo. A time-varying, non-axisymmetric potential creates a torque on the disk, and it is expected that symmetric, $`m`$=2 structures should form. Figure 5 shows an image of the fiducial model at a time of 655 Myr. The bisymmetric structure is apparent. To quantify the response of the disk to the halo potential, we performed a Fourier decomposition of the surface density of the disk into logarithmic spirals (Sellwood & Athanassoula 1986; Binney & Tremaine 1987) of the form $`\alpha u+m\varphi =ct`$, where $`u=\mathrm{ln}R`$ and $`\varphi `$ is the azimuthal angle. The amplitude $`A`$ of each Fourier component is the Fourier transform of the surface density $`S(u,\varphi )`$;

$$A_m(\alpha )=\frac{1}{(2\pi )^2}_{\mathrm{}}^{\mathrm{}}𝑑u_0^{2\pi }S(u,\varphi )\text{e}^{i(\alpha u+m\varphi )}𝑑\varphi .$$
(11)
In practice, since an $`N`$-body simulation is resolved by discrete particles, the integral signs can be replaced by a summation over all particles and equation (11) then becomes

$$A_m(\alpha )=\frac{1}{N}\underset{j=1}{\overset{N}{}}e^{i(\alpha u_j+m\varphi _j)}.$$
(12)
The parameter $`m`$ was set to 2 to look for bisymmetric structure, and coefficients were calculated for values of $`25\alpha 25`$. The decomposition of the fiducial model at 655 Myr is shown in Figure 6. The decomposition of the model with the spherical halo at the same time interval is shown for comparison. Power at negative values of $`\alpha `$ corresponds to leading spiral waves. Power at positive values of $`\alpha `$ corresponds to trailing spiral waves, which in the fiducial model have a distinct peak at $`\alpha =12`$. The pitch angle $`p`$ for this logarithmic spiral can be calculated from $`p=\mathrm{tan}^1(m/\alpha )=9.5^{}`$. Smaller values of $`\alpha `$ correspond to more loosely-bound spiral patterns, or arms with a larger pitch angle. The peak seen in Figure 6 centered near $`\alpha =0`$ is a linear bar structure.

The presence of a bulge in the fiducial model prevents a significant bar from forming. The peak of the power at $`\alpha _{max}=12`$ is significantly higher than the power at $`\alpha =0`$. An identical simulation to the fiducial model was performed, with the only change being the removal of the bulge component of the galaxy. The Fourier decomposition for this bulgeless model at $`t=655`$ Myr is also plotted on Figure 6. The amplitude of the trailing arms is slightly reduced and $`\alpha _{max}`$ is at smaller $`\alpha `$, but the bar component of the disk has been greatly enhanced. The Fourier decomposition for the simulation with a spherical halo is also plotted for comparison. The spectrum of this model shows no significant peaks and is at low amplitude.

### 3.3 Growth of the Spiral Modes

Because the disk consists solely of collisionless particles, it is expected that the spiral modes excited by the halo should damp away simply due to relaxation effects. A galactic disk made entirely of stars could not sustain spiral patterns indefinitely, although since the potential of such a system is much smoother than our simulations the relaxation time would be much longer. We quantify “arm strength” as the sum of the amplitudes from $`0\alpha 25`$. A strong bar would systematically bias this definition of arm strength, since $`\alpha =0`$ represents linear structures. However, the simulations for which this calculation was done did not contain strong bars. The amplitude of the Fourier spectrum at $`\alpha =0`$ in Figure 6 is not the dominant peak. Above $`\alpha =25`$, which corresponds to a pitch angle of $`4.6^{}`$, the Fourier spectrum is usually noise. Figure 7 plots arm strength against time for the fiducial simulation. It shows a rapid growth of the spiral modes, peaking around 0.5 Gyr, then exponentially decaying to an asymptotic value of around 0.5. The rms value of the arm strength for the model with the spherical halo is 0.2.

It is important to note that these simulations all begin with the disk being out of equilibrium with the triaxial halo. Several simulations were performed to test the response of the disk to a adiabatic conversion of the halo potential from spherical symmetry to the full triaxial form. Initially, the model is begun with the halo potential being spherical. The triaxial potential is phased in by the function

$$a_i=(1\xi )\frac{GM_h}{(r_i+r_0)^2}\xi \underset{n,l,m}{}A_{nlm}\mathrm{\Phi }_{nlm}$$
(13)
$$\xi =\left(\frac{t}{t_{ph}}\right)^2\left(32\frac{t}{t_{ph}}\right)$$
(14)
while $`t<t_{ph}`$. This function has zero derivatives at $`t=0`$ and $`t=t_{ph}`$, so that the change in potential is smooth at the limits. Simulations were run with $`t_{ph}=`$ 20 and 40 in $`N`$-body units, corresponding to physical times of 0.26 and 0.52 Gyr. These simulations, plotted in Figure 7, reach their peak arm strength near 2$`t_{ph}`$ and then proceed to decay along the same curve as the fiducial model. The peak power is reduced monotonically with $`t_{ph}`$, from 2.5 to 1.7 and 1.2. All the curves in Figure 7 begin to overlap at a time of 0.8 Gyr and decay to the same arm strength thereafter. The structure is very similar in all three simulations at $`t0.8`$ Gyr. Even though the peak arm strength is reduced as the triaxial potential is phased in, the final state of the disk is unchanged.

### 3.4 Increasing the Simulation Resolution

To explore the $`N`$-dependence of our simulations, we increased the resolution of our fiducial model by an order of magnitude, realizing the disk with $`10^6`$ particles and lowering the smoothing length to $`ϵ=0.027`$. In comparison to the $`N=10^5`$ simulation, the structures in the disk, shown in Figure 8, are extremely apparent and well-resolved. Inside three scale lengths, strong $`m=2`$ spiral arms are seen. These arms are encircled by two different structures at different radii. These are caustic-like structures, where the particle orbits from different radii have crossed each other, leading to a local density enhancement. The outer caustic is highly elliptical and is not continuous with azimuth. It has an axis ratio of 0.7 and a semimajor axis of 6.8$`h`$. The inner caustic, which nearly connects with the spiral pattern inside it, has an axis ratio of 0.9 and a semimajor axis of 3.9$`h`$. The inner Lindblad resonance for this model is shown on top of the image in Figure 8. Located just outside the inner caustic, at 5$`h`$, it connects the inner and outer caustics at the minor axis of the outer caustic.

Figure 9 shows the Fourier decomposition of the $`N=10^6`$ simulation in comparison to the fiducial model. The peaks in the power of both simulations are at the same locations, with a slightly higher peak at $`\alpha =12`$ and lower peak at $`\alpha =0`$ for the higher resolution simulation. Figure 9 also plots the arm strength as a function of time. The peak arm strength of the high resolution simulations is slightly higher than that of the fiducial simulation, and the two curves decay at nearly the same rate. These results demonstrate that the fiducial resolution is adequate for the analysis performed.

### 3.5 Altering Halo Parameters

Cold dark matter simulations do not predict a single shape to dark matter halos; rather a distribution of axis ratios and triaxialities. For the halo figure rotation, the fiducial value of $`\mathrm{\Omega }_h=0.06`$ is somewhat arbitrarily chosen and limits on this parameter are not defined. We ran simulations to test the effects of varying these parameters, first varying $`q_b`$ with $`\mathrm{\Omega }_h`$ held constant, then varying $`\mathrm{\Omega }_h`$ with $`q_b`$ held at 0.85.

Figure 10 shows the Fourier decompositions of five models, including the fiducial model, with different levels of flattening. The most extreme model has $`q_c=0.5`$ and $`q_b=0.75`$. The most moderate model has axis ratios of 0.9 and 0.95. (The value of $`q_c`$ for our models is always $`2q_b1`$.) All the models show a peak near $`\alpha =12`$, showing that the pitch angle of the spiral arms created remains nearly the same in these different models. As $`q_b`$ gets closer to unity, the response of the disk falls off; at $`q_b=0.95`$, the amplitude of the Fourier peak is a factor of 2.3 lower than the fiducial model, and the disk itself shows little visible evidence for spiral structure. As $`q_b`$ gets smaller, however, response of the disk around $`\alpha =0`$ increases dramatically. In the most flattened model, the amplitude at $`\alpha 3`$ is greater than that at $`\alpha =12`$.

Figure 11 shows how the integrated arm strength and the value of $`\alpha _{max}`$ each change with the rotation rate of the halo, $`\mathrm{\Omega }_h`$. As the halo rotation increases, the response of the disk is monotonically increasing when $`\mathrm{\Omega }_h`$ $`0.045`$. As $`\mathrm{\Omega }_h`$ increases to $`\mathrm{\Omega }_h`$ = 0.12, the value of $`\alpha _{max}`$ decreases from 12 to 7. This can be seen in Figure 12, which shows images of four of these models.

## 4 Exploring Parameter Space

The mass distribution of the Milky Way, although studied for several decades, is not known with certainty. The mass distribution of other galaxies are not known with much confidence at all, although it is known that there exists a range of halo and disk masses. Low surface brightness galaxies are thought to be dark-matter dominated, while many researchers espouse the maximal disk hypothesis presented by van Albada & Sancisi (1986) for galaxies with high surface brightness

To explore these parameters, we constructed several sets of models in which all three components of a disk galaxy–disk, bulge, and halo–were varied in systematic ways. First we explored models where the disk mass was reduced and the dynamics became more heavily influenced by the potential field of the halo. Second, we ran models where the parameters of the halo, namely the interior slope $`\gamma `$ and the scale radius $`r_0`$, were altered such that the mass distribution of the galaxy became disk-dominated at small radii. The total mass of the halo is constant, or nearly so, if one considers the truncation radius of the halo to be at a value of 100$`h`$, the value chosen to calculate the moment of inertia of the halo. We characterize each of these models by its mass ratio, $`(M_d/M_h)_R_{}`$, the ratio of the disk to halo mass within the solar radius, 8.5 kpc. Table 2 lists the general properties of the disk-dominated models, denoted D1–D5, where D5 is the most nearly maximal disk with a mass ratio of 4.53, and the properties of the halo-dominated models, denoted H1–H5, where H1 has the least massive disk with a mass ratio of 0.22. Figures 13 and 14 show the circular velocity profiles for the halo- and disk-dominated models respectively. In all models the bulge mass is set to 0.3$`M_d`$.

Third, we created a set of models in which the disk and halo mass distributions were held constant and the bulge mass and scale length were varied. The disk mass was $`M_d`$ = 0.8 and the halo parameters were the same as the fiducial model. The bulge masses were $`M_b`$ = 0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.3, and 2.0. The value of each $`r_b`$ was calculated such that the average density inside one bulge scale radius was the same as that of the fiducial bulge, namely

$$r_b=r_{b,0}\left(\frac{M_b}{M_{b,0}}\right)^{1/3}$$
(15)
where $`r_{b,0}`$=0.2 and $`M_{b,0}=0.3`$.

Our last set of models centered on models that had lower rotation velocities, values closer to 100 km s-1 than to the fiducial value of 220 km s-1. The motivation behind these models will be discussed in §6. The specific parameters used are as follows: The halo parameters are identical to the halo in model D5, which is the model with the lowest central concentration. The disk mass was 0.5. The bulge mass was varied, having values of 0, 0.075 and 0.15.

The results of these simulations will be discussed in the following two sections, focusing on the formation of bars and morphology of the spiral arms induced in the disk by the halo.

## 5 Bar Formation

Bars are simple to detect in our simulations, since they represent the $`\alpha =0`$ modes of the Fourier decompositions. Models with strong spiral arms, however, can also return high amplitudes at $`\alpha =0`$ compared to models with no structure at all. To prevent false detection of bars, we followed the procedure described in Efstathiou et al. (1982) to quantify the strength of the bar. The particles are binned in radial rings of linear separation and the amplitude and phase of the $`\alpha =0`$, $`m=2`$ component is calculated for each ring. If a bar is present, it is defined as existing within the rings that have phases which are coherent to within 10%. The bar strength is the maximum Fourier component within this coherence and is denoted by the symbol $`\delta _2^{max}`$. The rings are 0.1$`h`$ in width and contain several thousand particles each. To be considered a valid detection, the coherence must be at least 5 rings in length but does not have to begin at the center-most ring– some simulations with obvious bars do not begin having coherent phases until $`R0.30.5h`$. From this method we can determine the amplitude, length, and pattern speed of the bar. The binning of particles in cylindrical rings eliminates false detection of bars in simulations with strong spiral patterns because these models will have high values of $`\delta _2^{max}`$ but will not have coherent phases.

Figure 15 shows simulations at four different levels of bar strength: $`\delta _2^{max}`$ = 0.1, 0.2, 0.4, and 0.6. Values of $`\delta _2^{max}`$ less than or equal to 0.2 have weak bars that are generally $`1`$ disk scale length in size. At $`\delta _2^{max}`$ = 0.4, the bar is very distinct and longer. At $`\delta _2^{max}`$ = 0.6, the highest value seen in our simulations, the bar dominates the mass of the disk and can be nearly four scale lengths in size. When integrated with spherical halos, all of the models presented in this section developed no detectable bars.

Figure 16 shows the time evolution of bar strength in models with different values of halo rotation rate. For the fiducial value of $`\mathrm{\Omega }_h`$$`=0.06`$, the disk-dominated models are generally stable. Only model D3 develops a detectable bar, one which asymptotes to $`\delta _2^{max}`$$`0.15`$. When $`\mathrm{\Omega }_h`$ is increased, however, bar instabilities become much larger. At $`\mathrm{\Omega }_h`$$`=0.075`$, the models with the lowest mass ratios (D1 and D2) are stable, but the the other models show significant bar formation. The growth in the bar modes is very gradual at first, but at 3 Gyr models D3, D4, and D5 rapidly grow their bars, which peak at values of $`\delta _2^{max}`$ between 0.5 and 0.6 and slowly decline in amplitude. In the bottom panel of Figure 16, $`\mathrm{\Omega }_h`$ has been increased to 0.09. All the models begin to show bar structures at $``$ 2.5 Gyr and become bar-dominated by 3 Gyr. Models D3 and D4 make a rapid transformation to being bar-dominated, while $`\delta _2^{max}`$ for the other models grow more gradually. Regardless of the time it takes for the bar to grow, the values of $`\delta _2^{max}`$ for all the models is close to 0.4 once the bar has reached its peak strength. The values of $`\delta _2^{max}`$, the bar strength averaged over the last 0.8 Gyr of the simulation, do not show any correlation with mass ratio. For model D1, $`\delta _2^{max}`$ = 0.39, while for model D5 $`\delta _2^{max}`$ = 0.40. The halo-dominated models, including the fiducial model H5, are all stable against bar formation, even with a halo rotation rate of 0.09.

As stated above, the stability of the disks shown in Figure 16 is in part due to the mass of the bulge. We decreased the value of the bulge mass to $`0.15M_d`$ and re-ran our disk-dominated models to investigate bar formation in galaxies with less prominent bulges. Figure 17 shows the results. At $`\mathrm{\Omega }_h`$ = 0.06, models D3–D5 develop distinct bars while the other models remain stable. At $`\mathrm{\Omega }_h`$ = 0.075, all the models show some level of bar strength, with the average value of $`\delta _2^{max}`$ loosely correlating with increasing mass ratio. At $`\mathrm{\Omega }_h`$ = 0.09, all the models show rapid bar growth and become bar-dominated between 2 and 2.5 Gyr. The values of $`\delta _2^{max}`$ also correlate roughly with mass ratio. For model D1 $`\delta _2^{max}=0.30`$, increasing to $`\delta _2^{max}=0.41`$ for model D5. In these simulations, the growth of the bar begins earlier but is slower than the models with $`M_b=0.3M_d`$. The values of $`\delta _2^{max}`$ are slightly higher in the $`M_b=0.3M_d`$ simulations for those models which develop bars.

The bar pattern speeds do not appear to correlate with either $`\mathrm{\Omega }_h`$ or $`\delta _2^{max}`$. The simulations that develop strong bars all have bar pattern speeds near 0.5 radians per unit time (in the same units as $`\mathrm{\Omega }_h`$). The pattern speeds do not change on average once the bar has reached its peak amplitude. Since our simulations employ a static halo potential our results do not reflect the dynamical friction and angular momentum transfer that would surely occur in physical galaxies, reducing the pattern speed (Debattista & Sellwood 2000; Valenzuela & Klypin 2002).

## 6 Spiral Morphology

### 6.1 Spiral Arm Pitch Angle

Although in some simulations the peak in the Fourier spectrum is well defined and narrow, in others a simple search for the maximum value of $`|A_2|`$ can vary significantly in time if there is more than one peak in the spectrum. Pitch angles in spiral galaxies are known to vary with radius and are measured as an average value (Kennicutt 1981). Therefore, to quantify the morphology of the spiral structure in our simulations we calculate an average value of $`\alpha `$, summed from $`0\alpha 25`$ and weighted by the Fourier amplitude at each value of $`\alpha `$. The simulations for which this quantity, $`\alpha `$, were calculated did not contain bars as defined in §5, and so using $`\alpha =0`$ as the starting point of the weighted average does not bias the results to smaller values of $`\alpha `$. Only amplitudes greater than half the rms amplitude in that range are actually counted in the summation. The value of $`\alpha `$ varies in time about a mean value at a level of $``$ 5-20%. Fourier spectra with higher overall amplitude, usually found in models with higher values of $`\mathrm{\Omega }_h`$, have smaller rms fluctuations. The weighted average of $`\alpha `$ is again averaged over the time range 0.4–1.3 Gyr. The lower of these limits is chosen to allow the disk to reach equilibrium with the halo. As the simulations progress the arm strength decays, and so the integrations were stopped at 1.3 Gyr. In our notation, $`\alpha `$ is the weighted average of the Fourier spectrum for a given moment in time, and $`\alpha _t`$ is the average of this value over time range mentioned. The quantity $`\alpha _t`$ can then be converted into an average pitch angle, $`p`$.

The top panel in Figure 18 shows the evolution of $`p`$ with time for the halo-dominated models H1–H5 with $`\mathrm{\Omega }_h`$ = 0.09. Although there is some variation, the average value with time for each simulation is well-defined. The bottom panel in Figure 18 shows four realizations of the low-velocity model with $`M_b=0.3M_d`$, the last of which was resolved at $`3\times 10^5`$ particles rather than the fiducial value of $`10^5`$. The simulations were sampled with much higher frequency to investigate how $`\alpha `$ varies on small time scales. The results show little difference between all four realizations, showing that one realization at the fiducial resolution, although not ideal, is adequate for our results.

The clear trend in the top panel of Figure18 is that $`\alpha `$ decreases with decreasing mass ratio. Figure 19 shows the value of $`\alpha _t`$ plotted against mass ratio for all the models listed in Table 2. The results for $`\mathrm{\Omega }_h`$ = 0.09 have smaller error bars than those for $`\mathrm{\Omega }_h`$ = 0.06 because of the higher amplitude of the Fourier spectra for $`\mathrm{\Omega }_h`$ = 0.09. For both sets of simulations, the value of $`\alpha _t`$ is maximum at the fiducial mass ratio of 0.76 and decreases as the mass ratio gets bigger or smaller. As expected from the results of Figure 11, the data for $`\mathrm{\Omega }_h`$ = 0.09 are shifted to lower values of $`\alpha _t`$ as compared to the data for $`\mathrm{\Omega }_h`$ = 0.06. Figure 20 shows images of models H1–H4 at 655 Myr with $`\mathrm{\Omega }_h`$ = 0.09. The increased response of the disk with decreasing mass ratio is clearly evident. The spiral structure for model H2 can be traced for well over 2$`\pi `$ radians. The image of model H1, with the lowest mass disk, resembles a classic Kalnajs disk (e.g. Binney & Tremaine 1987). The caustic structures seen in the high-resolution run of the fiducial simulation are seen throughout the disk and are highly elliptical. Differential rotation has caused the oval caustics to overlap and create a spiral pattern.

### 6.2 Bulge to Disk Ratio and the Hubble Sequence

The third set of models described in §4, those with varying bulge mass, can be thought of as an ersatz Hubble sequence, at least along one dimension of the required criteria. The properties of the disk and halo were chosen as to produce significant amplitudes in the Fourier decomposition, reducing noise in the measured value of the pitch angle.

Figure 21 shows images of four models, $`M_b=0.1,0.3,0.9`$ and 2.0 with $`\mathrm{\Omega }_h`$ = 0.105. As the bulge mass varies from small to large, the spiral arms vary from loosely to tightly wound. The arms themselves vary in thickness; the larger the bulge mass, the narrower the spiral arms. For the model with $`M_b=0.1`$, the spiral arms can be traced to within one scale radius of the center. For $`M_b=2.0`$, the arms do not begin until outside two scale radii (although a high-resolution simulation could reveal structure inside this radius, as it did for the fiducial simulation). Figure 22 plots the average pitch angle $`p`$ against the logarithm of the bulge-to-disk ratio for each of these models. (Except for the model with $`M_b=0`$ of course, which had pitch angles of $`16.6^{}`$ and $`23.1^{}`$ for the halo rotation rates 0.075 and 0.105 respectively.) The plot shows a linear increase in $`p`$ with $`\mathrm{log}(M_b/M_d)`$. As the bulge mass increases, the spiral arms become tighter. The pitch angles for the simulations with $`\mathrm{\Omega }_h`$=0.105 are systematically higher than those with $`\mathrm{\Omega }_h`$=0.075. Plotted as well are data taken from Kennicutt (1981) for the pitch angles of galaxies along the Hubble sequence, from Sa-type to Sc-type. The bulge-to-disk ratios corresponding to these classifications were taken from Köppen & Arimoto (1990). The scatter in the observations, mentioned at the top of this section, is easily apparent.

The simulation data, however, follow the average trend of pitch angle with Hubble type. The slopes of the simulation data are flatter than that of the observations. Several explanations are viable for explaining this: The simulated Hubble sequence created here is artificial since it only varies one parameter, the bulge mass. There could be systematic variations of total galaxy mass and length scales along the Hubble sequence which could effect the spiral response of the disk. Only two values of $`\mathrm{\Omega }_h`$ were used, and the steeper slope of the observations could be reproduced if it were assumed that higher bulge mass correlated with lower halo rotation rates. Lastly, our definition of $`\alpha `$, since it is an average over a range of $`\alpha `$, would tend to return values that are closer to the center of that range, reducing the slope of the correlation. The pitch angles in Kennicutt (1981) were measured by obtaining the arm coordinates from H$`\alpha `$ or blue continuum observations and plotting these coordinates on the log $`r`$ vs $`\varphi `$ plane. A straight line was fitted to each arm and the results for the two main arms were averaged.

### 6.3 Pitch Angle and Rotational Velocity

As noted above, a better predictor than Hubble classification of spiral pitch angle is the rotational velocity of the galaxy (Kennicutt 1981). The relation found in the Kennicutt paper is plotted by the solid lines in Figure 23. These lines represent the mean and upper and lower bounds on his data. The slope and intercept of these lines were calculated by inspection, and are not meant to be quantitatively compared to our simulation results. Qualitatively, however, it is obvious that our simulations follow the same relation. Both sets of bulge simulations lie well within the bounds set by the observations. What is remarkable as well is that the simulations listed in Table 2 also show this trend. The simulations with the lower rotation velocity generally have higher pitch angles. Although the halo-dominated models do not return the same pitch angles as the disk-dominated models with the same rotation velocity, the trend is still preserved regardless of whether the mass dominating the kinematics is from the halo or the disk or the bulge. The scatter in our results lies within the scatter seen observationally.

The lowest rotation velocity in the models listed in Table 2 is D5, with $`v_{rot}=165`$ km s-1. To explore the response of systems with even lower velocities, we used the models described at the end of §3 with both low mass disks and halos. Three different bulge masses were used: $`M_b=`$ 0, $`0.15M_d`$ and $`0.3M_d`$. The first two of these models developed a strong and a weak bar respectively. The method of bar detection used in §4 allowed us to eliminate the particles that contribute to the bar from being used in calculating $`\alpha `$. Since the length of the bar is known, all particles with radii below the bar length were simply not counted in the summation in equation (12). The results are plotted in Figure 23, and they follow the observational relation.

A tighter correlation can be acheived by including what we know of the mass distribution and halo rotation. Figure 24 plots $`p`$ against the “normalized” rotational energy of the halo, a dimensionless quantity given by

$$^{}=_R_{}^{1/2}/(v_{rot}M_R_{}^{1/2}),$$
(16)
where $``$ is energy, $`v_{rot}`$ is the maximum rotation velocity of the disk also used in Figure 23, and $`M`$ is the total mass of the galaxy. The energy in the right hand side of equation (16) is the rotational energy of the halo, $`=I\mathrm{\Omega }_h^2`$, where $`I`$ is the moment of inertia of the halo. The subscript $`R_{}`$ on the energy and the total mass represents the radius at which these two quantities are calculated. The four windows in Figure 24 plot $`p`$ against $`^{}`$ calculated at four different values of $`R_{}`$. The correlation is apparent in all four panels, but is excellent for $`R_{}=3.25`$.

The choice of $`R_{}`$ significantly influences the slope of the correlation, with smaller values of $`R_{}`$ resulting in steeper slopes. The scatter changes only moderately with varying $`R_{}`$ for the first four sets of models in Figure 24, but effects the low-velocity models most. The quantity $`^{}`$ has a power-law dependence on $`R_{}`$, which is highest with the low-velocity models. Thus the values of $`^{}`$ for that set will be affected most by changes in $`R_{}`$.

## 7 Discussion

We have presented simulations of galaxy disks embedded in rigid, rotating, triaxial halos. We have constructed a suite of models in which the masses of the disk, bulge and halo have all been varied. With these models, the effects of changing the halo rotation rate and flattening have been investigated. Our results have focused on the formation of bars and spiral structure. Models with a mass ratio $`(M_d/M_h)_R_{}`$ $`>0.8`$ and $`\mathrm{\Omega }_h`$ $`=0.09`$ develop a strong bar, even if a significant bulge is present. For all models, the torque created by the halo induces spiral structure in the disk. The average pitch angle of this structure, $`p`$, is strongly anti-correlated with the maximum rotation speed of the model.

It is clear that estimations of the stability of disk galaxies cannot ignore the shape of the halo. Models with significant amounts of dark matter and a significant bulge (like D1 and D2, which have rotational velocities of over 200 km s-1 in our system of units) can become bar dominated with halo rotation rates of $`\mathrm{\Omega }_h`$ = 0.09 and greater. Maximal disk models, like D5, are very susceptible to bar instabilities if torqued by the halo. In all the simulations presented in §5, the only set of halo and bulge parameters for which models D3–D5 were stable against bar formation was $`M_b=0.3M_d`$ and $`\mathrm{\Omega }_h`$ = 0.06. Even with those parameters, model D3 still developed a low amplitude but persistent bar. A simple criterion for determining whether a disk galaxy will develop a bar, such as that determined by Efstathiou et al. (1982), is not readily apparent from the data. With more free parameters, $`M_b`$, $`q_b`$, and $`\mathrm{\Omega }_h`$, in our simulations, it may be that a simple criterion cannot quantify the stability of all our models used in our simulations.

The correlation of rotational velocity and pitch angle in our models provides a physical explanation of the same correlation found in observations of spiral galaxies. The pitch angle of the arms induced by an external torque is related to the total mass of the system, loosely through $`v_{max}`$, but more tightly to the dimensionless rotational energy parameter $`^{}`$. The differences seen in our results for different halo rotation rates and flattenings fit well with the amount of scatter seen in the correlation of $`v_{max}`$ with $`p`$ for real galaxies. This result does not have to be unique to the properties of a galaxy’s halo; it has been shown that encounters with satellite galaxies, which would also create a time-varying potential field in which to disrupt the disk, induce two-armed spiral structure in disk galaxies (e.g. Toomre & Toomre 1972). But encounters with smaller galaxies are by their nature transient and it remains that spiral structure is a ubiquitous feature of disk galaxies both with and without satellite perturbers.

In Kennicutt (1981), and references therein, the properties of the spiral arms are measured through their blue light, usually from HII regions, dust lanes, or the blue continuum, all of which are associated with star formation. Our collisionless simulations better represent the red stellar population which accounts for most of the mass in a galactic disk. Images of spiral galaxies in the near infrared have shown that spiral structure is evident in old stellar populations (Eskridge et al. 2002). The arms observed are smoother and wider than those seen in blue light, but generally follow the same morphology.

The Fourier decomposition method used in this paper has been applied to observational data of HII regions in spiral galaxies (Puerari & Dottori 1992; García-Goméz & Athanssoula 1993). The Fourier spectra obtained resemble in many ways the spectra of our simulated galaxies, often with several distinct peaks inside one spectrum. The pitch angle was calculated by these authors from the value of $`\alpha `$ for which the amplitude is highest. This definition of $`p`$ reproduces the general trend with Hubble type as described by Kennicutt (1981) but does not reflect the natural variations in pitch angle within a single galaxy seen both in the Fourier spectra of the observations and in Kennicutt’s calculations. Kennicutt noted that the variations of $`p`$ for a typical galaxy were significantly larger than the measurement error. An average value of the pitch angle better represents the multiple spiral components.

Hydrodynamic simulations by Bekki & Freeman (2002) have shown that the extended spiral structure seen in the blue compact dwarf galaxy NGC 2915 could be due to a rotating triaxial halo. Their simulations contained a disk fully composed of gas particles and a halo 100 times more massive than the disk, effectively making the disk massless. Our collisionless simulations show spiral structure in the inner regions of galactic disks. Since gas fractions in disk galaxies are usually no higher than 30%, the gas should be gravitationally coupled to the stars and follow the spiral patterns seen in our simulations. The dissipation and subsequent star formation in the gas would lead to self-consistent spiral structure which would outlast a collisionless disk and work to preserve its structure as well.

An open question is the correct value or range of values to use for $`\mathrm{\Omega }_h`$. The figure rotation of NGC 2915 has been calculated to be 8.0 km s-1 kpc-1 (Bureau et al. 1999), which is 0.107 in our $`N`$-body units. These authors also propose that halo figure rotation, a distinct phenomenon from halo angular momentum, is present in a significant fraction of dark matter halos in CDM simulations.

Regardless of the the details in determining $`\mathrm{\Omega }_h`$ or $`p`$, it has been clearly demonstrated that an external torqueing mechanism can reproduce the type of correlation observed between spiral pitch angle and rotational velocity. Even in simulations in which the magnitude of the applied torque is the same, i.e. the density profile, flattening, and rotation rate of that halo are all held constant, the pitch angle of the spiral structure varies according to the total mass of the system.

JLT would like to thank Richard Pogge for his helpful discussions, as well as Volker Springel and Stein Sigurdsson for help with the GADGET and SCF codes respectively. This work was performed at the Ohio Supercomputer Center under grant PAS0825.


format: text

text:
# References

Investigating Lorentz and CPT Symmetry with Antihydrogen111Invited talk at the XIVth Rencontres de Blois on Matter-Antimatter Asymmetry, Chateau de Blois, France, 16-22 June 2002
Neil Russell

Physics Department, Northern Michigan University, Marquette, MI 49855, U.S.A.
email: nrussell@nmu.edu
This talk discusses theoretical aspects of tests of CPT and Lorentz Symmetry that will in principle be possible with trapped antihydrogen. The framework is the standard-model extension, which admits minuscule violations of CPT and Lorentz symmetry in a general manner without giving up other features of the standard model of particle physics. Spectroscopic transitions in hydrogen and antihydrogen that exhibit leading-order effects are identified. Such comparisons of spectral frequencies in antimatter with the corresponding frequencies in regular matter will bound parameter combinations that are not accessible with regular-matter atoms alone.

1. The Standard-Model Extension
The laws of physics appear to be unchanged under the Lorentz and CPT transformations. These symmetries are built into the standard model of particle physics, and so calculations in the usual standard-model context cannot predict how violations of these symmetries might occur. However, the standard-model extension  allows symmetry-violating effects to be investigated in an explicit microscopic framework. In this proceedings, I give an overview of results in this context, mention the experimental findings, and then focus on hydrogen-antihydrogen symmetry tests investigated in collaboration with R. Bluhm and V.A. Kostelecký .

It is possible that a Lorentz- and CPT-symmetric fundamental theory underlying the standard model could violate these symmetries , perhaps through spontaneous symmetry breaking at lower energies . No satisfactory fundamental theory is known, but a candidate is string theory. A practical way to proceed is to perform calculations in the context of the standard-model extension . The standard-model extension adds to the conventional standard model all possible terms built from conventional fields that violate Lorentz symmetry while preserving observer coordinate independence. This involves a large set of parameters, which can be reduced by making simplifying assumptions. For example, it is often useful to require that the terms be gauge invariant, position independent, and renormalizable. Since the standard-model extension is expected to be the low-energy limit of an underlying theory, the symmetry violations would most likely be suppressed by ratios involving the low-energy mass and the Planck mass. Other theoretical issues about the standard-model extension include supersymmetry , noncommutative field theory , and anomaly cancellations .

The standard-model extension is applicable to all areas of physics. Applications of this framework include a possible mechanism for generating the baryon asymmetry in the universe . For the neutral mesons, the prediction of sidereal dependencies  in transition probabilities have been tested in the neutral kaons  and in the neutral-B systems . Other relevant aspects of neutral-meson physics include efforts to bound symmetry parameters in the neutral-D system  and studies of analogue models . In the photon sector, excellent bounds on Lorentz symmetry have been placed through analysis of light from distant cosmological sources . In the lepton sector, recent results on muonium and on anomaly-frequency comparisons of oppositely-charged muons at CERN and BNL have been obtained . Electron-positron comparisons using Penning traps have placed several bounds on standard-model extension parameters . An experiment with a spin-polarized torsion pendulum has obtained impressive results .

Clock-comparison experiments , comparing spectral lines of atomic transitions, are capable of Planck-scale tests of Lorentz symmetry . In the context of the standard-model extension, tests have been performed with Hg and Cs magnetometers , a noble-gas maser , and with an H maser . In such experiments, rotational symmetry is tested by monitoring frequency variations of a Zeeman hyperfine transition as the quantization axis of the clock changes direction with the Earth’s rotation. A practical way to do this is to operate two atomic clocks of different atomic species in the same laboratory, so as to avoid issues with signal propagation. To increase the sensitivity of clock-comparison experiments to other components of the standard-model extension couplings, and to address various other limitations, there is considerable interest in mounting experiments of this type on the International Space Station or other spacecraft .

Calculations in the standard-model extension context indicate that comparisons of matter and antimatter counterparts provide clean tests of CPT symmetry. It would be ideal to compare the narrow spectral lines of an atomic clock with the spectral lines of the corresponding antimatter clock. However, high-Z atoms like Rb or Cs used in current clocks are unlikely to be produced in antimatter form. For atoms of low atomic number, however, developments indicate that it may soon be possible to investigate the spectral lines of cold antimatter.

2. Antihydrogen
After experiments at CERN and Fermilab produced hot free antihydrogen ($`\overline{\mathrm{H}}`$) , the ATHENA and ATRAP collaborations  at CERN have worked towards the goals of obtaining cold and trapped antihydrogen . Numerous technological advances have been made along the way . The reported production of cold $`\overline{\mathrm{H}}`$ atoms  by the ATHENA collaboration marks a milestone on this path. When spectroscopy of $`\overline{\mathrm{H}}`$ atoms is performed, it will mark the first direct comparisons of neutral atoms with the corresponding antimatter atoms, and clean tests of CPT symmetry are likely to result.

From the experimental perspective, a comparison of the 1S to 2S transition in the two systems would be natural because of the parts in 1014 and 1012 precisions possible with hydrogen (H) beams  and with trapped H . In principle, the maximum possible precision is about a part in 1018, or about a 1-mHz resolution . With the same precision in $`\overline{\mathrm{H}}`$, an excellent test of CPT would be achieved.

From the theoretical perspective, there are many possible transitions in H and $`\overline{\mathrm{H}}`$ atoms, and the challenge is to identify those that are most likely to show CPT-violating differences. The point is that some will be suppressed by powers of the fine structure constant, or perhaps by the speed of the atoms in the reference frame. The standard-model extension provides a calculational framework of powerful generality for analyzing the various transitions in H and $`\overline{\mathrm{H}}`$. In the following sections, I consider the spectroscopy of free atoms, and then atoms confined in magnetic traps. The 1S to 2S transitions and the hyperfine transitions are considered.

3. Free Hydrogen and Antihydrogen
The effect of the standard-model extension on free atoms of H or $`\overline{\mathrm{H}}`$ can be studied by calculating the perturbative shifts in energy levels of the quantum stationary states, starting from the modified Dirac equation,

$$\left(i\gamma ^\mu D_\mu m_ea_\mu ^e\gamma ^\mu b_\mu ^e\gamma _5\gamma ^\mu \frac{1}{2}H_{\mu \nu }^e\sigma ^{\mu \nu }+ic_{\mu \nu }^e\gamma ^\mu D^\nu +id_{\mu \nu }^e\gamma _5\gamma ^\mu D^\nu \right)\psi =0.$$
(1)
In this expression, $`\psi `$ is a four-component electron field describing an electron of mass $`m_e`$ and charge $`q=|e|`$. We define $`iD_\mu i_\mu qA_\mu `$, and use units with $`\mathrm{}=c=1`$. The field interacts with the proton Coulomb potential $`A^\mu =(|e|/4\pi r,0)`$. The coefficients $`a_\mu ^e`$, $`b_\mu ^e`$, $`H_{\mu \nu }^e`$, $`c_{\mu \nu }^e`$ and $`d_{\mu \nu }^e`$ parameterize Lorentz violation, and the coefficients $`a_\mu ^e`$ and $`b_\mu ^e`$ also parameterize CPT violation. All are expected to be minuscule . Energy shifts in the H and $`\overline{\mathrm{H}}`$ systems could also be due to perturbative couplings $`a_\mu ^p`$, $`b_\mu ^p`$, $`H_{\mu \nu }^p`$, $`c_{\mu \nu }^p`$, and $`d_{\mu \nu }^p`$ in a modified Dirac equation for free protons.

In conventional relativistic quantum mechanics, the H and $`\overline{\mathrm{H}}`$ hamiltonians are identical: reversing the charges of the particles leaves the interaction potential unchanged, and the reduced mass is also unaffected. Thus, the eigenfunctions and associated energies are the same for H and $`\overline{\mathrm{H}}`$, and within conventional physics there should be no observable differences between the spectra of H and $`\overline{\mathrm{H}}`$. This argument carries over to all perturbative effects from conventional quantum electrodynamics.

If the Lorentz- and CPT-violating couplings are included in a perturbative calculation, the energy shifts in the spectrum of the H electron can differ from those for the $`\overline{\mathrm{H}}`$ positron. To make the comparison, the perturbative effect in $`\overline{\mathrm{H}}`$ is found from Eq. (1) by a method involving charge conjugation and appropriate field redefinitions . The energy perturbations arising from the proton and antiproton may also be included by using relativistic two-fermion techniques . In the uncoupled basis with electronic angular momentum $`J`$ and nuclear angular momentum $`I`$, and corresponding third components $`m_J`$, $`m_I`$, the perturbative Lorentz-violating energy corrections for the basis states $`|m_J,m_I`$ may be found. With this approach, the energy corrections for protons or antiprotons have the same mathematical form as those for electrons or positrons, except for the replacement of superscripts $`e`$ with $`p`$. The shifts in the 1S and 2S energy levels for free H are found to be matched at leading order:

$`\mathrm{\Delta }E^H(m_J,m_I)`$ $``$ $`(a_0^e+a_0^pc_{00}^em_ec_{00}^pm_p)`$ (4)
$`+(b_3^e+d_{30}^em_e+H_{12}^e)m_J/|m_J|`$
$`+(b_3^p+d_{30}^pm_p+H_{12}^p)m_I/|m_I|,`$
where $`m_p`$ is the proton mass. The same is true for $`\overline{\mathrm{H}}`$; the leading-order energy shifts $`\mathrm{\Delta }E^{\overline{H}}`$ in the 1S and 2S levels are identical. They are obtained from the expression (4) with the substitutions $`a_\mu ^ea_\mu ^e`$, $`d_{\mu \nu }^ed_{\mu \nu }^e`$, $`H_{\mu \nu }^eH_{\mu \nu }^e`$; $`a_\mu ^pa_\mu ^p`$, $`d_{\mu \nu }^pd_{\mu \nu }^p`$, $`H_{\mu \nu }^pH_{\mu \nu }^p`$.

To find the effect of the Lorentz-violation on the frequency spectrum of H and $`\overline{\mathrm{H}}`$, we must calculate appropriate energy differences. Since the electron and proton spins in H interact through the hyperfine interaction, the appropriate basis is the coupled basis $`|F,m_F`$, where $`F`$ is the total angular momentum, and $`m_F`$ is the third component. Similarly, this is true for the positron and antiproton spins in $`\overline{\mathrm{H}}`$. There are four allowed 1S–2S two-photon transitions in H and $`\overline{\mathrm{H}}`$, determined by the selection rules $`\mathrm{\Delta }F=0`$ and $`\mathrm{\Delta }m_F=0`$ . In each, the spins remain unchanged. It follows from (4) that the frequencies are unshifted at leading order since both energy levels in the transition are increased or decreased by identical amounts. So, the standard-model extension shows no leading-order evidence of Lorentz or CPT violation in 1S–2S spectroscopy of free H or $`\overline{\mathrm{H}}`$.

One way to proceed is to consider the next-largest energy-level shifts in the Lorentz-violating couplings for H and $`\overline{\mathrm{H}}`$. They are suppressed by a factor of $`\alpha ^2`$. However, the shifts differ in some cases, and so effects observable in principle could occur. One possibility is the $`m_F=1m_F^{}=1`$ line, which is shifted by

$$\delta \nu _{1S2S}^H\alpha ^2b_3^e/8\pi $$
(5)
relative to the unshifted $`m_F=0m_F^{}=0`$ line. Suppression by a factor at least of order $`\alpha ^25\times 10^5`$ would also occur in the proton-antiproton corrections to the 1S to 2S lines.

Because of these suppression effects in free H and $`\overline{\mathrm{H}}`$ lines, it is important to note that electron-position anomaly-frequency comparisons in Penning traps  are capable of producing tighter bounds on similar combinations of Lorentz and CPT violating parameters. Similar g-2 experiments with protons and antiprotons could in principle also place tighter bounds on proton-antiproton parameters than would be possible with 1S to 2S spectroscopy.

The advantages of Penning-trap experiments over spectroscopy in free H and antihydrogen may seem counterintuitive given that g-2 experiments compare gyromagnetic ratios with fractional resolutions of about $`2\times 10^{12}`$, as compared to the ultimate fractional resolution of about $`10^{18}`$ in H, and perhaps $`\overline{\mathrm{H}}`$, experiments. However, bounds on the coefficients of the Lorentz- and CPT-violating terms are determined by absolute, not relative, resolutions. In principle, 1S to 2S transitions could reach absolute resolutions of about one mHz, as compared with Penning-trap frequencies measured with about one-Hz resolution. However, the suppression in the former by about five orders of magnitude would effectively give the Penning-trap experiments a 100-fold advantage in placing bounds. These limitations of H and antihydrogen in testing Lorentz and CPT symmetry apply only to the free case where the transitions maintain the spin direction. Other considerations apply to H and antihydrogen confined by external fields.

4. Trapped hydrogen and antihydrogen
Given the technical difficulties of creating antihydrogen atoms, the advantages of trapping them appear clear. From the perspective of testing Lorentz and CPT symmetry, the external trapping fields offer the possibility of separating different spin states and hence of avoiding the suppression effects seen in the previous case of free H and $`\overline{\mathrm{H}}`$. We limit our discussion to the case of H or antihydrogen in a constant uniform magnetic field. This provides an approximation to the physical environment of the Ioffe-Pritchard trap , relevant in the CERN antihydrogen experiments.

In a uniform magnetic field $`B`$, the four eigenstates corresponding to any given principal quantum number $`n`$ are dependent on the magnetic field. In increasing energy, we denote each of the 1S and 2S hyperfine Zeeman levels by $`|a_n`$, $`|b_n`$, $`|c_n`$, $`|d_n`$, with $`n=1`$ or $`2`$, for both H and $`\overline{\mathrm{H}}`$. The states for H, expressed in terms of the uncoupled basis $`|m_J,m_I`$ are

$`|d_n`$ $`=`$ $`|\frac{1}{2},\frac{1}{2},`$
$`|c_n`$ $`=`$ $`\mathrm{sin}\theta _n|\frac{1}{2},\frac{1}{2}+\mathrm{cos}\theta _n|\frac{1}{2},\frac{1}{2},`$
$`|b_n`$ $`=`$ $`|\frac{1}{2},\frac{1}{2},`$
$`|a_n`$ $`=`$ $`\mathrm{cos}\theta _n|\frac{1}{2},\frac{1}{2}\mathrm{sin}\theta _n|\frac{1}{2},\frac{1}{2}.`$ (6)
The mixing angles $`\theta _n`$ are magnetic-field dependent, and differ for the 1S and 2S states:

$$\mathrm{tan}2\theta _n\frac{(51\mathrm{mT})}{n^3B}.$$
(7)
The energies of the $`|c_1`$ and $`|d_1`$ states increase with increasing $`|B|`$, and so are low-field seekers. In principle, they remain confined near the region of minimum field in the trap. Collisions within the trap, however, lead to spin exchange collisions $`|c_1+|c_1|b_1+|d_1`$ that tend to decrease the $`|c_1`$ population, and create a predominance of confined $`|d_1`$ states.

The transition between these unmixed-spin states would appear attractive since it avoids broadening due to magnetic-field dependence. A Lorentz or CPT test could be envisaged involving comparison of the 1S-2S transition $`|d_1|d_2`$ in H and in $`\overline{\mathrm{H}}`$. However, in both H and $`\overline{\mathrm{H}}`$, the spin configurations of the $`|d_1`$ and $`|d_2`$ states are the same, so there are again no leading-order shifts in these frequencies.

A comparison between the $`|c_1|c_2`$ transitions in H and antihydrogen is of more interest for testing Lorentz or CPT symmetry. The difference here is that the spin-mixing coefficients $`\mathrm{sin}\theta _n`$ and $`\mathrm{cos}\theta _n`$ depend on the principal quantum number $`n`$, and lead to an unsuppressed shift in the corresponding frequency. For H:

$$\delta \nu _c^H\kappa (b_3^eb_3^pd_{30}^em_e+d_{30}^pm_pH_{12}^e+H_{12}^p)/2\pi ,$$
(8)
where the spin-mixing function $`\kappa `$ is

$$\kappa \mathrm{cos}2\theta _2\mathrm{cos}2\theta _1.$$
(9)
It is positive, less than one, and has maximum value $`\kappa 0.67`$ occurring at $`B0.011`$ T. At this field value the atoms are optimally sensitive to any Lorentz- or CPT-violating effects.

The frequency $`\nu _c^H`$ depends on spatial components of Lorentz-violating couplings and would therefore vary sidereally due to the rotation of the Earth. This effect means it is possible in principle to detect Lorentz-violating signals with H alone. Such clock-comparison tests  involve comparing two different transitions, and have been done with a H maser  and other atomic systems . A limitation of these experiments is that they bound combinations of couplings that violate Lorentz symmetry, but do not isolate CPT-violating couplings alone. The analysis can be complicated because of the geometry of motion of the laboratory relative to a suitable inertial reference frame.

In the case of $`\overline{\mathrm{H}}`$, the calculation of $`\delta \nu _c^{\overline{H}}`$ for the same magnetic field direction and magnitude involves using the opposite positron and antiproton spins compared to the electron and proton spins in H. The result is identical to $`\delta \nu _c^H`$ in Eq. (8) except that the signs of $`b_3^e`$ and $`b_3^p`$ are reversed:

$$\delta \nu _c^{\overline{H}}\kappa (b_3^e+b_3^pd_{30}^em_e+d_{30}^pm_pH_{12}^e+H_{12}^p)/2\pi .$$
(10)
The result (10) can be compared with (8), to give the instantaneous difference

$$\mathrm{\Delta }\nu _{1S2S,c}\nu _c^H\nu _c^{\overline{H}}\kappa (b_3^eb_3^p)/\pi .$$
(11)
This expression involves couplings that are both CPT violating, and so provides a clean test of CPT symmetry. It should be noted that the H and antihydrogen measurements need to be made at effectively equal times to ensure that the magnetic field has the same orientation for both cases.

5. Hyperfine Transitions
Hyperfine transitions within the 1S level of H can be measured with accuracies exceeding 1 mHz in masers . So transitions of this type in trapped H and antihydrogen are interesting candidates for performing tests of Lorentz or CPT symmetry.

To find the energy shifts in the possible transitions within the 1S level of H or $`\overline{\mathrm{H}}`$, we may ignore the constant shift of $`a_0^e+a_0^pc_{00}^em_ec_{00}^pm_p`$ occurring in each of them, since this cannot contribute to the difference in energy levels. The terms of relevance are

$`\mathrm{\Delta }E_a^H`$ $``$ $`\widehat{\kappa }(b_3^eb_3^pd_{30}^em_e+d_{30}^pm_pH_{12}^e+H_{12}^p),`$
$`\mathrm{\Delta }E_b^H`$ $``$ $`b_3^e+b_3^pd_{30}^em_ed_{30}^pm_pH_{12}^eH_{12}^p,`$
$`\mathrm{\Delta }E_c^H`$ $``$ $`\mathrm{\Delta }E_a^H,`$ (12)
$`\mathrm{\Delta }E_d^H`$ $``$ $`\mathrm{\Delta }E_b^H,`$ (13)
where $`\widehat{\kappa }\mathrm{cos}2\theta _1`$. Conventionally, the triplet $`|b_1`$, $`|c_1`$, and $`|d_1`$ is degenerate at zero field. However, with the presence of Lorentz- and CPT-violating couplings, this degeneracy is broken even at zero field: equal and opposite energy shifts occur for $`|b_1`$ and $`|d_1`$. In fact, the difference in the $`|d_1|a_1`$ and $`|b_1|a_1`$ transitions is unsuppressed and varies sidereally,

$$|\mathrm{\Delta }\nu _{db}^H||b_3^e+b_3^pd_{30}^em_ed_{30}^pm_pH_{12}^eH_{12}^p|/\pi .$$
(14)
In a nonzero magnetic field, energies are shifted in all four hyperfine Zeeman levels. The function $`\widehat{\kappa }\mathrm{cos}2\theta _1`$ controls the shifts in $`|a_1`$ and $`|c_1`$. As $`B`$ increases, so does $`\widehat{\kappa }`$, reaching $`\widehat{\kappa }1`$ when $`B0.3`$ T.

The leading-order effects from CPT and Lorentz violation in the standard maser line, $`|c_1|a_1`$, are suppressed. This is because at the magnetic field of about $`B\text{ }<10^6`$ Tesla, the function $`\widehat{\kappa }`$ is approximately $`\widehat{\kappa }\text{ }<10^4`$. Two other candidate transitions for detecting unsuppressed Lorentz or CPT violation are in the field-dependent transitions $`|d_1|a_1`$ and $`|b_1|a_1`$. The difference between these two frequencies is $`\mathrm{\Delta }\nu _{db}^H`$, and could be used to study unsuppressed sidereal variations. One might, however, expect experimental limitations such as broadening of the line due to field dependence.

Another possibility would be to measure radio-frequency transitions between states within the triplet of hyperfine levels in H and $`\overline{\mathrm{H}}`$. This would also offer the possibility of working at a field-independent point where the limitations of magnetic-field inhomogeneities can be eliminated while maintaining sensitivity to Lorentz- and CPT-violating effects. An interesting option would be spectroscopy of the $`|d_1|c_1`$ transition at the field-independent point $`B0.65`$ Tesla. One might expect various experimental challenges requiring high homogeneity of the field and extremely low temperatures. Assuming these issues can be handled, frequency resolutions of order 1 mHz might be envisaged.

For this $`|d_1|c_1`$ transition at this field of $`B0.65`$ Tesla, the electron and proton spins in state $`|c_1`$ interact more strongly with the field than with each other; they are highly polarized with $`m_J=1/2`$ and $`m_I=1/2`$. So, the transition $`|d_1|c_1`$ is basically a proton spin-flip, and the frequency shifts are found to be

$$\delta \nu _{cd}^H(b_3^p+d_{30}^pm_p+H_{12}^p)/\pi ,\delta \nu _{cd}^{\overline{H}}(b_3^p+d_{30}^pm_p+H_{12}^p)/\pi $$
(15)
for H and $`\overline{\mathrm{H}}`$. One would thus expect these frequencies to exhibit sidereal variations. Another option would be to consider the instantaneous difference in these variations,

$$\mathrm{\Delta }\nu _{cd}\nu _{cd}^H\nu _{cd}^{\overline{H}}2b_3^p/\pi .$$
(16)
This comparison could provide a clean and unsuppressed test of the CPT-violating coupling $`b_3^p`$ for the proton. If, for example, a frequency resolution of 1 mHz were attained, this would correspond to an upper bound of about

$$|b_3^p|\text{ }<10^{27}GeV$$
(17)
Compared with estimated bounds in Penning traps, this is an improvement of about 1000 . It is also more than 10,000 times better than estimated bounds attainable from 1S-2S transitions.

Clock-comparison experiments are able to resolve spectral lines to about 1 $`\mu `$Hz . However, the coupling expressions they bound  are complex and isolation of $`b_3^p`$ is difficult, due to the structure of the nuclei involved. We also note that the experiments discussed here are sensitive only to spatial components of Lorentz-violating couplings. Sensitivity to timelike couplings, like $`b_0^e`$, would require a boost, but would magnify CPT- and Lorentz-violating effects .

6. Conclusion
The standard-model extension has been used to analyze various transitions in the H and antihydrogen systems, with the aim of understanding ways in which Lorentz and CPT symmetry violations may occur. In free H and $`\overline{\mathrm{H}}`$, signals for Lorentz and CPT violation are suppressed by at least one power of the fine-structure constant. For trapped H or $`\overline{\mathrm{H}}`$ atoms, 1S-2S transitions involving the mixed-spin $`|c`$ states, or the spin-flip $`|d_1|c_1`$ hyperfine transition could give rise to unsuppressed signals of Lorentz and CPT violation. Such signals would indicate qualitatively new physics existing in a fundamental theory at the Planck scale.

Acknowledgments
I thank Robert Bluhm and Alan Kostelecký for their collaboration. I am grateful for support from the organizers of the Rencontres de Blois XIV and from Northern Michigan University.

References


text:
# Analytical approach to electromagnetic processes in loosely bound nuclei: application to {^𝐵}8

## Abstract

In this paper we develop an analytical model in order to study electromagnetic processes involving loosely bound neutron–rich and proton–rich nuclei. We construct a model wave function, to describe loosely bound few–body systems, having the correct behaviour both at large and small distances. The continuum states are approximated by regular Coulomb functions. As a test case we consider the two–body Coulomb dissociation of $`{}_{}{}^{B}8`$ and, the inverse, radiative capture reaction. The difference between using a pure two–body model and the results obtained when incorporating many–body effects, is investigated. We conclude that the interpretation of experimental data is highly model dependent and stress the importance of measuring few–body channels.

journal: Phys. Lett. B

, and Electromagnetic processes such as radiative capture, photo–dissociation and Coulomb dissociation have always been, and still remain, an excellent tool to investigate nuclear structure. Unfortunately, for radioactive nuclei the radiative capture experiments are very difficult and studies of photo–dissociation are virtually impossible. However, with the advent of radioactive beam facilities, nuclear structure of dripline nuclei can be studied using Coulomb dissociation on heavy targets. The Coulomb breakup of these nuclei is of interest also in nuclear astrophysics, since it can be related to the corresponding radiative capture process at astrophysical energies . In this paper, we will present an analytical model for Coulomb dissociation (and consequently for radiative capture) of loosely bound dripline nuclei (see also our recent conference contribution ). We want to stress that our approach will be general in the sense that both neutron–rich, and proton–rich systems can be studied. The possibility to study reactions analytically is very appealing since it often allows a deeper physical understanding of the process. For some cases such studies are not only possible but might even give results that are directly comparable to experimental data and more advanced, numerical calculations. In particular, this is the case for electromagnetic transitions between a loosely bound nuclear state and a pure Coulomb continuum. “Loosely bound” implies that the nucleus will exhibit a large degree of clusterization. We will consider the case of two clusters, but will also discuss the three–cluster channel (see also similar approaches for one–neutron  and two–neutron  halo nuclei). After a general description of our model we will exemplify it with an application to $`{}_{}{}^{B}8`$. Our main focus will be the nuclear structure effects, and we will utilize the advanced three–body model of Grigorenko et al. , while the Coulomb dissociation is considered only to first order. As it was shown recently , higher order effects should not be significant at beam energies higher than around 70 MeV/$`A`$.

Our starting point for calculating electromagnetic cross sections will be the $`\mathrm{E}\lambda `$ strength function for a transition from an initial state to a final, continuum state with energy $`E`$

$$\frac{\mathrm{d}B(\mathrm{E}\lambda )}{\mathrm{d}E}=\frac{1}{2J_i+1}\underset{f}{}\tau _f|f||(\mathrm{E}\lambda )||i|^2\delta \left(E_fE\right),$$
(1)
where $`\mathrm{d}\tau _f`$ is the phase space element for final states, $`(\mathrm{E}\lambda )`$ is the electric multipole operator and $`|i`$, $`|f`$ are the initial and final states in the center of mass subsystem.

We will consider loosely bound systems of two clusters $`(c+x)`$ and, in particular, we will study transitions to the low–energy continuum in which excitations are manifested as relative motion between the clusters $`E=\mathrm{}^2k^2/2\mu _{cx}`$, where $`\mu _{cx}`$ is the reduced mass of the system. Introducing the intercluster distance, $`r`$, the corresponding cluster $`\mathrm{E}\lambda `$ operator (operating only on the relative motion of clusters) is

$$(\mathrm{E}\lambda ,\mu )=eZ(\lambda )r^\lambda Y_{\lambda \mu }(\widehat{r}),$$
(2)
with the effective multipole charge $`Z(\lambda )=\mu _{cx}^\lambda (Z_x/m_x^\lambda +(1)^\lambda Z_c/m_c^\lambda )`$.

The strength function is the key to study several reactions. The photo–dissociation $`(A+\gamma c+x)`$ cross section is given by

$$\sigma _\gamma ^{\mathrm{E}\lambda }(E_\gamma )=\frac{(2\pi )^3(\lambda +1)}{\lambda [(2\lambda +1)!!]^2}\left(\frac{E_\gamma }{\mathrm{}c}\right)^{2\lambda 1}\frac{\mathrm{d}B(\mathrm{E}\lambda )}{\mathrm{d}E},$$
(3)
where the photon energy, $`E_\gamma =E+E_0`$, is larger than the binding energy, $`E_0`$. The inverse radiative capture reaction can be studied using detailed balance

$$\sigma _{\mathrm{rc}}^{\mathrm{E}\lambda }(E)=\left(\frac{E_\gamma }{\mathrm{}ck}\right)^2\frac{2(2J_A+1)}{(2J_c+1)(2J_x+1)}\sigma _\gamma ^{\mathrm{E}\lambda }(E_\gamma ),$$
(4)
where $`J_i`$ is the spin of particle $`i`$. Using first order perturbation theory, and the method of virtual quanta , the energy spectrum for Coulomb dissociation on a high–$`Z`$ target can be written as a sum over multipole, $`\pi \lambda `$, photo–dissociation cross sections multiplied by the corresponding spectra of virtual photons, $`n_{\pi \lambda }(E_\gamma )`$,

$$\frac{\mathrm{d}\sigma _\mathrm{C}}{\mathrm{d}E}=\underset{\mathrm{E}\lambda }{}\frac{n_{\mathrm{E}\lambda }(E_\gamma )}{E_\gamma }\sigma _\gamma ^{\mathrm{E}\lambda }(E_\gamma )+\underset{\mathrm{M}\lambda }{}\frac{n_{\mathrm{M}\lambda }(E_\gamma )}{E_\gamma }\sigma _\gamma ^{\mathrm{M}\lambda }(E_\gamma ).$$
(5)
Note that since M$`\lambda `$ transitions are usually  strongly suppressed we will not consider them in this work. In order to study all these reactions analytically we will propose a model function to describe the (loosely) bound state of a two–body system. We will only be interested in direct transitions to a “clean” continuum, i.e., with all nuclear phase shifts equal to zero. Thus, the final state, with angular momentum $`l_f`$ between the clusters, will be described by a regular Coulomb function

$`\varphi _{l_f}(k,r)`$ $`=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{1}{k}}i^{l_f}e^{i\sigma _{l_f}}F_{l_f}(k,r),\mathrm{where}`$ (6)
$`F_{l_f}(k,r)`$ $`=C_{l_f}(k)e^{ikr}(kr)^{l_f+1}{}_{1}{}^{}F_{1}^{}(l_f+1+i\eta (k);2l_f+2;2ikr),`$ (7)
and $`\sigma _{l_f}`$ is the Coulomb phase, $`\eta (k)=Z_cZ_xe^2\mu _{cx}/\mathrm{}^2k`$ is the Sommerfeld parameter, $`C_{l_f}(k)=2^{l_f}e^{\pi \eta (k)/2}|\mathrm{\Gamma }(l_f+1+i\eta (k))|/(2l_f+1)!`$, and $`{}_{1}{}^{}F_{1}^{}(a;b;z)`$ is the confluent hypergeometric function .

The reduced matrix element introduced in the definition of the strength function, Eq. (1), contains a radial integral. With our approximation for the continuum state this integral takes the form

$$I_{l_f}(k)=_0^{\mathrm{}}dre^{ikr}r^{l_f+1}{}_{1}{}^{}F_{1}^{}(l_f+1i\eta (k);2l_f+2;2ikr)r^\lambda \varphi _b(r).$$
(8)
Here, $`\varphi _b(r)`$ is the two–body, relative motion wave function (WF) describing the initial state. For large $`r`$, and with angular momentum $`l_i`$ between the clusters, this radial function should be proportional to the Whittaker function $`W_{\eta _0,l_i+1/2}(2\kappa _0r)`$ (see Ref. ), where $`\eta _0=Z_cZ_xe^2\mu _{cx}/\mathrm{}^2\kappa _0`$ and $`E_0=\mathrm{}^2\kappa _0^2/2\mu _{cx}`$ is the binding energy.

In most studies on loosely bound systems, the Whittaker function has been used to describe the bound state for all $`r`$. However, this approximation is only motivated if the transition matrix element is dominated by contributions from very large $`r`$. This should be the case for reactions at very small energies. For real experimental energies ($`E100`$ keV), the WF of the bound state should be constructed in a more realistic way. Therefore we will introduce a model function that describes the bound state $`c+x`$ WF accurately for all distances. This is done by considering the behavior at small and large $`r`$. We have already pointed out the the WF should be described by a Whittaker function at large $`r`$. Furthermore, for a two–body system, consisting of point–like particles, it should behave as $`r^{l_i+1}`$ as $`r0`$. Both asymptotics are fulfilled using the following model function

$$\varphi _{b,\stackrel{~}{\gamma }}^{\mathrm{`}\mathrm{`}\mathrm{exact}^{\prime \prime }}(r)=\frac{1}{\sqrt{N_{\stackrel{~}{\gamma }}}}W_{\eta _0,l_i+1/2}(2\kappa _0r)\left(1e^{\kappa _1r}\right)^{2l_i+1},$$
(9)
where $`N_{\stackrel{~}{\gamma }}`$ is the normalization constant and $`\stackrel{~}{\gamma }`$ denotes the parameters $`\{\kappa _0,\kappa _1,\eta _0\}`$. The parameters $`\kappa _0`$ and $`\eta _0`$ are defined by the binding energy, charges and masses, while $`\kappa _1`$ can be fitted to give the correct size of the system. Using this WF, and solving the integral (8) numerically, it is possible to get very good estimates for the electromagnetic reaction cross sections. As a remark we want to point out that, for a one–neutron halo nucleus ($`\eta _0=0`$), the Whittaker function will transform into a modified, spherical Bessel function: $`W_{0,l_i+1/2}(2\kappa _0r)=\sqrt{2\kappa _0r/\pi }K_{l_i+1/2}(\kappa _0r)`$. In this case, the integral (8) can be solved exactly.

However, we are searching for a completely analytical model which will also enable us to incorporate many–body effects. Our model function (9) has to be modified accordingly. First, we note that the asymptotic form of the Whittaker function as $`r\mathrm{}`$ is

$$W_{\eta _0,l_i+1/2}(2\kappa _0r)e^{\kappa _0r}r^{\eta _0}\left[1+𝒪(1/r)\right].$$
(10)
Secondly, for two–body systems in which the clusters have an internal structure, the centrifugal barrier is effectively larger and the WF should behave as $`r^n`$, where $`n>l_i+1`$, as $`r0`$.

Motivated by this, we put forward the following model function

$`\varphi _{b,\gamma }(r)`$ $`={\displaystyle \frac{1}{\sqrt{N_\gamma }}}{\displaystyle \frac{e^{\kappa _0r}}{r^{\eta _0^{}}}}\left(1e^{\kappa _1r}\right)^p,\mathrm{with}\mathrm{norm}`$ (11)
$`N_\gamma `$ $`={\displaystyle \underset{m=0}{\overset{2p}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2p}{m}}\right)(1)^m(2\kappa _0+m\kappa _1)^{2\eta _0^{}1}\mathrm{\Gamma }(12\eta _0^{}),`$
where $`\gamma `$ denotes the parameters $`\{\kappa _0,\kappa _1,\eta _0^{},p\}`$. The parameter $`\kappa _0`$ is defined by the binding energy and effective mass. By putting $`\eta _0^{}=\eta _0`$ we would ensure to reproduce the tail of the WF at very large $`r`$. However, higher order terms in Eq. (10) remain important for $`r100`$ fm. Therefore, $`\eta _0^{}`$ and $`\kappa _1`$ are used as free parameters in a fit to the “exact” WF (9) in the interval of interest. In this way $`\eta _0^{}`$ will be an “effective” Sommerfeld parameter while $`\kappa _1`$ will still mainly be connected with the size. Finally, the integer $`p`$ is fixed by the small $`r`$ behaviour, $`r^n`$, and will thus allow us to take many–body effects into account.

It is now possible to make the integral (8) analytically (see Ref. )

$$\begin{array}{c}\hfill I_{l_f,\gamma }(k)=\frac{1}{\sqrt{N_\gamma }}\underset{m=0}{\overset{p}{}}\left(\genfrac{}{}{0pt}{}{p}{m}\right)(1)^m(m\kappa _1+\kappa _0+ik)^{(l_f+2+\lambda \eta _0^{})}\mathrm{\Gamma }(l_f+2+\lambda \eta _0^{})\\ \hfill \times {}_{2}{}^{}F_{1}^{}(l_f+2+\lambda \eta _0^{};l_f+1i\eta (k);2l_f+2;\frac{2ik}{m\kappa _1+\kappa _0+ik}).\end{array}$$
(12)
Many-body nuclear structure can further be taken into account by considering the possibility that the bound state WF contains several different two–body components

$$\varphi _b(r)=\underset{i}{}a_i\varphi _{b,\gamma _i}(r).$$
(13)
Note that pure many–body components will not contribute to two–body breakup, but will instead lead to $`_ia_i^2<1`$. Note also that the threshold for two–body breakup will be higher for components where one (or both) of the clusters is excited. Therefore, we define the continuum strength function separately for each component. Finally, we arrive at an analytical formula for the strength function

$$\frac{\mathrm{d}B(\mathrm{E}\lambda )}{\mathrm{d}E}|_i=\frac{e^2Z^2(\lambda )\mu _{cx}}{\mathrm{}^2}\frac{2\lambda +1}{2\pi ^2}\underset{l_f}{}a_i^2k^{2l_f+1}C_{l_f}^2(k)l_i0\lambda 0|l_f0^2|I_{l_f,\gamma _i}(k)|^2.$$
(14)
As an example we will apply our analytical approach to the 8B nucleus. Our key point here will be to show the applicability of our analytical model and to demonstrate the importance of many–body nuclear structure. The interest in 8B stems from its key role in the production of high-energy solar neutrinos. The probability of the reaction $`{}_{}{}^{B}7`$$`+p`$$`{}_{}{}^{B}8`$$`+\gamma `$ at solar energies strongly depends on the structure of $`{}_{}{}^{B}8`$ and, in particular, on the asymptotics of the valence proton WF. The reaction can be studied indirectly through Coulomb dissociation, using a radioactive 8B beam on a heavy target . We should also mention the recent progress in radiative capture measurements , where the cross section has been measured at energies around 200 keV with an accuracy of $`15`$ percent. However, in all cases theoretical models are needed to extrapolate the measured cross sections down to solar energies.

The low–lying $`{}_{}{}^{B}8`$ continuum can, with relatively good precision, be approximated as a pure Coulomb one. There are no negative parity states at low excitation energies  and the electromagnetic processes are, in all cases we are considering, dominated by $`\mathrm{E1}`$ transitions.

We start by treating $`{}_{}{}^{B}8`$ as a pure two–body ($`{}_{}{}^{B}7+p`$) system with binding energy $`E_0=137`$ keV and relative orbital momentum $`l_i=1`$. The single free parameter, $`\kappa _1`$, in our “exact” model function (9) is then fitted to an rms intercluster distance of $`r_{\mathrm{rms}}=4.57`$ fm (extracted from Ref. ). In order to get analytical results, we then introduce the model function (11). With $`p=2l_i+1=3`$ and $`\kappa _0`$ fixed from the binding energy, the remaining parameters $`\kappa _1`$ and $`\eta _0^{}`$ are fitted to the behavior of the “exact” model function (9), see Table 1. The resulting E1 strength function is showed as a dashed line in Fig. 1. This analytical approximation agrees, to a very high precision, with the numerical results obtained keeping the “exact” model function. The error is less than 2 % in the region of interest.

However, concerning the structure of the $`{}_{}{}^{B}8`$ ground state, one should keep in mind that the $`{}_{}{}^{B}7`$ core is in itself a weakly bound system with an excited $`1/2^{}`$ state at 429 keV. The common treatment of $`{}_{}{}^{B}8`$ as a two–body system is therefore highly questionable. We now want to investigate what effect the many–body structure of $`{}_{}{}^{B}8`$ might have on the strength function. We utilize a recent three–body calculation , where it was shown that, after projection on the two–body channel, there are three main components (adding up to 94 % of the total WF), and that the rest are pure three–body channels, see Table 1. For each of these two–body components we fit our parameters $`\kappa _1`$ and $`\eta _0^{}`$. The binding energy, $`E_0=137`$ keV, determines $`\kappa _0`$ for the two first components and $`E_0=566`$ keV for the third, $`{}_{}{}^{B}7`$ excited state, component. The best fit of the small $`r`$ behavior is obtained with $`p=5`$ which reflects the effectively larger centrifugal barrier in the three–body case. This centrifugal barrier will push the WF away from $`r=0`$ and will, thus, force it to become more narrow than the two–body WF. We therefore expect the distribution in momentum/energy space to be broader. This effect is clearly seen in Fig. 1 where the E1 strength function obtained using this three–body model111Note that this is not strictly a three–body model, but rather the two–body projection of a three–body WF. However, in the following we will consistently refer to it as three–body results. is shown as a solid line. This difference, seen in the strength function, should be even more pronounced in the energy spectrum. From these results one can conclude that the interpretation of energy spectra is highly model dependent.

In Fig. 2(a) we compare our analytical results for Coulomb dissociation, including both E1 and E2 transitions, to the experimental data from Davids et al. . This experiment is very appealing since the selection of scattering angles ($`\theta _{{}_{}{}^{B}8}1.77^{}`$) minimizes the contribution from nuclear scattering and the relatively high beam energy (82.7 MeV/$`A`$) justifies the use of first order perturbation theory. Concerning the shape of the energy spectrum we have an excellent agreement between the experimental data and our results obtained using the three–body model (see thin, dotted line), while the pure two–body calculation gives a too narrow peak. As to the absolute values, the three–body model gives about 20 % larger cross section than the experimental data. In Fig. 2(b) we compare the results from our two–body model (dashed line) to the results from a potential model calculation without final state interaction, from  (dash–dotted). Here, the intercluster distance for both models is $`r_{\mathrm{rms}}=4.23`$ fm which is lower than before ($`r_{\mathrm{rms}}=4.57`$ fm). The lesson from this figure is twofold. Firstly, it demonstrates the sensitivity to the intercluster distance and the importance of knowing this parameter. Secondly, we note that the difference between the two–body potential model calculation and our analytical results is not so large. The results deviate mainly in the peak amplitude.

In connection to the same experimental conditions we show in Fig. 3 the fraction of the cross section that can be attributed to E1 transitions. This information is of importance when studying the inverse radiative capture reaction and, in contrast to the conclusion drawn by Davids et al. , we claim that although the E2 contribution increases for low relative energies, it still does not dominate. Note that including M1 transitions will give raise to a sharp dip at 640 keV.

We have also compared our results for the radiative capture reaction to the experimental data by Hammache et al. . Our three–body result of $`\sigma _{\mathrm{rc}}(186\mathrm{keV})=16.1`$ nb agrees very well with the experimental value $`16.7\pm 2.1`$ nb. We also find that the E2 contribution to this cross section is around 0.03 %.

In conclusion, we have performed analytical studies of electromagnetic processes for loosely bound nuclei. This has been accomplished by using model radial functions that describe two–body WFs or the two–body projections of many–body WFs, accurately for all radii, and by only studying direct transitions to/from a pure Coulomb continuum. We have examined the difference between a pure two–body model and the results obtained incorporating many–body effects. From this we concluded that the interpretation of experimental data is highly model dependent. Comparisons have also been made to experimental results on $`{}_{}{}^{B}8`$ Coulomb dissociation and, the inverse, radiative capture reaction. We found that our three–body results coincide with the radiative capture cross section measured by Hammache et al.  while our Coulomb dissociation energy spectrum is about 20 % larger than the experimental data from Davids et al. . However, the shapes of the experimental and theoretical energy spectra are in excellent agreement. In contrast, our pure two–body model of $`{}_{}{}^{B}8`$, having the same intercluster distance $`r_{\mathrm{rms}}`$, does not agree with the experimental data. A final word of wisdom is that, in order to interpret these data correctly, it is very important to fix the spectroscopic factors of different two–body and many–body components. Therefore, we want to stress the usefulness of experiments where Coulomb dissociation is studied in complete kinematics. Examples of interesting channels in the $`{}_{}{}^{B}8`$ case is $`{}_{}{}^{B}8`$$``$$`{}_{}{}^{B}7`$($`1/2^{}`$)$`+p+\gamma `$ and $`{}_{}{}^{B}8`$$``$$`{}_{}{}^{H}3`$$`+\alpha +p`$. Recently, $`{}_{}{}^{B}7`$ fragments and $`\gamma `$–rays were measured in coincidence after breakup on a light target by Cortina–Gil et al.  and the excited core component of the WF was clearly observed.

N. B. S. is grateful for support from the Royal Swedish Academy of Science. The support from RFBR Grants No 00–15–96590, 02–02–16174 are also acknowledged.