format: text

text:
# Global Anisotropy Versus Small-Scale Fluctuations in Neutrino Flux in Core-Collapse Supernova Explosions

## 1 INTRODUCTION

For many years since the first work of Colgate & White (1966), numerical simulations of core-collapse supernova explosions have been exciting topics. Until the beginning of the 1990s, almost all the simulations included the assumption of spherical symmetry (e.g., Wilson, 1985). Such one-dimensional simulations, however, were unable to explain the observed explosion energy and often failed to produce explosions (see, e.g., Liebendoerfer et al. (2001)). Inclusion of convective motion via the mixing length theory cures this problem to some extent (Wilson & Mayle, 1993; Bruenn, Mezzacappa, & Dineva, 1995). However, spherical simulations based only on the Rayleigh-Taylor instability require considerably large initial fluctuations in density to explain the large-scale matter mixing. In addition, aspherical explosion is also supported by the observation of SN1987A, where asymmetric ejecta are clearly observed (e.g., Wang et al., 2002). These lead us to multidimensional simulations of supernova explosions.

At this time, the two- and three-dimensional simulations have been performed by several groups (Miller, Wilson, & Mayle, 1993; Herant et al., 1994; Burrows, Hayes, & Fryxell, 1995; Janka & Müller, 1996; Mezzacappa et al., 1998; Fryer & Heger, 2000; Shimizu et al., 2001; Fryer & Wallen, 2002; Kifonidis et al., 2003). In many multidimensional simulations, special attention is paid to the role of convection either near the surface of a nascent neutron star or in neutrino-heated regions above the neutrinosphere. It has been shown that large-scale mixing, caused by convection and convective overturn around the neutrino-heated region, increases the explosion energy and can trigger a successful explosion (Herant et al., 1994; Keil, Janka, & Müller, 1996; Janka & Müller, 1996).

Because supernova progenitors such as OB stars are generally observed to be fast rotators ($`200`$ km s<sup>-1</sup> at the surface, $`P1`$ day; see, e.g., Tassoul 1978; Fukuda 1982), the resulting proto$``$neutron star can have a large amount of angular momentum after the gravitational collapse. Centrifugal force then deforms the rotating core into an oblate form. This will cause asymmetric neutrino radiation, in which the flux along the pole is enhanced over that on the equatorial plane. Janka & Mönchmeyer (1989) first discussed the possibility of aspherical neutrino emission from a rapidly rotating inner core. They intended to evaluate the total neutrino energy outputs using the neutrino data of SN 1987A detected with the Kamiokande and Irvine-Michigan-Brookhaven experiments. It was argued in their paper that a neutrino flux along the pole might be up to a factor of 3 greater than that on the equatorial plane.

Inspired by this work, Shimizu and coworkers (Shimizu, Yamada, & Sato, 1994; Shimizu et al., 2001) proposed that the anisotropic neutrino radiation should play a critical role in the explosion mechanism itself, and carefully investigated the effects of anisotropic neutrino radiation on the explosion energy. They found that only a few percent enhancement in the neutrino emission along the pole is sufficient to increase the explosion energy by a large factor, and that this effect saturates around a certain degree of anisotropy. It should be noted here that the assumed rotational velocity of the inner core is very different between Janka & Mönchmeyer (1989) and Shimizu et al. (1994, 2001). Shimizu et al. (2001) concluded that the increase in the explosion energy due to anisotropic neutrino radiation occurs because cooling due to neutrino reemission is suppressed in anisotropic models. This is due to earlier shock revival and hence more efficient decrease of the matter temperature than those in spherical models. On the other hand, the neutrino heating itself is almost unchanged by the effect of anisotropic neutrino radiation. The neutrino heating dominates the cooling as a result, which increases the explosion energy, and leads to a successful explosion.

In Shimizu et al. (2001), the geometric effects of neutrino radiation have been rigorously treated outside the neutrinosphere for the first time, although its flux on the neutrinosphere was assumed. Only a global form of anisotropy was assumed there; the maximum peak in the neutrino flux distribution was located at the pole, and the minimum at the equatorial plane. However, Burrows et al. (1995) have suggested that the neutrino flux can fluctuate with angle and time. Such fluctuations are due to gravitational oscillation on the surface of the proto$``$neutron star and have a completely different origin from that of globally anisotropic neutrino radiation. Thus, it is interesting to investigate how the small-scale fluctuations affect the explosion mechanism and compare the results with those of the global anisotropy.

In this paper, we therefore introduce fluctuations in the neutrino flux in our numerical code by modifying the angular distribution of the neutrino flux. We aim to study the effects of these small-scale fluctuations on the shock position, the explosion energy, and the asymmetric explosion. Our numerical simulation is described in § 2, and the results are presented and discussed in § 3. Our conclusion is given in § 4.

## 2 NUMERICAL SIMULATION

Our simulation is performed by solving two-dimensional hydrodynamic equations in spherical coordinates. A generalized Roe’s method is employed to solve the hydrodynamic equations with general equations of state (EOSs). The details of our numerical technique, together with the EOS and the initial condition used, are described in the previous article (Shimizu et al., 2001). In our study, we have improved the numerical code of Shimizu et al. (2001); the cells in the $`\theta `$-direction were shifted by half of the cell size (Shimizu, 1995) in order to avoid a numerical error near the pole, although the error was not serious for the investigation of the explosion energy. The computational region is divided into 500 ($`r`$-direction) $`\times `$ 62 ($`\theta `$-direction) numerical cells.

In the present paper, the local neutrino flux is assumed to be

$$l_\nu (r,\theta )=\frac{7}{16}\sigma T_\nu ^4c_1\left(1+c_2\mathrm{cos}^2(n_\theta \theta )\right)\frac{1}{r^2},$$
(1)
where $`\sigma `$ is the Boltzmann constant, and $`T_\nu `$ is the temperature on the neutrinosphere. In equation (1), the parameter $`c_2`$ represents the magnitude of anisotropy, and $`n_\theta `$ the number of waves in the $`\theta `$-direction. The case of $`n_\theta =1`$ corresponds to the global anisotropy, namely, no fluctuation. We see in equation (1) that the neutrino fluxes in the $`x`$ (equatorial) and $`z`$ (polar) directions become $`l_xl_\nu (r,\theta =90^{})c_1`$ and $`l_zl_\nu (r,\theta =0^{})c_1(1+c_2),`$ respectively. The degree of anisotropy $`l_z/l_x`$ is then represented as

$$\frac{l_z}{l_x}=1+c_2.$$
(2)
Note that equation (2) is different from that defined by Shimizu et al. (2001), $`(l_z/l_x)_{\mathrm{Shimizu}}`$; for $`n_\theta =1`$ and sufficiently small $`c_2`$, we can relate them as $`c_2\left[(l_z/l_x)_{\mathrm{Shimizu}}^21\right]/2`$.

The value of $`c_1`$ is calculated from given $`c_2`$ and $`n_\theta `$ so as to adjust the total neutrino flux to that in the spherical model. The total neutrino luminosity is obtained by integrating equation (1) over the whole solid angle,

$$L_\nu =r^2l_\nu (r,\theta )𝑑\mathrm{\Omega }=\frac{7}{16}\sigma T_\nu ^4\mathrm{\hspace{0.33em}4}\pi c_1\left(1+c_2\frac{2n_\theta ^21}{4n_\theta ^21}\right),$$
(3)
which is equated to that of spherical explosion with the same $`T_\nu `$,

$$L_\nu ^{\mathrm{sp}}=\frac{7}{16}\sigma T_\nu ^4\mathrm{\hspace{0.33em}4}\pi R_{\mathrm{NS}}^2.$$
(4)
In the above, $`R_{\mathrm{NS}}`$ is the radius of a proto$``$neutron star and fixed to be 50 km. By comparing equations (3) with (4), we obtain

$$c_1=\frac{R_{\mathrm{NS}}^2}{1+c_2(2n_\theta ^21)/(4n_\theta ^21)}.$$
(5)
It should be noted here that the magnitude of fluctuations in the neutrino flux distribution for an observer far from the neutrinosphere (represented by $`c_2`$) and that on the neutrino-emitting surface (here we denote it as $`a`$) are different: the local neutrino flux is seen as equation (1) when we observe fluctuations on the surface of neutrino emission far from the neutrinosphere. It is preferable that we compare the results for the same value of $`a`$, since $`a`$ is more directly related to explosion dynamics. The value of $`c_2`$ should, therefore, be calculated from a given $`a`$, depending on $`n_\theta `$. Although it is difficult to calculate the exact relationship between $`c_2`$ and $`a`$, we can estimate it as follows. First, we assume that the strength of the neutrino flux on the neutrinosphere is represented by a profile of step functions:

