text:
# Triaxial Modeling of Halo Density Profiles with High-resolution N-body Simulations.

## 1 Introduction

The density profiles of dark matter halos have attracted a lot of attention recently since Navarro, Frenk, & White (1996, 1997, NFW hereafter) discovered the unexpected scaling behavior in their simulated halos. Subsequent independent higher-resolution simulations (e.g., Fukushige & Makino, 1997, 2001; Moore et al., 1998; Jing, 2000; Jing & Suto, 2000) confirmed the validity of the NFW modeling, in particular the presence of the central cusp, although the inner slope of the cusp seems somewhat steeper than they originally claimed. Those previous models, however, have been based on the spherical average of the density profiles. Actually it is also surprising that the fairly accurate scaling relation applies after the spherical average despite the fact that the departure from the spherical symmetry is quite visible in almost all simulated halos (e.g., Fig.1 of Jing & Suto, 2000).

A more realistic modeling of dark matter halos beyond the spherical approximation is important in understanding various observed properties of galaxy clusters and non-linear clustering (especially the high-order clustering statistics) of dark matter in general. In particular, the non-sphericity of dark halos is supposed to play a central role in the X-ray morphologies of clusters (Jing, et al., 1995; Buote & Xu, 1997), in the cosmological parameter determination via the Sunyaev-Zel’dovich effect (Birkinshaw, Hughes & Arnaud, 1991; Inagaki, Suginohara & Suto, 1995; Yoshikawa, Itoh & Suto, 1998) and in the prediction of the cluster weak lensing and the gravitational arc statistics(Bartelmann et al., 1998; Meneghetti et al., 2000, 2001; Molikawa & Hattori, 2001; Oguri, Taruya, & Suto, 2001; Keeton & Madau, 2001). Nevertheless useful analytical modeling of the non-sphericity is almost impossible, and numerical simulations are the only practical means to provide statistical information.

While the non-sphericity of the dark matter halos is a poorly studied topic, some seminal studies do exist which attempt to detect and characterize the non-spherical signature (e.g., Barnes & Efstathiou, 1987; Warren et al., 1992; Dubinski, 1994; Jing, et al., 1995; Thomas et al., 1998; Yoshida et al., 2000; Meneghetti et al., 2001; Bullock, 2001). Nevertheless there is no systematic and statistical study to model and characterize the density profiles of simulated halos. This is exactly what we will present in the rest of the paper. In particular, much higher mass and spatial resolutions of our current N-body simulations enable us to characterize the statistics of the halo non-sphericity with an unprecedented precision.

This paper is organized as follows; two different sets of N-body simulations that we extensively analyze here are described in §2. In §3, we discuss how to define the iso-density surfaces of dark mater halos from simulation data and then argue that they are well approximated by a sequence of the concentric triaxial model. Section 4 characterizes the statistical distribution of the triaxial model parameters. Finally §5 is devoted to summary and discussion.

## 2 Simulations for dark matter halos

We use two different simulations for the current purpose. The first is our new set of cosmological N-body simulations with $`N=512^3`$ particles in a $`100h^1`$Mpc box, and the other is a set of high-resolution halo simulation runs. We describe the two simulations in the next subsections in order.

### 2.1 Cosmological simulations

The first set of simulations is our new runs with $`N=512^3`$ particles in a $`100h^1`$Mpc box. These runs have been carried out in 2001 with our Particle-Particle-Particle-Mesh ($`\mathrm{P}^3\mathrm{M}`$) code on the vector-parallel machine VPP5000 at the National Astronomical Observatory of Japan. The code adopts the standard $`\mathrm{P}^3\mathrm{M}`$ algorithm (Hockney & Eastwood, 1981; Efstathiou et al., 1985), is vectorized (Jing & Suto, 1998), and has been recently parallelized. A mesh of $`1200^3`$ grid points is used for the Particle-Mesh (PM) force computation with the optimized Green function (Hockney & Eastwood, 1981). The short-range force is compensated for the PM force calculation at the separation less than $`ϵ=2.7H`$, where $`H`$ is the mesh cell size (Efstathiou et al., 1985). The linked-list technique has been used for computing the short-range Particle-Particle (PP) interaction with $`448^3`$ linked-list cells. The computer has a total of 64 processors, and we use $`N_{\mathrm{CPU}}=8`$ to 32 processors, upon their availability, to run our code. The most important advantage of the machine for our work is that each processor has a big memory of $`16\mathrm{G}\mathrm{B}`$, sufficient for storing all the information of the code. The PM computation can be easily parallelized, and it is crucial to parallelize the PP computation which dominates the CPU computation time for a strongly clustered simulation like our present case. We sliced the simulation box in one direction (e.g. z-axis) with the thickness chosen to be the cell size of the linked-list cell. Those 448 slices in total are sorted in the descending order according to the number of particles they contain. We distribute the PP force computation among the different processors in a simple way; $`n`$-th processor ($`n=1`$ to $`N_{\mathrm{CPU}}`$) is assigned the force computation for those slices with indices of $`jN_{\mathrm{CPU}}+n`$ where $`j`$ runs from 0 to $`j_{\mathrm{max}}448/N_{\mathrm{CPU}}1`$. The PP interaction of the particles in the same slice and in the adjacent lower slice is considered, so the interaction for each pair of particles is computed only once. With this computation partition, we find that the load-balance problem, which becomes progressively serious for $`\mathrm{P}^3\mathrm{M}`$ simulations in the later strongly clustered regime, can be overcome to a satisfactory degree; even at the final stage of our simulation runs, the CPU time for the PP part is nearly inversely proportional to the number of the processors used. This implies that the code has achieved a good parallelization efficiency.

We consider two representative cold dark matter (CDM) models; a low-density flat cosmological model (LCDM) with $`\mathrm{\Omega }_0=0.3`$ and $`\lambda _0=0.7`$, and the Einstein-de Sitter model with $`\mathrm{\Omega }_0=1`$ (SCDM). The primordial density fluctuation is assumed to obey the Gaussian statistics, and the power spectrum is given by the Harrison-Zel’dovich type. The linear transfer function for the dark matter power spectrum is taken from Bardeen et al. (1986). The shape and the normalization of the linear power spectrum are specified by the shape parameter, $`\mathrm{\Gamma }=\mathrm{\Omega }_0h`$, and $`\sigma _8`$ respectively, where $`h`$ is the Hubble constant in $`100\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`\sigma _8`$ is the rms linear density fluctuation within the sphere of the radius $`8h^1\mathrm{Mpc}`$. Table 1 summarizes the physical and simulation parameters used for these simulations. We adopted $`\sigma _8=0.9`$ for LCDM and 0.55 for SCDM, both of which are slightly smaller than those in our previous simulations (Jing & Suto, 1998), but seem more consistent with recent observations (e.g., Seljak 2002; Lahav et al. 2002). With the adopted values for those physical parameters, the LCDM model satisfies almost all current observations while the SCDM model is known to have many difficulties. Therefore we mainly analyze the LCDM model for our purpose, and sometimes use the SCDM simulation just for comparison.

The boxsize of our cosmological simulations is $`100h^1\mathrm{Mpc}`$, so the particle mass is $`m_p=6.2\times 10^8h^1M_{}`$ and $`2.1\times 10^9h^1M_{}`$, respectively, for the LCDM and SCDM simulations (Table 1). The force resolution is $`\eta =20h^1\mathrm{kpc}`$ for the linear density softening form (Efstathiou et al. 1985; this roughly corresponds to $`\eta /3`$ for the Plummer-type softening length). The simulations are evolved by 1200 time steps from the initial redshift $`z_i=72`$. Two realizations are computed for each model. One additional LCDM simulation (LCDMa) uses a smaller force softening $`\eta =10h^1\mathrm{kpc}`$ and is evolved with 5000 time steps in order to check the possible effect of the force softening on the final dark matter distribution especially at small scales. As far as the shape of the virialized halos is concerned, we made sure that both simulations (LCDM and LCDMa) yield almost identical results. In what follows, therefore, we do not distinguish LCDM and LCDMa, and simply refer to them as LCDM.

### 2.2 Identification of dark halos in the cosmological simulations

The Friends-Of-Friends (FOF) method is a widely used algorithm to identify dark matter clumps in N-body data. The mean overdensity within the clumps is approximately proportional to $`b^3`$, where $`b`$ is the bonding length. It has been shown that the FOF clumps with $`b=0.2\overline{d}`$, where $`\overline{d}L/N^{1/3}`$ is the mean separation of particles, approximately correspond to the virialized dark matter halos of the mean overdensity 180 (e.g., Davis et al., 1985; Lacey & Cole, 1994). On the other hand, a large fraction of the FOF clumps identified with $`b=0.2\overline{d}`$ are known to form a system of multiple virialized halos that are bridged via thin filaments (e.g., Suto, Cen & Ostriker 1992; Suginohara & Suto 1992; Jing & Fang 1994, hereafter JF94). JF94 proposed to compute the overdensity around the local potential minima within each FOF clump to separate the virialized halos. While this can effectively achieve the goal, it is time-consuming to find the local potential minima (because there may be multiple minima within a single FOF clump).

Here we propose to use an alternative method which works faster. The thin bridges connecting the halos identified with $`b=0.2\overline{d}`$ can be effectively eliminated by reducing $`b`$. By trial and test, we found that the thin bridges almost disappear if we adopt $`b=0.1\overline{d}`$. With this recipe, however, the resulting FOF clumps have a smaller size and a higher overdensity than those defined according to the spherical collapse model. Therefore our scheme should be interpreted to identify first the central parts or the substructures of the entire halo. Next, for each FOF clump of $`b=0.1\overline{d}`$, we compute the gravitational potential of every member particle. The position of the particle of the minimum potential is defined as the center of the hosting halo. Then the spherical overdensity is computed around the halo center with increasing the radius, and the virial radius $`r_{\mathrm{vir}}`$ is found when the overdensity reaches the value predicted in the spherical collapse model. Here we use the fitting formula of Bryan & Norman (1998) for the spatially-flat ($`\mathrm{\Omega }(z)+\lambda (z)=1`$) models:

$$\mathrm{\Delta }_{\mathrm{vir}}(z)\frac{3M_{\mathrm{vir}}}{4\pi r_{\mathrm{vir}}^3\rho _{crit}}=18\pi ^2+82\left[\mathrm{\Omega }(z)1\right]39\left[\mathrm{\Omega }(z)1\right]^2,$$
(1)
where $`\rho _{crit}`$ is the critical density of the Universe. Since our choice $`b=0.1\overline{d}`$ preferentially selects smaller clumps than those predicted in the spherical model, some fraction of such clumps turn out to be substructures within the virial radius of a larger halo defined in the above equation. If the virial spheres of more than one halos overlap, we simply retain the most massive clump and throw away the others from the final halo list.

### 2.3 High-resolution halo simulations

Our cosmological simulations which we described above have a sufficient spatial resolution to discuss the statistics concerning the halo shapes and the concentration of the density profile (§4 and 5) as was conducted by Jing (2000) in the framework of the spherical approximation. Actually except for a delicate problem of determining the slope of the central cusp at $`r0.01r_{vir}`$, a larger simulation volume is more important than the higher resolution for the current purpose. Nevertheless we also use our higher-resolution halo simulations (Jing & Suto 2000; hereafter simply referred to halo simulations) to demonstrate that our triaxial modeling indeed provides a better description for halo profiles than the conventional spherical modeling (§3).

These halos are simulated with about a million particles within their virial radii (see Table 1 of Jing & Suto 2000). For mass scales of clusters, groups, and galaxies, there are four halos, respectively, and thus twelve halos in total. They are simulated in the LCDM model except the fact that the fluctuation amplitude, $`\sigma _8=1`$ (Kitayama & Suto 1997), is a bit larger than our current choice $`\sigma _8=0.9`$. Another advantage of the halo simulations is that those halos are simulated with almost the equal number of particles independently of the mass of the halos, and thus the resolution relative to the virial radius and the halo mass is kept constant. This is not the case for the cosmological simulations in which massive halos would have a better resolution in terms of the number of particles involved. Thus the possible artificial effect due to the variable resolution is suppressed in the halo simulations.

After Jing & Suto (2000) was published, we have completed runs of additional two halos with a galactic mass and with a group mass. Those new halos are referred to GX5 and GR5, respectively, according to our previous convention. While we add these two, we also eliminate two previous halos from the list of halos that we examine below; GR2 which shows a clear bi-modal structure, and GX1 which is seriously disrupted at $`z0.5`$ due to the tidal force of a nearby massive object. This is because the major purpose of analyzing the halo simulation catalogs is to check the validity of the triaxial modeling for typical halos. The fraction of those atypical halos is properly taken into account in the statistics drawn from the cosmological simulations. Thus the above replacement does not bias our conclusion.

## 3 Modeling the non-spherical density profiles of dark matter halos

In this section, we propose that a non-sphericity in the density profiles of dark halos is well described by a triaxial model on the basis of the detailed analysis of the halo simulations. In fact, we demonstrate that the triaxial modeling significantly improves the fit to the simulated profiles, at least for relatively relaxed halos, compared to the conventional spherical model. The statistical description including the probability distribution functions for axis ratios and the concentration parameters will be discussed in the next section using the cosmological simulations.

### 3.1 Defining the iso-density surfaces inside individual halos

The shapes of dark halos have been previously studied by many authors (e.g., Barnes & Efstathiou, 1987; Warren et al., 1992; Jing, et al., 1995; Thomas et al., 1998), and it is already well known that they exhibit a significant amount of departure from spherical distribution. Those previous studies first compute the inertial tensor for each halo, and then compute the distribution of the axial ratios and the correlation of the direction of the principal axes. While this is a well-defined method to characterize the shape of halos in principle, we do not employ this for two reasons.

First, this method assumes that we know in advance which particles belong to each halo. In reality this is not the case since we usually attempt to determine the member particles of a halo and its shape simultaneously. This is serious because the inertia tensor is sensitive to the outer boundary of the halo where the membership of particles is also difficult to define. Previous studies get around the problem by applying the procedure iteratively; first, all particles within a certain spherical radius from the center of halo are included to compute the inertial tensor and the resulting ellipsoidal configuration. Next, those particles outside the ellipsoid are thrown away from the member particles of the halo, and the inertia tensor is re-calculated. This procedure is repeated until the solution converges. While this method seems to work well in previous low-resolution N-body data, we were not able to obtain a good convergence in the case of our high-resolution halos. This is ascribed to the fact that our high-resolution halos retain a significant amount of substructures which have been artificially erased due to the overmerging effect in previous lower-resolution simulations. The iteration procedure is not stable in the presence of significant substructures especially at the boundary region of halos, since the inertial tensor is quite sensitive to them.

Second, our main interest here is not simply to define the overall shape of halos, but to characterize the density profile. Therefore we would like to have a sequence of iso-density surfaces with different overdensities. The ellipsoidal surface obtained from the above procedure, even if it converges, is not related to those iso-density surfaces, and thus not so useful after all for our purpose here.

With the above problems of the previous method in mind, we propose another approach to find the iso-density surfaces. This begins with the computation of a local density at each particle’s position. We adopt the smoothing kernel widely employed in the Smoothed Particle Hydrodynamics (SPH) method (e.g., Hernquist & Katz, 1989):

$$W(r,h_i)=\frac{1}{\pi h_i^3}=\{\begin{array}{cc}1\frac{3}{2}\left(\frac{r}{h_i}\right)^2+\frac{3}{4}\left(\frac{r}{h_i}\right)^3\hfill & (rh_i)\hfill \\ \frac{1}{4}\left(2\frac{r}{h_i}\right)^3\hfill & (h_i<r<2h_i),\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}$$
(2)
where $`h_i`$ is the smoothing length for the $`i`$-th particle. We use 32 nearest neighbor particles to compute the local density $`\rho _i`$, and $`h_i`$ is set to be a half of the radius of the sphere that contains those 32 neighbors. Using $`\rho _i`$, we construct the iso-density surfaces corresponding to the 5 different thresholds:

$`\rho _\mathrm{s}^{(n)}`$ $`=`$ $`A^{(n)}\rho _{\mathrm{crit}},`$ (3)
$`A^{(n)}`$ $`=`$ $`100\times 5^{n1}(n=15).`$ (4)
In practice, we collect all particles satisfying $`0.97\rho _\mathrm{s}^{(n)}<\rho _i<1.03\rho _\mathrm{s}^{(n)}`$ to define the $`n`$-th isodensity surface. The typical sizes (the mean radii) of those surfaces are 0.6, 0.4, 0.25 0.12 and 0.06 times the virial radius of the halo, respectively. Note that $`\rho _\mathrm{s}^{(n)}`$ is the local density, and thus the mean density of the halo inside the corresponding radius of $`\rho _\mathrm{s}^{(n)}`$ is generally much higher.

Actually a straightforward application of equation (3) results in many small distinct regions with the identical density threshold inside an individual halo. This is again due to the presence of the strong substructures in the halo. Since we are interested in the isodensity surfaces which represent the overall density profile of the parent halo, we have to eliminate those small regions corresponding to the substructures. For this purpose, we again use the FOF technique but with a different bonding length from that we used when identifying the virialized halos. After some trial and error, we find that an adaptive (i.e., dependent on each isodensity value) bonding length of $`b_n=3(\rho _\mathrm{s}^{(n)}/m_p)^{1/3}`$ works well (c.f., Suto, Cen & Ostriker, 1992).

### 3.2 Triaxial model fits to the iso-density surfaces

Figure 1 plots typical examples of the projected particle distributions within the isodensity surfaces for four different halos (CL3, GR1, GR5 and GX3) after particles in strong substructures are eliminated as described above. Those plots clearly suggest that the isodensity surfaces are typically approximated as triaxial ellipsoids. So we performed the following triaxial fit to the iso-density surfaces with 5 different thresholds separately:

$`R^2(\rho _\mathrm{s})={\displaystyle \frac{X^2}{a^2(\rho _\mathrm{s})}}+{\displaystyle \frac{Y^2}{b^2(\rho _\mathrm{s})}}+{\displaystyle \frac{Z^2}{c^2(\rho _\mathrm{s})}}.`$ (5)
The origin of the coordinates is always set at the center of mass of each surface, and the principal vectors $`𝐚`$, $`𝐛`$ and $`𝐜`$ ($`abc`$) are computed by diagonalizing the inertial tensor of particles in the surface (Fig.2). The projected views of the corresponding fitted ellipsoids are shown at the bottom panels in Figure 1, which implies that the ellipsoid fitting is a good approximation (at least visually).

Figure 3 plots the dependence of the axis ratios, $`a/c`$ and $`b/c`$, on the isodensity threshold $`\rho _\mathrm{s}`$. Naturally each halo exhibits different behavior which may reflect the different merging history and/or tidal force field. Nevertheless, several systematic dependences are quite visible. The halos of cluster mass generally have smaller axial ratios than those of galactic mass, implying that the halos of the galactic mass are rounder on average than those of cluster mass. This mass-dependence will be quantified with a large sample of halos from the cosmological simulations in §4.

On the other hand, we also note that the axial ratios decrease with increasing the density; the iso-density surfaces become more elongated in the central region than in the outer region. The mean (with the error bar of the mean) of the axial ratios computed from the twelve halos are plotted in the right panels of Figure 3 (the symbols). The solid lines show the single power-law fit for the mean axis ratios:

$`{\displaystyle \frac{a}{c}}`$ $`=`$ $`0.56\left({\displaystyle \frac{\rho _\mathrm{s}/\rho _{\mathrm{crit}}}{2500}}\right)^{0.052}`$ (6)
$`{\displaystyle \frac{b}{c}}`$ $`=`$ $`0.71\left({\displaystyle \frac{\rho _\mathrm{s}/\rho _{\mathrm{crit}}}{2500}}\right)^{0.040}.`$ (7)
Figure 4 shows the degree of the alignment of the axis directions among isodensity surfaces at different densities (radii). We define $`\theta _{11}`$ as the angle between the major axis of the isodensity surfaces and that of the $`A^{(3)}=2500`$ isodensity surface as shown in Figure 2. Similarly, $`\theta _{22}`$ is defined with respect to their middle axes. According to our definition, $`\mathrm{cos}\theta _{11}=\mathrm{cos}\theta _{22}=1`$ at $`\rho _s/\rho _{\mathrm{crit}}=A^{(3)}=2500`$.

We find that the major axes align pretty well within a halo; for about $`70\%`$ of the halos $`\mathrm{cos}\theta _{11}`$ at different radii is larger than 0.7. For about half of the sample, $`\mathrm{cos}\theta _{11}`$ is larger than 0.9. In a few cases (3 out of 12 halos), however, the alignment of the major axes is poor. When we check these halos individually (e.g. GR1), it turns out that $`b/c`$ for the two halos is quite close to unity, indicating they are oblate halos with $`bc`$ and thus the direction of the major axis is difficult to measure (if $`b=c`$, the direction of the major axis is arbitrary within a plane). Thus the apparent mis-alignment of their major axes is not meaningful. Only for the remaining one halo (GX3; Figure 1), the major axes of the outer and the innermost isodensity surfaces are indeed perpendicular to that at the middle. This is the real case that the major axes are significantly mis-aligned.

The alignments of the middle axes show similar behavior: for most of the halos the degree of the alignment is satisfactory. For those which show significant misalignment of the middle axes, their $`a/b`$ or $`b/c`$ ratio is usually quite close to unity and the direction of the middle axes (and the minor or major axes) can be poorly determined at best. Only in the case like GX3, no simple ellipsoid description can be found, but this is fairly exceptional. The alignment seems slightly better for cluster-sized halos, but this would be simply because galactic halos are more spherical and thus the direction of the major axis is less accurate than that for cluster-sized halos.

### 3.3 Triaxial versus spherical modeling of dark halos

In the last subsection, we have seen that the isodensity ellipsoids at different radii are approximately aligned, and the axial ratios of the ellipsoids are nearly constant. These facts suggest the possibility that the internal density distribution within a halo can be approximated by a sequence of the concentric ellipsoids of a constant axis ratio. To show this to be an improved description over the conventional spherical description, we compute the quadrupole of the particle distribution within a spherical shell ($`Q_s`$) or an ellipsoid shell ($`Q_e`$). For a spherical shell, the positions of particles inside the shell can be described by

$$\{\begin{array}{c}x=r\mathrm{sin}\theta \mathrm{cos}\varphi \hfill \\ y=r\mathrm{sin}\theta \mathrm{sin}\varphi \hfill \\ z=r\mathrm{cos}\theta \hfill \end{array}$$
(8)
with $`r`$ being the (conventional) spherical radius. Similarly, the positions of the particles in an ellipsoidal shell can be described by

$$\{\begin{array}{c}X=R\left(\frac{a}{c}\right)\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\mathrm{\Phi }\hfill \\ Y=R\left(\frac{b}{c}\right)\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\mathrm{\Phi }\hfill \\ Z=R\mathrm{cos}\mathrm{\Theta }\hfill \end{array},$$
(9)
where $`X`$, $`Y`$ and $`Z`$ axes are the principal vectors of the ellipsoidal shell, and $`a/c`$ and $`b/c`$ are the axis ratios. In the rest of the paper, we preferentially use the capital $`R`$ to refer to the length of the major axis defined in the triaxial model.

Then the quadrupole moments of the iso-density surfaces in the spherical and triaxial models, $`Q_s`$ and $`Q_e`$, are computed as

$`Q_s`$ $``$ $`{\displaystyle \frac{1}{5N_p}}{\displaystyle \underset{m=2}{\overset{+2}{}}}\left|{\displaystyle \underset{j}{}}Y_{2m}(\theta _j,\varphi _j)\right|^21,`$ (10)
$`Q_e`$ $``$ $`{\displaystyle \frac{1}{5N_p}}{\displaystyle \underset{m=2}{\overset{+2}{}}}\left|{\displaystyle \underset{j}{}}Y_{2m}(\mathrm{\Theta }_j,\mathrm{\Phi }_j)\right|^21,`$ (11)
where the summation over $`j`$ runs for all particles ($`N_p`$) in the iso-density surface, and $`Y_{lm}`$ is the spherical harmonics. If the spherical (triaxial) model is exact, $`Q_s`$ $`(Q_e)`$ vanishes. Using these measures, we will show the extent to which the triaxial model indeed provides a significantly improved description for the simulated halos.

In practice, we compute $`Q_s(r)`$ and $`Q_e(R)`$ for 5 shells of each halo at $`r=R=0.65r_{\mathrm{vir}}`$, $`0.35r_{\mathrm{vir}}`$, $`0.2r_{\mathrm{vir}}`$,$`0.12r_{\mathrm{vir}}`$, and $`0.065r_{\mathrm{vir}}`$ with the shell thickness $`\mathrm{\Delta }r/r=\mathrm{\Delta }R/R=\mathrm{ln}10\times 0.1=0.23`$. Those shells are centered at the potential minimum of the halo. In the triaxial model, we assume that the shells have the same axis ratios and the same principal axis directions as measured from the isodensity surface at $`A^{(3)}=2500`$. Thus those shells do not necessarily correspond to the iso-density surfaces that we have discussed. Actually this treatment is important because otherwise the triaxial model (with more degrees of freedom) should always provide a better fit. Also this approximation is most likely what one would like to apply statistically to halos of visible objects, which would yield a practical and fair comparison between the spherical and triaxial models.

In the top and the middle panels of Figure 5, we present the quadruple moments $`Q_s`$ and $`Q_e`$ for the twelve halos. The quadruple $`Q_s`$ for the spherical modeling increases nearly monotonically with the radius. The $`Q_e`$ for our ellipsoidal modeling (in its simplified version as described above) stays flat at $`R<0.3r_{\mathrm{vir}}`$ but increases with the radius at the larger radius. The ratio of the two quadruples is shown in the bottom panel of Figure 5 which indicates that our triaxial model, even simplified, fits the simulated halo profiles much better than the spherical model. For 10 out of 12 halos, the ratio, $`Q_e(R)/Q_s(r)`$, is much smaller than 1 at all scales ($`r=R`$). Even for the remaining two halos (GX3 and GR5), the ratio exceeds unity a bit only at the largest radius, and the triaxial description shows a significant improvement over the spherical model. The ratio $`Q_e(R)/Q_s(r)`$ seems to approach unity as $`r`$ becomes closer to $`r_{\mathrm{vir}}`$. The reason might be that the subclustering is more prominent in the outskirt region than in the central region, since the subhalos are tidally stripped when they fall into the central region (see Fig. 1 of Jing & Suto (2000)). The strong subclustering makes both the spherical model and the ellipsoidal model difficult to accurately describe the complicated density distribution at $`r_{vir}`$. But the figure also clearly shows that our triaxial model works significantly better than the spherical model for $`r<0.7r_{vir}r_{200}`$, i.e. almost the entire halo (the definition for $`r_{200}`$ will be given shortly).

### 3.4 Density profiles in the triaxial model

The next important task is to describe the density profiles in the triaxial model generalizing the previous results in the spherical approximation (NFW; Moore et al., 1998; Jing & Suto, 2000; Klypin, et al., 2001). In the same spirit of the previous subsection, we do not perform the fit to the iso-density surfaces that we identified, but rather compute the mean density $`\rho (R)`$ at the simplified triaxial shells (i.e., the same axis ratios and the axis directions for the entire halo as those measured from its isodensity surface at $`A^{(3)}=2500`$) within a thickness of $`\mathrm{\Delta }R/R=0.12`$.

Figure 6 plots the density profiles measured in this way for individual halos as a function of $`R`$. As in the spherical case, we adopt the following form:

$$\frac{\rho (R)}{\rho _{\mathrm{crit}}}=\frac{\delta _\mathrm{c}}{\left(R/R_0\right)^\alpha \left(1+R/R_0\right)^{3\alpha }},$$
(12)
where $`R_0`$ is a scale radius and $`\delta _\mathrm{c}`$ is a characteristic density. Again following the definition of $`r_{200}`$ in the spherical model (within which the mean matter density is $`200\rho _{\mathrm{crit}}`$), we define a radius $`R_e`$ so that the mean matter density within the ellipsoid of the major axis radius $`R_e`$ is $`\mathrm{\Delta }_e\rho _{\mathrm{crit}}`$ with

$$\mathrm{\Delta }_e=5\mathrm{\Delta }_{\mathrm{vir}}\left(\frac{c^2}{ab}\right)^{0.75}.$$
(13)
The non-trivial dependence of $`\mathrm{\Delta }_e`$ on the axis ratios in the above equation is chosen so that $`R_e`$ becomes a fixed fraction of the virial radius $`r_{\mathrm{vir}}`$ (see Fig.7 below).

The best-fits to equation (12) for each halo are shown in Figure 6 for $`\alpha =1.5`$ (solid lines) and for $`\alpha =1.0`$ (dotted lines). Up to the resolution limit of the halo simulations ($`R/R_e0.02`$), equation (12) yields a good fit both for $`\alpha =1`$ and for $`\alpha =1.5`$. If comparing the fits to the simulation data more carefully, however, $`\alpha =1`$ works better for the halos of cluster mass and $`\alpha =1.5`$ better for those of galactic mass, which is consistent with the finding of Jing & Suto (2000) in the spherical model (but see Fukushige & Makino, 2001, for a different point of view).

We also introduce a concentration parameter in our triaxial model:

$$c_\mathrm{e}\frac{R_\mathrm{e}}{R_0},$$
(14)
which is plotted in the upper panel of Figure 7 adopting $`\alpha =1.0`$ (crosses) and $`\alpha =1.5`$ (filled circles) in the fit. In what follows we will not address the issue related to the inner slope of the density profiles, and adopt $`\alpha =1`$. It should be noted, however, that our statistical results presented in the next section can be readily applied to the $`\alpha =1.5`$ case since the ratio $`c_e(\alpha =1.5)/c_e(\alpha =1)`$ is always close to 1/2.

Before moving to the statistical analysis of halos in the cosmological simulations, we note that the value of $`R_e`$ and thus that of $`c_e`$ are dependent on our specific definition of $`\mathrm{\Delta }_e`$ (eq.). As the middle and bottom panels in Figure 7 indicate, both $`R_e/r_{\mathrm{vir}}`$ and $`c_e/c_{\mathrm{vir}}`$ (where $`c_{\mathrm{vir}}`$ is the ratio of the virial halo radius to the scale radius $`r_s`$ in the spherical model) remain constant ($`0.45`$) independently of the mass of the halos when we adopt equation (13) for $`\mathrm{\Delta }_e`$. This property is quite useful in applying our results for a variety of theoretical predictions, since for a halo of given virial mass $`M_{vir}`$:

$$M_{vir}=\frac{4\pi }{3}r_{vir}^3\mathrm{\Delta }_{\mathrm{vir}}\rho _{\mathrm{crit}},$$
(15)
the radius $`R_e`$ in our triaxial model is easily computed. It is also known that the $`c_{\mathrm{vir}}`$ is a function of the halo mass (NFW; Eke, Navarro & Steinmetz 2001) with the scatter described by the lognormal distribution function (Jing 2000; Bullock et al. 2001). Therefore, once the shape of a halo at a given mass is specified, the density profile of the halo is completely fixed. The statistical distribution function of the halo shape is discussed in the next section.

## 4 Statistics of triaxial density profiles

High-resolution halo simulations, like those used in the last section, is well suited for studying the detailed internal structures of individual halos, but the number of such halos is too small for a statistical description. Therefore we switch to the halo catalogs constructed from our cosmological simulations in order to study the probability distribution of the shape of halos. As emphasized in §2, the cosmological simulations employ $`N=512^3`$ particles in a 100$`h^1`$Mpc box and thus the mass resolution is even better than that of individual halo simulations of the original NFW paper, for instance.

We consider halos which contain more than $`10^4`$ particles within the virial radius. The lower mass limits are $`6.2\times 10^{12}M_{}`$ and $`2\times 10^{13}M_{}`$ in the LCDM and SCDM models, respectively. We also consider three epochs at redshifts $`z=0`$, 0.5 and 1.0 to examine the time-dependence. At these redshifts, we have 2494, 2160, and 1534 halos in the LCDM model, and 1806, 879, and 263 halos in the SCDM model, respectively.

### 4.1 Probability distribution of axis ratios

Following the prescription presented in the last section, we determine the halo shapes at the iso-density surfaces with $`A^{(3)}=2500`$. Since the typical radius of the surfaces is about $`0.3r_{\mathrm{vir}}`$, they are well resolved in our cosmological simulations; the force softening length is typically smaller by one order of magnitude.

Left panels of Figures 8 and 9 present the ratio $`a/c`$ of the minor axis $`a`$ to the major axis $`c`$ for halos from the cosmological simulations in the LCDM and SCDM models, respectively; solid, dotted and dashed histograms indicate the results for $`10^4N_{\mathrm{halo}}<2\times 10^4`$, $`2\times 10^4N_{\mathrm{halo}}<6\times 10^4`$, $`6\times 10^4N_{\mathrm{halo}}`$, where $`N_{\mathrm{halo}}`$ is the number of particles within the virial radius of each halo; in those Figures we use $`M_4N_{\mathrm{halo}}/10^4`$, and thus $`M_4=1`$ corresponds to $`M_{\mathrm{vir}}=6.2\times 10^{12}h^1M_{}`$, and $`2.1\times 10^{13}h^1M_{}`$ for our LCDM and SCDM models). Top, middle, and bottom panels show the results at $`z=0`$, 0.5, and 1.0.

Two systematic trends are visible; the ratio is slightly larger for less massive halos, and decreases at higher redshifts. This motivates us to attempt the following empirical scaling for the axis ratio $`a/c`$:

$$\stackrel{~}{r}_{\mathrm{ac}}\left(\frac{a}{c}\right)_{\mathrm{sc}}=\left(\frac{a}{c}\right)\left(\frac{M_{\mathrm{vir}}}{M_{}}\right)^{0.07[\mathrm{\Omega }(z)]^{0.7}},$$
(16)
where $`M_{}`$ is the characteristic non-linear mass at $`z`$ so that the rms top-hat smoothed overdensity at the scale $`\sigma (M_{},z)`$ is $`\delta _\mathrm{c}=1.68`$. The $`M_{}`$ at $`z=0`$, 0.5, and 1.0 are $`9.4\times 10^{12}h^1M_{}`$, $`2.0\times 10^{12}h^1M_{}`$, and $`3.8\times 10^{11}h^1M_{}`$ respectively for LCDM, and $`8.5\times 10^{12}h^1M_{}`$, $`9.7\times 10^{11}h^1M_{}`$, and $`1.4\times 10^{11}h^1M_{}`$ respectively for SCDM.

Such scaled axis ratios $`\stackrel{~}{r}_{\mathrm{ac}}`$ show a fairly universal distribution almost independently of the mass and the epoch (histograms in the the right panels of Figs. 8 and 9). The universal probability distribution function of the ratio $`\stackrel{~}{r}_{\mathrm{ac}}`$ is well fitted to the following Gaussian:

$$p(\stackrel{~}{r}_{\mathrm{ac}})d\stackrel{~}{r}_{\mathrm{ac}}=\frac{1}{\sqrt{2\pi }\sigma _s}\mathrm{exp}\left(\frac{(\stackrel{~}{r}_{\mathrm{ac}}0.54)^2}{2\sigma _s^2}\right)d\stackrel{~}{r}_{\mathrm{ac}}$$
(17)
with $`\sigma _s=0.113`$.

Next we decompose the joint probability distribution function of the axis ratios as

$`p(a/c,b/c)d(a/c)d(b/c)`$ $`=`$ $`p(a/c)d(a/c)p(b/c|a/c)d(b/c)`$ (18)
$`=`$ $`p(a/c)d(a/c)p(a/b|a/c)d(a/b)`$ (19)
in terms of the conditional probability distribution functions, $`p(b/c|a/c)`$ and $`p(a/b|a/c)`$. The second equality holds because once $`a/c`$ is fixed, the distribution of $`a/b`$ is uniquely determined from that of $`b/c`$. Since we have shown that the distribution function $`p(a/c)`$ is well approximated by equations (16) and (17), we compute the conditional probability distribution $`p(a/b|a/c)`$. Figures 10 and 11 plot the results for the LCDM and SCDM models, respectively. Different panels correspond to $`p(a/b|a/c)`$ for different ranges of $`a/c`$. Solid, dotted, and dashed histograms indicate $`p(a/b|a/c)`$ at $`z=0`$, 0.5, and 1.0, respectively.

The conditional functions appear to be insensitive to the redshift. In both cosmological models, they are accurately fitted to

$$p(a/b|a/c)=\frac{3}{2(1r_{\mathrm{min}})}\left[1\left(\frac{2a/b1r_{\mathrm{min}}}{1r_{\mathrm{min}}}\right)^2\right],$$
(20)
for $`a/br_{\mathrm{min}}`$, where $`r_{\mathrm{min}}=a/c`$ for $`a/c0.5`$ and $`r_{\mathrm{min}}=0.5`$ for $`a/c<0.5`$. $`p(a/b|a/c)=0`$ for $`a/br_{\mathrm{min}}`$.

### 4.2 Probability distribution of the concentration parameter

We apply the triaxial density profile (eq.) obtained in the halo simulations to the halo catalogs in the cosmological simulations. Considering the resolution limits, we adopt $`\alpha =1`$ and use the data points at $`\eta <R_e<r_{\mathrm{vir}}`$ in the fit, where $`\eta `$ the force softening length (see §2). Since we do not address the innermost structures of the halos and rather focus on the value of the concentration parameter $`c_e`$, this catalog has a sufficient resolution to yield an unbiased estimate (e.g., Jing 2000; Bullock et al. 2001; Eke, Navarro & Steinmetz 2001 for discussion). As already found in the spherical model (Jing 2000), the distribution of $`c_e`$ in the triaxial model has a significant scatter even if the range of the halo mass is fairly specified reflecting the dependence of the merging history of the individual halo.

The resulting probability distribution functions for $`c_e`$ are presented in Figure 12, which are well fitted by the lognormal distribution:

$$p(c_e)dc_e=\frac{1}{\sqrt{2\pi }\sigma _{c_e}}\mathrm{exp}\left[\frac{(\mathrm{ln}c_e\mathrm{ln}\overline{c}_e)^2}{2\sigma _{c_e}^2}\right]d\mathrm{ln}c_e$$
(21)
with a dispersion of $`\sigma _{c_e}0.3`$ both in the SCDM and LCDM models. The dispersion is slightly larger than the value estimated in the spherical model ($`0.2`$) for equilibrium halos, but is comparable to the value for all halos put together (Jing 2000). It should be noted here that despite the fact that the triaxial model is superior in describing the density distribution of halos, the scatter in $`c_e`$ is comparable to the scatter of concentrations in spherical profile fits, which probably means that the scatter originates from the halo merger histories rather than non-sphericity of the halos.

The probability distribution (eq.) is completed by specifying the mean of the concentration parameter $`\overline{c}_e`$. The result from our simulations is plotted in Figure 13 as a function of the halo mass at $`z=0`$, 0.5 and 1.0. NFW proposed a semi-analytic fitting formula for the concentration $`c_{\mathrm{vir}}`$ in the spherical model. <sup>1</sup><sup>1</sup>1Originally NFW defined the concentration parameter as $`c_{200}r_{200}/r_s`$, where $`r_{200}`$ is the radius within which the mean overdensity is $`200\rho _{\mathrm{crit}}`$. Their recipe, however, can be easily generalized to $`c_{\mathrm{vir}}`$, since $`r_{200}/r_v`$ is almost constant for a given cosmology. More recently Bullock et al. (2001) have shown that in their LCDM model (the parameters are similar to those of our LCDM model here) $`c_{\mathrm{vir}}`$ of a given mass decreases with $`z`$ as $`(1+z)^1`$. The redshift dependence is stronger than that predicted in the NFW recipe. Thus Bullock et al. (2001) have proposed another recipe which successfully describes the concentration $`c_{\mathrm{vir}}`$. Since we have already shown that the ratio, $`c_e/c_{\mathrm{vir}}`$, is almost constant (Fig.7), it is interesting to see if the formula of Bullock et al. (2001) also describes the behavior of $`c_e`$ in our triaxial model.

In the LCDM model, we find that the redshift dependence of $`c_e`$ for a given mass is approximately $`(1+z)^1`$ in good agreement with their result. In the SCDM model, however, our result of $`c_e`$ shows a stronger redshift dependence than their prediction. This is also likely that the fitting formula of Bullock et al. (2001) was designed to fit spherical concentrations, $`c_{\mathrm{vir}}`$, and is not applicable to the non-spherical case due to the evolution of $`c_e/c_{\mathrm{vir}}`$.

Following NFW and Bullock et al. (2001), we propose a new fitting formula for $`\overline{c}_e`$ in the triaxial model:

$$\overline{c}_e(M,z)=A_e\sqrt{\frac{\mathrm{\Omega }(z)}{\mathrm{\Omega }(z_c)}}\left(\frac{1+z_c}{1+z}\right)^{\frac{3}{2}}.$$
(22)
In the above, $`z_c`$ is the collapse redshift of the halo of mass $`M`$ (NFW):

$$\mathrm{erfc}\frac{\delta _c(z_c)\delta _c}{\sqrt{2[\sigma ^2(fM)\sigma ^2(M)]}}=\frac{1}{2},$$
(23)
where $`\sigma (M)`$ is the rms top-hat mass variance at $`z=0`$, $`\delta _c=1.68`$, $`\delta _c(z)=1.68/D(z)`$, $`D(z)`$ is the linear growth factor, and $`f=0.01`$. Solid lines in Figure 13 indicate the predictions of equation (22), implying that the formula describes our simulation results very accurately. In those plots, we adopt $`A_e=1.1`$ and 1.0 for the LCDM and SCDM models, respectively.

We also made sure that the formula also agrees well with our halo simulations in the LCDM model, while the results appear 10 to 20 % higher (i.e., $`A_e=1.21.3`$) than those of the cosmological simulations ($`A_e=1.1`$). Considering both the typical 30% scatter in $`c_e`$ and the limited number of the high-resolution halos (12 in total), the above level of difference may not be interpreted so seriously at this point. In fact, the difference may be attributed partly to the fact that halos with significant substructures (like GR2) have been eliminated in the high-resolution halo samples (§2), while we have not attempted such a selection in the cosmological simulations. Indeed Jing (2000) has noted that halos in equilibrium are systematically more centrally concentrated than those with significant substructures. We also note that most previous studies including NFW and Eke, Navarro & Steinmetz (2001) have preferentially selected isolated halos in re-simulating with higher resolution, which would have less substructures and therefore have slightly higher concentration than average. If one is interested in halos in nearly equilibrium, the best-fit value of $`A_e`$ should become 1.3. Since $`c_e/c_{vir}`$ remains constant(Figs.7 and 15), the fitting formula (eq.) can also be used for predicting $`c_{vir}`$ in CDM models.

Finally we have checked if the fitted values of $`R_e`$ and $`c_e`$ are dependent on the shapes of halos. Figure 14 presents the ratio of $`R_e`$ to the virial radius $`r_{\mathrm{vir}}`$ as a function of the axis ratio $`a/b`$. Clearly $`R_e/r_{\mathrm{vir}}`$ is independent of $`a/b`$ and of the redshift (or equivalently the halo mass in unit of $`M_{}`$, see also Fig.7), and approximately given by 0.45. Similarly, we find that $`R_e/r_{\mathrm{vir}}`$ is independent of $`b/c`$ and $`a/c`$. On the other hand, the concentration parameter $`c_e`$ is slightly dependent on the halo shape. Figure 15 indicates that halos with smaller $`a/c`$ are less centrally concentrated.

In terms of the scaled axis ratio $`(a/c)_{\mathrm{sc}}`$ (eq.), the ratio of the mean concentration $`c_e`$ for a given $`\stackrel{~}{r}_{ac}=(a/c)_{\mathrm{sc}}`$ and the overall average $`\overline{c}_e(M,z)`$ (eq.) is well approximated by

$$\frac{c_e[\stackrel{~}{r}_{ac},M,z]}{\overline{c}_e(M,z)}=1.35\mathrm{exp}\left[\left(\frac{0.3}{\stackrel{~}{r}_{ac}}\right)^2\right].$$
(24)
This fit is plotted in the solid line in Figure 15, which is in good agreement with the simulation data for different halo masses and both in the LCDM and SCDM models.

In this section, we used the halos identified from cosmological simulations which have $`N>10^4`$ particles. The smallest halos are resolved much more poorly than the largest halos and the high-resolution halos (§3) that consist of $`10^6`$ particles. In order to make sure that our results are robust to the mass resolution, we repeat the same analysis of the axis ratios and the density profile by randomly selecting $`N=10^4`$ particles from each of the twelve high-resolution halos of the previous section. The results are compared in Figure 16 with those obtained in the last section where we consider all particles. Both the axis ratios and the concentration of the randomly selected sample agree well with those of the original halo sample; the typical dispersion of the axis ratios between the two samples is $`10\%`$, and the concentration $`c_e`$ of the randomly selected sample is slightly lower ($`8\%`$). This comparison indicates that the mass resolution does not affect our results in this section significantly.

## 5 Summary and Discussion

This paper has presented a triaxial modeling of the dark matter halo density profiles extensively on the basis of the combined analysis of the high-resolution halo simulations (12 halos with $`N10^6`$ particles within their virial radius) and the large cosmological simulations (5 realizations with $`N=512^3`$ particles in a $`100h^1`$Mpc boxsize). In particular, we found that the universal density profile discovered by NFW in the spherical model can be also generalized to our triaxial model description. Our triaxial density profile is specified by the concentration parameter $`c_e`$ and the scaling radius $`R_0`$ (or the virial radius $`R_e`$ in the triaxial modeling) as well as the axis ratios $`a/c`$ and $`a/b`$.

We have obtained several fitting formulae for those parameters which are of practical importance in exploring the theoretical and observational consequences of our triaxial model (in doing so we have adopted $`\alpha =1`$ since the precise value of the inner slope is difficult to reliably determine even with the resolution of the current simulations);

* the mass and redshift dependence of the axis ratio, or equivalently the definition of the scaled axis ratio $`\stackrel{~}{r}_{ac}(a/c)_{\mathrm{sc}}`$ : equation(16)
* the probability distribution of the axis ratio $`p(\stackrel{~}{r}_{ac})`$ : equation(17)
* the conditional probability distribution of the axis ratios $`p(a/b|a/c)`$ : equation(20)
* the mean value of the concentration parameter $`\overline{c}_e(M,z)`$ : equation(22)
* the dependence of the concentration parameter on the axis ratio $`\stackrel{~}{r}_{ac}`$ : equation(24)
* the probability distribution of the concentration parameter $`p(c_e)`$ : equation(21)

Since $`c_e/c_{vir}`$ remains constant(Figs.7 and 15), the fitting formula (eq.) can also be used for predicting $`c_{vir}`$ in CDM models.

We have focused on the triaxial modeling and characterization of dark halos in the present paper, and plan to show specific applications elsewhere. Nevertheless it would be worthwhile to mention several important examples of the current model.

The results of the paper are applicable in a fairly direct manner to the following three areas. (i) the weak and strong lens statistics. The comparison with the weak lensing observations provides information of the degree of triaxiality of observed clusters, mainly at outer regions. In addition, the frequency of the lensing arc is known to be sensitive to the non-sphericity of the halo mass profile especially in the central regions (e.g., Bartelmann et al., 1998; Meneghetti et al., 2000, 2001; Molikawa & Hattori, 2001; Oguri, 2002). (ii) predictions of the non-linear clustering of dark matter based on the halo model (e.g., Mo, Jing, & Börner, 1997; Ma & Fry, 2000; Hamana et al., 2001; Kang et al., 2002). The high-order statistics of clustering, e.g. the three-point correlation and the bispectrum, should be quite sensitive to the non-sphericity. (iii) Dynamics of galactic satellites. Recently this has been argued to be very sensitive to the non-sphericity of the host halo (e.g., Ibata et al., 2001). The combination of those three approaches would even yield a direct test of the cold dark matter paradigm (Spergel & Steinhardt, 2000; Yoshida et al., 2000).

Of course the non-sphericity of dark matter halos is critical to understanding that of the gas density profile of clusters of galaxies. Since gas dynamics is characterized by the isotropic pressure tensor and does not directly follow the dark matter distribution in halos. In fact, most hydrodynamic simulations of galaxy clusters suggest that the gas distribution is generally rounder than that of dark matter. Nevertheless we would like to mention a couple of important examples where the non-sphericity in the gas distribution has crucial and observable consequences. (iv) the gas and temperature profiles of X-ray clusters. Almost all previous analytical models for the X-ray profiles of galaxy clusters have adopted the spherical approximation perhaps due to the lack of any specific model for the non-sphericity. Since our triaxial model specifies the gravitational potential of the hosting halos, one may compute the gas or temperature profiles, with an additional assumption of the hydrostatic equilibrium for instance, as performed in the NFW model (e.g., Makino, Sasaki & Suto, 1998; Suto, Sasaki & Makino, 1998; Komatsu & Seljak, 2001). If combined with the observed surface brightness distribution of clusters, one may in principle solve for the gas and temperature profiles simultaneously for a given non-spherical distribution of dark matter (Silk & White, 1978; Yoshikawa & Suto, 1999; Zaroubi et al., 1998). (v) the systematic bias and statistical distribution of the Hubble constant estimated via the Sunyaev – Zel’dovich effect. In view of the on-going observational projects, it is of vital importance to re-evaluate the reliability of the estimates taking account of the non-sphericity effect of the clusters. With the above modeling of the gas and temperature profiles for individual clusters, one may discuss the statistical properties of the estimates of the Hubble constant(e.g., Fox & Pen, 2002), combining the extensive fitting formula for the probability distribution functions of the triaxial model parameters and the halo mass function (e.g., Sheth & Tormen, 1999; Jenkins et al., 2001).

We thank an anonymous referee for helpful comments on the earlier manuscript. Numerical simulations presented in this paper were carried out at ADAC (the Astronomical Data Analysis Center) of the National Astronomical Observatory, Japan, and at KEK (High Energy Accelerator Research Organization, Japan). Y.P.J. was supported in part by the One-Hundred-Talent Program, by NKBRSF (G19990754) and by NSFC (No.10125314), and Y.S was supported in part by the Grant-in-Aid from Monbu-Kagakusho (07CE2002, 12640231), and by the Supercomputer Project (No.00-63) of KEK.


format: text

text:
# Renormalization group approach to the interacting boson fermion systems

## Abstract

We study a pseudogap region of the mixed boson fermion system using a recent formulation of the renormalization group technique through the set of infinitesimal unitary transformations. Renormalization of fermion energies gives rise to a depletion of the low energy states (pseudogap) for temperatures $`T^{}>T>T_c`$ which foreshadows appearance of the pairwise correlations with their long range phase coherence being missed. With a help of the flow equations for boson and fermion operators we analyze spectral weights and finite life times of these quasiparticles caused by interactions.

Many experimental measurements of e.g. the angle resolved photoemission (ARPES), tunneling spectroscopy (STM), specific heat as well as the optical and magnetic characteristics (see for an overview ) point at a presence of pseudogap in the normal phase spectrum of underdoped cuprate superconductors. It is believed that a proper description of this phenomenon would be crucial for understanding of the HTSC mechanism.

Model of itinerant fermions coupled via the charge exchange interaction to boson particles , for example of bipolaronic origin, has been shown to posses such a pseudogap in the normal phase fermion spectrum. Its occurrence has been obtained on a basis of the selfconsistent perturbation theory , diagrammatic expansion  and also by solving the dynamical mean field equations for this model within the non-crossing approximation (NCA) . Unfortunately, these techniques have not been able to study a superconducting phase of the model, so evolution of the pseudogap into the true superconducting gap could not be investigated. Only very recently we have reported some prior results for the boson fermion (BF) model using the flow equation method  where superconducting gap and normal phase pseudogap could be treated on equal footing. Both gaps originate from the same source (from the exchange interaction between fermions and bosons) but they are differently scaled against a magnitude of this interaction and the total concentration of charge carriers. Still the most intriguing question concerning the evolution from the pseudo- to the superconducting gap is poorly understood due to technical difficulties (it is a very tough problem to solve the flow equations for $`2+ϵ`$ dimensional system). Some attempts in this direction are currently under consideration.

In this paper we address another important issue concerning life time effects of the quasiparticles. The high resolution ARPES data  reveal that normal phase of underdoped cuprates has the marginal Fermi liquid type properties where the quasiparticle weight and life time are very sensitive to temperature. We want to check how the life times of fermions and bosons depend on temperature $`T`$ and momentum using the BF model scenario.

Flow equation for any arbitrary operator, say $`\widehat{A}(l)`$, is given by

$`{\displaystyle \frac{d\widehat{A}(l)}{dl}}=[\widehat{\eta }(l),\widehat{A}(l)]`$ (1)
where $`\widehat{\eta }`$ is a generating operator of the continuous canonical transformation. The initial BF model Hamiltonian consists of the free part

$`\widehat{H}_0={\displaystyle \underset{k,\sigma }{}}(\epsilon _k^\sigma \mu )\widehat{c}_{k\sigma }^{}\widehat{c}_{k\sigma }+{\displaystyle \underset{q}{}}(E_q2\mu )\widehat{b}_q^{}\widehat{b}_q`$ (2)
and the interaction

$`\widehat{H}_{int}={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{k,p}{}}\left(v_{k,p}\widehat{b}_{k+p}^{}\widehat{c}_k\widehat{c}_p+v_{k,p}^{}\widehat{b}_{k+p}\widehat{c}_p^{}\widehat{c}_k^{}\right).`$ (3)
In order to decouple fermion from boson subsystems we have previously chosen $`\eta `$ in a form

$`\widehat{\eta }(l)={\displaystyle \underset{k,q}{}}(\alpha _{k,q}(l)v_{k,q}(l)\widehat{b}^{}(l)\widehat{c}_{k,}(l)\widehat{c}_{q,}(l)h.c.)`$ (4)
with $`\alpha _{k,q}=\epsilon _k^{}+\epsilon _q^{}E_{k+q}`$. Using the general flow equation (1) we transformed the Hamilton operator $`\widehat{H}(l)`$ so that in the limit $`l\mathrm{}`$ the charge exchange interaction got eliminated $`v_{k,q}(\mathrm{})=0`$.

In a course of the continuous canonical transformation the initial boson and fermion operators change according to the following flow equations

$`{\displaystyle \frac{d\widehat{b}_q(l)}{dl}}`$ $`=`$ $`{\displaystyle \underset{k}{}}\alpha _{k,qk}(l)v_{k,qk}(l)\widehat{c}_{k,}(l)\widehat{c}_{qk,}(l)`$ (5)
$`{\displaystyle \frac{d\widehat{c}_{k,}(l)}{dl}}`$ $`=`$ $`{\displaystyle \underset{q}{}}\alpha _{q,k}(l)v_{k,q}(l)\widehat{b}_{q+k}(l)\widehat{c}_{q,}^{}(l).`$ (6)
Similar equations can be derived for their hermitean conjugates $`\widehat{b}_q^{}(l)`$, $`\widehat{c}_{k,}^{}(l)`$. By inspecting (5) we notice that initial boson annihilation operator is coupled to the electron pair operator $`\widehat{c}_{k,}(l)\widehat{c}_{qk,}(l)`$. Thereof one should in a next step determine a corresponding flow equation for the finite momentum pair operator, however this process will never end (like equations of motion for the Greens functions). There would be the higher and higher order product of operators getting involved into the differential flow equations.

Following the other works using such continuous canonical transformation - we terminate the flow equations by postulating some (physical) Ansatz. In a context of the BF model it seems natural to postulate the following decompositions

$`\widehat{b}_q(l)`$ $`=`$ $`X_q(l)\widehat{b}_q+{\displaystyle \underset{k}{}}Y_{q,k}(l)\widehat{c}_{k,}\widehat{c}_{qk,},`$ (7)
$`\widehat{c}_{k,\sigma }(l)`$ $`=`$ $`A_k(l)\widehat{c}_{k,\sigma }+{\displaystyle \underset{q}{}}B_{k,q}(l)\widehat{b}_{q+k}\widehat{c}_{q,\sigma }^{},`$ (8)
where $`\widehat{b}_q\widehat{b}_q(l=0)`$ and $`\widehat{c}_{k,\sigma }\widehat{c}_{k,\sigma }(l=0)`$. The initial $`l=0`$ values are

$`\begin{array}{cc}A_k(0)=1,\hfill & B_{k,q}(0)=0,\hfill \\ X_q(0)=1,\hfill & Y_{k,q}(0)=0.\hfill \end{array}`$ (11)
Unknown parameters $`Y_{q,k}(l)`$, $`B_{k,q}(l)`$ can be found directly from the flow equations (5,6)

$`{\displaystyle \frac{dB_{k,q}(l)}{dl}}`$ $`=`$ $`\alpha _{q,k}(l)v_{q,k}(l)X_{q+k}(l)A_q(l)`$ (12)
$`{\displaystyle \frac{dY_{q,k}(l)}{dl}}`$ $`=`$ $`\alpha _{k,qk}(l)v_{k,qk}(l)A_k(l)A_{qk}(l).`$ (13)
There is no unambiguous way for determining the other missing coefficients $`A_k(l)`$, $`X_q(l)`$ because with the Ansatz (7,8) the flow equations still do not close. In what follows bellow we determine these parameters using the constraint

$`[b_q(l),b_q^{}^{}(l)]`$ $`=`$ $`\delta _{q,q^{}},`$ (14)
$`\{c_{k,\sigma }(l),c_{k^{},\sigma ^{}}^{}(l)\}`$ $`=`$ $`\delta _{k,k^{}}\delta _{\sigma ,\sigma ^{}}.`$ (15)
These statistical commutation/anticommutation relations yield

$`A_k(l)^2+{\displaystyle \underset{q}{}}B_{k,q}(l)^2\left(f_{q,}b_{k+q}\right)`$ $`=`$ $`1,`$ (16)
$`X_q(l)^2+{\displaystyle \underset{k}{}}Y_{q,k}(l)^2\left(1f_{k,}f_{qk,}\right)`$ $`=`$ $`1.`$ (17)
Here $`f_{k,\sigma }`$, $`b_q`$ denote the Fermi and Bose distribution functions. They appear after introducing the normal ordered forms for operators left of the commutation (14) and anticommutation (15) when substituting (7,8).

We solved numerically a set of the coupled equations (12,13,16,17) simultaneously with the additional flow equations for $`\epsilon _k^\sigma (l)`$,$`E_k(l)`$ and $`v_{k,q}(l)`$ given in the Ref. . Initially, for $`l=0`$, we assumed: a) the one dimensional tight binding dispersion $`\epsilon _k^\sigma (l=0)=2t\mathrm{cos}k`$, b) localized bosons $`E_q(l=0)=\mathrm{\Delta }_B`$ and c) local exchange interaction $`v_{k,q}(l=0)=v`$. For the sake of comparison with the previous studies of this model we used $`D=4t1`$, $`\mathrm{\Delta }_B=0.6`$, $`v=0.1`$ and total concentration of carries fixed to be 1.

In figure 1 we plot the momentum dependence of $`A_k`$ coefficient for $`l=\mathrm{}`$ (when fermion subsystem is decoupled from the boson one). Considerable variation of the parameter $`A_k`$ from its initial value (being $`1`$) can be observed for such momenta which are located near $`k^{}0.89`$. As discussed in the Ref.  this situation corresponds to the resonant scattering processes when $`\epsilon _k+\epsilon _k=\mathrm{\Delta }_B`$. It is worth mentioning that for $`n_{tot}=1`$ the Fermi momentum approaches this $`k^{}`$ from bellow when temperature decreases. We could think of the coefficient $`A_k`$ as a measure of the spectral weight for fermions. The missing part $`|1A_k|`$ is transfered towards some composite objects (mixtures of fermions and bosons or fermion pairs). Taking into account the fact that for a decreasing temperature $`k_F`$ shifts closer and closer to $`k^{}`$, we shall thus observe a diminishing spectral weight of fermions. This is in agreement with interpretation of the ARPES experimental results .

For a completeness we illustrate also the resulting ($`l=\mathrm{}`$) values of the boson coefficients given in (7). Parameter $`X_q`$ shows only a small deviation from its initial value $`1`$. This time such a change is most effective for boson momenta close to $`2k^{}`$. It is exactly the same region of the Brillouin zone where the effective boson dispersion shows a kink (see figure 3).

The missing part of the boson spectral weight is transferred to the fermion pairs, as given by the Ansatz (7). Bottom figure 2 shows that $`_kY_{q,k}`$ (which roughly measures the fermion pairs spectral weight) changes a lot for $`q2k^{}`$. Again, if we recall that for a decreasing temperature $`k_Fk^{}`$ we can conclude that there is an increasing probability of finding more and more fermion pairs, even though our system is not in the superconducting phase. This is one of possible ways for observing the precursor effect.

In conclusion we studied evolution of fermion and boson operators in the BF model. During a continuous transformation which eliminated the exchange interaction we observe transformation of the initial fermion and boson particles into the composite objects. In particular we find that fermion (boson) spectral weights at $`k_F`$ ($`2k_F`$) are reduced when temperature decreases. This effect is in an agreement with the recent interpretation of the ARPES measurements for underdoped Bi2212 cuprates. To clarify a problem of the quasiparticles life time we need a more detailed study of the correlation functions. Such an analysis is currently under investigation and the results will be published elsewhere.


text:
# Explicit 𝑆𝑈(5) monopole solutions

## I Introduction

Monopoles are topological solitons in three space dimensions in Yang-Mills-Higgs gauge theories. They arise in theories in which the gauge group $`G`$ is spontaneously broken to a residual symmetry group $`H`$ by the action of a Higgs field which attains a nonzero vacuum expectation value (VEV), and in which the vacuum manifold $`G/H`$ has non-trivial second homotopy group. Bogomolnyi-Prasad-Sommerfield monopoles Bog-76 ; PS-75  occur in models with the Higgs field in the adjoint representation, in the limiting case in which the Higgs potential vanishes, but the nonzero VEV is kept as a boundary condition. They are static solutions of a first-order partial differential equation, the Bogomolnyi equation, and achieve the lower bound on monopole energy which was derived by Bogomolnyi Bog-76 . The simplicity of the Bogomolnyi equation, relative to the full second-order field equations, makes it easier both to derive exact solutions, and to study the structure of the moduli space of all solutions, and a number of methods have been developed to study them (see Sut  for a review).

In this paper, we describe in detail how BPS monopoles in several (possibly different) theories may be embedded simultaneously into a theory with a larger gauge group. This generalizes a construction used by a number of authors to embed $`SU(2)`$ monopoles into larger theories (see e.g. Wein ). We apply this embedding scheme to present explicit fields for several monopoles in the $`SU(5)[SU(3)\times SU(2)\times U(1)]/_6`$ model. It is probably well known to experts that an embedding scheme for BPS monopoles such as is described here exists, but it has not been described in the literature in this level of generality; nor has such a scheme been applied to construct monopole solutions in a phenomenologically interesting model such as the $`SU(5)`$ model considered here.

The construction here is in a sense complementary to that described by Lepora in Lepora . Lepora discusses embeddings of ’t Hooft-Polyakov monopoles into larger gauge groups away from the BPS limit; by working in the BPS limit we are able here to present explicit solutions. The embedding can be treated more simply in the BPS limit as well, by taking advantage of the fact that the Bogomolnyi equation is almost linear, in a sense discussed in the next section. Working with the Bogomolnyi equation allows us to avoid considering topological issues, and also makes trivial the matter of showing that an embedding does indeed produce a solution.

## II The general construction

We consider a model with gauge group $`G`$ with Lie algebra $`𝔤`$. The Higgs field $`\mathrm{\Phi }`$ is a $`𝔤`$-valued scalar field on $`^3`$ and the gauge potential $`A`$ is a $`𝔤`$-valued one-form on $`^3`$. The field strength is defined by $`F=dA+AA`$, and the magnetic field is $`B=F`$ where $``$ is the duality operator on forms. The covariant derivative of the Higgs field is $`D\mathrm{\Phi }=d\mathrm{\Phi }+[A,\mathrm{\Phi }]`$. For simplicity of notation, we take the coupling strength to be unity, and the generators of $`G`$ to be antihermitian. Also for simplicity, we work in a fixed gauge.

The Bogomolnyi equation is

$$D\mathrm{\Phi }=\pm B.$$
(1)
For our purposes, a monopole (or $`𝔤`$-monopole) is a pair $`(\mathrm{\Phi },A)`$ which satisfies the Bogomolnyi equation, along with the boundary conditions

$`\underset{r\mathrm{}}{lim}\mathrm{\Phi }(0,0,r)`$ $`=`$ $`\stackrel{~}{\mathrm{\Phi }}`$
$`\underset{r\mathrm{}}{lim}rB𝐫`$ $`=`$ $`M`$
for fixed constants $`\stackrel{~}{\mathrm{\Phi }},M𝔤`$, where $`𝐫`$ is the radial vector field on $`^3`$ and $`r=𝐫`$. $`M`$ is called the charge of the monopole, and we call $`\stackrel{~}{\mathrm{\Phi }}`$ the Higgs field at $`(0,0,\mathrm{})`$. Note that by possibly making the replacement $`\mathrm{\Phi }\mathrm{\Phi }`$, we are free to consider only the case where the $``$ sign occurs in (1).

We first wish to note that the Bogomolnyi equation is translation-invariant, so a monopole may be arbitrarily recentered. Second, note that the Bogomolnyi equation is invariant under the rescaling

$$(\mathrm{\Phi }(𝐫),A(𝐫))(\lambda \mathrm{\Phi }(\lambda 𝐫),\lambda A(\lambda 𝐫))$$
(2)
for any constant $`\lambda >0`$. This means that we may arbitrarily rescale the Higgs field at $`(0,0,\mathrm{})`$. Note also that the rescaling (2) leaves the charge $`M`$ unchanged.

Suppose now that we have $`n`$ pairwise commuting Lie subalgebras $`𝔤_1,𝔤_2,\mathrm{},𝔤_n`$ of $`𝔤`$ (i.e., for every $`g𝔤_i`$ and $`h𝔤_j`$, we have $`[g,h]=0`$ whenever $`ij`$) and for each $`i=1,2,\mathrm{},n`$ a $`𝔤_i`$-monopole $`(\mathrm{\Phi }_i,A_i)`$ (that is, $`(\mathrm{\Phi }_i,A_i)`$ satisfies the Bogomolnyi equation (1) for $`𝔤_i`$). Let $`\stackrel{~}{\mathrm{\Phi }}_i`$ denote the Higgs field at $`(0,0,\mathrm{})`$ for $`(\mathrm{\Phi }_i,A_i)`$ and $`M_i`$ its charge. (Note that throughout this section, subscripts will always refer to these subalgebras, not components.) Suppose that we also have a constant $`C𝔤`$ which commutes with each of the $`𝔤_i`$. We first make the definitions

$`A`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}A_i`$ (3)
$`\mathrm{\Phi }`$ $`=`$ $`C+{\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{\Phi }_i.`$ (4)
Now since $`𝔤_i`$ and $`𝔤_j`$ commute for $`ij`$, $`A_iA_j=0`$ and $`[A_i,\mathrm{\Phi }_j]=0`$ in that case. This implies that