$$\begin{array}{ccc}\text{for }n_\theta =1\hfill & & \\ & f(\theta )\hfill & \{\begin{array}{c}1+a\left(1\mathrm{cos}\theta >1/2\right)\hfill \\ 1\left(1/2>\mathrm{cos}\theta 0\right),\hfill \end{array}\hfill \\ \text{for }n_\theta =3\hfill & & \\ & f(\theta )\hfill & \{\begin{array}{c}1+a(1\mathrm{cos}\theta >3/4,\mathrm{\hspace{0.33em}1}/2>\mathrm{cos}\theta >1/4)\hfill \\ 1(3/4>\mathrm{cos}\theta >1/2,\mathrm{\hspace{0.33em}1}/4>\mathrm{cos}\theta 0)),\hfill \end{array}\hfill \\ \text{for }n_\theta =5\hfill & & \\ & f(\theta )\hfill & \{\begin{array}{c}1+a(1\mathrm{cos}\theta >5/6,\mathrm{\hspace{0.33em}2}/3>\mathrm{cos}\theta >1/2,\mathrm{\hspace{0.33em}1}/3>\mathrm{cos}\theta >1/6)\hfill \\ 1(5/6>\mathrm{cos}\theta >2/3,\mathrm{\hspace{0.33em}1}/2>\mathrm{cos}\theta >1/3,\mathrm{\hspace{0.33em}1}/6>\mathrm{cos}\theta 0).\hfill \end{array}\hfill \end{array}$$
(6)
We continue similarly for larger values of $`n_\theta `$. The parameter $`a`$ in equation (6) represents the magnitude of fluctuations in the neutrino flux on the neutrino emitting surface; the bright regions ($`f1+a`$) correspond to those where rising convective motion occurs, while the dark regions ($`f1`$) correspond to those of sinking convection. We have assumed that the areas of bright and dark regions are the same. Note that the neutrino flux function is always bright ($`f1+a`$) when $`\theta =0`$ (along the pole) and dark ($`f1`$) in the case of $`\theta =\pi /2`$ (on the equatorial plane), and that the fluctuations are essentially added to the global anisotropic model ($`n_\theta =1`$).

The neutrino flux observed far from the neutrinosphere is obtained by averaging all contributions from the flux on the surface of the neutrinosphere, $`l(\mathrm{\Theta })=𝑑\varphi 𝑑\theta f(\theta )\mathrm{sin}\theta \mathrm{cos}(\theta \mathrm{\Theta })`$. Here $`\mathrm{\Theta }`$ is the inclination angle between the line of sight of an observer and the polar axis of the proto$``$neutron star. This means that fully geometric effects from an anisotropically radiating surface are included in neutrino radiation field. For example, there is a contribution from fluxes along the pole axis to those on the equatorial plane, which will tend to reduce an efficiency of anisotropy. The ratio of the local neutrino flux along the polar axis ($`l_z`$) to that on the equatorial plane ($`l_x`$) for an observer far from the neutrinosphere is described as

$`{\displaystyle \frac{l_z}{l_x}}`$ $`=`$ $`{\displaystyle \frac{{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle _0^{\pi /2}}f(\theta )\mathrm{sin}\theta \mathrm{cos}\theta d\theta }{{\displaystyle _0^\pi }𝑑\varphi {\displaystyle _0^\pi }f(\theta )\mathrm{sin}\theta \mathrm{cos}\left(\theta \pi /2\right)𝑑\theta }}`$ (10)
$`=`$ $`\{\begin{array}{c}\frac{1+0.750a}{1+0.391a}\text{for }n_\theta =1,\hfill \\ \frac{1+0.625a}{1+0.438a}\text{for }n_\theta =3,\hfill \\ \frac{1+0.583a}{1+0.457a}\text{for }n_\theta =5,\mathrm{}\hfill \end{array}`$
By comparing equation (10) with equation (2), we finally obtain the value of $`c_2`$ for each value of $`n_\theta `$.

In the following, we examine two model series: $`a=0.31`$ (model series A) and $`a=0.71`$ (model series B). These values of $`a`$ are chosen in such a way that the value of $`l_z/l_x`$ for the global model ($`n_\theta =1`$) becomes 1.10 and 1.20, respectively. The values of $`c_2`$ for each fluctuation model ($`n_\theta =3,5`$) are accordingly calculated. These are summarized in Table 1. The neutrino temperature on the neutrino-emitting surface $`T_\nu `$ is assumed to be 4.65 and 4.70 MeV. Note that the value of 1+$`a`$, 1.71, for the model series B roughly corresponds to the variation in the neutrino flux obtained by Burrows et al. (1995), which is at a factor of 1.6.

## 3 RESULTS AND DISCUSSION

Figures 1 and 2 show the color-scale maps of the dimensionless entropy (Shimizu et al., 2001) distribution with the velocity fields for the model of global anisotropy ($`n_\theta =1`$; models A1-T470 and B1-T470). At $`t=82`$ ms after the shock stall, the shock front reaches $`r430`$ km on the equatorial plane and $`r530`$ km at the pole for the model A1-T470. The shock front is prolate, since the neutrino heating along the pole is more intensive than that on the equatorial plane, resulting in a jetlike explosion. At a later stage ($`t=244`$ ms), the shock wave is around a few thousand kilometers with large distortion. The entropy distribution of the model B1-T470 has a similar profile except that the shock front is more extended.

In Figures 3 and 4, the results for models with fluctuations (models A3-T470, A5-T470, B3-T470, and B5-T470) are depicted. When $`n_\theta =3`$, we find that the degree of asymmetry is smaller than that of the global anisotropy. In particular, the shock position for these models is less extended than the globally anisotropic one. This trend becomes more obvious for the model of $`c_2=+0.035`$ and $`n_\theta =5`$ (model A5-T470), where the shock front is almost spherical and its radius is only about 1300 km. In the case of the model B5-T470, the shock front is distorted because of a strong hydrodynamic flow along the pole, although this does not affect the explosion energy (see discussion on the energy later in this section).

The profile of the explosion energy is found to be closer to that of spherical explosion as the mode number of fluctuations increases. Figure 5 shows the evolution of the explosion energy, as well as the thermal, kinetic, and gravitational energies for the models of $`T_\nu =4.70`$ MeV. The difference between the globally anisotropic model and the models with fluctuations is prominent: the energy gain for the case of $`n_\theta =1`$ is the highest among others at all stages of the explosion. It is also seen that the explosion energy decreases as the mode number of fluctuations in the neutrino flux increases and finally approaches that of spherical explosion. For series B, we have obtained similar results. We observe that asymmetry in explosive motion is more enhanced than that for each model of series A.

In Figure 5, we compare the results of the explosion energy for the two model series. No significant difference is found between the two, although the result of the model A1-T470 becomes larger than that of the model B1-T470 at the later stages of the explosion. It has been shown (Shimizu et al., 2001) that the explosion energy increases as the degree of anisotropy becomes larger for not so large a degree of anisotropy, and finally saturates at $`(l_z/l_x)_{\mathrm{Shimizu}}1.2`$. Our new result shows that the final explosion energy of the more anisotropic model B1-T470 is smaller than that of the model A1-T470 (see Fig. 5). The difference may be attributed to the fact that the assumed forms of the local neutrino flux are different (compare eq.  here with eq.  in Shimizu et al. 2001). In the present paper, we have assumed the form of the neutrino fluxes that has more sharply concentrated flux on the pole. Therefore, the neutrino heating and rising convection are focused on the pole, and those in the equatorial direction are extremely reduced. The shock wave of the model B1-T470 appears to be too weak on the equatorial plane at $`t\begin{array}{c}>\hfill \\ \hfill \end{array}300`$ms, which causes a energy loss (Shimizu et al., 2001). Such features are clearly seen in the entropy distribution of the model B1-T470.

The effectiveness of global anisotropy becomes more pronounced as $`T_\nu `$ is decreased. Figure 6 shows the same energy evolution as Figure 5, except for $`T_\nu =4.65`$ MeV. The difference between the model with global anisotropy and those with fluctuations is extremely remarkable. The globally anisotropic models succeed, while all the other fluctuated and spherical models fail to explode except for the model B3-T465. Note that difference in the neutrino temperature is only 1%, which indicates a sensitivity of the supernova problem on the neutrino luminosity and energy. Note also that an increase of the explosion energies of the models with fluctuations at $`t500`$ ms is physically meaningless, because we have not taken into account a decay of the neutrino luminosity when $`t\begin{array}{c}>\hfill \\ \hfill \end{array}500`$ ms (e.g., Wilson & Mayle 1988; see also Shimizu et al. 2001).

We found that there are remarkable differences in the explosion energy depending on the mode of the fluctuations and that larger number of modes in the fluctuations makes the result closer to that of spherical explosion, irrespective of the model series A or B. Any small-scale fluctuations on the neutrinosphere are greatly averaged out when the neutrino emission is observed far enough from the neutrino-emitting surface. Moreover, we found that a certain broad space is needed to be heated by neutrinos to revive the stalled shock wave rigorously and that the global anisotropy ($`n_\theta =1`$) is the most effective to increase the explosion energy. Burrows et al. (1995) suggested that the neutrino flux can fluctuate not only with angle but with time. Such time fluctuations are expected to reduce further the efficiency of anisotropy, which needs to be confirmed in the future.

## 4 CONCLUSION

We have investigated the effects of small-scale fluctuations in the neutrino flux on the core-collapse supernova explosion. In order to examine the effect of the degree of anisotropy itself on the explosion, we have studied two model series parametrically. We specified the neutrino radiation field taking its geometric effects into account in each model. We found that the global anisotropy ($`n_\theta =1`$) and the local fluctuations ($`n_\theta >1`$) in the neutrino flux have quite different effects on the explosion mechanism, that is, the shock dynamics, the explosion energy, and the explosion asymmetry. Since the small-scale fluctuations are averaged out for radiative and hydrodynamic reasons, the results including fluctuations become closer to that of spherical explosion. Consequently, the global anisotropy is the most effective mechanism of increasing the explosion energy. Note here that the explosion energy could differ substantially depending on the neutrino temperature (the difference between 4.70 and 4.65 MeV is only 1%). This indicates that the supernova problem is very sensitive to the neutrino energy and luminosity. However, the total luminosity cannot be simply increased to explain the observed explosion energy because such treatment leads to the problem of Ni overproduction, especially in the case of essentially spherical models. We therefore conclude that globally anisotropic neutrino radiation is of great importance in actual supernova explosions. This supports the claim made by Shimizu et al. (2001).

The global anisotropy can originate from rotation of a proto$``$neutron star or a hot spot on the neutrino-emitting region, while the small-scale fluctuations are considered to be resulted from gravitational oscillation or uniform convection. It will be very interesting if any evidence of anisotropic neutrino radiation is observed at facilities like Super Kamiokande (Hirata et al., 1987; Suzuki, 1998) and SNO (Poon et al., 2001), together with detailed optical observations (e.g., Wang et al., 2002).

We are grateful to the anonymous referee for useful comments and suggestions that improved this paper.


text:
# On nonuniqueness for Calderón’s inverse problem

## 1 Introduction

An anisotropic conductivity on a domain $`\mathrm{\Omega }𝐑^n`$ is defined by a symmetric, positive semi-definite matrix-valued function, $`\sigma =(\sigma ^{ij}(x))`$. In the absence of sources or sinks, an electrical potential $`u`$ satisfies

$`(\sigma )u=_j\sigma ^{jk}(x)_ku`$ $`=`$ $`0\text{ in }\mathrm{\Omega },`$ (1)
$`u|_\mathrm{\Omega }`$ $`=`$ $`f,`$
where $`f`$ is the prescribed voltage on the boundary. Above and hereafter we use the Einstein summation convention where there is no danger of confusion. The resulting voltage-to-current (or Dirichlet-to-Neumann) map is then defined by

$`\mathrm{\Lambda }_\sigma (f)=Bu|_\mathrm{\Omega }`$ (2)
where

$`Bu=\nu _j\sigma ^{jk}_ku,`$ (3)
$`u`$ is the solution of (1) and $`\nu =(\nu _1,\mathrm{},\nu _n)`$ is the unit normal vector of $`\mathrm{\Omega }`$.

Applying the divergence theorem, we have

$`Q_\sigma (f)=:{\displaystyle _\mathrm{\Omega }}\sigma ^{jk}(x){\displaystyle \frac{u}{x^j}}{\displaystyle \frac{u}{x^k}}dx={\displaystyle _\mathrm{\Omega }}\mathrm{\Lambda }_\sigma (f)fdS,`$ (4)
where $`u`$ solves (1) and $`dS`$ denotes surface measure on $`\mathrm{\Omega }`$. $`Q_\sigma (f)`$ represents the power needed to maintain the potential $`f`$ on $`\mathrm{\Omega }`$. By (4), knowing $`Q_\sigma `$ is equivalent with knowing $`\mathrm{\Lambda }_\sigma `$. If $`F:\mathrm{\Omega }\mathrm{\Omega },F=(F^1,\mathrm{},F^n)`$, is a diffeomorphism with $`F|_\mathrm{\Omega }=\text{Identity}`$, then by making the change of variables $`y=F(x)`$ and setting $`u=vF^1`$ in the first integral in (4), we obtain

$`\mathrm{\Lambda }_{F_{}\sigma }=\mathrm{\Lambda }_\sigma ,`$
where

$`(F_{}\sigma )^{jk}(y)={\displaystyle \frac{1}{\text{det}[\frac{F^j}{x^k}(x)]}}{\displaystyle \underset{p,q=1}{\overset{n}{}}}{\displaystyle \frac{F^j}{x^p}}(x){\displaystyle \frac{F^k}{x^q}}(x)\sigma ^{pq}(x)|_{x=F^1(y)}`$ (5)
is the push-forward of the conductivity $`\sigma `$ by $`F`$. Thus, there is a large (infinite-dimensional) class of conductivities which give rise to the same electrical measurements at the boundary. The uniqueness question, first proposed by Calderón \[C\], we wish to address is whether two conductivities with the same Dirichlet-to-Neumann map must be such pushforwards of each other. By a direct construction, we will show that the answer is no. Let $`D\mathrm{\Omega }`$ be a smooth subdomain. We will construct in $`\mathrm{\Omega }`$ a conductivity $`\widehat{\sigma }`$ for which the boundary measurements coincide with those made for homogeneous conductivity $`\gamma =1`$ in $`\mathrm{\Omega }`$. We note that the conductivity is singular in the sense that some components of the conductivity tensor go to zero, i.e., correspond to perfectly insulating directions and some components can go to infinity, i.e., correspond to perfectly conducting directions, as one approaches $`D`$. The physical meaning of these counterexamples is evident: For instance in medical imaging, e.g., in Electrical Impedance Tomography (EIT) certain anisotropic structures can form barriers which, including their interior, appear in measurements to be a homogeneous medium.

To construct these counterexamples, we need to consider a variant of (1), which comes from the Laplace-Beltrami operator on a compact Riemannian manifold with boundary. Let us assume now that $`(M,g)`$ is an $`n`$-dimensional Riemannian manifold with smooth boundary $`M`$. The metric $`g`$ is assumed to be symmetric and positive definite. The invariant object analogous to the conductivity equation (1) is the Laplace-Beltrami operator, which is given by

$`\mathrm{\Delta }_gu=G^{1/2}_j(G^{1/2}g^{jk}_ku)`$ (6)
where $`G=\text{det}(g_{jk}),`$ $`[g_{jk}]=[g^{jk}]^1`$. The Dirichlet-to-Neumann map is defined by solving the Dirichlet problem

$`\mathrm{\Delta }_gu=0\text{in}M,`$ (7)
$`u|_M=f.`$
The operator analogous to $`\mathrm{\Lambda }_\sigma `$ is then

$`\mathrm{\Lambda }_g(f)=G^{1/2}\nu _jg^{jk}{\displaystyle \frac{u}{x_k}}|_M,`$ (8)
with $`\nu =(\nu _1,\mathrm{},\nu _n)`$ the outward unit normal to $`M`$. In dimension three or higher, the conductivity matrix and the Riemannian metric are related by

$`\sigma ^{jk}=\text{det}(g)^{1/2}g^{jk},\text{or}g^{jk}=\text{det}(\sigma )^{2/(n2)}\sigma ^{jk}.`$ (9)
Moreover, $`\mathrm{\Lambda }_g=\mathrm{\Lambda }_\sigma `$, and $`\mathrm{\Lambda }_{\psi ^{}g}=\mathrm{\Lambda }_g`$, where $`\psi ^{}g`$ denotes the pullback of the metric $`g`$ by a diffeomorphism of $`M`$ fixing $`M`$ \[LeU\].

In dimension two, (9) is not valid; in this case, the conductivity equation can be reformulated as

$`\text{Div}_g(\beta \text{Grad}_gu)=0\text{in}M,`$ (10)
$`u|_M=f`$
where $`\beta `$ is the scalar function $`\beta =|\text{det}\sigma |^{1/2}`$, $`g=(g_{jk})`$ is equal to $`(\sigma _{jk})`$, and $`\text{Div}_g`$ and $`\text{Grad}_g`$ are the divergence and gradient operators with respect to the Riemannian metric $`g`$. Thus we see that, in two dimensions, Laplace-Beltrami operators correspond only to those conductivity equations for which $`\text{det}(\sigma )=1`$.

For domains in two dimensions, Sylvester\[Sy\] showed, using isothermal coordinates, that one can reduce the anisotropic problem to the isotropic one; combining this with the isotropic result of Nachman\[Na\], one obtains

###### Theorem 1

If $`\sigma `$ and $`\stackrel{~}{\sigma }`$ are two $`C^3`$-smooth anisotropic conductivities in $`\mathrm{\Omega }𝐑^2`$ for which $`\mathrm{\Lambda }_\sigma =\mathrm{\Lambda }_{\stackrel{~}{\sigma }}`$, then there is a diffeomorphism $`F:\mathrm{\Omega }\mathrm{\Omega }`$, $`F|_\mathrm{\Omega }=Id`$ such that $`\stackrel{~}{\sigma }=F_{}\sigma .`$

In dimensions three and higher, the following result is known (see \[LU\], \[LTU\], and \[LeU\]):

###### Theorem 2

If $`n3`$ and $`(M,M)`$ is a $`C^\omega `$ manifold with connected, $`C^\omega `$ boundary, and $`g,\stackrel{~}{g}`$ are $`C^\omega `$ metrics on $`M`$ such that $`\mathrm{\Lambda }_g=\mathrm{\Lambda }_{\stackrel{~}{g}}`$, then there exists a $`C^\omega `$ diffeomorphism $`F:MM`$ such that $`F|_D=Id`$.

## 2 Counterexamples

Returning now to domains $`\mathrm{\Omega }𝐑^n,n3`$, let $`D\mathrm{\Omega }`$ be an open subset with smooth boundary and $`g=g_{ij}`$ be a metric on $`\mathrm{\Omega }`$. Let $`yD`$ be such that there is a diffeomorphism $`F:\mathrm{\Omega }\{y\}\mathrm{\Omega }\overline{D}`$, and let $`\stackrel{~}{g}=F_{}g`$ on $`\mathrm{\Omega }\overline{D}`$. To obtain a conductivity on all of $`\mathrm{\Omega }`$, first extend $`\stackrel{~}{g}`$ to a bounded metric inside $`D`$ and denote this new metric on $`\mathrm{\Omega }`$ by $`\widehat{g}`$. We make this continuation so that the conductivity jumps on $`D`$ and that $`\widehat{g}c>0`$ is smooth inside $`D`$, but is otherwise arbitrary. Let $`\widehat{\sigma }`$ be the conductivity corresponding to $`\widehat{g}`$ by (9) in $`\mathrm{\Omega }`$. We say that $`v`$ is a solution of the conductivity equation if

$`\widehat{\sigma }v(x)=0\text{ in the sense of distributions in }\mathrm{\Omega },`$ (11)
$`v|_\mathrm{\Omega }=f_0,`$
$`vL^{\mathrm{}}(\mathrm{\Omega })`$
where $`f_0H^{\frac{1}{2}}(\mathrm{\Omega })`$. That the equation in the sense of distributions means that $`vH^1(\mathrm{\Omega })`$ and $`\sigma vH(\mathrm{\Omega };)=\{wL^2(\mathrm{\Omega }):wL^2(\mathrm{\Omega })\}`$. If this problem has a unique solution, we define the Dirichlet-to-Neumann map

$`\mathrm{\Lambda }_{\widehat{\sigma }}f_0=n\widehat{\sigma }v|_\mathrm{\Omega }.`$
Our aim is to prove the following result, first announced in \[GLU\]:

###### Theorem 3

Let $`\mathrm{\Omega }𝐑^n`$, $`n3`$, and $`g=g_{ij}`$ a metric on $`\mathrm{\Omega }`$. Let $`D\mathrm{\Omega }`$ be such there is a $`C^{\mathrm{}}`$-diffeomorphism $`F:\mathrm{\Omega }\{y\}\mathrm{\Omega }\overline{D}`$ satisfying $`F|_\mathrm{\Omega }=Id`$ and that

$`dF(x)c_0I,\text{det}(dF(x))c_1\text{dist}_{_{𝐑^n}}(x,y)^1`$ (12)
where $`dF`$ is the Jacobian matrix in Euclidean coordinates of $`𝐑^n`$ and $`c_0,c_1>0`$. Let $`\stackrel{~}{g}=F_{}g`$ and $`\widehat{g}`$ be an extension of $`\stackrel{~}{g}`$ into $`D`$ such that it is positive definite in $`D^{int}`$. Finally, let $`\gamma `$ and $`\widehat{\sigma }`$ be the conductivities corresponding to $`g`$ and $`\widehat{g}`$ by (9). Then the boundary value problem for the conductivity equation with conductivity $`\widehat{\sigma }`$ is uniquely solvable and

$`\mathrm{\Lambda }_{\widehat{\sigma }}=\mathrm{\Lambda }_\gamma .`$
Note that here is no diffeomorphism $`H:\mathrm{\Omega }\mathrm{\Omega }`$ such that $`\widehat{\sigma }=H_{}\gamma `$, so the Riemannian manifolds corresponding to $`\widehat{\sigma }`$ and $`\gamma `$ cannot be the same. Also, $`\widehat{\sigma }`$ can be changed in arbitrary way inside $`D`$ without changing boundary measurements. The proof of Theorem 3 will be given below.

Example. Let $`\mathrm{\Omega }=B(0,2)𝐑^3`$ be the ball with center 0 and radius 2. Consider $`y=0D=B(0,1)`$ and the map $`F:\mathrm{\Omega }\{0\}\mathrm{\Omega }\overline{D}`$ given by

$`F:x({\displaystyle \frac{|x|}{2}}+1){\displaystyle \frac{x}{|x|}}.`$ (13)
Let $`\gamma =1`$ be the homogeneous conductivity in $`\mathrm{\Omega }`$ and define $`\sigma =F_{}\gamma `$. Now the metric tensor $`g`$ and the corresponding conductivity $`\sigma _g`$ are related by $`\sigma _g=|\text{det}g|^{1/2}g^{jk}`$. Let $`g`$ be the metric corresponding to $`\gamma `$ and $`\stackrel{~}{g}`$ be the metric corresponding to $`\sigma `$. Consider these in the standard spherical coordinates on $`\mathrm{\Omega }\{0\}`$, $`(r,\varphi ,\theta )(r\mathrm{sin}\theta \mathrm{cos}\varphi ,r\mathrm{sin}\theta \mathrm{sin}\varphi ,r\mathrm{cos}\theta )𝐑^3`$. With respect to these coordinates, we see that the metric $`g`$ and conductivity $`\gamma `$ correspond to the matrices

$`g=\left(\begin{array}{ccc}1& 0& 0\\ 0& r^2& 0\\ 0& 0& r^2\mathrm{sin}^2\theta \end{array}\right),\gamma =\left(\begin{array}{ccc}r^2\mathrm{sin}\theta & 0& 0\\ 0& \mathrm{sin}\theta & 0\\ 0& 0& (\mathrm{sin}\theta )^1\end{array}\right)`$
and $`\stackrel{~}{g}`$ and $`\sigma `$ correspond in the annulus $`1<r<2`$ to the matrices

$`\stackrel{~}{g}=\left(\begin{array}{ccc}4& 0& 0\\ 0& 4(r1)^2& 0\\ 0& 0& 4(r1)^2\mathrm{sin}^2\theta \end{array}\right),`$
$`\sigma =\left(\begin{array}{ccc}2(r1)^2\mathrm{sin}\theta & 0& 0\\ 0& 2\mathrm{sin}\theta & 0\\ 0& 0& 2(\mathrm{sin}\theta )^1\end{array}\right).`$
Let $`\widehat{\sigma }`$ be a continuation of $`\sigma `$ that is $`C^{\mathrm{}}`$-smooth in $`D`$. Then these metric are as in for Theorem 3 and in particular $`F`$ satisfies (12).

To start to prove Theorem 3, we start with the following result:

###### Proposition 1

Let $`\mathrm{\Omega }𝐑^n`$, $`n3`$, and $`g=g_{ij}`$ a metric on $`\mathrm{\Omega }`$. Let $`u`$ satisfy

$`\mathrm{\Delta }_gu(x)=0\text{in }\mathrm{\Omega },`$
$`u|_\mathrm{\Omega }=f_0C^{\mathrm{}}(\mathrm{\Omega }).`$
Let $`D\mathrm{\Omega }`$ be such that there is a diffeomorphism $`F:\mathrm{\Omega }\{y\}\mathrm{\Omega }\overline{D}`$ satisfying $`F|_\mathrm{\Omega }=Id`$. Let $`\stackrel{~}{g}=F_{}g`$ and $`v`$ be a function satisfying

$`\mathrm{\Delta }_{\stackrel{~}{g}}v(x)=0\text{in }\mathrm{\Omega }\overline{D},`$
$`u|_\mathrm{\Omega }=f_0,`$
$`uL^{\mathrm{}}(\mathrm{\Omega }\overline{D}).`$
Then $`u`$ and $`F^{}v`$ coincide and have the same Cauchy data on $`\mathrm{\Omega }`$,

$`_\nu u|_M=_{\stackrel{~}{\nu }}F^{}v|_M`$ (17)
where $`\nu `$ is unit normal vector in metric $`g`$ and $`\stackrel{~}{\nu }`$ is unit normal vector in metric $`\stackrel{~}{g}`$. Moreover, for the constant $`c_0:=u(y)`$ we have

$`\underset{xD}{lim}v(x)=c_0.`$
Proof. Let $`g=(g_{ij})`$ be a Riemannian metric tensor defined on $`\mathrm{\Omega }𝐑^n`$.

First, we continue the $`C^{\mathrm{}}`$-metric $`g`$ to a metric $`𝔤`$ in $`𝐑^n`$ such that $`𝔤_{ij}(x)=g_{ij}(x)`$ for $`x\mathrm{\Omega }`$ is such that for any $`y𝐑^m`$ there is a positive Green’s function $`G(x,y)`$ satisfying

$`\mathrm{\Delta }_𝔤G(,y)=\delta _y\text{in }𝐑^n.`$
There are several easy ways to obtain this continuation. For instance, we can continue $`g`$ to a metric $`𝔤`$ such that outside some ball $`B(0,R)`$ the metric is hyperbolic. This implies that the manifold $`(𝐑^n,𝔤)`$ is non-parabolic and has a positive non-constant super-harmonic function (see \[Gr\]). By \[LiT\] there thus exists a positive Green’s function $`G(x,y)`$.

Next we consider the probability that Brownian motion $`B_t^x`$ on manifold $`(𝐑^n,𝔤)`$ sent from the point $`x`$ at time $`t=0`$ enters an open set $`U`$: Let

$`e_U(x)=P(\{\text{there is }t>0\text{ such that }B_t^xU\})`$
Let $`A=𝐑^nU`$. By Hunt’s theorem (see, \[H\], or \[Gr\], Prop. 4.4) we have

$`e_U(x)=s_A(x)`$
where $`s_A(x)`$ is the super-harmonic potential of $`A`$, that is, $`s_A(x)`$ is infimum of all bounded super-harmonic functions $`h`$ in $`𝐑^n`$ such that $`h|_U=1`$ and $`h0`$ in $`𝐑^n`$.

Let $`U=B(y,r)`$ and

$`m(y,r)=\underset{xB(y,r)}{inf}G(x,y).`$
The function

$`h_{y,r}(x)={\displaystyle \frac{1}{m(y,r)}}\mathrm{min}(G(x,y),m(y,r))`$
is positive super-harmonic function and satisfies $`h_{y,r}|_{B(y,r)}=1`$. Since $`e_U(x)h_{y,r}(x)`$ and $`lim_{r0}m(y,r)=\mathrm{}`$, we see that

$`\underset{r0}{lim}P(\{\text{there is }t>0\text{ such that }B_t^xB(y,r)\})=0.`$ (18)
In particular, taking limit $`r0`$ we see that the probability that $`B_t^x=y`$ for some $`t>0`$ is zero.

Now, let $`\stackrel{~}{u}`$ be any solution of

$`\mathrm{\Delta }_g\stackrel{~}{u}(x)=0\text{in }\mathrm{\Omega }\{y\}`$ (19)
$`\stackrel{~}{u}|_\mathrm{\Omega }=f_0,`$
$`\stackrel{~}{u}L^{\mathrm{}}(\mathrm{\Omega }\{y\}).`$
Denote $`f_r=\stackrel{~}{u}|_{B(y,r)}`$ and let $`\tau (r,x)(0,\mathrm{}]`$ be the first time when $`B_t^x\mathrm{\Omega }B(y,r)`$. Similarly, let $`\tau (0,x)`$ be the hitting time to $`\mathrm{\Omega }`$.

Let next $`\chi _{0,r}^x`$ be a random variable defined by $`\chi _{0,r}^x=1`$ if $`B_{\tau (r,x)}^x\mathrm{\Omega }`$ and zero otherwise, and let $`\chi _{1,r}^x=1\chi _{0,r}^x`$. Then by Kakutani’s formula,

$`\stackrel{~}{u}(x)=𝔼(\chi _{0,r}^xf_0(B_{\tau (r,x)}^x))+𝔼(\chi _{1,r}^xf_r(B_{\tau (r,x)}^x)).`$
Letting $`r0`$ and using the fact that $`f_r_{\mathrm{}}\stackrel{~}{u}_{\mathrm{}}`$ are uniformly bounded and (18) we see $`\stackrel{~}{u}(x)=𝔼(f_0(B_{\tau (0,x)}^x))`$. Thus $`\stackrel{~}{u}(x)=u(x)`$ for $`x\mathrm{\Omega }\{y\}`$ where

$`\mathrm{\Delta }_gu(x)=0\text{in }\mathrm{\Omega },`$ (20)
$`u|_\mathrm{\Omega }=f_0.`$
Thus we have shown that boundary value problem (19) is solvable, and that the solution is unique. Next we change this problem to an equivalent problem. Let $`F:\mathrm{\Omega }\{y\}\mathrm{\Omega }\overline{D}`$ be a diffeomorphism that is identity at boundary $`\mathrm{\Omega }`$. Define the metric $`\stackrel{~}{g}=F_{}g`$; then the boundary value problem (20) is equivalent to

$`\mathrm{\Delta }_{\stackrel{~}{g}}v(x)=0\text{in }\mathrm{\Omega }\overline{D}`$ (21)
$`v|_\mathrm{\Omega }=f_0,`$
$`vL^{\mathrm{}}(\mathrm{\Omega }\overline{D}).`$
and solutions of the problems (20) and (21) are related by $`v(x)=u(F(x))`$ for $`x\mathrm{\Omega }\{y\}`$.

Clearly, the Cauchy data of the equations (20) and (21) coincide in the sense of (17). Moreover,

$`\underset{xD}{lim}v(x)=c_0:=u(y)`$
This concludes the proof of Prop. 1. $`\mathrm{}`$

We remark that the application of Brownian motion above is not essential but makes the proof perhaps more intuitive. Alternatively, in the proof of Prop. 1 one can use properties of $`L^p`$-Sobolev spaces. Indeed, for $`\stackrel{~}{u}`$ satisfying $`\mathrm{\Delta }_g\stackrel{~}{u}(x)=0`$ in $`\mathrm{\Omega }\{y\}`$, $`\stackrel{~}{u}|_\mathrm{\Omega }=f_0`$, and $`\stackrel{~}{u}L^{\mathrm{}}(\mathrm{\Omega }\{y\})`$, we can consider the extension $`\stackrel{~}{u}L^{\mathrm{}}(\mathrm{\Omega })`$ that satisfies

$`\mathrm{\Delta }_g\stackrel{~}{u}(x)=F\text{in }\mathrm{\Omega },`$
where $`F`$ is a distribution is supported in $`y`$. Then $`F`$ has to be a finite sum of derivatives of the Dirac delta distribution supported at $`y`$. Now, $`FW^{2,p}(\mathrm{\Omega })`$ for all $`1<p<\mathrm{}`$, and since $`\delta _yW^{2,p}(\mathrm{\Omega })`$ for $`p>\frac{n}{n2}`$ and the same is true for the derivatives of the delta distribution, we see that that $`F=0`$. This implies that $`\stackrel{~}{u}=u`$.

We now turn to the proof of Theorem 3. From now on, we assume that $`F:\mathrm{\Omega }\{y\}\mathrm{\Omega }\overline{D}`$, $`F(x)=(F^1(x),\mathrm{},F^n(x))`$ is such that condition (12) is satisfied. First, continue $`\stackrel{~}{g}`$ into $`D`$ so that that it is positive definite in $`D^{int}`$. Next, extend $`v`$ inside $`D`$ to a function

$`h(x)=v(x)\text{ for }x\mathrm{\Omega }\overline{D},`$
$`h(x)=c_0\text{ for }x\overline{D}.`$
Our aim is to show that $`h`$ is the solution of (11). Since any solution of (21) is constant on the boundary of $`D`$ and $`\stackrel{~}{g}`$ is positive definite in $`D`$, we see that the solution has to be constant inside $`D`$. Thus $`h`$ is the unique solution if it is a solution.

Now we are ready to show that $`v`$ is a solution also in the sense of distributions. First we note that when $`y(x)=F^1(x)`$ we have $`v(x)=u(F^1(x))`$ and

$`{\displaystyle \frac{v}{x^j}}(x)={\displaystyle \frac{u}{y^k}}(y(x)){\displaystyle \frac{y^k}{x^j}}(x)`$
and since $`uH^1(\mathrm{\Omega })`$ and $`\frac{y^k}{x^j}L^{\mathrm{}}(\mathrm{\Omega }\overline{D})`$ we have $`vH^1(\mathrm{\Omega }\overline{D})`$. Also, as $`hC(\mathrm{\Omega }D^{int})`$ and the trace $`vv|_D`$ is a continuous map $`C(\mathrm{\Omega }D^{int})H^1(\mathrm{\Omega }\overline{D})L^2(D)`$, we see that trace $`v|_D`$ is well defined and is the constant function having value $`c_0`$. Since restrictions of $`h`$ to $`D`$ and $`\mathrm{\Omega }\overline{D}`$ are in $`H^1(D)`$ and $`H^1(\mathrm{\Omega }\overline{D})`$, correspondingly, and the trace from both sides of $`D`$ coincide, we see that $`hH^1(\mathrm{\Omega })`$.

Next we show that $`\sigma vL^2(\mathrm{\Omega }\overline{D})`$. Let $`e_1,e_2,\mathrm{},e_n`$ be the standard Euclidean coordinate vectors in $`𝐑^n`$. Then for $`x\mathrm{\Omega }\overline{D}`$

$`\stackrel{~}{g}^{jk}_kv(x)=(e_j,_{\stackrel{~}{g}}v(x))_{\stackrel{~}{g}}=(d(F^1)e_j,_gu(F^1(x)))_g`$
and we see from (12) that this inner product is uniformly bounded. By (5) and (12), $`|\text{det}(\stackrel{~}{g})(x)|c_3\text{dist}_{𝐑^n}(F^1(x),y)^2`$ and thus in $`\mathrm{\Omega }\overline{D}`$ the functions $`V_j(x)=|\text{det}\stackrel{~}{g}(x)|^{1/2}\delta _{jk}\stackrel{~}{g}^{ki}_iv`$ satisfy

$`|V_j(x)|c_4\text{dist}_{𝐑^n}(F^1(x),y).`$
This implies that $`\stackrel{}{V}=(V_1(x),\mathrm{},V_n(x))`$ is in $`H(\mathrm{\Omega }\overline{D};)C(\overline{\mathrm{\Omega }}D^{int})`$. Moreover, by \[Gr\] the normal trace $`\stackrel{}{W}n\stackrel{}{W}|_D`$ is a continuous map $`H(\mathrm{\Omega }\overline{D};)C(\overline{\mathrm{\Omega }}D^{int})H^{1/2}(D)`$, we see that $`n\stackrel{}{V}|_D=0`$. Since the normal traces of $`\stackrel{}{V}`$ coincide from both sides of $`D`$ we see that $`\stackrel{}{V}H(\mathrm{\Omega };)`$ and

$`{\displaystyle \underset{j,k=1}{\overset{n}{}}}_k(|\text{det}\stackrel{~}{g}(x)|^{1/2}\stackrel{~}{g}^{ki}_ih)=\stackrel{}{V}=0\text{in }\mathrm{\Omega }`$
in the sense of distributions. This means that $`h`$ satisfies the conductivity equation in $`\mathrm{\Omega }`$ in the sense of distributions. This proves Theorem 3.

A similar approach can be used to construct a different type of counterexample. In contrast to the previous one, here the conductivity is not bounded above near the singular surface. Consider again the sets $`\mathrm{\Omega }=B(0,2)`$ and $`D=B(0,1)`$ in $`𝐑^3`$. In spherical coordinates we define the metric $`\widehat{g}`$ and the conductivity $`\widehat{\sigma }`$ in the domain $`\{1<r<2\}`$ by the matrices

$`\widehat{g}=\left(\begin{array}{ccc}(r1)^2& 0& 0\\ 0& \rho ^2& 0\\ 0& 0& \rho ^2\mathrm{sin}^2\theta \end{array}\right),`$
$`\widehat{\sigma }=\left(\begin{array}{ccc}(r1)\rho ^2\mathrm{sin}\theta & 0& 0\\ 0& (r1)^1\mathrm{sin}\theta & 0\\ 0& 0& (r1)^1\mathrm{sin}^1\theta \end{array}\right)`$
where $`\rho >0`$ is a constant. Thus, $`(\mathrm{\Omega }\overline{D},\widehat{g})`$ is isometric to the product $`𝐑_+\times S_\rho ^2`$ with the standard metric where $`S_\rho ^2`$ is a 2-sphere with radius $`\rho `$. Now extend the conductivity $`\widehat{\sigma }`$ to $`\mathrm{\Omega }`$ so that it is positive definite in $`D`$. It can be shown that, in the domain $`\mathrm{\Omega }\overline{D}`$, the equation

$`\widehat{\sigma }v(x)=0\text{ in }\mathrm{\Omega }\overline{D}`$
$`v|_\mathrm{\Omega }=f_0,`$
$`vL^{\mathrm{}}(\mathrm{\Omega }\overline{D})`$
has a unique solution. Also, when $`\rho `$ is small enough, we can extend the definition of $`v(x)`$ to the whole domain $`\mathrm{\Omega }`$ by defining it to have the constant value $`c_0=lim_{xD}v(x)`$ in $`D`$. The function $`v(x)`$ thus obtained is a solution of the boundary value problem

$`\widehat{\sigma }v(x)=0\text{ in }\mathrm{\Omega },`$
$`v|_\mathrm{\Omega }=f_0,`$
$`vL^{\mathrm{}}(\mathrm{\Omega })`$
in the sense of distributions.

In this case, the boundary measurements do not give information about the metric inside $`D`$. To see this, consider $`(\mathrm{\Omega }\overline{D},\widehat{g})`$ as $`𝐑_+\times S_\rho ^2`$ with coordinates $`(t,\varphi ,\theta )`$, where $`t=t(r)=\mathrm{log}(r1)`$, i.e., $`t𝐑_+`$ and $`\varphi ,\theta `$ are coordinates on $`S^2`$. We see that for sets $`B(0,1+R^1)\mathrm{\Omega }`$, $`R>1`$ in original coordinates we can use the superharmonic potentials

$`h(r,\varphi ,\theta )=\mathrm{min}(1,{\displaystyle \frac{\mathrm{log}(r1)}{\mathrm{log}(R1)}})=\mathrm{min}(1,{\displaystyle \frac{t(r)}{t(R)}})`$
to see that the Brownian motion sent from $`x\mathrm{\Omega }\overline{D}`$ does not enter $`D`$ a.s.. By writing the solution $`u`$ in terms of spherical harmonics, we have

$`u(t,\varphi ,\theta )={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}a_{n,m}e^{\lambda _nt}Y_m^n(\varphi ,\theta )`$
where $`\lambda _n=\rho ^1\sqrt{n(n+1)}`$. Using the original coordinates $`(r,\varphi ,\theta )`$ of $`\mathrm{\Omega }\overline{D}`$ we see

$`u(r,\varphi ,\theta )={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}a_{n,m}(r1)^{\lambda _n}Y_m^n(\varphi ,\theta ).`$
Thus, if $`\rho `$ is small enough, the solution goes to constant and its derivative to zero when $`r1^+`$ so fast that the solutions of the conductivity equation satisfy the equation in sense of distributions.


text:
# A schematic model for QCD. I: Low energy meson states.

## 1 Introduction

QCD is the favored theory of the strong interactions. At low energy, however, the description of the hadronic spectrum based on QCD becomes difficult due to the non linear structure of the theory. This non-perturbative regimen, in contrast to lattice gauge calculations, can be explored by means of schematic models. The use of such models is common to other fields of physics where the many-body structure of the theory may be explored by introducing effective degrees of freedom and their couplings. Nuclear structure physics is one of these examples and it is as complicated and involved as the low-energy domain of QCD. Like in the case of nuclear structure physics, QCD descriptions based on simple models may help in the understanding of basic concepts and procedures. The Lipkin model  is one of the most famous schematic models, and it helped substantially to appreciate the importance of pairing two-body interactions as well as the importance of collectivity in building the low-energy part of the nuclear spectrum. An extended version of the Lipkin model was applied to the description of pion condensates in nuclei . A variety of many body techniques have been tested with Lipkin-type models . In  some realistic, less schematic, nuclear interactions, suitable to describe various nuclear properties, were investigated in this way. In Ref.  a Lipkin model was applied to describe a system of many quarks. As seen in these examples, the predictions of schematic models can be also rather rich in their complexities. This fact was shown, for a simple model of many gluon systems, in Ref. .

Until now, the only formalism which can handle QCD from first principles is the lattice gauge theory . Particularly, in many gluon systems, a good description is obtained without considering finite volume effects . The problem with lattice gauge calculations, to treat QCD at low energies, is that only the lowest state, and in some cases also the next to the lowest state, for a given spin, charge conjugation and parity, can be calculated. Lattice theory is numerically quite involved, and the inclusion of quarks and antiquarks brings in additional problems which are difficult to solve. Effective models of the hadrons, like the MIT model , can help to shed some light into the structure of QCD at low energy. In  the spectrum of gluons, as obtained in , was reproduced and the sequence of levels explained by simple assumptions. In other works  many body methods were used to describe the spectrum of QCD at low energy. After these considerations, it is obvious that it would be nice to have a model which: i) must be able to describe the basic structure of QCD at low and high energy, and ii) must be solbable exactly. Probably, such a model does not exist, due to the complicated structure of QCD. Nevertheless, one can try to construct a model which comes as near as possible to QCD.