$$B=\underset{i=1}{\overset{n}{}}B_i\text{ and }D\mathrm{\Phi }=\underset{i=1}{\overset{n}{}}D\mathrm{\Phi }_i.$$

Therefore $`(\mathrm{\Phi },A)`$ satisfies the Bogomolnyi equation (1) if each of the $`(\mathrm{\Phi }_i,A_i)`$ does. In this case the Higgs field at $`(0,0,\mathrm{})`$ is

$$\stackrel{~}{\mathrm{\Phi }}=C+\underset{i=1}{\overset{n}{}}\stackrel{~}{\mathrm{\Phi }}_i$$

and the monopole has charge

$$M=\underset{i=1}{\overset{n}{}}M_i.$$

In more physical terms, what we have is a monopole which is a superposition of simpler monopoles which live in mutually noninteracting sectors of the gauge theory. The essential point that makes this work is that although the Bogomolnyi equation (1) is nonlinear, both $`B`$ and $`D\mathrm{\Phi }`$, and therefore the Bogomolnyi equation, are additive with respect to sums of fields which commute with each other.

In practice, one begins with the Higgs field at $`(0,0,\mathrm{})`$, $`\stackrel{~}{\mathrm{\Phi }}`$, and charge $`M`$ for a monopole and wishes to find corresponding fields $`(\mathrm{\Phi },A)`$. To apply the above construction one must first identify appropriate subalgebras of $`𝔤`$ for which monopole solutions are known and such that their charges add up to $`M`$. In order to match the Higgs field at $`(0,0,\mathrm{})`$, it may be necessary to rescale the known solutions as in (2) and add a constant $`C`$ which commutes with the subalgebras of $`𝔤`$. We will illustrate this process in the next section. Note that the above construction does not imply that a monopole corresponding to a particular Higgs field at $`(0,0,\mathrm{})`$ and charge must arise from this kind of embedding.

## III BPS monopoles in the $`SU(5)`$ model

Consider now the symmetry-breaking pattern $`SU(5)[SU(3)\times SU(2)\times U(1)]/_6`$. This symmetry-breaking is achieved by a Higgs field with VEV

$$\stackrel{~}{\mathrm{\Phi }}=i\mathrm{diag}(2,2,2,3,3).$$
(5)
The fundamental monopole in this model has charge

$$M=\frac{i}{2}\mathrm{diag}(0,0,1,1,0).$$

To construct the fields for this monopole, we use the following copy of $`𝔰𝔲(2)`$ contained in $`𝔰𝔲(5)`$:

$$\left(\begin{array}{ccccc}0& & & & \\ & 0& & & \\ & & \alpha & \beta & \\ & & \overline{\beta }& \alpha & \\ & & & & 0\end{array}\right),$$
(6)
and the $`𝔰𝔲(2)`$-monopole solution found by Prasad and Sommerfield in PS-75 . The Prasad-Sommerfield solution is given by

$`\mathrm{\Phi }`$ $`=`$ $`f(r){\displaystyle \underset{i=1}{\overset{3}{}}}x_i\tau _i,`$ (7)
$`A`$ $`=`$ $`g(r){\displaystyle \underset{i,j,k=1}{\overset{3}{}}}ϵ_{ijk}x_i\tau _jdx^k,`$ (8)
where