The purpose of the present work is to present a model which fulfills the above requirements. In Ref.  the most simple version of such a model was presented. Like the models mentioned at the beginning, it is based on a Lipkin type model and it consists of two levels for the description of quarks and antiquarks. These quarks are coupled to a boson level which describes gluon pairs coupled to spin zero. The other gluon states are treated as spectators. The basic ingredients of the model are the quark-antiquark pairs coupled to flavor singlet and spin zero and gluon pairs with spin zero. In Ref.  it is shown that the model is able to describe the appearance of a quantum phase transition at zero temperature, when the interaction is turned on, and a phase transition to the non interacting case at non-vanishing temperature. In Ref.  the basic features of the model were discussed. Only flavor singlet and spin zero mesons were taken into account. The appearance of a Goldstone boson was obtained for a sufficiently strong interaction. This state consisted of a meson with negative parity. Also, it possesses a very collective nature, i.e. it is a superposition of many particles (quarks, antiquarks and gluons) states. The behavior of the model at high temperature was discussed, together with some consequences for the Quark-Gluon Plasma (QGP) in .

In the first part of the paper we will introduce the general form of the model for the description of meson states. The discussion will concentrate on the behavior of the model at low energy, corresponding to the zero temperature regime of QCD. The study of the high energy behavior and the transition to the Quark Gluon Plasma  will be presented in the forthcoming paper of the series . In section 2 the basic ingredients of the model are introduced, with the proposition of a Hamilton operator. Because of the difficulties to treat fermion pairs exactly, we shall diagonalize them in a boson mapping scheme. The basis used to deal with the bosonic images of the fermion pairs, and the corresponding matrix elements of the proposed Hamiltonian, are given in the same section 2. There we show how to assign charge conjugation and G-parity symmetries to the states belonging to the basis. In section 3 the model is applied to the description of the low energy meson spectrum. Conclusions are drawn in section 4.