$`f(r)`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}(r\mathrm{coth}r1),`$ (9)
$`g(r)`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}(r\mathrm{csch}r1),`$ (10)
and $`\tau _j=\frac{i}{2}\sigma _j`$. The Prasad-Sommerfield monopole (7, 8) has $`\mathrm{\Phi }(0,0,\mathrm{})=\tau _3`$ and charge $`M=\tau _3`$.

Embedding the Prasad-Sommerfield monopole into $`𝔰𝔲(5)`$ as indicated in (6) gives a solution with the desired charge; however to get the Higgs field at $`(0,0,\mathrm{})`$ to match (5), we must rescale as in (2) and add some constant matrix $`C`$ to the Higgs field which commutes with the copy of $`𝔰𝔲(2)`$ given by (6). To see how to achieve this, we write $`C=i\mathrm{diag}(a,b,c,c,ab2c)`$, and let $`\lambda `$ be the rescaling factor. When we rescale the solution as in (2), the Higgs field at $`(0,0,\mathrm{})`$ gets multiplied by $`\lambda `$. Thus to match the Higgs field at $`(0,0,\mathrm{})`$ we must solve the equation

$$(2,2,2,3,3)=(0,0,\lambda /2,\lambda /2,0)+(a,b,c,c,ab2c),$$

which has the solution $`\lambda =5`$, $`C=i\mathrm{diag}(2,2,1/2,1/2,3)`$. Rescaling the Prasad-Sommerfield monopole according to (2) with $`\lambda =5`$ and adding $`C`$ to the rescaled Higgs field, we obtain the fields

$`\mathrm{\Phi }`$ $`=`$ $`\left(\begin{array}{cccc}2i& & & \\ & 2i& & \\ & & \frac{1}{2}iI_2+_{j=1}^325f(5r)x_j\tau _j& \\ & & & 3i\end{array}\right)`$ (15)
$`A`$ $`=`$ $`{\displaystyle \underset{i,j,k=1}{\overset{3}{}}}\left(\begin{array}{cccc}0& & & \\ & 0& & \\ & & 25g(5r)ϵ_{ijk}x_i\tau _j& \\ & & & 0\end{array}\right)dx^k,`$ (20)
where $`f(r)`$ and $`g(r)`$ are as given by (9) and (10), and $`I_2`$ is the $`2\times 2`$ identity matrix. This solution has been used by Pogosian and Vachaspati in PV  as a starting point for numerical study of the interaction of monopoles and domain walls in the $`SU(5)`$ model.

We could of course have equally well used another copy of $`𝔰𝔲(2)`$ for this embedding, such as

$$\left(\begin{array}{ccccc}0& & & & \\ & \alpha & & & \beta \\ & & 0& & \\ & & & 0& \\ & \overline{\beta }& & & \alpha \end{array}\right).$$
(21)
It is obvious however that this would result in an $`𝔰𝔲(5)`$ monopole which is gauge-equivalent to the one constructed above. The same would be true of any other embedding of the same kind, which embeds $`𝔰𝔲(2)`$ matrices as $`2\times 2`$ submatrices of an $`𝔰𝔲(5)`$ matrix. To determine whether $`𝔰𝔲(2)`$ monopoles may be embedded into $`𝔰𝔲(5)`$ in any other way, we note first that any such embedding determines a five-dimensional representation of $`𝔰𝔲(2)`$. Such a representation must be a direct sum of irreducible representations. Since the VEV for the $`SU(5)`$ model has only two distinct eigenvalues, the same must be true of the Higgs field at $`(0,0,\mathrm{})`$ for the embedded $`𝔰𝔲(2)`$ monopole. But the number of distinct eigenvalues for any generator for an irreducible representation of $`𝔰𝔲(2)`$ is equal to the dimension of the representation. So the only possible embeddings of an $`𝔰𝔲(2)`$ monopole into $`𝔰𝔲(5)`$ must define a representation of $`𝔰𝔲(2)`$ which is a direct sum of two-dimensional representations and trivial representations. Since at most two two-dimensional representation can fit into $`𝔰𝔲(5)`$, there are only two possibilities. If the representation of $`𝔰𝔲(2)`$ contains only one irreducible two-dimensional representation, we obtain the embedding discussed above. If it contains two irreducible two-dimensional representations, then we obtain a special case of the next embedding to be discussed (namely, the case in which the two embedded $`𝔰𝔲(2)`$ monopoles have the same center).

The next monopole of interest in the $`SU(5)`$ model has charge $`M=\frac{i}{2}\mathrm{diag}(0,1,1,1,1)`$. In topological terms, this monopole has winding two. In the same way we obtained the fields of the fundamental monopole above, we could obtain the fields for this monopole by rescaling and embedding the fields for a monopole in a theory with gauge group $`SU(4)`$ and a Higgs field with VEV $`\stackrel{~}{\mathrm{\Phi }}=\frac{i}{2}\mathrm{diag}(1,1,1,1)`$, i.e., with $`SU(4)[SU(2)\times SU(2)\times U(1)]/[_2\times _2]`$ symmetry breaking. Unfortunately, the general explicit solutions for this model are not known. We can, however, obtain a family of such solutions by instead using the following two commuting copies of $`𝔰𝔲(2)`$ in $`𝔰𝔲(5)`$:

$$\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& \alpha _1& 0& 0& \beta _1\\ 0& 0& \alpha _2& \beta _2& 0\\ 0& 0& \overline{\beta }_2& \alpha _2& 0\\ 0& \overline{\beta }_1& 0& 0& \alpha _1\end{array}\right),$$
(22)
embedding a rescaled Prasad-Sommerfield monopole into each of the positions indicated by (22) exactly as was done in (15, 20). In this case, to find rescaling factors $`\lambda _1`$ and $`\lambda _2`$ (one for each embedded $`𝔰𝔲(2)`$ monopole) and constant matrix $`C=i\mathrm{diag}(2a2b,a,b,b,a)`$ to make the Higgs field at $`(0,0,\mathrm{})`$ match, we solve

$$(2,2,2,3,3)=(0,\lambda _1/2,0,0,\lambda _1/2)+(0,0,\lambda _2/2,\lambda _2/2,0)+(2a2b,a,b,b,a),$$