## 2 The Model

As indicated in Ref.  the fermion sector is described by a Lipkin type model , consisting of two levels, one at energy $`\omega _f`$ and the other at energy $`+\omega _f`$ (see Fig. 1). This is the Dirac picture for fermions, where antiquarks are regarded as holes in the lower level. Alternatively, one also describes quarks and antiquarks in the level at $`+\omega _f`$. The quarks and antiquarks are coupled to a 1.6 GeV level, which is occupied by gluon pairs with spin zero , a value which was obtained in the description of a many gluon system and which was adjusted to lattice gauge calculations . In consequence, we shall take the energy of the gluon pair state as an externally fixed data. The value $`\omega _f`$ is fixed to one third of the nucleon mass (0.33GeV). There are further gluon states  which do not interact with the quarks and antiquarks. These states will be treated as spectators and should be taken into account in the final spectrum.

The degeneracy of each fermion level is 2$`\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ refers to color, flavor and eventually other degrees of freedom. If only spin ($`n_s`$), flavor ($`n_f`$) and color ($`n_c`$) degrees of freedom are considered, the value of $`\mathrm{\Omega }`$ is given by the product $`2\mathrm{\Omega }`$ $`=`$ $`n_cn_fn_s`$. For flavor (0,0) and spin 0 pairs, only, the model has similarities to the one of Ref. , but with a different interaction. In Ref.  the pion condensate in nuclei was the dominant phenomena. In the present case, nucleons are replaced by quarks and the pions by gluons. The model has some similarities to Ref. , which is also a Lipkin-type model. There, only quarks were considered and the interaction conserves their number.

For zero temperature and no interactions the lower level is filled by fermions. The creation (annihilation) operators of these fermions are $`𝒄_{\alpha (1,0)f\sigma i}^{}`$ ($`𝒄^{\alpha (1,0)f\sigma i}`$), in co- and contra-variant notation for the indices. The symbol $`(1,0)f`$ refers to the flavor part, where $`(1,0)`$ is the SU(3)-flavor notation and $`f`$ is a short hand notation for the hypercharge $`Y`$, the isospin $`T`$ and its third component $`T_z`$. The index $`\sigma `$ represents the two spin components $`\pm \frac{1}{2}`$, the index $`i`$ $`=`$ 1 or 2, stands for the upper or lower level and the index $`\alpha `$ represents all remaining degrees of freedom, which are at least 3 because of the color degree of freedom (when only color is taken into account, instead of $`\alpha `$ we will use the index $`c`$). Lowering and raising the indices of the operators introduces a phase, which depends on the convention used , and a change of the indices to their conjugate values, i.e., the quantum numbers $`(1,0)YTT_z\sigma `$ change to $`(0,1)YTT_z\sigma `$.

The operators, defined above, contain the relevant degrees of freedom of QCD, i.e. color, spin and flavor. These basic degrees of freedom appear at all energies, no matter how the resulting particles are defined, i.e., either in the perturbative or in the non-perturbative regime. In the non-perturbative regime one usually denotes them as constituent or effective particles. This is mainly due to the difference in the spatial properties, while color, spin and flavor have the same meaning as in QCD. Here, the quarks and antiquarks are constituent particles at low energy and have little in common (except for the quantum numbers mentioned) with the ones at high energy. We shall show that a model which contains these basic degrees of freedom and which takes into account the dynamic coupling with gluons can describe the main characteristics of QCD at low energy.

The quark and antiquark creation and annihilation operators are given in terms of the operators $`𝒄`$ and $`𝒄^{}`$

$`𝒂_{\alpha f\sigma }^{}`$ $`=`$ $`𝒄_{\alpha f\sigma 1}^{},𝒅_{\alpha f\sigma }=𝒄_{\alpha f\sigma 2}^{}`$
$`𝒂^{\alpha f\sigma }`$ $`=`$ $`𝒄^{\alpha f\sigma 1},𝒅^{\alpha f\sigma }=𝒄^{\alpha f\sigma 2},`$ (1)
which corresponds to the Dirac picture of particles and antiparticles: quarks are described by fermions in the upper level and antiquarks by holes in the lower level.

The gluon sector of the model space is described by bosons which represent pairs of gluons coupled to spin zero. The energy of a boson state is fixed at the value $`\omega _b=1.6`$ GeV , as mentioned before.

The quark-antiquark pairs of the model are given by

$`𝑪_{f_1\sigma _11}^{f_2\sigma _22}`$ $`=`$ $`𝑩_{f_1\sigma _1}^{f_2\sigma _2}={\displaystyle \underset{\alpha }{}}𝒄_{\alpha f_1\sigma _11}^{}𝒄^{\alpha f_2\sigma _22}={\displaystyle \underset{\alpha }{}}𝒂_{\alpha f_1\sigma _1}^{}𝒅^{\alpha f_2\sigma _2}`$
$`𝑪_{f_1\sigma _12}^{f_2\sigma _21}`$ $`=`$ $`𝑩_{f_1\sigma _1}^{f_2\sigma _2}={\displaystyle \underset{\alpha }{}}𝒄_{\alpha f_1\sigma _12}^{}𝒄^{\alpha f_2\sigma _21}={\displaystyle \underset{\alpha }{}}𝒅_{\alpha f_1\sigma _1}𝒂^{\alpha f_2\sigma _2}`$
$`𝑪_{f_1\sigma _11}^{f_2\sigma _21}`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}𝒄_{\alpha f_1\sigma _11}^{}𝒄^{\alpha f_2\sigma _21}={\displaystyle \underset{\alpha }{}}𝒂_{\alpha f_1\sigma _1}^{}𝒂^{\alpha f_2\sigma _2}`$
$`𝑪_{f_1\sigma _12}^{f_2\sigma _22}`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}𝒄_{\alpha f_1\sigma _12}^{}𝒄^{\alpha f_2\sigma _22}={\displaystyle \underset{\alpha }{}}𝒅_{\alpha f_1\sigma _1}𝒅^{\alpha f_2\sigma _2}.`$ (2)
The first two equations describe the creation and annihilation of quark-antiquark pairs. The pairs can be coupled to definite flavor $`(\lambda ,\lambda )`$ $`=`$ $`(0,0)`$ or $`(1,1)`$ and spin $`S`$ $`=`$ 0 or 1. We shall write, in this coupling scheme, $`𝑩_{(\lambda ,\lambda )f,SM}^{}`$, where $`f`$ is the flavor, $`S`$ is the spin and $`M`$ is the spin-projection. The operators $`𝑩_{(\lambda ,\lambda )f,SM}`$ annihilate the vacuum $`|0>`$, which can be taken as the configuration where the lower state is completely filled and the upper one is empty. Note, that the vacuum state is not uniquely defined . All states, which contain only quarks in the upper level and where the lower level is completely filled (so that antiquarks are not activated), regardless of color, as for example the three quark baryon states, are annihilated by $`𝑩_{(\lambda ,\lambda )f,SM}`$. This property derives from the fact that the operators $`𝑩_{(\lambda ,\lambda )f,SM}`$ contain an antiquark annihilation operator which anticommutes with all the quark creation operators. Therefore, the Hilbert space of the model may be divided into sectors, each one with a different vacuum state having a given baryon number. The one with the baryon number zero is the real particle vacuum.

### 2.1 Group Theory of the Fermion Part

From now on, we restrict to $`2\mathrm{\Omega }=n_sn_cn_f`$ $`=`$ 18 with $`n_s=2`$, $`n_c=3`$ and $`n_f=3`$, for the spin, color and flavor degrees of freedom, respectively. The largest group, whose generators are $`c_{c_1f_1\sigma _1i}^{}c^{c_2f_2\sigma _2j}`$ ($`c_i=1,2,3`$, $`f_i=1,2,3`$, $`\sigma _i=1,2`$ and $`i,j=1,2`$), is the $`U(4\mathrm{\Omega })`$ group. One possible group chain for the classification of the states, which include the flavor ($`SU_f(3)`$) and the spin ($`SU_S(2)`$) groups, is given by

$`[1^N]`$ $`[h]=[h_1h_2h_3]`$ $`[h^T]`$
$`U(4\mathrm{\Omega })`$ $`U({\displaystyle \frac{\mathrm{\Omega }}{3}})`$ $`U(12)`$ (3)
$``$ $``$
$`(\lambda _C,\mu _C)SU_C(3)(\lambda _f,\mu _f)`$ $`SU_f(3)SU_S(2)S,M,`$
where the irreducible representation (irrep) of $`U(4\mathrm{\Omega })`$ is the completely antisymmetric one and $`N`$ is the number of particles involved. The upper index in $`[h^T]`$ refers to the transposed Young diagram of $`[h]`$, where the columns and rows are interchanged . Due to the antisymmetric irrep $`[1^N]`$ of $`U(4\mathrm{\Omega })`$ the irreps of $`U(\mathrm{\Omega }/3)`$ and $`U(12)`$ are complementary and the irrep of $`U(\mathrm{\Omega }/3)`$, which is for $`\mathrm{\Omega }=9`$ the color group, has maximally three rows . In the group chain (3) no multiplicity labels are indicated. There is a multiplicity $`\rho _f`$ for $`(\lambda _f,\mu _f)`$ and $`\rho _S`$ for the spin $`S`$. The color labels $`(\lambda _C,\mu _C)`$ are related to the $`h_i`$ via $`\lambda _C=h_1h_2`$ and $`\mu _C=h_2h_3`$. The complete state is given by

$`|N,(\lambda _C,\mu _C),\rho _f(\lambda _f,\mu _f)YTT_z,\rho _SSM>,`$ (4)
where $`Y`$ is the hypercharge, $`T`$ is the isospin and $`T_z`$ its third component. For meson-like states, the color quantum numbers to be considered are $`(\lambda _C,\mu _C)=(0,0)`$. These states will be located in an elementary volume of about $`\frac{4\pi }{3}`$fm3, corresponding to a sphere of radius 1 fm.

To obtain the values of $`h_i`$ one has to consider all possible partitions of $`N=h_1+h_2+h_3`$, which fixes the color. For colorless states we have $`h_1=h_2=h_3=h`$. Each partition of $`N`$ appears only once. The irrep $`[hhh]`$ of $`U(\frac{\mathrm{\Omega }}{3})=U(3)`$ ($`\mathrm{\Omega }=9`$) fixes the irrep of $`U(12)`$, as indicated in (3). For the reduction of the irrep of $`U(12)`$ we have written a computer code , which is available to the interested reader. As an example, let us consider the $`U(12)`$ irrep $`[3^60^6]`$ and the two $`U(4)`$ irreps $`[9^20^2]`$ and $`[9720]`$, where the first one contains the state where the lower level is completely filled and the upper one empty, and the second irrep is the next highest one which contains flavor $`(0,0)`$ . The first is accompanied by flavor $`(0,0)_1`$ and the second one by $`(0,0)_1`$ and $`(1,1)_1`$, where the subindex denotes the multiplicity. The spin content of $`[9^20^2]`$ is given by $`0_{55}`$, $`1_{45}`$, $`2_{36}`$, $`3_{28}`$, $`4_{21}`$, $`5_{15}`$, $`6_{10}`$, $`7_6`$,$`8_3`$ and $`9_1`$. The spin content of $`[9720]`$ is $`0_{81}`$, $`1_{171}`$, $`2_{189}`$, $`3_{135}`$, $`4_{90}`$, $`5_{54}`$, $`6_{27}`$ and $`7_9`$. The lowest dimensional irrep is $`[5^24^2]`$ with the spin content $`0_1`$ and $`1_1`$.

### 2.2 The Boson Mapping

The explicit construction of the basis states, Eq. (4) and the calculation of the matrix elements, become very involved, which is in conflict with the idea to develop a simple model. A way out of it is to use a boson mapping of the pair operators $`𝑩^{}`$ and $`𝑩`$ and work in the boson model space.

The basic ingredients of the model are the pair operators, given in (2). They can be mapped onto boson operators

$`𝑩_{f_1\sigma _1}^{f_2\sigma _2}`$ $``$ $`𝒃_{f_1\sigma _1}^{f_2\sigma _2}`$
$`𝑩_{f_1\sigma _1}^{f_2\sigma _2}`$ $``$ $`𝒃_{f_1\sigma _1}^{f_2\sigma _2}.`$ (5)
where the operators on the right satisfy the normal boson commutation relations.

$`[𝒃_{f_1\sigma _1}^{f_2\sigma _2},𝒃_{f_3\sigma _3}^{f_4\sigma _4}]=\delta _{f_3f_2}\delta _{f_4f_1}\delta _{\sigma _3\sigma _2}\delta _{\sigma _4\sigma _1}.`$ (6)
The exact boson mapping is quite involved, but it can be obtained in general . For the sake of this work, it is worth to show that the mapping can be performed, indeed. We shall: i) work from the beginning in the boson space, ii) define a Hamiltonian which corresponds to the Hamiltonian acting in the fermion space, and iii) select a basis. The advantage of working in the boson space is the simplification in getting the matrix elements (see below). The price to pay is related to the appearance of non-physical states , as we shall discuss later on.

In order to choose a basis in the boson space, we profit from the fact that the basic degrees of freedom are given by the boson creation operators $`b_{(\lambda ,\lambda )fSM}^{}`$ $`=`$ $`b_{\lambda fSM}^{}`$, with $`\lambda `$ $`=`$ 0 or 1 and $`S`$ $`=`$ 0 or 1. This gives four possible combinations of $`[\lambda ,S]`$: , ,  and . Consequently, the total Hilbert space is the direct product of a one, three, eight and 24 dimensional harmonic oscillators .

For each harmonic oscillator we can define a seniority basis

$`𝒩_{N_{\lambda S}\nu _{\lambda S}}(𝒃_{\lambda S}^{}𝒃_{\lambda S}^{})^{\frac{N_{\lambda S}\nu _{\lambda S}}{2}}|\nu _{\lambda S}\alpha _{\lambda S}>,`$ (7)
where $`N_{\lambda S}`$ is the number of bosons of type $`[\lambda ,S]`$, $`\nu _{\lambda S}`$ the corresponding seniority and $`𝒩_{N_{\lambda S}\nu _{\lambda S}}`$ is a normalization constant. The seniority is the number of $`𝒃_{\lambda S}`$-bosons not coupled into pairs. The $`\alpha _{\lambda S}`$ contain all other quantum numbers for a particular harmonic oscillator. The dot in the factor refers to the scalar product.

The choice of a seniority basis is particularly useful for the calculation of the matrix elements of the interaction, which contains expressions of the form $`(𝒃_{\lambda S}^{}𝒃_{\lambda S}^{})`$, $`(𝒃_{\lambda S}𝒃_{\lambda S})`$ and $`(𝒃_{\lambda S}^{}𝒃_{\lambda S})`$, where the latter is just the number operator of the bosons of the type $`[\lambda ,S]`$. The exact structure of $`|\nu _{\lambda S}\alpha _{\lambda S}>`$ is not needed, except for the knowledge of the quantum numbers $`\alpha _{\lambda S}`$.

For the one dimensional harmonic oscillator () the seniority can take the values 0 or 1. The state is of the form $`(𝒃_{00}^{})^{N_{00}}|0>`$ $`=`$ $`(𝒃_{00}^{}𝒃_{00}^{})^{\frac{N_{00}\nu _{00}}{2}}`$ $`(𝒃_{00}^{})^{\nu _{00}}|0>`$. For the three dimensional harmonic oscillator the seniority is equal to the spin $`S_{\lambda S}`$. The explicit expression of the state is given in Ref. . The eight dimensional oscillator contributes to flavor only and it is discussed in Appendix A. The 24 dimensional oscillator can be found in Ref. , where the color part in Ref.  has to be interpreted here as the flavor part. In Ref.  only singlet states are listed, but the procedure to obtain non-singlet flavor states is outlined.

The parity of each state is given by $`P=(1)^N`$, where $`N`$ $`=`$ $`_{\lambda ,S}N_{\lambda S}`$ is the total number of bosons. Each boson stems from a particle-antiparticle pair, which carries negative parity.

In order to obtain the property under charge conjugation, one has to apply the charge conjugation operator $`𝑪`$ to the pair creation operator $`𝑩_{(\lambda ,\lambda )f,SM}^{}`$. The result is (see Appendix B)

$`𝑪𝑩_{\lambda f,SM}^{}𝑪^1`$ $`=`$ $`(1)^S𝑩_{\overline{\lambda }\overline{f},SM}^{},`$ (8)
where $`\overline{\lambda }=(\lambda ,\mu )`$, $`\overline{f}`$ $`=`$ $`Y`$, $`T`$, $`T_z`$. From this it is clear that only states with $`Y=0`$, $`T_z=0`$, and $`\mu =0`$, can have a definite $`C`$-parity. In Eq. (8) we make use of the application of the operator $`𝑪`$ which interchanges quark and antiquark operators ($`a^{}d^{}`$) and inverts the magnetic quantum numbers ($`Y_iY_i`$ and $`T_{iz}T_{iz}`$) of flavor and of color only 111After here, for the sake of notation, we shall indicate charge conjugate states with a bar on the index f.

For products of two pair creation operator we obtain

$`𝑪\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]_{f,M}^{\rho (\lambda ,\mu )}𝑪^1`$ $`=`$ $`(1)^{S_1+S_2\lambda \mu +\rho _{max}\rho }\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]_{\overline{𝐟},𝐌}^{\rho (\mu ,\lambda )},`$ (9)
where $`\rho `$ is the multiplicity label of $`(\lambda ,\mu )`$ in the product $`(\lambda _1,\lambda _1)(\lambda _2,\lambda _2)`$. The symbol $`\rho _{max}`$ denotes the maximal value of $`\rho `$. The phase convention of Ref.  was used. The symbol $``$ denotes the combined product in $`SU_f(3)`$ and $`SU_S(2)`$.

In analogy, the action of the charge conjugation on a product of three pair operators can be obtained:

$`𝑪\left[\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]^{\rho _{12}(\lambda _{12},\mu _{12}),S_{12}}𝑩_{\lambda _3,S_3}^{}\right]_{f,M}^{\rho (\lambda ,\mu ),S}𝑪^1=`$ (10)
$`(1)^{S_1+S_2+S_3\lambda \mu +\rho _{12,max}\rho _{12}+\rho _{max}\rho }\left[\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]^{\rho _{12}(\mu _{12},\lambda _{12}),S_{12}}𝑩_{\lambda _3,S_3}^{}\right]_{\overline{𝐟},𝐌}^{\rho (\mu ,\lambda ),S}`$ $`,`$
where $`\rho _{12}`$ is the multiplicity of $`(\lambda _{12},\mu _{12})`$ in the product of $`(\lambda _1,\lambda _1)(\lambda _2,\lambda _2)`$, $`\rho `$ is the multiplicity of the total irrep in the last coupling of the above equation, and $`\rho _{12,max}`$ is the maximal value of $`\rho _{12}`$.

The procedure outlines here can be used in a recursive way for more involved coupling schemes. For our purpose it is sufficient to go up to three pairs, which will be the dominant structure at low energy.

For the G-parity the additional rotation in the isospin space has to be applied, which changes $`T_{i,z}`$ to $`T_{i,z}`$ . For a polynomial in the pair operators this gives an additional phase $`(1)^T`$, where $`T`$ is the total isospin .

The same phase properties under C- and G-parity transformation have to be valid for the mapped boson operators $`𝒃_{\lambda f,SM}^{}`$.

In a seniority basis, the matrix elements are easily obtained, and they are written

$`<N_{\lambda S}+2\nu _{\lambda S}\alpha _{\lambda S}|(𝒃_{\lambda S}^{}𝒃_{\lambda S}^{})|N_{\lambda S}\nu _{\lambda S}\alpha _{\lambda S}>`$ $`=`$ $`\sqrt{(N_{\lambda S}\nu _{\lambda S}+2)(N_{\lambda S}+\nu _{\lambda S}+d_{\lambda S})}`$
$`<N_{\lambda S}2\nu _{\lambda S}\alpha _{\lambda S}|(𝒃_{\lambda S}𝒃_{\lambda S})|N_{\lambda S}\nu _{\lambda S}\alpha _{\lambda S}>`$ $`=`$ $`\sqrt{(N_{\lambda S}\nu _{\lambda S})(N_{\lambda S}+\nu _{\lambda S}+d_{\lambda S}2)}`$
$`<N_{\lambda S}\nu _{\lambda S}\alpha _{\lambda S}|(𝒃_{\lambda S}^{}𝒃_{\lambda S})|N_{\lambda S}\nu _{\lambda S}\alpha _{\lambda S}>`$ $`=`$ $`N_{\lambda S},`$
where $`d_{\lambda S}`$ is 1, 3, 8 or 24 for the case of the one, three, eight or 24 dimensional harmonic oscillator. As a short hand notation we will use instead of $`(𝒃_{\lambda S}^{}𝒃_{\lambda S}^{})`$ the expression $`(𝒃_{\lambda S}^{})^2`$, and similarly for the other products, $`(𝒃_{\lambda S})^2`$ and $`𝒃_{\lambda S}^{}𝒃_{\lambda S}`$.

As a Hamiltonian, invariant under rotation, charge conjugation and G-parity, we propose

$`𝑯`$ $`=`$ $`2\omega _f𝒏_𝒇+\omega _b𝐧_𝐛+`$ (12)
$`{\displaystyle \underset{\lambda S}{}}V_{\lambda S}\{[(𝒃_{\lambda S}^{})^2+2𝒃_{\lambda S}^{}𝒃_{\lambda S}+(𝒃_{\lambda S})^2](1{\displaystyle \frac{𝒏_𝒇}{2\mathrm{\Omega }}})𝒃+`$
$`𝒃^{}(1{\displaystyle \frac{𝒏_𝒇}{2\mathrm{\Omega }}})[(𝒃_{\lambda S}^{})^2+2𝒃_{\lambda S}^{}𝒃_{\lambda S}+(𝒃_{\lambda S})^2]\}.`$
Due to symmetry arguments, the interaction strength $`V_{\lambda S}`$ is the same for the two last lines in Eq. (12). The term $`(𝒃_{\lambda S}^{})^2`$ ($`(𝒃_{\lambda S})^2`$) describes the creation (annihilation) of two quark-antiquark pairs with the simultaneous creation or annihilation of a gluon pair. The term $`𝒃_{\lambda S}^{}𝒃_{\lambda S}`$, in Eq. (12), describes the scattering of a fermion pair with the emission or annihilation of a gluon pair. All processes can be depicted by a Feynman graph and all graphs can be obtained from any other one by an appropriate interchange of lines. Because the strength $`V_{\lambda S}`$ should be, basically, invariant under the exchange of lines, we shall use the same interaction strength for all channels, as a first approximation. The terms which appear in Eq. (12) originate in the normal product of $`:(𝒃_{fS}^{}+𝒃_{fS})^2:`$, where the square implies a scalar product. The factor $`(1\frac{𝒏_𝒇}{2\mathrm{\Omega }})`$ represents a cutoff which can be traced back to the boson mapping of the fermion pairs with flavor (0,0) and spin 0. This term simulates the effect of an exact boson mapping  and it is responsible for the disappearance of the interaction when the number of pairs reaches $`2\mathrm{\Omega }`$. In other words, this cutoff term simulates the Pauli principle which does not allow more than $`2\mathrm{\Omega }`$ pairs.

The Hamiltonian (12) is the most simple form we can think of and it contains only four parameters (the values of $`V_{\lambda S}`$). The value of $`\omega _f`$ is fixed to 0.33 GeV, which is about $`\frac{1}{3}`$ of the mass of a nucleon. The most notorious difficulty, associated to the use of the boson mapping, lies in the Hilbert space of the boson operators. It is larger than the Hilbert space of the fermion pairs. In some situations one can identify the source of the spurious dimensions. If, for example, only flavor (0,0) and spin 0 pairs are taken into account, the relevant group structure is $`U(4\mathrm{\Omega })`$ $``$ $`U(2\mathrm{\Omega })U(2)`$. The irrep of $`U(4\mathrm{\Omega })`$ has to be antisymmetric, which implies that the irreps of $`U(2\mathrm{\Omega })`$ and $`U(2)`$ have to be complementary, i.e. if $`U(2)`$ is given by a Young diagram with two rows, the one of $`U(2\mathrm{\Omega })`$ has to be the adjoint, which is obtained by interchanging rows and columns . The upper limit, up to which no spurious states appear, is $`2\mathrm{\Omega }`$ because $`U(2\mathrm{\Omega })`$ allows $`2\mathrm{\Omega }`$ rows in the Young diagram. This is also the maximum number of pairs allowed, i.e. for this case no un-physical states occur. If flavor values (0,0) and (1,1) and spin 0 are used, only, we have $`U(4\mathrm{\Omega })`$ $``$ $`U(\frac{2\mathrm{\Omega }}{3})U(3)`$ and up to $`\frac{2\mathrm{\Omega }}{3}`$ pairs there is no problem with respect to the appearance of un-physical states. This implies that states with explicit flavor will present un-physical states only for large number of bosons. If flavor (0,0) and spin 0 and 1 are considered, we have $`U(4\mathrm{\Omega })`$ $``$ $`U(\mathrm{\Omega })U(4)`$, which has as an upper limit the number $`\mathrm{\Omega }`$ up to which no spurious states appear. Finally, for all pairs, i.e. flavor (0,0), (1,1) and spin 0, 1 the relevant group chain is $`U(4\mathrm{\Omega })`$ $``$ $`U(\frac{\mathrm{\Omega }}{3})U(12)`$ and the upper limit is $`\frac{\mathrm{\Omega }}{3}`$. This gives us a hint about the group sequence where un-physical states may appear. The upper limit up to which all states are physical is lowered in the sequence where bosons with more degrees of freedom appear. There are states which can be described either by one or the other type of bosons or even by a combination of bosons. For example, when both levels, the lower and the upper ones, are filled there is only one allowed state, which is the flavor singlet with spin zero. However, all types of boson pairs can describe it, e.g. when the number of bosons coupled to flavor singlet and spin zero is equal to $`2\mathrm{\Omega }`$. In view of these considerations, and for the one, three, eight and twenty-four dimensional harmonic oscillator basis, we have introduced the limits $`2\mathrm{\Omega }`$, $`\mathrm{\Omega }`$, $`\frac{2\mathrm{\Omega }}{3}`$ and $`\frac{\mathrm{\Omega }}{3}`$, respectively. The higher non physical states do not play an essential role because, as shown below, the dominant contribution at low energy comes from configurations with a small number of quark-antiquark pairs . By working with these dimensional cut-off values the influence of non-physical states is minimized. Also, for each case, the total number of bosons is restricted to $`2\mathrm{\Omega }`$. For a reasonable interaction strength, however, the dominant contribution comes from a small number of bosons. In such cases, the number of un-physical states is small and they do not influence much the result. The dimensional cut-off in the Hilbert space has the advantage that most un-physical states are excluded. In principle, one can eliminate the spurious states by applying another, more involved, procedure. For that one has to reduce the irrep of $`U(12)`$ to the flavor and spin groups, as done in the last section. This gives us the allowed content of flavor and spin for a given irrep of $`U(12)`$. The matching condition, i.e. by comparing for a given number of pairs the spin and flavor content on the boson side to the one on the fermion side, eliminates un-physical states. If on the boson and fermion sides, for a given flavor and spin, the number of states are equal, all states in the model space are taken into account. This is the case of low lying basis states. A simple counting procedure can be used for other situations. If there are, for a given flavor and spin, more states in the boson space than in the fermion space, one can not decide easily which combination is allowed. However, one can reduce the number of states of the model space to the same dimension as the one of the fermion space. As a rule one can first eliminate the states which contain most of the bosons with the largest degree of freedom, i.e. which are of the type , and in this way proceed, if necessary, until only states with flavor (0,0) and spin zero bosons are left. At least, the proposed procedure eliminates most of the spurious states. The error made can be absorbed in the parameters of the model, a general practice in dealing with phenomenological models, because in the end the correct number of degrees of freedom (dimension of the Hilbert space) dominates in a successful description of the spectrum. The idea of the proposed procedure is not new and it was used in another context by J. Cseh et al. .

## 3 Application to the meson spectrum

The Hamiltonian (12) commutes with the isospin operator and it does not depend explicitly on the hypercharge $`Y`$. As a consequence, all states which belong to the same flavor irrep are degenerate. In principle, we can add terms proportional to $`𝑻^2`$, $`𝒀`$ and $`𝒀^2`$ in order to lift the degeneracy. These terms will add new parameters to the four already present ($`\omega _b`$ is fixed, as in Ref. ). Also, a flavor mixing term could be added, as suggested by the $`\eta `$-$`\eta ^{}`$ mixing . In order to simplify the discussion, we shall first ignore these additional interaction terms.

In fitting the spectrum of the mesons we will use, as an experimental input for each multiplet, only the state with $`T=0`$ and $`Y=0`$. For an octet all states have the same energy as the isospin singlet and hypercharge zero state. Because later on we shall take into account flavor mixing interactions too, the position of the singlet and octet state are not fixed at the measured energy values but at the values obtained when the flavor mixing interactions are switched off . The mixing angle is introduced for two multiplets: the (1,1) $`0^{}`$ and (0,0) $`0^{}`$ irrep, containing the pions and the $`\eta `$, $`\eta ^{}`$, and the (1,1) $`1^{}`$, (0,0) $`1^{}`$, containing the $`\omega `$, $`\varphi `$ and $`\rho `$ particles. The mixing angles are, respectively, $`23.7^0`$ and $`35.3^0`$, . For other multiplets one assumes that the mixing angle is zero, because of missing data, and because of the smallness of the energy splitting between members, as compared to the energy splitting within the multiplet which contains the pions or the $`\rho `$ mesons. The uncorrected masses for, e.g., the octet (before flavor mixing) are $`m_8=615`$ MeV in the first case (see notation of Ref. ) and $`m_8=940`$ MeV for the second case (see also Table 1).

In Table 1 we show the states used in the fit. Their flavor, spin and parity are indicated together with the experimental values. In total, to perform the fit, we have considered thirteen states with spin zero and one in the four parameter fitting procedure. All other states are predicted, particularly those with spin 2 and 3.

In Figure 2 we give the spectrum for spin 0 and 1 meson states without any particle number changing interaction. On each side of a level the flavor quantum numbers and its degeneracy are indicated. This serves to illustrate that the multiplicity at energies lower than 2 GeV is already very large. This is a consequence of the various manners in which the same set of quantum numbers can be obtained, for a given configuration, when many quarks, antiquarks and gluons are considered. This is known as the multiplicity problem. The result of the best fit values, obtained after the interaction is turned on, is given in Figures 3-6 for spin 0, 1, 2 and 3 respectively. Only states which correspond to non exotic parity/charge conjugation quantum numbers are shown. Most of them appear above 2 GeV and some can be deduced from the gluon spectrum published in Ref. . In Figures 3-6 each theoretical spectrum is compared to the experimental one. On the right hand side of each level the theoretical interpretation in terms of flavor and the multiplicity is indicated. On the left hand side of each spin ($`J^{PC}`$) the experimental information is given, taken from the particle data group . The energy of these states, appearing in the summary table of Ref. , is given in boxes and the experimental error is reflected by the size of the box. If the error is very small, the box is replaced by a line. States which are not in the summary table of  are indicated by dashed boxes (lines). Only states which correspond to isospin singlet and hypercharge zero, after having corrected for the isospin mixing, are listed.

Note that very few states have a multiplicity. Most states were pushed upwards due to the interaction. This is an effect of the interaction, because it changes the number of particles and relates the quark-antiquark sector with the gluons. Models with a particle conserving interaction will always present the multiplicity problem. Thus, the particle mixing interaction is essential to remove the multiplicity problem.

The spin 2 and 3 states where not adjusted. The theoretical results seems to agree with the data, because the states are predicted at the correct energy domain. The density of states for a given flavor also seems to be reproduced.

In Table 2 we show the quark-antiquark pair and gluon pair contents of some selected states. The total number of quark-antiquark pairs is denoted by $`<n_q>`$, where the symbol $`<\mathrm{}>`$ indicates the expectation value of this number. The quantity $`<n_{ij}>`$ gives the average number of boson pairs of the type \[i,j\], while $`<n_g>`$ is the expectation value of the number of gluon pairs with spin zero. The total number of gluons is twice $`<n_g>`$.

The structure of the calculated ground state (physical vacuum) is an interesting piece of information about the model predictive power. The calculated value of the ground state energy is equal to -0.726 GeV. The physical vacuum state contains about 3.1 quark-antiquark pairs of type , and the other configurations contribute with 0.06 pairs. The dominance of the $`[1,0]`$ quark-antiquark pairs is consistent with the strengths of the interactions. The parameters, obtained in the fit, are : $`m_f=0.33`$ GeV, $`V_{0,0}=0.0337`$ GeV, $`V_{0,1}=0.0422`$ GeV, $`V_{1,0}=0.1573`$ GeV and $`V_{1,1}=0.0177`$ GeV. The channel  is clearly the strongest one. Also, the number of gluons is noticeable, i.e 1.7 pairs, which correspond to more than three gluons contained in the elementary hadron volume. This indicates a collective behavior of QCD states at low energy. Indeed, the pion state (corresponding to $`\eta (541)`$ in Table 2) contains about 2.7 pairs of the type  while the rest contributes with about 0.06 pairs, as in the ground state. The number of gluon pairs in the pion, 1.2 pairs or 2.4 gluons, is similar to the number of gluon pairs in the ground state. They constitute about 30 per cent of the particle content. In brief, the calculated states contain a large number of quark-antiquark pairs and gluons. Roughly speaking, no single state can be described approximately by a pure quark-antiquark pair. Note that in this respect theory and experiment do agree, in spite of the simplicity of the model. We think that the complex structure of the meson spectrum is described qualitatively by our model, as well as several other features, like the position of the first states with spin 2 and 3, and the density of states with flavor (0,0) or (1,1) for each spin and parity, charge conjugation number.

Concerning baryons, we have to include them yet in the model. Up to now they are described as spectators, i.e. without having an explicit coupling to the quark, antiquark, and gluon sea. For that, further interaction terms should be introduced. For example, one can introduce the interaction

$`𝒏_{D,(0,1)0}(𝒃^{}+𝒃)+𝒏_{D,(2,0)1}(𝒃^{}+𝒃),`$ (13)
where $`𝒏_{D,(\lambda ,\mu )S}`$ is the number operator of a Di-quark coupled to flavor $`(\lambda ,\mu )`$ and spin $`S`$. This is analogous to the above ansatz of the Hamiltonian. The product of two pair creation operators of Di-quarks can not appear because this would mix the baryon number. The interaction in (13) is a direct extension from the pair operator interaction of the former Hamiltonian. The first term in (13) acts only on states like the nucleon octet and the last one on particles like the baryon decouplet. The residual interaction will mix the number of gluons, and, mainly, increase the gluon content in the baryons, from a gluon content of about 30 percent to, say, 50 percent. The inclusion of baryons will be reported in another work  .

## 4 Conclusions

In the present paper we have advanced a schematic model of QCD, based on a Lipkin-type model for fermions and an interaction to one gluon pair states. We have discussed the low energy structure of the model. The Hamiltonian is composed by a diagonal, particle conserving part, and an interaction which couples the quark-antiquark pairs to the gluons and changes the number of particles. The model contains only 4 parameters which were adjusted to reproduce 13 observed meson states with spin 0 and 1. After fixing these parameters, we have predicted the remaining part of the spectrum. The complex structure of the meson spectrum was qualitatively reproduced by our results.

Due to the schematic nature of the model one cannot expect to be able to reproduce all details of the low energy meson spectrum. However, the results are in qualitative agreement with data, a fact that shows the validity of the model as a toy model for QCD. Baryons where not considered but the extension to this sector was briefly indicated. The baryons would correspond to states where three extra quarks are added in the upper level. The corresponding operators will commute with the boson pair operators and an interaction of baryon states with the quark-antiquark sea should be included in the model.

We have found that the inclusion of particle mixing interactions turns out to be essential in order to remove the multiplicity problem encountered in other models, when states with many quarks and antiquarks are considered. This particle changing interaction also introduces ground state correlations resulting in many quarks, antiquarks and gluons configurations in the states. It produces a large contribution of the gluons and the total spin is not a simple product of a quark-antiquark state but of many quarks, antiquarks and gluons. This illustrates the fact that, even at low energy, the structure of the hadron states is by no means as simple as suggested by earlier particle conservation models. It also shows that phenomenological potentials, which simulate the presence of gluons in a pure quark model, cannot resolve the problem of multiplicity.

## Acknowledgment

We are pleased to thank Prof. M. Moshinsky for his renewed interest in our work. We acknowledge financial support through the CONACyT-CONICET agreement under the project name Algebraic Methods in Nuclear and Subnuclear Physics and from CONACyT project number 32729-E. (S.J.) acknowledges financial support from the Deutscher Akademischer Austauschdienst (DAAD) and SRE, (S.L) acknowledges financial support from DGEP-UNAM. Financial support from DGAPA, project number IN119002, is acknowledged.

## Appendix A: The eight dimensional oscillator

The reduction of the eight dimensional oscillator group, $`U(8)`$, to the flavor group $`SU_f(3)`$ is discussed in Ref. . As an intermediary group, between $`U(8)`$ and $`SU_f(3)`$ one can use the $`SO(8)`$ group. Though, in Ref.  only the reduction to flavor singlet groups is listed, the general procedure is outlined. Programs are available on request  and the procedure has been published elsewhere .

The generators of the $`U(8)`$ group are given by $`𝒃_{(1,1)f_1,00}^{}𝒃_{(1,1)f_2,00}`$, where the zeros refer to zero spin and its projection. Therefore, these bosons can only contribute to the flavor content. In Table 3 we list the flavor content of up to four bosons of the type $`𝒃_{(1,1)f,00}^{}`$. As one can see, the multiplicity raises especially for the (1,1) flavor irrep. With the help of the $`SO(8)`$ group one can further reduce the multiplicity. For our purpose this is not necessary.

## Appendix B: Parity, Charge Conjugation and G-Parity

The charge conjugation operator acts as follows on the quark and antiquark creation operators

$`𝑪𝒂_{cf\sigma }^{}𝑪^1`$ $`=`$ $`𝒅_{\overline{c}\overline{f}\sigma }^{}`$
$`𝑪𝒅_{cf\sigma }^{}𝑪^1`$ $`=`$ $`𝒂_{\overline{c}\overline{f}\sigma }^{},`$ (14)
where the $`𝒂^{}`$ transforms in color and flavor as a (1,0) $`SU(3)`$ irrep, while the $`𝒅^{}`$ transform as (0,1). If one of the color or flavor index is raised then the $`𝒅^{}`$ transform as (1,0). The ”$`\overline{c}`$” and ”$`\overline{f}`$” refer to the reflection in the magnetic quantum numbers of $`SU(3)`$. I.e. $`f`$ stands for $`Y`$, $`T`$ and $`T_z`$ and $`\overline{f}`$ for $`Y`$, $`T`$ and $`T_z`$ and similar for ”$`\overline{c}`$”.

With this, the action of the charge conjugation operator on a quark-antiquark pair is given by

$`𝑪𝑩_{\lambda f,SM}^{}𝑪^1=𝑪{\displaystyle \underset{cf_1f_2\sigma _1\sigma _2}{}}𝒂_{cf_1\sigma _1}^{}𝒅_{f_2\sigma _2}^c<(1,0)f_1,(0,1)f_2|(\lambda _1,\lambda _1)f>_1`$
$`({\displaystyle \frac{1}{2}}\sigma _1,{\displaystyle \frac{1}{2}}\sigma _2|SM)𝑪^1`$
$`={\displaystyle \underset{cf_1f_2\sigma _1\sigma _2}{}}𝒅_{\overline{c}\overline{f}_1\sigma _1}^{}𝒂_{\overline{f}_2\sigma _2}^{\overline{c}}<(1,0)f_1,(0,1)f_2|(\lambda _1,\lambda _1)f>_1({\displaystyle \frac{1}{2}}\sigma _1,{\displaystyle \frac{1}{2}}\sigma _2|SM)`$
$`={\displaystyle \underset{cf_1f_2\sigma _1\sigma _2}{}}𝒂_{\overline{f}_2\sigma _2}^c𝒅_{c\overline{f}_1\sigma _1}^{}<(1,0)\overline{f}_1,(0,1)\overline{f}_2|(\lambda _1,\lambda _1)\overline{f}>_1({\displaystyle \frac{1}{2}}\sigma _1,{\displaystyle \frac{1}{2}}\sigma _2|SM)`$
$`=(1)^{2\lambda _1+1S}{\displaystyle \underset{cf_1f_2\sigma _1\sigma _2}{}}𝒂_{\overline{f}_2\sigma _2}^c𝒅_{c\overline{f}_1\sigma _1}^{}<(1,0)\overline{f}_2,(0,1)𝒇_1|(\lambda _1,\lambda _1)\overline{f}>_1({\displaystyle \frac{1}{2}}\sigma _1,{\displaystyle \frac{1}{2}}\sigma _2|SM)`$
$`=(1)^S𝑩_{\overline{\lambda }\overline{f},S_1M}^{},`$ (15)
where we made use of the properties of the $`SU(2)`$  and $`SU(3)`$  Clebsch-Gordan coefficients. The subindex 1 in the SU(3) Clebsch-Gordan coefficient indicates a multiplicity of one .

For the product of two pair operators we have

$`𝑪\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]_{f,M}^{\rho (\lambda ,\mu )}𝑪^1`$
$`=𝑪{\displaystyle \underset{f_1f_2M_1M_2}{}}𝑩_{\lambda _1f_1,S_1M_1}^{}𝑩_{\lambda _2f_2,S_2M_2}^{}<(\lambda _1,\lambda _1)f_1,(\lambda _2,\lambda _2)f_2|(\lambda ,\mu )f>_\rho `$
$`(S_1M_1,S_2M_2|SM)𝑪^1`$
$`=`$ $`{\displaystyle \underset{f_1f_2M_1M_2}{}}(1)^{S_1+S_2}𝑩_{\lambda _1\overline{f}_1,S_1M_1}^{}𝑩_{\lambda _2\overline{f}_2,S_2M_2}^{}<(\lambda _1,\lambda _1)f_1,(\lambda _2,\lambda _2)f_2|(\lambda ,\mu )f>_\rho `$
$`(S_1M_1,S_2M_2|SM)`$
$`=`$ $`(1)^{S_1+S_2\lambda \mu +\rho _{max}\rho }{\displaystyle \underset{f_1f_2M_1M_2}{}}𝑩_{\lambda _1\overline{f}_1,S_1M_1}^{}𝑩_{\lambda _2\overline{f}_2,S_2M_2}^{}<(\lambda _1,\lambda _1)\overline{f}_1,(\lambda _2,\lambda _2)\overline{f}_2|(\lambda ,\mu )\overline{f}>_\rho `$
$`(S_1M_1,S_2M_2|SM)`$
$`=`$ $`(1)^{S_1+S_2\lambda \mu +\rho _{max}\rho }\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]_{\overline{𝐟},𝐌}^{(\mu ,\lambda )},`$ (16)
where $`\rho `$ is the multiplicity index of $`(\lambda ,\mu )`$ in the product $`(\lambda _1,\lambda _1)(\lambda _2,\lambda _2)`$.

For the product of three pairs we have

$`𝑪\left[\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]^{(\lambda _{12},\mu _{12}),S_{12}}𝑩_{\lambda _3,S_3}^{}\right]_{f,M}^{(\lambda ,\mu ),S}𝑪^1`$
$`={\displaystyle \underset{f_{12}f_3M_{12}M_3}{}}𝑪\left[𝑩_{\lambda _1S_1}^{}𝑩_{\lambda _2S_1}^{}\right]_{f_{12}M_{12}}^{(\lambda _{12},\mu _{12})S_{12}}𝑪^1𝑪𝑩_{\lambda _3f_3S_3M_3}^{}𝑪^1`$
$`<(\lambda _{12},\mu _{12})f_{12},(\lambda _3,\lambda _3)f_3|(\lambda ,\mu )f>_\rho (S_{12}M_{12},S_3M_3|SM)𝑪^1`$
$`=`$ $`{\displaystyle \underset{f_{12}f_3M_{12}M_3}{}}(1)^{S_1+S_2\lambda _{12}\mu _{12}+\rho _{12,max}\rho _{12}}\left[𝑩_{\lambda _1S_1}^{}𝑩_{\lambda _2S_1}^{}\right]_{\overline{f}_{12}M_{12}}^{(\mu _{12},\lambda _{12})S_{12}}(1)^{S_3}𝑩_{\lambda _3\overline{f}_3S_3M_3}^{}`$
$`<(\lambda _{12},\mu _{12})f_{12},(\lambda _3,\lambda _3)f_3|(\lambda ,\mu )f>_\rho (S_{12}M_{12},S_3M_3|SM)𝑪^1`$
$`=`$ $`(1)^{S_1+S_2+S_3\lambda _{12}\mu _{12}+\rho _{12,max}\rho _{12}}{\displaystyle \underset{f_{12}f_3M_{12}M_3}{}}\left[𝑩_{\lambda _1S_1}^{}𝑩_{\lambda _2S_1}^{}\right]_{\overline{f}_{12}M_{12}}^{(\mu _{12},\lambda _{12})S_{12}}𝑩_{\lambda _3\overline{f}_3S_3M_3}^{}`$
$`<(\mu _{12},\lambda _{12})\overline{f}_{12},(\lambda _3,\lambda _3)\overline{f}_3|(\mu ,\lambda )\overline{f}>_\rho (1)^{\lambda _{12}+\mu _{12}\lambda \mu +\rho _{max}\rho }(S_{12}M_{12},S_3M_3|SM)`$
$`=`$ $`(1)^{S_1+S_2+S_3\lambda \mu +\rho _{12,max}\rho _{12}+\rho _{max}\rho }\left[\left[𝑩_{\lambda _1,S_1}^{}𝑩_{\lambda _2,S_2}^{}\right]^{(\mu _{12},\lambda _{12}),S_{12}}𝑩_{\lambda _3,S_3}^{}\right]_{\overline{𝐟},𝐌}^{(\mu ,\lambda ),S},`$ (17)
with the use of the notation of Ref.  for the $`SU(3)`$ Clebsch-Gordan coefficients and their symmetry properties.