and obtain $`\lambda _1=\lambda _2=5`$, $`C=i\mathrm{diag}(2,1/2,1/2,1/2,1/2)`$. Since each of the $`𝔰𝔲(2)`$-monopoles may be recentered independently, and is spherically symmetric, this construction results in a family of $`𝔰𝔲(5)`$-monopoles which are symmetric about the axis connecting the centers of the $`𝔰𝔲(2)`$-monopoles. If we assume the $`𝔰𝔲(2)`$ monopole which is embedded as in (6) is centered at the origin and the $`𝔰𝔲(2)`$ monopole which is embedded as in (21) is centered at $`𝐜=(c_1,c_2,c_3)`$, then we obtain the fields

$`\mathrm{\Phi }`$ $`=`$ $`C+{\displaystyle \underset{i=1}{\overset{3}{}}}25[f(5r)x_i\tau _i^{}+f(5𝐫𝐜)(x_ic_i)\tau _i^{\prime \prime }],`$ (23)
$`A`$ $`=`$ $`{\displaystyle \underset{i,j,k=1}{\overset{3}{}}}25[g(5r)ϵ_{ijk}x_i\tau _j^{}+g(5𝐫𝐜)ϵ_{ijk}(x_ic_i)\tau _j^{\prime \prime }]dx^k,`$ (24)
where $`C=i\mathrm{diag}(2,1/2,1/2,1/2,1/2)`$, $`\tau _i^{}`$ is $`\tau _i`$ embedded in $`𝔰𝔲(5)`$ as in (6), and $`\tau _i^{\prime \prime }`$ is $`\tau _i`$ embedded in $`𝔰𝔲(5)`$ as in (21).

Finally, we wish to consider monopoles with charge $`\frac{i}{2}\mathrm{diag}(1,1,1,1,2)`$, i.e., monopoles with winding three. In this case, we will construct the monopole by embedding a rescaled Prasad-Sommerfield monopole exactly as in (15,20), and an $`𝔰𝔲(3)`$ monopole with VEV $`i\mathrm{diag}(1,1,2)`$ and charge $`\frac{i}{2}\mathrm{diag}(1,1,2)`$ into the following $`𝔰𝔲(3)`$ subalgebra of $`𝔰𝔲(5)`$:

$$\left(\begin{array}{ccccc}\alpha _1& \beta & 0& 0& \gamma \\ \overline{\beta }& \alpha _2& 0& 0& \delta \\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ \overline{\gamma }& \overline{\delta }& 0& 0& \alpha _1\alpha _2\end{array}\right),$$
(25)
(An argument similar to the one outlined above for the fundamental monopole implies that an embedding of an $`𝔰𝔲(3)`$ monopole must be of essentially this type.) The general solution for such an $`𝔰𝔲(3)`$ monopole is not known; however the Higgs field for solutions with spherical symmetry has been found by Dancer Dancer , using the Nahm method Nahm . The Nahm data found by Dancer could in principle also be used to compute the gauge fields of such monopoles, although this calculation has not been carried out, but the Higgs field is sufficient to find the energy density of a monopole, which is proportional to $`\mathrm{\Delta }\mathrm{tr}(\mathrm{\Phi }^2)`$, where $`\mathrm{\Delta }`$ is the Laplacian on $`^3`$ Ward . Dancer’s solutions are given in the form

$$\mathrm{\Phi }=i\mathrm{diag}(\mathrm{\Phi }_{11}(r),\mathrm{\Phi }_{22}(r),\mathrm{\Phi }_{33}(r)),$$
(26)
where

$`\mathrm{\Phi }_{11}(r)`$ $`=`$ $`{\displaystyle \frac{e^{6r}\left(\frac{3}{2r}+\frac{1}{r^2}\frac{1}{4r^3}\right)+\frac{1}{4r^3}+\frac{1}{2r^2}}{e^{6r}\left(\frac{3}{2r}\frac{1}{4r^2}\right)+\frac{1}{4r^2}}},`$
$`\mathrm{\Phi }_{22}(r)`$ $`=`$ $`(\mathrm{\Phi }_{11}(r)+\mathrm{\Phi }_{33}(r)),`$ (27)
$`\mathrm{\Phi }_{33}(r)`$ $`=`$ $`\mathrm{\Phi }_{11}(r).`$
To match our $`𝔰𝔲(5)`$ Higgs field at $`(0,0,\mathrm{})`$, we must rescale Dancer’s solutions by a factor of $`\lambda =5/3`$ (and also multiply by $`1`$ since Dancer uses the opposite sign convention in the Bogomolnyi equation), and add to $`\mathrm{\Phi }`$ the constant matrix $`C=i\mathrm{diag}(1/3,1/3,1/2,1/2,1/3)`$. We thus obtain the Higgs field

$$\mathrm{\Phi }=\left(\begin{array}{cccc}\frac{5}{3}i\mathrm{\Phi }_{11}(\frac{5}{3}r)+\frac{1}{3}i& & & \\ & \frac{5}{3}i\mathrm{\Phi }_{22}(\frac{5}{3}r)+\frac{1}{3}i& & \\ & & \frac{1}{2}iI_2+_{j=1}^325f(5r)x_j\tau _j& \\ & & & \frac{5}{3}i\mathrm{\Phi }_{33}(\frac{5}{3}r)+\frac{1}{3}i\end{array}\right),$$
(28)
where $`\mathrm{\Phi }_{11}`$, $`\mathrm{\Phi }_{22}`$, and $`\mathrm{\Phi }_{33}`$ are from (27), $`f(r)`$ is given in (9), and $`I_2`$ is again the $`2\times 2`$ identity matrix. Again, the $`𝔰𝔲(2)`$ and $`𝔰𝔲(3)`$ monopole may be independently recentered, so that we obtain a family of axially symmetric monopoles.

We have restricted our attention above to embeddings for which explicit solutions can be given based on solutions in the literature, and which furthermore are likely to be of interest as a first approximation to solutions outside the BPS limit (as in PV  for example). One could also consider, for example, $`𝔰𝔲(5)`$ monopoles with charge $`\frac{i}{2}\mathrm{diag}(0,0,k,k,0)`$ which arise by embedding an $`𝔰𝔲(2)`$ monopole with charge $`k\tau _3`$ as in (6), where $`k>1`$. It is known that there are such $`𝔰𝔲(2)`$ monopoles in the BPS limit, but outside the BPS limit such monopoles are unstable. (In the BPS case, there is a long-range Higgs attraction between fundamental monopoles which exactly cancels the magnetic repulsion.) For a review of cases in which exact solutions are known, see Sut .

###### Acknowledgements.

The author would like to thank Prof. T. Vachaspati for guidance and advice throughout the work that led to this paper, and for suggesting that this paper be written. The author would also like to thank Prof. E. J. Benveniste for several helpful discussions.


text:
# Triangular and 𝑌-shaped hadrons in QCD

## Abstract

Gauge invariant extended configurations are considered for the three fundamental (quarks) or adjoint (gluons) particles. For quarks it is shown that the $`Y`$-shaped configuration is the only possible. For adjoint sources both the $`Y`$-shaped and triangular configurations may realize. The corresponding static potentials are calculated in the Method of Field Correlators and in the case of baryon shown to be consistent with the lattice simulations. For adjoint sources the potentials of $`Y`$-shaped and $`\mathrm{\Delta }`$-shaped configurations turn out to be close to each other, which leads to almost degenerate masses of $`3^{}`$ 3g glueballs and odderon trajectories.

1. To make conclusions on the structure of gluonic fluxes confining color charges in physical states, one has to start from the considering of the space-extended gauge invariant wave function of the hadron. It is easy then to show that in the case of the static charges the Green function of the hadron reduces to the Wilson loop, which in the case of the baryon has the $`Y`$-type shape and consists from the three contours formed by the quark trajectories and joined at the point of the string junction . The Wilson loop of the tree adjoint sources is less known than the 3q one. In the paper we will show that it can have both $`Y`$-type and $`\mathrm{\Delta }`$-type shape.

Using the formalism of the Method of the Field Correlators (MFC), we compute the static potentials, corresponding to the Wilson loops of the hadrons. In the case of the baryon the static potential was used long ago in many dynamical calculations . Recently it has been computed in lattice gauge theory in a number of papers . We show in the paper that our potential is in a full agreement with the lattice studies.

In the case of adjoint sources we find that the $`Y`$-type and $`\mathrm{\Delta }`$-type potentials remain near each other at the characteristic hadronic size. Using that we estimate masses of lowest $`3g`$ glueballs, lying on the corresponding odderon trajectories, and show that they are close to each other, implying that there are two possible odderon trajectories with not much different Regge slopes. A short discussion of physical implications of these results concludes the paper.

To avoid confusion, we should stress that the term ”$`\mathrm{\Delta }`$ configuration” used in  in the context of the static baryon potential, refers to the perimeter behaviour of the potential and not to the gauge invariant configurations as well as to the structures of fluxes discussed in the present paper.

2. Hadron building in SU(3) starts with listing elementary building blocks: quarks $`q^\alpha ,\alpha =1,2,3`$, gluons (or adjoint static sources) $`g^a,a=1,\mathrm{}8`$, parallel transporters (PT) in fundamental representation $`\mathrm{\Phi }_\alpha ^\beta (x,y)=(P\mathrm{exp}igA_\mu (z)𝑑z_\mu )_\alpha ^\beta `$, adjoint parallel transporters $`\mathrm{\Phi }_{ab}(x,y)`$, generators $`t_\alpha ^{(a)\beta }`$ symmetric symbols $`\delta _\alpha ^\beta ,\delta _{ab},d^{abc}`$, and antisymmetric ones, $`e_{\alpha \beta \gamma }`$ and $`f^{abc}`$. Note that we always use Greek indices for fundamental representation and Latin ones for the adjoint.

To construct a real extended (not point-like) hadron one uses all listed elements, PT included, and forms a white (gauge-invariant) combination. It is convenient to form an extended quark (antiquark) operator

$$q^\alpha (x,Y)q^\beta (x)\mathrm{\Phi }_\beta ^\alpha (x,Y);$$

$$\overline{q}_\alpha (x,Y)=\overline{q}_\beta (x)\mathrm{\Phi }_\alpha ^\beta (x,Y).$$
(1)
In this way one has for the Y-shaped baryon:

$$B_Y(x,y,z,Y)=e_{\alpha \beta \gamma }q^\alpha (x,Y)q^\beta (y,Y)q^\gamma (z,Y).$$
(2)
One can also define quark operator with two lower indices: $`e_{\alpha \beta \gamma }q^\alpha (x)q_{\beta \gamma }(x)`$. However an attempt to create a gauge-invariant combination from 3 operators $`q_{\beta \gamma }(x)`$ and 3 PT to construct a $`\mathrm{\Delta }`$-type configuration fails: the structure

$$B_\mathrm{\Delta }(x,y,z)=q_{\alpha \beta }(x)\mathrm{\Phi }_\gamma ^\beta (x,y)q_{\gamma \delta }(y)\mathrm{\Phi }_\epsilon ^\delta (y,z)q_{\epsilon \rho }(z)\mathrm{\Phi }_\alpha ^\rho (z,x)$$
(3)
is not gauge invariant, that can be checked directly, substituting in (3) $`q^\alpha (x)U_\beta ^\alpha (x)q^\beta (x)`$. One can try all combinations, but it is impossible to form a continuous chain of indices to represent the $`\mathrm{\Delta }`$-type structure using as operators $`q_\alpha `$ as $`q_{\alpha \beta }`$. Thus one can conclude that the $`Y`$-shaped configuration is the only possible gauge-invariant configuration of wave function for baryons.

One may wonder, what is the relation between the spacial structures of wave function of hadrons and their gluonic flux? The answer from the flux-tube models is that these structures coincide. In realistic lattice calculations one is to use the Wilson loop, which describes the gauge invariant state of hadron generated at some initial and annihilated at final moment of time. As is well known, for static charges the Wilson loop consists of two wave functions considered above, joined by parallel transporters. Let us imagine now some evolution of the fluxes in baryon which would lead to the emerging of the $`\mathrm{\Delta }`$-type configuration of fluxes in some intermediate time. First of all we should note that the cross-section of the Wilson loop by the time-like hypersurface will recover a gauge-invariant 3q state, i.e. the $`Y`$-type configuration. To have a $`\mathrm{\Delta }`$-shape for fluxes one should admit that the fluxes have no relation with the wave function, which is unprobable.

Consider now the adjoint source $`g^a(x)t_\alpha ^{(a)\beta }G_\alpha ^\beta (x)`$. We do not specify here the Lorentz structure of $`g^a(x)`$, but only impose condition that it should gauge transform homogeneously, $`G_\alpha ^\beta U_\beta ^{}^{+\beta }G_\alpha ^{}^\beta ^{}U_\alpha ^\alpha ^{}`$. Therefore $`g^a(x)`$ can be either the field strength $`F_{\mu \nu }^a(x)`$, or valence gluon field $`a_\mu ^a(x)`$ in the background-field perturbation theory . It is easy to construct a $`\mathrm{\Delta }`$-type configuration for 3 such sources;

$$G_\mathrm{\Delta }(x,y,z)=G_\alpha ^\beta (x)\mathrm{\Phi }_\beta ^\gamma (x,y)G_\gamma ^\delta (y)\mathrm{\Phi }_\delta ^\epsilon (y,z)G_\epsilon ^\rho (z)\mathrm{\Phi }_\rho ^\alpha (z,x).$$
(4)
It is clear that in (4) all repeated indices form gauge-invariant combinations, and $`G_\mathrm{\Delta }(x,y,z)`$ is a gauge-invariant $`\mathrm{\Delta }`$-type configuration, which was used previously for the $`3g`$ glueball in .

But one can persuade oneself that (4) is not the only $`3g`$ gauge-invariant configuration. Consider adjoint sources and adjoint PT (here distinguishing upper and lower indices is not necessary) and form as in (1) an extended gluon operator:

$$g_a(x,Y)g^b(x)\mathrm{\Phi }_{ab}(x,Y)$$
(5)
and an $`Y`$-shaped configuration

$$G_Y^{(f)}(x,y,z,Y)=f^{abc}g_a(x,Y)g_b(y,Y)g_c(z,Y).$$
(6)
In the same way one constructs $`G_Y^{(d)}`$replacing $`f`$ by $`d`$ in (6).It is clear that $`G_Y`$ is gauge-invariant and should be considered on the same grounds as $`G_\mathrm{\Delta }`$.

At this point it is necessary to clarify how (2), (6) generate Green’s functions and Wilson loops.

To this end consider initial and final states made of (2), (4), (6) and for simplicity of arguments take all fundamental and adjoint sources to be static, i.e. propagating only in Euclidean time.

Then the Green’s function for the object will be

$$𝒢_i(\overline{X},X)=\mathrm{\Psi }_i^+(\overline{X})\mathrm{\Psi }_i(X)$$
(7)
where $`\mathrm{\Psi }_i=G_\mathrm{\Delta },G_Y,B_Y;X=x,y,z`$ for $`G_\mathrm{\Delta }`$ and $`x,y,z,Y`$ otherwise. Now it is important that vacuum average in (7) produces a product of Green’s functions for quarks or for valence gluons in the external vacuum gluonic field, which is proportional to the corresponding PT, fundamental – for quarks and adjoint – for gluons. Namely,

$$\overline{q}_\beta (\overline{x})q^\alpha (x)\mathrm{\Phi }_\beta ^\alpha (\overline{x},x),$$

$$g^a(\overline{x})g^b(x)\mathrm{\Phi }_{ab}(\overline{x},x).$$
(8)
(This statement is well known for static sources, for relativistic quarks and gluons this follows directly from the exact Fock-Feynman-Schwinger representation (FFSR), see  and for a review ).

As a result one obtains a gauge-invariant Wilson-loop combination for each Green’s function (7). In particular for $`B_Y`$ (2) one has a familiar 3-lobe Wilson loop $`W_Y`$:

$$W_Y(\overline{X},X)=\mathrm{tr}_Y\underset{i=1}{\overset{3}{}}W_i(C_i),$$
(9)
where $`\mathrm{tr}_Y=\frac{1}{6}e_{\alpha \beta \gamma }e_{\alpha ^{}\beta ^{}\gamma ^{}}`$, and the contour $`C_i`$ in the open loop $`W_i`$ passes from $`Y`$ to $`\overline{Y}`$ through points $`x,\overline{x}(i=1)`$, $`y,\overline{y}(i=2)`$, or $`z,\overline{z}(i=3)`$, as shown in Fig.1.

This situation is well-known and was exploited in numerous applications. Relatively less known are the Wilson-loop configurations for $`G_Y`$ and $`G_\mathrm{\Delta }`$. In the first case the structure is the same with the replacement of fundamental lines and symbols by the adjoint ones: $`e_{\alpha \beta \gamma }f^{abc}`$ or $`e_{\alpha \beta \gamma }d^{abc}`$, $`\mathrm{\Phi }_\alpha ^\beta \mathrm{\Phi }_{ab}`$, so that the whole structure in (9) is the same with this replacement. Contrary to the baryon case, we can contract adjoint indices in two ways, using antisymmetric symbol $`f^{abc}`$ or symmetric one $`d^{abc}`$. The proper choice is related to the Bose-statistics of the gluon system which ensures the full coordinate-spin function to be symmetric.

In the case of $`G_\mathrm{\Delta }`$ using (4) and (8) one can write the resulting structure symbolically as follows

$$𝒢_\mathrm{\Delta }(\overline{X},X)=\mathrm{\Delta }_{a^{}b^{}c^{}}(\overline{x},\overline{y},\overline{z})\mathrm{\Phi }_{a^{}a}(\overline{x},x)\mathrm{\Phi }_{b^{}b}(\overline{y},y)\mathrm{\Phi }_{c^{}c}(\overline{z},z)\mathrm{\Delta }_{abc}(x,y,z)$$
(10)
where we have denoted

$$\mathrm{\Delta }_{abc}(x,y,z)=t_\alpha ^{(a)\beta }\mathrm{\Phi }_\beta ^\gamma (x,y)t_\gamma ^{(b)\delta }\mathrm{\Phi }_\delta ^\epsilon (y,z)t_\epsilon ^{(c)\rho }\mathrm{\Phi }_\rho ^\alpha (z,x).$$
(11)
To understand better the structure of (10), one can use the large $`N_c`$ approximation, in which case one has

$$\mathrm{\Psi }_{\alpha \alpha ^{}}^{\beta \beta ^{}}t_\alpha ^{}^{(a^{})\beta ^{}}\mathrm{\Phi }_{a^{}a}(\overline{x},x)t_\alpha ^{(a)\beta }\frac{1}{2}\mathrm{\Phi }_\alpha ^\beta ^{}(x,\overline{x})\mathrm{\Phi }_\alpha ^{}^\beta (\overline{x},x).$$
(12)
As a result in this approximation $`𝒢_\mathrm{\Delta }`$ appears to be a product of 3 fundamental closed loops, properly oriented with respect to each other

$$𝒢_\mathrm{\Delta }(\overline{X},X)W(\overline{x},\overline{y}|x,y)W(\overline{y},\overline{z}|y,z)W(\overline{z},\overline{x}|z,x)W_\mathrm{\Delta }(\overline{X},X),$$
(13)
it is displayed in Fig. 2.

3. Static potentials for configurations (2), (4), (6) can be computed using Field Correlator Method (FCM) , through the equation

$$V=\underset{T\mathrm{}}{lim}\frac{1}{T}\mathrm{ln}W,$$
(14)
where $`T`$ is the time extension of the Wilson loop.

For the baryon in the case of 3 quarks at the vertices of an equilateral triangle, at the distance $`R`$ from the string junction Y, the static baryon potential reads as

$$V^{(B)}(R)=3V^{(M)}(R)+V^{(\mathrm{nd})}(R),$$
(15)
where

$$V^{(M)}(R)=\frac{2\sigma }{\pi }\left\{R_0^{R/T_g}𝑑xxK_1(x)T_g\left(2\frac{R^2}{T_g^2}K_2\left(\frac{R}{T_g}\right)\right)\right\}$$
(16)
is the mesonic confining potential with the asymptotic slope $`\sigma 0.18`$ GeV2 and the gluonic correlation length $`T_g=0.12÷0.2`$ fm , and the nondiagonal part of the potential,

$$V^{(\mathrm{nd})}(R)=\frac{2}{\sqrt{3}}\sigma T_g\frac{3\sqrt{3}}{2\pi }\frac{\sigma R^2}{T_g}_{\frac{\pi }{6}}^{\frac{\pi }{3}}\frac{d\phi }{\mathrm{cos}\phi }K_2\left(\frac{\sqrt{3}R}{2T_g\mathrm{cos}\phi }\right),$$
(17)
appears due to the interference of the gluonic fields on different lobes of the Wilson loop. Note the difference in the overall factor $`1/2`$ with the previous calculations , where it was erroneously omitted. Let us denote $`L3R`$ the total length of the string. In Fig. 3 from  the dependence of lattice nonperturbative baryon potential from  on $`L`$ along with the MFC potential (15)-(17) is shown. One can see that our potential is in the complete agreement with the lattice results. In the asymptotic region $`L\begin{array}{c}>\hfill \\ \hfill \end{array}1.5`$ fm the potential has a linear form

$$V^{(B)}(R)\sigma L+\left(\frac{2}{\sqrt{3}}\frac{12}{\pi }\right)\sigma T_g.$$
(18)
The dotted tangent in Fig. 3 demonstrates that in the range 0.3 fm$`\begin{array}{c}<\hfill \\ \hfill \end{array}L\begin{array}{c}<\hfill \\ \hfill \end{array}1.5`$ fm the lattice data can be described by the linear potential with the slope some 10% less than $`\sigma `$.

The potential written so far contains only the nonperturbative confining part. To obtain the total potential we should add to it the perturbative color-coulombic potential

$$V_{(\mathrm{fund})}^{\mathrm{pert}}(r)=\frac{3}{2}\frac{C_2(\mathrm{fund})\alpha _s}{r},$$
(19)
where $`r=\sqrt{3}R`$ is the interquark distance in the equilateral triangle and $`C_2(\mathrm{fund})=4/3`$.

In Fig. 4 lattice data from the last ref. of  and the potential $`V_B(r)+V_{(\mathrm{fund})}^{\mathrm{pert}}(r)`$ are shown. One can see that our results are in the complete agreement with this independent set of the lattice data as well.

In a similar way one can write the static potential for the adjoint sources, neglecting the nondiagonal term, which is different in symmetric and antisymmetric states:

$$V_Y^{(G)}(R)=\frac{C_2(\mathrm{adj})}{C_2(\mathrm{fund})}V_Y^{(B)}(R)=\frac{9}{4}V_Y^{(B)}(R).$$
(20)
Consider the $`\mathrm{\Delta }`$-configuration in the approximation (12). In this case $`V_\mathrm{\Delta }^{(G)}(R)`$ reduces to the sum of the mesonic potentials corresponding to area laws for all three loops minus the nondiagonal interference term, and one obtains just as in

$$V_\mathrm{\Delta }^{(G)}(r)=3V^{(M)}(r)+V^{(\mathrm{nd})}(r).$$
(21)
Along with the adjoint perturbative potential

$$V_{(\mathrm{adj})}^{\mathrm{pert}}(r)=\frac{3}{2}\frac{C_2(\mathrm{adj})\alpha _s}{r},$$
(22)
where $`C_2(\mathrm{adj})=3`$, we plot both $`V_Y^{(G)}`$ and $`V_\mathrm{\Delta }^{(G)}`$ in Fig. 5 without the interference terms. We see from the figure that the curves intersect at $`r0.5`$ fm, and are very close to each other.

4. To summarize our results, we have considered possible gauge-invariant configurations of 3 fundamental or adjoint sources and corresponding Wilson loops, which have $`Y`$-type shape for fundamental charges and may have both $`Y`$-type and $`\mathrm{\Delta }`$-type shapes for the adjoint ones. We have shown that the static baryon potential obtained in the MFC is in the complete agreement with the lattice data.

For adjoint sources it was demonstrated that two possible configurations yield static potentials differing only a little. This in turn implies that 3-gluon glueballs  may be of two distinct types, with no direct transitions between them (quark-containing hadrons must be involved as intermediate states). The mass of the $`\mathrm{\Delta }`$-shaped $`3^{}`$ glueball was found in  to be $`M_\mathrm{\Delta }^{(3g)}=3.51`$ GeV for $`\sigma _f=0.18`$ GeV2 (or 4.03 GeV for $`\sigma _f=0.238`$ GeV2 to be compared with lattice one calculated in  4.13$`\pm `$0.29 GeV). The mass of the $`Y`$-shaped glueball can easily be computed from the baryon mass calculated in , multiplying it by $`\sqrt{9/4}=3/2`$. In this way one obtains $`M_Y^{(3g)}=3.47`$ GeV ($`\sigma _f=0.18`$ GeV2). The slope of the corresponding odderon trajectory is almost the same and corresponds to $`ggg`$-configuration. Thus one obtains the $`\mathrm{\Delta }`$-odderon (slope)-1 to be twice the standard Regge slope, while for $`Y`$-odderon it is $`\frac{9}{4}`$ of the standard slope. In both cases the intercept comes out as in  to be rather low ($`1.8`$ for the $`Y`$-shape and $`2.4`$ for the $`\mathrm{\Delta }`$-shape) implying very small odderon contribution to reactions under investigation  in agreement with measurements. We plan to perform more accurate calculations of glueball potentials and spectra taking into account the string-string interference in subsequent publications.

Both authors acknowledge partial support from grants RFFI 00-15-96786, 00-02-17836 and INTAS 00-00110, 00-00366. One of the authors (Yu.S.) is grateful for useful discussions to D. Richards and R. Edwards; he was supported by DOE contract DE-AC05-84ER40150 under which SURA operates the Thomas Jefferson National Accelerator Facility.