text:
# Rosseland and Planck Mean Opacities for Protoplanetary Discs

## 1 Introduction

Recently, significant progress toward the understanding of the possible composition and properties of dust grains and gas species in many astrophysical environments has been achieved. For instance, (sub)millimetre observations of molecular lines provided basic information about chemical composition and dynamical properties of the gas in discs around pre-main-sequence stars and young stellar objects (e.g., Bujarrabal et al. Bea97 (1997); Olofsson, Liseau, & Brandeker OLB01 (2001); Thi et al. Tea01 (2001); Piétu et al. PDK03 (2003)). The infrared-to-millimetre continuum observations of such environments constrain the properties of dust grains and can be used to estimate masses and thermal structure of the objects (e.g., Bouwman et al. Bea00 (2000); Boogert, Hogerheijde, & Blake BHB02 (2002); Tuthill et al. Tea02 (2002)). Finally, experimental studies on the formation and spectra of various gaseous species (e.g., Butler et al. BLP01 (2001); Sanz, McCarthy, & Thaddeus SMT02 (2002)) as well as the composition and properties of meteoritic, cometary, and interplanetary dust together with their laboratory analogues (e.g., Chihara et al. Cea02 (2002); Mutschke et al. MPFD02 (2002)) form a basis for theoretical investigations.

On the other hand, the increase of computer power and new numerical methods stimulate the development of more sophisticated hydrodynamical models (e.g., Kley, D’Angelo, & Henning KDH01 (2001); D’Angelo, Henning, & Kley AHK02 (2002); Struck, Cohanim, & Willson SCW02 (2002)). Many of such simulations need an accurate treatment of the energy transport within the dusty medium (see e.g., Klahr, Henning, & Kley KHK99 (1999); Niccolini, Woitke, & Lopez nwl03 (2003)), which requires a detailed description of the radiative properties of matter. Consequently, the adopted opacity model is an important issue.

In this paper, we deal with physical conditions typical for protostellar nebulae and protoplanetary discs around low-mass young stellar objects. Virtually everywhere within the medium dust grains are the main opacity source, as they absorb radiation much more efficiently compared to the gas and because the temperature in these regions is low enough to prevent their destruction. However, for hotter domains ($`T\mathrm{1\hspace{0.17em}500}`$ K), where even the most stable dust materials cannot survive, it is necessary to take absorption and scattering due to molecular species into account.

Recently, several extensive models describing the properties and evolution of dust grains in protostellar cores and protoplanetary discs were proposed by Henning & Stognienko (HS96 (1996)), Schmitt, Henning, & Mucha (Sea97 (1997)), and Gail (G01 (2001)G02 (2002)).

Henning & Stognienko (HS96 (1996)) studied the influence of particle aggregation on the dust opacity in the early evolutionary phases of protoplanetary discs. They concluded that distribution of iron within the particles affects their optical properties in a great respect. Schmitt et al. (Sea97 (1997)) for the first time investigated collisional coagulation of dust grains in protostellar and protoplanetary accretion discs coupled with hydrodynamical evolution of these objects. They reported significant alteration of the thermal disc structure caused by the modification of the opacity due to dust growth. Gail (G01 (2001, 2002)) considered annealing and combustion processes leading to the destruction of silicate and carbon dust grains consistently with the evolution of a steady-state accretion disc. They found that the modification of the dust composition in the inner regions due to these processes and its consequent transport toward outer disc domains affect the opacity and, eventually, the entire disc structure.

A number of papers deal with the computation of Rosseland or/and Planck mean gas opacities in atmospheres of cool stars, protostars, and stellar winds. Alexander & Ferguson (AF94 (1994)) computed a set of opacity tables for temperatures between $`700`$ K and $`\mathrm{12\hspace{0.17em}500}`$ K for several compositions. They considered a condensation sequence111which is not to be confused with a formation sequence (see also Woitke woi2000 (2000)) of refractory materials based on chemical equilibrium calculations and took into account absorption and scattering properties of these solids as well as various gas species. Note that the actual formation process is not modelled by such calculations. Helling et al. (H00 (2000)) calculated gas mean opacities for wide ranges of density, temperature and various chemical compositions based on up-to-date spectral line lists of the Copenhagen SCAN data base and studied the importance of the molecular opacity for the dynamics of the stellar winds of cool pulsating stars. An extension of this work will be used to construct our opacity model for protoplanetary accretion discs.

In Appendix A, we give a brief overview of the most common opacity models and studies and where they have been applied. It can be clearly seen that there is a lack of papers which focus on calculations of both Rosseland and Planck mean opacities of grain and gas species for temperatures between several K and few thousands K in a wide range of densities based on both the best estimates of the dust composition and properties and recent improvements in molecular line lists.

The first goal of this paper is to define such a model. In addition, we study how Rosseland and Planck mean opacities depend on the properties of dust grains and compare them with other opacity tables. Furthermore, we investigate how different opacity models affect the hydrodynamical structure of accretion discs. Our opacity model is freely available in the Internet: http://www.astro.uni-jena.de/Laboratory/labindex.html.

The paper is organised as follows. We introduce the opacity model in Sect. 2. The influence of the grain properties on resulting Rosseland mean opacities is described in Sect. 3.1. The Rosseland and Planck mean opacities are compared to other recent opacity tables in Sect. 3.2. We study how different opacity tables affect the hydrodynamical structure of active accretion discs in Sect. 3.3. A summary follows in Sect. 4. Finally, an overview of the opacity models cited in the paper is given in Appendix A.

## 2 The opacity model

### 2.1 Dust opacities

In order to calculate frequency-dependent, Rosseland and Planck mean opacities of dust grains, we partly follow the schemes, proposed in Pollack et al. (Pea94 (1994), hereafter PHB), Henning & Stognienko (HS96 (1996), HS), and Voshchinnikov & Mathis (VM99 (1999), VM).

#### 2.1.1 Composition model and grain sizes

In this paper we follow the dust composition model for accretion discs by Pollack et al. (Pea94 (1994); Sect. 2.3) which is based on an analysis of a wide range of theoretical, laboratory, and observational dust data. This composite model has been frequently used, for instance, in the evolutionary modelling of accretion discs around young stellar objects or estimates on their mass from millimetre observations (e.g., Greaves et al. Gea98 (1998); Jura & Werner JW99 (1999); D’Alessio, Calvet, & Hartmann DCH01 (2001)).

The main dust constituents include amorphous pyroxene (\[Fe, Mg\]SiO3), olivine (\[Fe,Mg\]2SiO4), volatile and refractory organics (CHON material), amorphous water ice, troilite (FeS), and iron222See: http://www.astro.uni-jena.de/Laboratory/labindex.html. Following HS, we vary the relative iron content in the silicates considering “iron-rich” (IRS) silicates with Fe/(Fe+Mg)=0.4, “normal” silicates (NRM) with Fe/(Fe+Mg)=0.3, and “iron-poor” (IPS) silicates with Fe/(Fe+Mg)=0. However, the absolute amount of metallic iron in all these models is kept constant, which leads to the absence of solid iron in the first case and enhanced mass fraction of Fe in the third case. Such a variety of silicate models allows us to study the influence of iron content within the grain constituents on the extinction properties of dust. Another reason is that the exact mineralogical composition of the silicates in the protostellar clouds and protoplanetary discs is poorly constrained and can be different for various environments. Compared to HS, we re-estimated the absolute abundances of the silicates (FexMg1-xSiO3, Fe2xMg2x-2SiO4), iron, and troilite (FeS) in the cases of the “iron-poor” and “iron-rich” models from the Fe-stocheometric fractions keeping constant the total amount of iron. For the mass fractions of all dust constituents and their densities we follow Table $`2`$ in PHB. Note the difference between the iron mass fractions in the different dust models.

The sublimation temperatures of the grain constituents are adopted from PHB (see Table $`3`$ therein). We suppose that destruction of dust materials occurs in a narrow range of temperatures ($`1030`$ K). Given that the evaporation of the silicates and iron happens at approximately the same conditions, we do not distinguish between their evaporation temperatures and assume that they evaporate in one wide temperature range, $`100`$ K. Thus, we account for six principal temperature regions:

1. $`T155`$ K – all dust constituents are present;
2. $`165\mathrm{K}<T270`$ K – no ice;
3. $`280\mathrm{K}<T410`$ K – no ice and volatile organics;
4. $`440\mathrm{K}<T675`$ K – silicates, iron, and troilite are present;
5. $`685\mathrm{K}<T\mathrm{1\hspace{0.17em}500}`$ K – silicates and iron are present;
6. $`\mathrm{1\hspace{0.17em}500}\mathrm{K}<T\mathrm{10\hspace{0.17em}000}`$ K – gas-dominated opacities;

We take into account the dependence of the evaporation temperatures of ice, silicates, and iron on gas density (the values shown above are given for a density $`10^{10}`$ g cm-3). Note that if one considers a condensation sequence of these materials (hysteresis behaviour), it results in higher values of the condensation temperatures compared to the evaporation temperatures.

We assume that for the fifth temperature region the absolute amount of solid iron raises due to the destruction of troilite. The corresponding mass fractions of metallic iron are $`6.15\times 10^4`$, $`2.42\times 10^4`$, and $`1.29\times 10^3`$ for the NRM, IRS, and IPS silicate models, respectively.

As for the size distribution, we apply a modified MRN (Mathis et al. MRN (1977)) function suggested by Pollack, McKay, & Christofferson (P85 (1985)). The modification consists of the inclusion of large ($`0.5\mu `$m – $`5\mu `$m) dust grains. Such particle growth is expected to proceed efficiently at the early phases of the protoplanetary disc evolution due to the coagulation of small dust grains (see, for instance, Mizuno, Markiewicz, & Voelk MMV88 (1988)). We do not consider other size distributions since the overall effect of the particle sizes on the dust opacity is well studied (e.g., Pollack et al. P85 (1985); Beckwith, Henning, & Nakagawa BHN00 (2000)). It further allows us to compare directly our results with other works.

#### 2.1.2 Grain structure and topology

It becomes evident from theoretical investigations and laboratory experiments that the dust agglomeration is an efficient process in dense and relatively cold environments, like protostellar cores or protoplanetary discs (e.g., Kesselman K80 (1980); Nuth & Berg NB94 (1994); Ossenkopf & Henning OH94 (1994); Kempf, Pfalzner, & Henning KPH99 (1999); Wurm & Blum WB98 (1998)WB00 (2000); Blum et al. Bea02 (2002)). Agglomeration leads to the formation of irregular particles consisting of hundreds or thousands of tiny subgrains. Usually, dust aggregates of two extreme kinds are considered, depending on the assumed coagulation processes, namely, PCA (particle-cluster aggregation) and CCA (cluster-cluster aggregation). As the laboratory and theoretical studies reveal, the PCA aggregates are sphere-like particles having a fractal dimension of about $`3`$. They have a compact “core” and a more rarefied “mantle”. The CCA process results in the formation of very filamentary grains with complicated structure. They have fractal dimension of roughly $`2`$ (Stognienko, Henning, & Ossenkopf SHO95 (1995)).

During the evolution of parent objects, like molecular clouds or protostellar discs, chemical and physical processes can further modify the properties of dust grains. For instance, accretion of volatile materials on dust surfaces and their subsequent chemical processing are efficient in outer regions of protoplanetary discs and in protostellar clouds (e.g. Greenberg G67 (1967); Brown, Charnley, & Millar BCM88 (1988); Hartquist & Williams HW90 (1990); Hasegawa & Herbst HH93 (1993); Willacy, Rawlings, & Williams WRW94 (1994); Aikawa et al. Aea99 (1999)). This results in well-defined “core-mantle” or, more probably, “onion-like” grain structure. In protostellar discs, dust can be transported by the accretion flow toward hotter regions, where their volatile mantle materials evaporate, and sputtering, annealing, combustion, and crystallisation processes may change the structure, composition, and sizes of the grains (Gail G01 (2001, 2002)). Therefore, it seems obvious that the real astronomical grains must have a very complicated structure and topology.

Unfortunately, modern computational methods and facilities allow only the consideration of somewhat simplified (but still reasonable) types of dust grains. In the present study, we focus on the following particle types:

1. Homogeneous and composite aggregates;
2. Homogeneous, composite, and porous composite spherical particles;
3. Multishell and porous multishell spherical particles.

An aggregate dust particle is assumed to be a cluster of small spherical subgrains sticking together. A particle is called “homogeneous” if it consists of only one dust component. On the contrary, “composite” means that a particle incorporates a fine mixture of various materials (heterogeneous particle). In addition to these two extreme cases of the chemical dust structure, we consider a model of “multishell” grains, where each particle includes all constituents distributed within concentric spherical shells. To study the influence of the porosity on the extinction properties of dust grains, we fill composite and multishell spherical particles with vacuum. It is reasonable to assume that the optical behaviour of these porous multishell and porous composite particles may resemble that of more realistic kinds of dust grains.

#### 2.1.3 Computational methods

The aggregate model and the numerical method to compute the optical properties of coagulated particles are adopted from HS. The aggregates are assumed to be in the form of PCA (50%) and CCA (50%) particles consisting of $`0.01\mu `$m spherical subgrains. The spectral representation of inhomogeneous media (Bergman B78 (1978)) and the discrete multipole method (DMM) are elaborated to calculate effective dielectric functions of aggregates (Stognienko et al. SHO95 (1995)). At first, we use the DMM to calculate the spectral function of the aggregated particles of a special topology, when their subgrains touch each other only at one point. Then, we account for the interaction strength between the subgrains (percolation), which varies with the size and type of the aggregates, by an analytical expression (see Equation $`13`$ in HS). Finally, the effective dielectric functions of the aggregates are obtained by the spectral representation (HS, Eq. $`5`$). In the case of the composite aggregates, we compute the optical constants of the composite material by the Bruggeman effective medium theory (Bruggeman B35 (1935)), generalised to many components by Ossenkopf (O91 (1991)). The optical properties of the dust aggregates are calculated with the usual Mie theory. It should be noted that this numerical approach is valid in the static limit only, which means that the scale of inhomogeneities within the particles must be small compared to the wavelength. Given the $`0.01\mu `$m size of the subgrains and the shortest wavelength of $`0.1\mu `$m we considered, this condition is fulfilled.

We use the numerical approach of VM to model (porous) composite and (porous) multishell dust particles. In this method, the composite grains are represented as spheres with many concentric shells, where each shell includes several layers of randomly distributed dust materials. The multishell grains are modelled exactly as the composite ones but each shell includes only one layer of a dust constituent. Then a generalised multilayered Mie theory can be applied to calculate their optical properties. As it has been shown by VM, a convergence in the optical behaviour of the multishell particles is achieved if the number of shells exceeds $`3`$ and dust materials within shells are randomly ordered. In our calculations, we found that this number must be at least $`20`$, because the highly absorbing material troilite is used, which induces interference within the shells and prevents fast convergence. Thus, in our case, a typical composite grain is modelled as a spherical particle with about hundred shells. On the contrary, a multishell grain is represented by a spherical particle with only a few shells.

We modify somewhat the dust model for the case of multishell and composite spherical particles. We mixed the silicates and iron into one material with the Bruggeman mixing rule. A similar mixture of silicates, sulphides, and metals (GEMS, Glass with Embedded Metals and Sulphur) is found to be a common component of interplanetary dust particles (Rietmeijer & Nuth RN00 (2000)). However, the main reason for this change is the convergence failure of the applied numerical method for the case of multishell grains with iron layers. This is due to numerical uncertainty which arises in the calculation procedure for the Mie scattering coefficients in the case of highly absorbing materials, like iron (for further explanation, see Gurvich, Shiloah, & Kleiman GSK01 (2001)).

We assume that each dust component has a total volume fraction in a particle according to its mass fraction and density, as specified in Table 1. For example, for the first temperature region and in the case of IPS silicate mineralogy, the mixture of iron and silicates occupies $`8.9`$%, troilite – $`1.6`$%, refractory organics – $`23.4`$%, volatile organics – $`6`$%, and water ice – $`60`$% of the entire particle volume, respectively. These values are similar in the case of NRM and IRS models. Thus, if the temperature is low, organics and ice are the dominant components of the dust grains.

Unlike to the case of the composite particles, in the case of multishell spherical grains it is assumed that the distribution of dust materials is not random but follows their evaporation sequence. Thus, for the first temperature region in the protoplanetary disc, multishell spherical grains consist of a refractory core made of a mixture of silicates and iron and subsequent shells of troilite, refractory organics, volatile organics, and water ice. For higher temperature ranges, the number of shells is smaller since some materials are evaporated. In total, the number of shells in the case of multishell spherical particles varies from $`2`$ to $`5`$. For the fifth temperature range ($`T>700`$ K), where all troilite is converted to solid iron, we let this iron form an additional layer on the grain surface. It is an extreme and probably physically unjustified case, but this allows us to study the influence of the formation of a highly absorbing surface layer on resulting opacities.

The porosity of particles is treated in a simple manner by the addition of $`50\%`$ of vacuum (by volume) inside. For the case of the porous composite spheres, we consider vacuum to be one of the grain constituents, which forms additional ”empty” layers. In the case of the multishell grain model, we mix each shell with vacuum in the same way which is applied to create the composite spherical particles. That is, we subdivide each individual shell in many layers and fill some of them with vacuum according to the requested porosity degree.

With the two computational approaches described above, we calculated the ensemble averaged absorption and scattering cross section as well as albedo and the mean cosine of the scattering angle for all kinds of dust grains. Applying Eqs. $`1`$-$`5`$ from Pollack et al. (P85 (1985); see Table 1 of the present paper), the dust monochromatic opacity and, consequently, the Rosseland and Planck mean opacities were obtained for temperatures below roughly $`\mathrm{1\hspace{0.17em}500}`$ K and a density range between $`10^{18}`$ g cm-3 and $`10^7`$ g cm-3. A convenient analytical representation of the Rosseland and Planck mean opacities for every temperature region is provided as a $`5`$-order polynomial fit. This representation allows to calculate the Rosseland and Planck mean dust opacities accurately ($`1`$%) and quickly for any given temperature and density values within the model applicability range, which is important for computationally extensive hydrodynamical simulations. The corresponding fit coefficients can be found in the code333http://www.astro.uni-jena.de/Laboratory/labindex.html.

### 2.2 Gas opacity

The opacity of the very inner part of the protoplanetary accretion disc is dominated by various gaseous species. Here, the temperature is too high for dust to be present.

Compared to the calculation of dust opacities, the calculation of accurate Rosseland and Planck mean gas opacities is more challenging due to the large variation in frequency, temperature, and density of the absorption coefficient of numerous molecules, atoms, and ions. In addition, the body of data to be handled easily amounts to several millions of absorption lines per molecule.

Missing data for absorption lines are critical for the calculation of Rosseland mean gas opacity since it is dominated by transparent spectral regions due to the harmonic nature of the averaging process. Therefore, each Rosseland mean is always only a lower limit of the correct value. The opposite is true for the case of the Planck mean opacity – missing data for weak lines or bands cause an overestimation of the strong lines. Therefore, a Planck mean is always an upper limit of the case of ideally complete data.

The dust opacity model for protoplanetary accretion discs outlined in the previous sections is supplied by a new table of gas opacities assembled on the basis of Helling (H\_PhD99 (1999); Copenhagen SCAN data base) and Schnabel (S\_MT01 (2001); HITRAN data base). The gas opacity model is outlined in Helling et al. (H00 (2000)) and only a short summary is given here. The Rosseland and the Planck mean opacities are calculated from opacity sampled lines lists. The data for the line absorption coefficients used in Helling et al. (H00 (2000)) (CO - Goorvitch & Chackerian CO (1994); TiO - Jørgensen TiO (1994); SiO - Langhoff & Bauschilder SiO (1993); H2O - Jørgensen & Jensen H2O (1993); CH - Jørgensen et al. CH (1996); CN, C2 \- Jørgensen & Larsson CN (1990); C3 \- Jørgensen C3 (1989); HCN, C2H2 \- Jørgensen C2H2 (1990)) were supplemented by data from the HITRAN 96 database (CH4, NH3, HNO3, H2CO, CO2, N2O, O3, SO2, NO2, HO2, H2, O2, NO, OH, N2). The opacity sampling of the latter was carried out by Schnabel (S\_MT01 (2001)). The set of continuum opacities and scattering includes continuum absorption from H i (Karzas & Latter Hi (1961)), H- (John H- (1988)), H$`+`$H (Doyle H&H (1968)), H$`{}_{}{}^{}{}_{2}{}^{}`$ (Somerville H2- (1964)), H$`{}_{}{}^{+}{}_{2}{}^{}`$ (Mihalas H2+ (1965)), He- (Carbon et al. He- (1969)), He i, C i, Mg i, Al i, Si i (all from Peach P70 (1970)) as well as Thompson scattering on free electrons and Rayleigh scattering. Collision-induced absorption has been considered for H2-H2 and H2-He according to Borysow et al. (bjz97 (1997)).

The number densities of the ions, atoms and molecules are computed from an updated chemical equilibrium routine, including 14 elements and 155 molecules based on the JANAF table (electronic version of Chase et al. Ch85 (1985); for more detail see Helling et al. H00 (2000)). The element abundances are chosen mainly according to Anders & Grevesse (AG89 (1989)), but have been updated for various elements (see Helling et al. H00 (2000)).

Using this approach, the Rosseland and Planck mean gas opacities were computed for temperature ranges between $`500`$ K and $`\mathrm{10\hspace{0.17em}000}`$ K and for gas densities between $`10^{18}`$ g cm-3 and $`10^7`$ g cm-3. In contrast to the dust opacity, no simplified analytical expression can be found for the gas mean opacities because of their sensitive dependence on temperature and density. Thus, we apply a second-order interpolation in order to calculate the gas opacities for any given temperature and density value from tabulated values.

### 2.3 Opacity table

In order to assemble the opacity table, we take into account either only dust opacity data for low temperatures, $`T\mathrm{1\hspace{0.17em}500}`$ K, or only gas opacity data for higher temperatures. As has been shown by many authors (e.g., Lenzuni, Gail, & Henning L95 (1995)), it is an accurate approach because dust dominates the absorption properties of matter whenever it is present. It has been shown that already 30% of total SiO content in the solid phase is enough to exceed even the Planck mean gas opacity (Helling H\_PhD99 (1999)). This fraction will decrease with increasing impurity of the grains. However, for the dust-to-gas transitional region ($`T`$ from $`\mathrm{1\hspace{0.17em}400}`$ K to $`\mathrm{1\hspace{0.17em}600}`$ K), where the last dust grain population evaporates, it is necessary to calculate the opacity of dust and gas simultaneously. In this narrow temperature range the resulting Rosseland and Planck mean opacities are certainly going down by few orders of magnitude, so one can apply a simple linear interpolation to estimate their values with a good accuracy. Still, we note that our table may give only approximate opacity data for that temperature range.

## 3 Results

In this section, we present the opacities of all types of dust grains and discuss the differences between them. Second, we compare the Rosseland and Planck mean opacities with other recently published opacity models. Finally, we study the influence of the adopted opacity model on the accretion disc structure.

### 3.1 Opacity and the dust models

The Rosseland mean opacities $`\kappa _\mathrm{R}`$ computed for all dust models are presented in Fig. 1. We compared two silicate models, namely, the IPS (left panel) and the IRS (right panel). Shown are the following dust models: multishell spheres, composite aggregates, composite spheres, homogeneous aggregates, homogeneous spheres, porous multishell grains, and porous composite spheres.

The most prominent trends in behaviour of the Rosseland mean dust opacities can be summarised as follows:

1. There is a significant difference in the calculated dust opacity values between the aggregates, (porous) composite and (porous) multishell spherical particles for the first ($`T150`$ K) and fifth ($`T700`$ K) temperature regions;
2. For intermediate temperatures, $`150\mathrm{K}T700`$ K, the resulting Rosseland mean opacities do not show profound variations with the applied dust models, except for the case of the IPS homogeneous aggregates;
3. The discrepancy between opacity curves are smaller for the case of the IRS silicate model compared to the IPS;
4. Addition of vacuum inside compact composite and multishell spherical particles leads to significant modification of their opacities.

In the first temperature region, organics and water ice are the dominant dust materials according to our compositional model (see Table 1). At that temperatures ($`T150`$ K), the main contribution to the Rosseland mean opacity comes from monochromatic opacities in the wavelength range about $`30\mu `$m – $`400\mu `$m. In this spectral range, iron and troilite have higher values of the real ($`n_\lambda `$) and imaginary ($`k_\lambda `$) indices of refraction compared to the other materials considered in our model. For instance, at $`\lambda =100\mu `$m, iron has $`n_{100\mu \mathrm{m}}^{\mathrm{Fe}}=95.02`$, $`k_{100\mu \mathrm{m}}^{\mathrm{Fe}}=181.95`$, troilite has $`n_{100\mu \mathrm{m}}^{\mathrm{FeS}}=8.5`$, $`k_{100\mu \mathrm{m}}^{\mathrm{FeS}}=0.73`$, whereas organics and water ice have $`n_{100\mu \mathrm{m}}^{Org}=2.14`$, $`k_{100\mu \mathrm{m}}^{Org}=0.15`$ and $`n_{100\mu \mathrm{m}}^{Ice}=1.82`$, $`k_{100\mu \mathrm{m}}^{Ice}=0.05`$, respectively. Furthermore, for troilite and especially for iron these values rapidly increase with wavelength. Thus, the resulting optical properties of a dust grain become sensitive to the absolute amount of Fe and FeS and its distribution inside the particle. In such a case, different theories to calculate optical properties of an aggregate particle may give different results and should be used with caution (see Stognienko et al. SHO95 (1995), Fig. 6, 8).

We considered aggregated particles of two kinds, namely, composite and homogeneous aggregates. The optical constants of the composite material do not show a peculiar behaviour at long wavelengths. Moreover, the refractive index of this composite changes only little if one is switching from the IPS to the IRS silicate model. One reason is that the amount of solid Fe is small compared to other constituent materials. Another reason is that the total iron abundance is kept constant in all compositional models. Hence the resulting optical properties of the composite aggregates are not very sensitive to the actual topology of the particles and adopted silicate model. Indeed, as it can be clearly seen in the Figure, the dust opacity values of the composite aggregates (triangles) are close to the $`\kappa _\mathrm{R}`$ of the homogeneous spheres (dot-dashed line). The maximum deviation of these opacity curves is achieved in the case of the IPS silicates for temperatures higher than $`700`$ K. Here, the absolute amount of metallic iron is increased due to conversion of FeS to Fe. In addition, the Rosseland mean opacities of the composite aggregates do not depend much on the adopted silicate model (compare triangles on the left and right panel).

On the contrary, the $`\kappa _\mathrm{R}`$ values for the homogeneous aggregates do demonstrate a strong variations with the compositional model. In the case of the ”iron-rich” silicates (right panel) the Rosseland mean values of the homogeneous aggregates (dashed line) lie much closer to the $`\kappa _\mathrm{R}`$ of the composite aggregates (triangles) then for the ”iron-poor” compositional model (left panel). As it has been shown already by HS, this is caused by the presence of bare iron aggregates in the case of the IPS homogeneous aggregate model. Due to the extremely high absorptivity induced by the strong interactions between individual aggregate sub-grains, the optical properties of such iron clusters determine the overall behaviour of the resulting opacities. Note, that in the case of the ”iron-rich” silicate composition, all iron is locked inside silicates and the absolute amount of troilite is reduced by a factor of two compared to the IPS and NRM models. In the absence of a population of highly absorbing grains, the Rosseland mean values of the homogeneous aggregates are close to that of the composite aggregates. However, at temperatures higher than $`700`$ K, troilite is converted to solid iron which increases the $`\kappa _\mathrm{R}`$ values of the homogeneous aggregates in respect to that of the composites (compare dashed lines with triangles on the right panel of Fig. 1).

The situation is different for the case of the (porous) composite and (porous) multishell spherical particles. As it was mentioned in Sect. 2.1.3, we changed the compositional model in this case and locked all solid iron in silicates using the Bruggeman rule of the effective medium theory. Thus, the only material with a high absorptivity at long wavelengths which remains in all compositional models is troilite. The metallic iron is another highly-absorbing dust component, but it is present in the fifth temperature region only, at $`T>700`$ K.

Voshchinnikov & Mathis (VM99 (1999)) have shown that conventional EMTs are rather approximate in the case of small composite spherical particles, $`x=2\pi a/\lambda 1`$, when one of the grain constituents has a large refractive index (see Fig. 3b therein). Given a typical size $`a=0.1\mu `$m of our dust grains, a typical wavelength of $`\lambda =100\mu `$m for the first temperature region, and the high refractive index of troilite at that wavelength, this condition is fulfilled. Thus, we adopted the approach of VM to model a composite grain as a spherical particle with many concentric layers and applied it for all temperature regions. Note that there is always an interference between the layers, which makes the optical properties of such composite particles to be different from that of composite grains with well mixed dust constituents. This is especially true if one of the dust components has a high absorptivity. In our compositional model for layered spherical particles, this is troilite in the first temperature region while in the fifth region it is iron. Therefore, we may expect to see the difference between the Rosseland mean opacity values of the composite spheres and composite aggregates, particularly for $`T150`$ K and $`T>700`$ K.

As it can be seen in Fig. 1, indeed the $`\kappa _\mathrm{R}`$ values in the case of the composite spherical particles are much higher than for the composite aggregates. For instance, at $`T=10`$ K this difference can reach factors of 10 and 20 for the IPS and IRS compositions, respectively (compare circles and triangles on the left and right panels). For higher temperatures, the Rosseland mean opacity curves of the composite spheres and aggregates are close to each other till $`T700`$ K is reached. In this temperature region, we assumed that iron forms a layer on the surface of the composite spherical grains. Such a layer ”screens out” all underlying materials and totally dominates the optical behaviour of the entire particle. Due to this fact, the dust opacity values of the composite spheres in the fifth temperature region is nearly the same for both the IPS and IRS compositions.

The multishell spherical particles have a restricted number of layers compared to the composite spheres, namely from 2 to 5, depending on the temperature region. The troilite layer is assumed to be the first layer after the silicate core and thus troilite is ”hidden” inside. It prevents a strong interference between the consequent particle shells as it is the case for the composite spheres. Then one may expect that the $`\kappa _\mathrm{R}`$ values of the multishell spheres should be lower than that of the composite spherical particles, especially for the IPS model. As it can be clearly seen in Fig. 1, this is true for $`T150`$ K, whereas for higher temperatures both opacity curves almost coincide (compare solid line with circles). Hence, the actual distribution of dust constituents within a multilayered spherical particle is not that important for the relevant Rosseland mean opacities at $`T150`$ K.

The addition of vacuum inside the compact composite and multishell spherical grains leads to a significant increase of the corresponding $`\kappa _\mathrm{R}`$ values for the first and fifth temperature regions (compare solid line with dotted line and circles with pluses). The first reason is that the density of the porous grains becomes lower than the density of the compact particles. Second, for the porous spheres the relevant extinction efficiencies are higher compared to that of the compact spherical particles if some of the dust constituents have a particularly high absorptivity, like troilite in the first and iron in the fifth temperature region, respectively. This is due to the coherence between the particle layers. The interference is more intense for the case of the composite sphere since it has more concentric layers and a nearly homogeneous distribution of the dust constituents from the centre to the surface compared to the multishell spherical particle. Note that in the fifth temperature region, both composite and multishell spheres have a similar topology, namely, a silicate core covered by iron mantle. Therefore, it is naturally to expect that they have a similar behaviour of the resulting Rosseland mean opacities.

As it can be clearly seen in Fig. 1, the $`\kappa _\mathrm{R}`$ values of the porous composite spheres are higher than that of the compact composite spherical particles at $`T<150`$ K and $`T>700`$ K by a factor of 2 for both the IPS and IRS compositions (compare pluses and circles). For other temperatures, the corresponding Rosseland mean curves lie close to each other. The situation is similar for the case of the multishell and porous multishell spherical particles (compare solid line with dotted line). As we expected, in the fifth temperature region the Rosseland mean opacities for the case of the porous composite and porous multishell spheres have almost the same values.

### 3.2 Comparison to other studies

In Fig. 2, the Rosseland (left panel) and Planck (right panel) mean opacities composed of the NRM composite aggregate for the low temperature range and gas opacities for the high temperature range are compared with other models. We plotted these values for a wide temperature range, $`T`$ from $`10`$ K to $`10^5`$ K and for gas densities which scale as $`10^{19}\times T^3`$ g cm-3. It allows us to make a comparison between the models in a wide temperature as well as density range simultaneously.

As it is clearly seen in Fig. 2 (left panel), the discrepancy between the Rosseland mean values provided by various models is negligible at high temperatures ($`T\mathrm{3\hspace{0.17em}000}`$ K). On the contrary, the Planck mean opacity values differ by a few order of magnitude in this temperature range (compare curves on the right panel). The reason is that the Rosseland mean is much less sensitive to differences in the material data than the Planck mean due to the nature of the averaging process. However, the Planck mean heavily depends on the adopted values of the band and line strengths which vary for different line lists and on the adopted chemical equilibrium constants444Note that the chemical equilibrium constants, $`K_p`$, used by different authors can cause differences in the resulting opacity values because it affects the number density of species. The same effect will be caused by the neglect of the metal ions in the chemical equilibrium calculations (see, e.g., Helling et al. H00 (2000)). It is, however, not very straightforward to decide which line list is the most correct (for a discussion see Jørgensen ugj03 (2003)). We nevertheless dare to demonstrate the difficulties arising from comparing Planck mean opacities calculated by different authors because it has – to our knowledge – not been pointed out clearly in the literature. Only comparisons of Rosseland mean opacities are presented, e.g. by Alexander & Ferguson (AF94 (1994)).

Therefore, the $`\kappa _\mathrm{P}`$ of the Opacity Project (OP, crosses) are much larger than all the other opacity models since it combines more atomic opacity sources (see e.g. Table $`3.3`$ in Helling H\_PhD99 (1999)). This model does not contain molecules for temperatures $`\mathrm{1\hspace{0.17em}500}`$ K – $`\mathrm{5\hspace{0.17em}000}`$ K, which become important absorbers in this temperature range. The difference between the Planck mean opacity values in the case of the AF model (dashed line) and our model (solid line) may be caused by different molecular line data and a different set of chemical species adopted. The same is true for the $`\kappa _\mathrm{P}`$ values of Sharp (Sharp92 (1992), diamonds), which lay somewhat in between the Planck mean opacities of our model (solid line) and the values of Alexander & Ferguson (dashed line).

For temperatures lower than about $`\mathrm{1\hspace{0.17em}500}`$ K, dust grains are the main opacity source. As it has been shown by Pollack et al. (Pea94 (1994)), in this case the difference between the Rosseland and Planck mean opacities computed for the same model is small, $`30`$% (see Fig. 4b therein). The reason is that both opacities are dominated by continuum absorption and scattering rather than absorption lines in this case. In what follows, we focus on the low-temperature Rosseland mean opacities only ($`T\mathrm{3\hspace{0.17em}000}`$ K).

The Rosseland mean opacity values $`\kappa _\mathrm{R}`$ computed by the model of Bell & Lin (BL94 (1994), BL) strongly deviate from the $`\kappa _\mathrm{R}`$ calculated by other models. For example, this difference can reach a factor of hundred for $`T\mathrm{1\hspace{0.17em}500}`$ K – $`\mathrm{1\hspace{0.17em}800}`$ K (compare dot-dashed and solid lines in Fig. 2). The reason is that the BL data are a modification of the old analytical Lin & Papaloizou (L80 (1980)) opacities, which are based on the opacity tables of Alexander (A75 (1975)) and Cox & Stewart (1970a ; 1970b ) supplemented by the data from Alexander, Augason, & Johnson (AAJ89 (1989)) for $`T<\mathrm{3\hspace{0.17em}000}`$ K. The model of Bell & Lin includes dust grains of two types, namely, homogeneous icy and metallic spherical grains. Since we used the more advanced dust compositional model of PHB and a different set of evaporation temperatures, this causes lower $`\kappa _\mathrm{R}`$ values in the case of the BL model compared to the other models for $`T\mathrm{1\hspace{0.17em}500}`$ K. For higher temperatures, $`\mathrm{1\hspace{0.17em}500}`$ K – $`\mathrm{3\hspace{0.17em}000}`$ K, this deviation is huge, $`100`$ times. As it has been shown by AF, the reason is that Bell & Lin truncated monochromatic opacities of water, which is one of the main absorbers at such temperatures, at a too short wavelength in their calculations. The missing opacity data affect the resulting Rosseland mean opacity. They do, so far, not draw any conclusions regarding Planck mean opacities.

In overall, the opacity curves of all other models do not show such a strong difference between each other (compare triangles, dashed line, solid line, and line with open squares). Our dust model differs from the model of Pollack et al. by taking into account an aggregate nature of cosmic dust grains and a new set of optical constants, but dust size distribution, composition, and evaporation temperatures are the same. However, the difference between the Rosseland mean opacity values of these two model can reach about a factor of two (see also Fig. 5a in Henning & Stognienko HS96 (1996)).

The $`\kappa _\mathrm{R}`$ values of the opacity model of Alexander (A75 (1975), line with open squares) is lower for the dust-dominated temperature region ($`T\mathrm{1\hspace{0.17em}500}`$ K) and higher for the dust-to-gas transitional region ($`\mathrm{1\hspace{0.17em}500}\mathrm{K}T\mathrm{1\hspace{0.17em}700}`$ K) compared to our model (solid line) by factors of 5 and 100, respectively. The reason is, as mentioned above, that we used an approximation to compute $`\kappa _\mathrm{R}`$ in that temperature region, where the last dust grains get evaporated, which is not a very accurate approach. On the contrary, the opacity model of Alexander assumes the presence (but no evaporation!) of small $`0.1\mu `$m spherical silicate grains in a rather approximate way, assuming that all dust is homogeneously condensed when the gas becomes saturated. Since here a phase-transition takes place, a supersaturated gas would be needed which results in higher molecular abundances than AF derive from their equilibrium consideration (for a discussion see e.g. Woitke & Helling wh2003 (2003)). Moreover, the condensation begins at higher temperature than the value we assumed for the evaporation of the last grain constituents because of the dust hysteresis. Finally, the model of AF neglects the presence of other refractory materials, like iron, in the dust-dominated temperature range, which makes the relevant Rosseland mean opacity values lower than provided by our model.

On the contrary, the model of AF does consider several dust materials, namely, iron, silicate, carbon, SiC assumed to be present as small ellipsoidal grain particles. Therefore, the corresponding Rosseland mean dust opacities nearly coincide with our values (compare dashed and solid lines, respectively). The same is true even for the dust-to-gas transitional region, where our opacity model gives rather approximate opacity values.

### 3.3 Opacity and disc structure

We compare the thermodynamical structure of a typical protoplanetary disc around a low-mass star computed with two different opacity tables in Fig. 3. The $`1+1`$D model of an active steady-state accretion disc of Ilgner (I03 (2003)) was used with the following input parameters: $`M_{}=1M_{\mathrm{}}`$, $`\dot{M}=10^7M_{\mathrm{}}`$ yr-1, and $`\alpha =0.01`$. Here, $`M_{}`$ is the stellar mass, $`\dot{M}`$ is the mass accretion rate, and $`\alpha `$ is the parameter describing the kinematic viscosity. Note that this model incorporates a star as a gravitational center only and does not take into account the effect of the stellar radiation on the disc structure.

The thermal structure shown on the left panel was obtained with the Rosseland mean opacity table of Bell & Lin (BL94 (1994)). On the right panel, we present the same disc structure but for the case of the IPS homogeneous aggregate model (IPSHA) supplemented by the gas mean opacity. We choose these two opacity models as the overall difference between them is the largest among the different models (compare solid and dash-dotted lines in Fig. 2).

It can be clearly seen that the higher values of the Rosseland mean opacity in the case of the IPSHA model leads to a hotter and more extended disc structure. For instance, the vertical scale height of the disc at $`20`$ AU is equal to $`1.3`$ AU for the former and $`2.1`$ AU for the latter opacity models, respectively. Consequently, there is also a variation of the density structure between the models, namely, the disc density is higher for the BL opacity model compared to that of the IPSHA model. The temperature difference is also prominent. For example, the midplane temperature of $`100`$ K, which roughly corresponds to the ice melting point, is reached at $`6.5`$ AU for the model of Bell & Lin, whereas in the case of the IPSHA model it is at $`8`$ AU.

To confirm our findings, we did a similar comparison with another code. We used a full 2D hydrodynamical code designed to simulate the interaction of the protoplanetary disc with a protoplanet (D’Angelo DA01 (2001)). The parameters of the model were as follows: $`M_{}=1M_{\mathrm{}}`$, $`M_\mathrm{d}=0.01M_{\mathrm{}}`$, $`\nu =10^{15}`$ cm2 s-1, $`\mu =2.39`$, and $`\gamma =1.4`$, where $`M_\mathrm{d}`$ is the total disc mass, $`\nu `$ is the kinematic viscosity, $`\mu `$ is the mean atomic weight of the gas, and $`\gamma `$ is the adiabatic exponent.

The Bell & Lin opacity model together with IPS and IRS homogeneous aggregate models were chosen for the comparison. The midplane temperature for all three opacity models is shown in Fig. 4. It can be clearly seen that the difference in the temperature values between all models can reach about $`50`$% for the disc radii $`r_\mathrm{d}2`$ AU. Note that it disappears at larger distances, $`r_\mathrm{d}10`$ AU. As expected, in the case of the BL model the midplane temperature is the lowest almost everywhere (dashed line), whereas for the IPS homogeneous aggregates it is the highest (solid line) and the IRS temperature values lie between them (dotted line). The reason is the same as for the case of the $`1+1`$D disc model, namely, lower opacity values of the BL model compared to both the IPS and IRS opacities and higher values of IPS opacity in comparison with those of BL and IRS models. It is interesting that at $`T130`$ K ($`r2`$ AU) all the temperature curves are very close to each other. This is due to the adopted ice evaporation temperatures, which is a little lower in our case compared to the model of Bell & Lin. That leads to nearly the same opacity values for all three models at a restricted set of temperatures around $`130`$ K (compare solid and dot-dashed lines in Fig. 3). For $`r10`$ AU the temperature curves are close to each other because the corresponding opacity models have rather similar low-temperature opacity values.

Thus, we showed that the difference in the Rosseland opacity tables applied in hydrodynamical calculations leads to deviations in the resulting disc structure. As it has been shown by Markwick et al. (Mea02 (2002)), the temperature distribution is a crucial factor for the chemical evolution in the inner parts of accretion discs. Therefore, proper opacity modelling is an important issue in order to follow the dynamical and especially chemical evolution of protoplanetary accretion discs.

## 4 Summary

We have compiled the Rosseland and Planck mean dust and gas opacities for the temperature range $`T[5,10^4]`$ K and for gas densities $`\rho [10^{18},10^7]`$ g cm-3 which is appropriate for the conditions in protoplanetary discs. The absorption and scattering due to dust grains of different compositions, shapes and topological structures and the absorption provided by over $`30`$ atoms, molecules, ions were taken into consideration. The corresponding well-documented numerical code together with representative data and figures are electronically available: http://www.astro.uni-jena.de/Laboratory/labindex.html.

It has been found that the topological distribution of highly absorbing materials, such as iron and troilite, within dust grains dominates the relevant optical properties to a large degree. Particularly, the composite and porous composite spherical grains modelled in a special way have remarkably high Rosseland mean opacities at $`T150`$ K and $`T700`$ K even in comparison with the composite aggregated particles. We have shown that at intermediate temperatures the Rosseland mean opacities of distinct grain models are close to each other. It has been demonstrated that the difference between the opacity values of various dust models is smaller in the case of the compositional model with a smaller amount of solid iron and troilite. We found that porous composite and porous multishell spherical particles show higher opacity values in comparison with their compact analogues.

We performed a comprehensive comparison of our results with other recent opacity models. We found a significant difference between the opacity models in the case of the Planck mean and a good agreement between them for the case of the Rosseland mean at $`T\mathrm{1\hspace{0.17em}500}`$ K, where gas species are the main opacity sources. For lower temperatures, where opacities are dominated by dust grains, there is a discrepancy (a factor of $`3`$ at most) in both the Rosseland and Planck mean values for all considered models.

We demonstrated that differences in the Rosseland mean opacity values provided by distinct opacity models affect the hydrodynamical structure of active steady-state accretion discs. Namely, higher values of the Rosseland mean opacity lead to a hotter and more extended disc structure in the case of $`1+1`$D and $`2`$D disc modelling.

###### Acknowledgements.

DS was supported by the German *Deutsche Forschungsgemeinschaft, DFG* project “Research Group Laboratory Astrophysics” (He 1935/17-1). The work of MI was supported by the Max Planck Society. We are grateful to D. Alexander for valuable comments and discussions. For the calculations of the optical properties of multishell spherical particles, we used a code by N.Voshchinnikov (http://www.astro.spbu.ru/staff/ilin2/ilin.html)

## Appendix A An overview of the opacity models

The goal of this compilation is to provide essential information about opacity models cited in the text. We show for which elemental composition they are developed, what kind of opacities they supply, and in what temperature and density ranges they work. In addition, key references to the studies, where these opacity models have been applied, are given and the aims of the investigations are briefly mentioned.

Primarily, we distinguish between two kinds of opacity models. The first kind is designed for stellar evolutions, where it is more convenient to use a special parameter, $`R=\rho /T_6^3,T_6=T/10^6`$ K instead of gas density $`\rho `$ (see discussion in Rogers & Iglesias RI92 (1992)). Thus, opacity data of such models are assembled in rectangular tables in R–$`T_6`$ space. However, in hydrodynamical simulations of accretion discs, it is more convenient to have opacity tables composed in $`\rho `$$`T`$ parameter space, for which we refer as to the second type of the models. The direct conversion of R to $`\rho `$ in the opacity models of the first kind leads to trapezoidal tables, where opacity values for different temperatures have different density intervals. We mark such models with “a” and show the maximum limits of density and temperature values for them.

In addition, we point out if opacity data are available on-line in the Internet or via E-mail (“b”). By default, the papers are supposed to contain the opacity data in a tabular form. Otherwise, when only opacity plots are available, we label the corresponding papers with “c”. The analytical opacity models are marked with “d”.


format: text

text:
# 1 Introduction

## 1 Introduction

The study of vertex models have led to important developments in the field of exactly solvable models in two dimensions . Their transfer matrices are in general constructed from local Boltzmann weights $`_{𝒜i}(\lambda )`$ where $`\lambda `$ denotes a spectral parameter. This operator can be viewed as a matrix on the space of states $`𝒜`$ representing, for instance, the horizontal degrees of freedom of the vertex model on the square lattice. Its matrix elements are operators on $`{\displaystyle \underset{i=1}{\overset{L}{}}}V_i`$ where $`V_i`$ represents the space of states of the vertical degrees of freedom at each site $`i`$ of a chain of size $`L`$. The corresponding transfer matrix can be expressed in terms of an ordered product of $`_{𝒜i}(\lambda )`$ operator over the space $`𝒜`$ denominated monodromy operator $`𝒯_𝒜(\lambda )`$ ,

$$𝒯_𝒜(\lambda )=_{𝒜L}(\lambda )_{𝒜L1}(\lambda )\mathrm{}_{𝒜1}(\lambda )$$
(1)
In terms of the monodromy matrix, a sufficient condition for integrability is the Yang-Baxter algebra  which reads

$$R(\lambda ,\mu )𝒯_𝒜(\lambda )𝒯_𝒜(\mu )=𝒯_𝒜(\mu )𝒯_𝒜(\lambda )R(\lambda ,\mu )$$
(2)
where $`R(\lambda ,\mu )`$ is an invertible matrix over complex numbers acting on the tensor product $`𝒜𝒜`$ space.

The Yang-Baxter algebra is invariant by the transformation $`𝒯_𝒜(\lambda )𝒢_𝒜𝒯_𝒜(\lambda )`$ provided that the group of c-numbers matrices $`𝒢_𝒜`$ satisfies the following property

$$[R(\lambda ,\mu ),𝒢_𝒜𝒢_𝒜]=0$$
(3)
An immediate consequence of this symmetry is the possibility to define the operator

$$T(\lambda )=Tr_𝒜[𝒢_𝒜𝒯_𝒜(\lambda )]$$
(4)
which gives origin to generalized families of commuting transfer matrices.

When the matrix $`𝒢_𝒜`$ is non singular a quantum spin chain can be associated with the transfer matrix (4). For sake of simplicity, consider the usual situation in which the spaces $`𝒜`$ and $`V_i`$ are isomorphic and that the $`_{𝒜i}(\lambda )`$ is proportional to the exchange operator $`P_{𝒜i}`$ at certain special point say $`\lambda =0`$. The corresponding one-dimensional Hamiltonian is obtained as a logarithmic derivative of the transfer matrix at point $`\lambda =0`$, which reads

$$=\underset{i=1}{\overset{L1}{}}P_{ii+1}\frac{d_{ii+1}(\lambda )}{d\lambda }_{\lambda =0}+𝒢_L^1P_{L1}\frac{d_{L1}(\lambda )}{d\lambda }_{\lambda =0}𝒢_L$$
(5)
Clearly, the admissible $`𝒢_𝒜`$ matrices play the role of more general toroidal boundary conditions than the particular periodic case when $`𝒢_𝒜`$ is the identity matrix, the simplest possibility satisfying relation (3). From the point of view of a vertex model, such general twisted boundary conditions correspond to the introduction of a line of defects along the infinite direction on the cylinder. Though boundary conditions are not expected to influence the infinite volume properties it can change the finite-size corrections of massless systems in a strip of width $`L`$ which contains fundamental informations concerning the underlying conformal field theories . For instance, in statistical mechanics boundary conditions provide the means to relate the critical behaviour of a variety of different lattice systems such as the Heisenberg spin chain, the Ashkin-Teller and the Potts models . In this sense, it is highly desirable to study integrable models with as much general boundary conditions as possible.

If the boundary matrix $`𝒢_𝒜`$ is diagonal the corresponding transfer matrix (4) can be diagonalized with very little difference from the periodic case because it does not change in a drastic way the properties of the monodromy matrix elements. The same does not occur when $`𝒢_𝒜`$ is non-diagonal, starting from the fact that the reference state of the periodic case, essential to implement Bethe ansatz approaches, is a priori no longer of utility due to the breaking of the original bulk symmetry by the boundary terms. In fact, progress towards solving commuting transfer matrices with general twists by Bethe ansatz techniques are modest as compared with the literature known for the periodic case, specially for solvable vertex models based on Lie algebras, e.g. refs.. To our knowledge, the six vertex model and its higher spin descendants  are the only solvable vertex systems analyzed so far with non-diagonal boundary conditions. Even in these cases, the functional relation method used in refs. gives the transfer matrix eigenvalues but not information on the corresponding eigenvectors. The latter is certainly an important step in the program of solving integrable systems.

The purpose of this paper is the formulation of the quantum inverse scattering method for the simplest multistate generalization of the six vertex model having $`N`$ independent degrees of freedom on each lattice bond. This turns out to be the isotropic $`SU(N)`$ vertex model whose origin goes back to the work by Uimin  and Sutherland  on generalized integrable Heisenberg chains with higher symmetry. Its corresponding $`_{𝒜i}(\lambda )`$ operators can be written as

$$_{𝒜i}(\lambda )=\lambda I_{𝒜i}+P_{𝒜i}$$
(6)
where as usual $`I_{𝒜i}`$ is the identity matrix on the space $`𝒜V_i`$.

The interesting feature of this system is that the admissible symmetries constitute of arbitrary $`N\times N`$ $`𝒢_𝒜`$ matrices due to the standard property $`P_{12}A_1B_2=B_2A_1P_{12}`$. Therefore this provides us a rich variety of possible diagonal and non-diagonal boundary conditions. In next section, we present the details of the solution of the eigenvalue problem for the transfer matrix in the simplest $`N=2`$ case. Interesting enough, we find that the Bethe ansatz solution depends on the eigenvalue problem related to the boundary $`𝒢_𝒜`$ matrix. In section 3 we generalize these results for arbitrary values of $`N`$ by using the nested Bethe ansatz approach. Our conclusions are presented in section 4 as well as a discussion of singular boundaries for the model $`N=2`$. In Appendix A we discuss briefly the completeness of the Hilbert space for $`N=2`$ and finite $`L`$.

## 2 Algebraic Bethe ansatz for Heisenberg model

The purpose of this section is to determine the eigenvalues and the eigenvectors of the following transfer matrix

$$T(\lambda )=Tr_𝒜[𝒢_𝒜_{𝒜L}(\lambda )\mathrm{}_{𝒜1}(\lambda )]$$
(7)
The operator $`_{𝒜i}(\lambda )`$ is the elementary Boltzmann weights of the isotropic six vertex model which can be written as

$$_{𝒜i}(\lambda )=\left(\begin{array}{cc}\frac{[a(\lambda )+b(\lambda )]}{2}I_i+\frac{[a(\lambda )b(\lambda )]}{2}\sigma _i^z& \sigma _i^{}\\ \sigma _i^+& \frac{[a(\lambda )+b(\lambda )]}{2}I_i\frac{[a(\lambda )b(\lambda )]}{2}\sigma _i^z\end{array}\right)$$
(8)
where $`\sigma _{𝒜,i}^\pm `$ and $`\sigma _{𝒜,i}^z`$ are Pauli matrices acting on the vertical space of states and the weights are $`a(\lambda )=\lambda +1`$ and $`b(\lambda )=\lambda `$. The boundary matrix $`𝒢_𝒜`$ is an arbitrary $`2\times 2`$ matrix over the complex numbers whose matrix elements are denoted by

$$𝒢_𝒜=\left(\begin{array}{cc}g_{11}& g_{12}\\ g_{21}& g_{22}\end{array}\right)$$
(9)
An essential ingredient of the quantum inverse scattering is the existence of a reference state such that the action of the monodromy operator in this state gives as a result a triangular matrix. Though each of the operators $`_{𝒜i}(\lambda )`$ when acting on the trivial spin up $`\left(\begin{array}{c}1\\ 0\end{array}\right)_i`$ or spin down $`\left(\begin{array}{c}0\\ 1\end{array}\right)_i`$ states becomes triangular, such property is not extended to the total monodromy because the off diagonal elements of $`𝒢_𝒜`$ are in general non-null. Therefore the standard ferromagnetic pseudovacuum is not useful when both $`g_{12}`$ and $`g_{21}`$ are different from zero. In order to find an appropriate reference state we have to introduce a set of gauge transformations similar to that used by Baxter  in the solution of the eight vertex model. We replace the local operators $`_{𝒜i}(\lambda )`$ by new matrices $`\stackrel{~}{}_{𝒜j}(\lambda )`$ such that

$$\stackrel{~}{}_{𝒜j}(\lambda )=M_{j+1}^1_{𝒜j}(\lambda )M_j$$
(10)
where $`M_j`$ are arbitrary invertible 2 $`\times `$ 2 c-number matrices acting on the space $`𝒜`$. After performing this gauge transformations the transfer matrix (7) becomes

$$T(\lambda )=Tr_𝒜[M_1^1𝒢_𝒜M_{L+1}\stackrel{~}{𝒯}_𝒜(\lambda )]$$
(11)
where $`\stackrel{~}{𝒯}_𝒜(\lambda )=\stackrel{~}{}_{𝒜L}(\lambda )\mathrm{}\stackrel{~}{}_{𝒜1}(\lambda )`$.

The next step is to look for gauge transformations $`M_j`$ such that $`\stackrel{~}{}_{𝒜j}(\lambda )`$ is annihilated for instance by its lower left element for arbitrary values of the spectral parameter. Representing the matrices $`M_j`$ by

$$M_j=\left(\begin{array}{cc}x_j& r_j\\ y_j& s_j\end{array}\right)$$
(12)
we can conclude  that such annihilation property occurs when the ratio $`\frac{x_j}{y_j}`$ is a constant for $`j=1,\mathrm{},L+1`$. As a consequence of that we can choose the local reference state $`|0_j`$ as

$$|0_j=\left(\begin{array}{c}\frac{x_j}{y_j}\\ 1\end{array}\right)_j$$
(13)
following that the action of the operator $`\stackrel{~}{}_{𝒜j}(\lambda )`$ in this state is given by

$$\stackrel{~}{}_{𝒜j}(\lambda )|0_j=\left(\begin{array}{cc}a(\lambda )\frac{y_j}{y_{j+1}}|0_j& \mathrm{\#}\\ 0& b(\lambda )\frac{y_{j+1}}{y_j}\frac{det[M_j]}{det[M_{j+1}]}|0_j\end{array}\right)$$
(14)
where the symbol $`\mathrm{\#}`$ represents general non-null values.

The remaining freedom that we have on the matrix elements of $`M_j`$ is now used to choose matrices $`M_1`$ and $`M_{L+1}`$ in such way that they transform the boundary matrix $`𝒢_𝒜`$ into a diagonal matrix. More precisely, by imposing that

$$M_1^1𝒢_𝒜M_{L+1}=\left(\begin{array}{cc}\stackrel{~}{g}_1& 0\\ 0& \stackrel{~}{g}_2\end{array}\right)$$
(15)
it follows that the constrains for the first column elements are

$`g_{11}x_{L+1}+g_{12}y_{L+1}`$ $`=`$ $`\stackrel{~}{g}_1x_1`$
$`g_{21}x_{L+1}+g_{22}y_{L+1}`$ $`=`$ $`\stackrel{~}{g}_1y_1`$ (16)
while for the second column elements we have

$`g_{11}r_{L+1}+g_{12}s_{L+1}`$ $`=`$ $`\stackrel{~}{g}_2r_1`$
$`g_{21}r_{L+1}+g_{22}s_{L+1}`$ $`=`$ $`\stackrel{~}{g}_2s_1`$ (17)
At this point we emphasize our assumption that we are dealing with a non-singular boundary matrix. While we have an enormous freedom to choose the second column elements the same does not occur for the first ones because the ratio $`\frac{x_j}{y_j}`$ needs to be kept fixed to preserve triangularity of $`\stackrel{~}{}_{𝒜j}(\lambda )`$. This latter fact together with relation (16) impose a restriction to this ratio which is precisely the same satisfied by ratio of the components of the eigenvectors of the boundary matrix $`𝒢_𝒜`$. Therefore, we have two possibilities for the ratio $`p^{(\pm )}=\frac{x_j}{y_j}`$ which are

$$p^{(\pm )}=\frac{(g_{11}g_{22})\pm \sqrt{(g_{11}g_{22})^2+4g_{12}g_{21}}}{2g_{21}}$$
(18)
Putting now all these informations together it is possible to build up two appropriate global reference states $`|0^{(\pm )}`$ by the tensor product

$$|0^{(\pm )}=\underset{j=1}{\overset{L}{}}\left(\begin{array}{c}p^{(\pm )}\\ 1\end{array}\right)_j$$
(19)
At this point the state (19) preserves at least the desirable triangular property of the total monodromy $`M_1^1𝒢_𝒜M_{L+1}\stackrel{~}{𝒯}_𝒜(\lambda )`$. Below we shall see that they are indeed eigenstates of the transfer matrix (7) independent of further choices of the elements of the gauge matrices $`M_j`$. Further progress is made by recasting the Yang-Baxter algebra for the gauge transformed monodromy $`\stackrel{~}{𝒯}_𝒜(\lambda )`$ in the form of commutation relations for the creation and annihilation fields. In order to do that it is convenient to represent $`\stackrel{~}{𝒯}_𝒜(\lambda )`$ by the following 2 $`\times `$ 2 matrix

$$\stackrel{~}{𝒯}_𝒜(\lambda )=\left(\begin{array}{cc}\stackrel{~}{A}(\lambda )& \stackrel{~}{B}(\lambda )\\ \stackrel{~}{C}(\lambda )& \stackrel{~}{D}(\lambda )\end{array}\right)$$
(20)
As a consequence of the triangular property (14) we are able to derive important relations for the diagonal elements of the transformed monodromy matrix

$`\stackrel{~}{A}(\lambda )|0^{(\pm )}`$ $`=`$ $`[a(\lambda )]^L{\displaystyle \frac{y_1}{y_{L+1}}}|0^{(\pm )}`$
$`\stackrel{~}{D}(\lambda )|0^{(\pm )}`$ $`=`$ $`[b(\lambda )]^L{\displaystyle \frac{y_{L+1}}{y_1}}{\displaystyle \frac{det[M_1]}{det[M_{L+1}]}}|0^{(\pm )}`$ (21)
besides the annihilation property

$$\stackrel{~}{C}(\lambda )|0^{(\pm )}=0$$
(22)
Now taking into account that gauge matrices $`M_j`$ are themselves symmetries allowed by the property (3) it is not difficult to show that the gauge transformed monodromy $`\stackrel{~}{𝒯}_𝒜(\lambda )`$ matrix satisfies the same Yang-Baxter algebra as the original monodromy matrix $`𝒯_𝒜(\lambda )`$. In other words, we have that $`\stackrel{~}{𝒯}_𝒜(\lambda )`$ satisfies the relation

$$R(\lambda ,\mu )\stackrel{~}{𝒯_𝒜}(\lambda )\stackrel{~}{𝒯_𝒜}(\mu )=\stackrel{~}{𝒯_𝒜}(\mu )\stackrel{~}{𝒯_𝒜}(\lambda )R(\lambda ,\mu )$$
(23)
where in our case the $`R`$-matrix is given by

$$R(\lambda ,\mu )=\left(\begin{array}{cccc}a(\lambda \mu )& 0& 0& 0\\ 0& 1& b(\lambda \mu )& 0\\ 0& b(\lambda \mu )& 1& 0\\ 0& 0& 0& a(\lambda \mu )\end{array}\right)$$
(24)
This means that we have the same set of commutation rules of the periodic six vertex model  however now for the gauged matrix elements. Out of sixteen possible relations three of them are of great use, namely

$`\stackrel{~}{A}(\lambda )\stackrel{~}{B}(\mu )`$ $`=`$ $`{\displaystyle \frac{a(\mu \lambda )}{b(\mu \lambda )}}\stackrel{~}{B}(\mu )\stackrel{~}{A}(\lambda ){\displaystyle \frac{1}{b(\mu \lambda )}}\stackrel{~}{B}(\lambda )\stackrel{~}{A}(\mu )`$ (25)
$`\stackrel{~}{D}(\lambda )\stackrel{~}{B}(\mu )`$ $`=`$ $`{\displaystyle \frac{a(\lambda \mu )}{b(\lambda \mu )}}\stackrel{~}{B}(\mu )\stackrel{~}{D}(\lambda ){\displaystyle \frac{1}{b(\lambda \mu )}}\stackrel{~}{B}(\lambda )\stackrel{~}{D}(\mu )`$ (26)
$`[\stackrel{~}{B}(\lambda ),\stackrel{~}{B}(\mu )]`$ $`=`$ $`0`$ (27)
The fields $`\stackrel{~}{B}(\lambda )`$ are then interpreted as a kind of creation operators over the pseudovacuum $`|0^{(\pm )}`$ and a natural ansatz for the eigenvectors $`|\varphi ^{(\pm )}`$ of the transfer matrix $`T(\lambda )`$ is

$$|\varphi ^{(\pm )}=\underset{j=1}{\overset{n_\pm }{}}\stackrel{~}{B}(\lambda _j^{(\pm )})|0^{(\pm )}$$
(28)
The eigenvalue problem $`T(\lambda )|\varphi ^{(\pm )}=\mathrm{\Lambda }^{(\pm )}(\lambda )|\varphi ^{(\pm )}`$ now becomes

$$\left[\stackrel{~}{g}_1\stackrel{~}{A}(\lambda )+\stackrel{~}{g}_2\stackrel{~}{D}(\lambda )\right]|\varphi ^{(\pm )}=\mathrm{\Lambda }^{(\pm )}(\lambda )|\varphi ^{(\pm )}$$
(29)
and it can be solved in the same way as the periodic six vertex model , i.e by taking the fields $`\stackrel{~}{A}(\lambda )`$ and $`\stackrel{~}{D}(\lambda )`$ through the creation operators $`\stackrel{~}{B}(\lambda )`$ with the help of the commutation rules (25-27). One peculiarity here, however, is that the calculations involving the action of the diagonal fields $`\stackrel{~}{A}(\lambda )`$ and $`\stackrel{~}{D}(\lambda )`$ over the reference state $`|0^{(\pm )}`$ requires extra simplifications to eliminate unnecessary dependence of the gauge matrices elements. They are carried out by using the help of Eqs.(16 \- 18) and our final result for the eigenvalues $`\mathrm{\Lambda }^{(\pm )}(\lambda )`$ are

$$\mathrm{\Lambda }^{(\pm )}(\lambda )=g^{(\pm )}[a(\lambda )]^L\underset{i=1}{\overset{n_\pm }{}}\frac{\lambda _i^{(\pm )}\lambda +\frac{1}{2}}{\lambda _i^{(\pm )}\lambda \frac{1}{2}}+g^{()}[b(\lambda )]^L\underset{i=1}{\overset{n_\pm }{}}\frac{\lambda \lambda _i^{(\pm )}+\frac{3}{2}}{\lambda \lambda _i^{(\pm )}+\frac{1}{2}}$$
(30)
provided that the rapidities $`\lambda _i^{(\pm )}`$ satisfy the following Bethe ansatz equations

$$\left[\frac{\lambda _i^{(\pm )}+\frac{1}{2}}{\lambda _i^{(\pm )}\frac{1}{2}}\right]^L=\frac{g^{()}}{g^{(\pm )}}\underset{\stackrel{j=1}{ji}}{\overset{n_\pm }{}}\frac{\lambda _i^{(\pm )}\lambda _j^{(\pm )}+1}{\lambda _i^{(\pm )}\lambda _j^{(\pm )}1}$$
(31)
where we have performed the convenient shift $`\lambda _i^{(\pm )}\lambda _i^{(\pm )}\frac{1}{2}`$. The phase factors $`g^{(\pm )}`$ are just the eigenvalues of the matrix $`𝒢_𝒜`$

$$g^{(\pm )}=\frac{(g_{11}+g_{22})\pm \sqrt{(g_{11}g_{22})^2+4g_{12}g_{21}}}{2}$$
(32)
Rather remarkably, we see that the final form of the eigenvalues and Bethe ansatz equations resemble much that of the isotropic six vertex model with diagonal boundary if we replace the diagonal twists by the eigenvalues of the non-diagonal boundary $`𝒢_𝒜`$ matrix. The eigenvectors can also be thought as multi-particle states in which the integers $`n_\pm L`$ play the role of particle number sectors. We emphasize, however, that the corresponding basic creation fields are much more sophisticated operators than that of the periodic six vertex model . It is tempting to think that the two possible ways we have at our disposal to build up the Hilbert space is related to the remaining $`Z_2`$ symmetry allowed by boundary terms. One expects therefore that it should be possible to obtain the eigenvalues of the transfer matrix either from the $`|0^{(+)}`$ or $`|0^{()}`$ pseudovacuums. Indeed, we have verified this fact by numerically solving the equations for some values of $`L`$ and comparing them to exact diagonalization of the transfer matrix (7). We note, however, that a given eigenvalue of the transfer matrix is in general obtained at different particle sectors $`n_\pm `$ over the $`|0^{(\pm )}`$ reference states. For example, the eigenvalue $`g^+[a(\lambda )]^L+g^{}[b(\lambda )]^L`$ can be obtained either from the zero-particle state $`|0^{(+)}`$ or as a $`L`$-particle state over the pseudovacuum $`|0^{()}`$. In Appendix A, we present details of our study for $`L=2`$ in which Eqs.(30-31) can be solved by analytical means. Our numerical results up to $`L=4`$ suggest that two possible branches of the Bethe ansatz solutions (31) produce the complete spectrum of the transfer matrix (7). It would be interesting to further investigate the completeness of the Bethe ansatz (31) by adapting the recent arguments developed by Baxter  to the case of non-diagonal twists.

We now can derive similar results for the spin chain that commutes with the transfer matrix (7). The corresponding spin-$`1/2`$ XXX Hamiltonian follows from expression (5) and it is given by

$$=𝒥\underset{j=1}{\overset{L}{}}\left(\sigma _j^+\sigma _{j+1}^{}+\sigma _j^{}\sigma _{j+1}^++\frac{\sigma _j^z\sigma _{j+1}^z}{2}\right)$$
(33)
with the following boundary condition

$$\left(\begin{array}{c}\sigma _{L+1}^+\\ \sigma _{L+1}^{}\\ \sigma _{L+1}^z\end{array}\right)=\frac{1}{g_{11}g_{22}g_{12}g_{21}}\left(\begin{array}{ccc}g_{11}^2& g_{21}^2& g_{11}g_{21}\\ g_{12}^2& g_{22}^2& g_{12}g_{22}\\ 2g_{11}g_{12}& 2g_{21}g_{22}& g_{11}g_{22}+g_{12}g_{21}\end{array}\right)\left(\begin{array}{c}\sigma _1^+\\ \sigma _1^{}\\ \sigma _1^z\end{array}\right)$$
(34)
Its eigenvalues $`E^{(\pm )}=\frac{dLog[\mathrm{\Lambda }^{(\pm )}(\lambda )]}{d\lambda }_{\lambda =0}`$ are

$$E^{(\pm )}=𝒥\underset{j=1}{\overset{n_\pm }{}}\frac{1}{\lambda _j^{\pm 2}\frac{1}{4}}+\frac{𝒥L}{2}$$
(35)
where $`\lambda _i^{(\pm )}`$ satisfy the Bethe ansatz equations (31).

Our final comment concerns with the comparison between our results (30-31) and that of refs. in the isotropic limit case when the trigonometric weights become rational functions. We see that they are in accordance for the common non-diagonal boundary $`𝒢_𝒜=\left(\begin{array}{cc}0& g_{12}\\ g_{21}& 0\end{array}\right)`$ apart from the fact that our numbers of roots $`n_\pm `$ can vary up to $`L`$ while that of refs. are fixed at $`L`$. This implies that for non-diagonal boundary conditions the complete solution of the isotropic limit does not follows directly from that found for the anisotropic six vertex model . This means that even for this particular non-diagonal boundary the results (30-31) are novel in the literature.

## 3 Nested Bethe ansatz for $`SU(N)`$ model

The purpose of this section is to generalize the results of the previous section for general $`N`$. We wish to diagonalize the transfer matrix (7), where now the operator $`_{𝒜i}(\lambda )`$ is

$$_{𝒜i}(\lambda )=a(\lambda )\underset{\alpha =1}{\overset{N}{}}\widehat{e}_{\alpha \alpha }^{(𝒜)}\widehat{e}_{\alpha \alpha }^{(i)}+b(\lambda )\underset{\stackrel{\alpha ,\beta =1}{\alpha \beta }}{\overset{N}{}}\widehat{e}_{\alpha \alpha }^{(𝒜)}\widehat{e}_{\beta \beta }^{(i)}+\underset{\stackrel{\alpha ,\beta =1}{\alpha \beta }}{\overset{N}{}}\widehat{e}_{\alpha \beta }^{(𝒜)}\widehat{e}_{\beta \alpha }^{(i)}$$
(36)
where $`\widehat{e}_{ij}^{(V)}`$ are the standard Weyl matrices whose elements acting on the space $`V`$ are $`[\widehat{e}_{ij}^{(V)}]_{kl}=\delta _{ik}\delta _{jl}`$. In this basis the boundary matrix $`𝒢_𝒜`$ is generally represented by

$$𝒢_𝒜=\underset{\alpha ,\beta =1}{\overset{N}{}}g_{\alpha \beta }\widehat{e}_{\alpha \beta }^{(𝒜)}$$
(37)
As before we have to seek for suitable references states by imposing the gauge transformation (10) for each operator (36) and require that they are up triangular when acting on such pseudovacuum. Denoting the gauge matrices by $`M_j={\displaystyle \underset{\alpha ,\beta =1}{\overset{N}{}}}m_j(\alpha ,\beta )\widehat{e}_{\alpha \beta }^{(𝒜)}`$ we find that such triangular property is fully achieved when the following ratios relations are satisfied

$$p_{\alpha ,\beta }=\frac{m_j(\alpha ,\beta )}{m_j(N,\beta )}=\frac{m_{j+1}(\alpha ,\beta )}{m_{j+1}(N,\beta )},\alpha ,\beta =1,\mathrm{},N1$$
(38)
for each $`j=1,\mathrm{},L+1`$. In terms of these ratios the local reference state $`|0_j`$ assume the form

$$|0_j=\left(\begin{array}{c}p_{1,1}\\ p_{2,1}\\ \mathrm{}\\ p_{N1,1}\\ 1\end{array}\right)_j$$
(39)
Other important ingredient is the action of the gauge transformed operator $`\stackrel{~}{}_{𝒜j}(\lambda )`$ over the local state of reference. This now can be represented by the following $`N\times N`$ matrix on the space $`𝒜`$

$$\stackrel{~}{}_{𝒜j}|0_j=\left(\begin{array}{ccccc}a(\lambda )\frac{f_1^j}{f_1^{j+1}}|0_j& \mathrm{\#}& \mathrm{\#}& \mathrm{}& \mathrm{\#}\\ 0& b(\lambda )\frac{f_2^j}{f_2^{j+1}}|0_j& \mathrm{\#}& \mathrm{}& \mathrm{\#}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& b(\lambda )\frac{f_N^j}{f_N^{j+1}}|0_j\end{array}\right)_{N\times N}$$
(40)
where the variables $`f_\alpha ^j`$ are given by

$$f_\alpha ^j=\{\begin{array}{cc}m_j(N,\alpha ),\alpha =1,\mathrm{},N1\hfill & \\ \left(\underset{i=1}{\overset{N1}{}}\frac{1}{m_j(N,i)}\right)det[M_j],\alpha =N\hfill & \end{array}$$
(41)
Similarly to the previous section we can take advantage of the remaining freedom of the elements of the gauge matrices to transform $`M_1^1𝒢_𝒜M_{L+1}`$ into a diagonal matrix. By imposing this condition the matrices elements of $`M_1`$ and $`M_{L+1}`$ become related by the expression

$$\underset{\gamma =1}{\overset{N}{}}g_{\alpha \gamma }m_{L+1}(\gamma ,\beta )=m_1(\alpha ,\beta )\stackrel{~}{g}_\beta ,\alpha ,\beta =1,\mathrm{},N$$
(42)
where $`\stackrel{~}{g}_\alpha `$ represent the diagonal elements of the transformed boundary matrix.

Equations (38) and (42) together impose constrains to the possible values ratios $`p_{\alpha ,\beta }`$ which turns out to be same conditions satisfied by the ratio of the components of the eigenvectors of boundary matrix $`𝒢_𝒜`$. This means that we have $`N`$ possible choices for $`p_{\alpha ,1}^{(l)}`$ $`l=1,\mathrm{},N`$ and consequently from Eq.(39) $`N`$ kind of suitable local references states $`|0_j^{(l)}`$. A natural ansatz for the $`N`$ possible choices of global reference states are

$$|0^{(l)}=\underset{j=1}{\overset{L}{}}|0_j^{(l)}l=1,\mathrm{},N$$
(43)
The next step is to write a suitable representation for the gauge transformed monodromy matrix in the auxiliary space $`𝒜`$. The triangular property (40) suggests us to seek for the structure used in nested Bethe ansatz diagonalization of the periodic $`SU(N)`$ vertex models  which is

$$\stackrel{~}{𝒯}(\lambda )=\left(\begin{array}{cccc}\stackrel{~}{A}(\lambda )& \stackrel{~}{B}_1(\lambda )& \mathrm{}& \stackrel{~}{B}_{N1}(\lambda )\\ \stackrel{~}{C}_1(\lambda )& \stackrel{~}{D}_{11}(\lambda )& \mathrm{}& \stackrel{~}{D}_{1N1}(\lambda )\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \stackrel{~}{C}_{N1}(\lambda )& \stackrel{~}{D}_{N11}(\lambda )& \mathrm{}& \stackrel{~}{D}_{N1N1}(\lambda )\end{array}\right)_{N\times N}$$
(44)
The triangularity property (40) implies that the fields $`\stackrel{~}{B}_i(\lambda )`$ play the role of creations operators, $`\stackrel{~}{C}_i(\lambda )`$ are annihilation fields while the diagonal operator $`\stackrel{~}{A}(\lambda )`$ and $`\stackrel{~}{D}_{ii}(\lambda )`$ acts on the reference state $`|0^{(l)}`$ as

$`\stackrel{~}{A}(\lambda )|0^{(l)}`$ $`=`$ $`[a(\lambda )]^L{\displaystyle \frac{f_1^1}{f_1^{L+1}}}|0^{(l)}`$ (45)
$`\stackrel{~}{D}_{ii}(\lambda )|0^{(l)}`$ $`=`$ $`[b(\lambda )]^L{\displaystyle \frac{f_{i+1}^1}{f_{i+1}^{L+1}}}|0^{(l)},i=1,\mathrm{},N1`$ (46)
To construct other eigenvectors we shall use the commutation relations between the gauge transformed monodromy matrix elements. The arguments used in section 2 allows to conclude that these commutation rules are the same as that already known for the periodic $`SU(N)`$ models . The most useful relations for subsequent derivations are

$`\stackrel{~}{A}(\lambda )\stackrel{~}{B}_i(\mu )={\displaystyle \frac{a(\mu \lambda )}{b(\mu \lambda )}}\stackrel{~}{B}_i(\mu )\stackrel{~}{A}(\lambda ){\displaystyle \frac{1}{b(\mu \lambda )}}\stackrel{~}{B}_i(\lambda )\stackrel{~}{A}(\mu )`$ (47)
$`\stackrel{~}{D}_{ij}(\lambda )\stackrel{~}{B}_k(\mu )={\displaystyle \frac{1}{b(\lambda \mu )}}\stackrel{~}{B}_p(\mu )\stackrel{~}{D}_{iq}(\lambda )r^{(1)}(\lambda \mu )_{pq}^{jk}{\displaystyle \frac{1}{b(\lambda \mu )}}\stackrel{~}{B}_j(\lambda )\stackrel{~}{D}_{ik}(\mu )`$ (48)
$$\stackrel{~}{B}_i(\lambda )\stackrel{~}{B}_j(\mu )=\stackrel{~}{B}_p(\mu )\stackrel{~}{B}_q(\lambda )r^{(1)}(\lambda \mu )_{pq}^{ij}$$
(49)
where $`r^{(1)}(\lambda )_{pq}^{ij}`$ are the elements of the $`R`$-matrix associated to the $`SU(N1)`$ vertex model.

In terms of the gauge transformed fields, the eigenvalue problem for the transfer matrix $`T(\lambda )`$ becomes

$$\left[\stackrel{~}{g}_1\stackrel{~}{A}(\lambda )+\underset{i=1}{\overset{N1}{}}\stackrel{~}{g}_{i+1}\stackrel{~}{D}_{ii}(\lambda )\right]|\varphi ^{(l)}=\mathrm{\Lambda }^{(l)}(\lambda )|\varphi ^{(l)}$$
(50)
where $`|\varphi ^{(l)}`$ denotes the eigenvectors. Previous experience with these models  suggests us to suppose that eigenvectors can be written in terms of the following linear combination

$$|\varphi ^{(l)}=\stackrel{~}{B}_{a_1}(\lambda _1^{(1,l)})\mathrm{}\stackrel{~}{B}_{a_{m_1^l}}(\lambda _{m_1^l}^{(1,l)})^{a_{m_1^l}\mathrm{}a_1}|0^{(l)}$$
(51)
where sum over repeated indices $`a_n=1,\mathrm{}N1`$ is assumed. At this stage the components $`^{a_{m_1^l}\mathrm{}a_1}`$ are thought as coefficients of an arbitrary linear combination that are going to be determined a posteriori.

By carrying on the fields $`\stackrel{~}{A}(\lambda )`$ and $`\stackrel{~}{D}_{ii}(\lambda )`$ over the multi-particle state (51) we generate terms that are proportional to $`|\varphi ^{(l)}`$ and those that are not the so-called unwanted terms. The first ones will contribute directly to the eigenvalue $`\mathrm{\Lambda }^{(l)}(\lambda )`$ and are obtained by keeping only the first term of the commutation rules (47-48). These calculations are by now standard in the literature and here we present only the main results of the action of the transfer matrix on the eigenvector $`|\varphi ^{(l)}`$ which is

$`T(\lambda )|\varphi ^{(l)}`$ $`=`$ $`\stackrel{~}{g}_1{\displaystyle \frac{f_1^1}{f_1^{L+1}}}[a(\lambda )]^L{\displaystyle \underset{j=1}{\overset{m_1^l}{}}}{\displaystyle \frac{a(\lambda _j^{(1,l)}\lambda )}{b(\lambda _j^{(1,l)}\lambda )}}|\varphi ^{(l)}`$ (52)
$`+[b(\lambda )]^L{\displaystyle \underset{j=1}{\overset{m_1^l}{}}}{\displaystyle \frac{1}{b(\lambda \lambda _j^{(1,l)})}}\stackrel{~}{B}_{b_1}(\lambda _1^{(1,l)})\mathrm{}\stackrel{~}{B}_{b_{m_1^l}}(\lambda _{m_1^l}^{(1,l)})T^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})_{b_1\mathrm{}b_{m_1^l}}^{a_1\mathrm{}a_{m_1^l}}^{a_{m_1^l}\mathrm{}a_1}|0^{(l)}`$
$`+\mathrm{unwanted}\mathrm{terms}`$
All the pieces entering the above expression can be summarized as follows. The terms $`T^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})_{b_1\mathrm{}b_{m_1^l}}^{a_1\mathrm{}a_{m_1^l}}`$ are transfer matrix elements of an auxiliary inhomogeneous problem related to the $`SU(N1)`$ vertex model with twisted boundaries $`\stackrel{~}{𝒢}`$ defined by

$$T^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})_{b_1\mathrm{}b_{m_1^l}}^{a_1\mathrm{}a_{m_1^l}}=r^{(1)}(\lambda \lambda _1^{(1,l)})_{b_1d_1}^{aa_1}r^{(1)}(\lambda \lambda _2^{(1,l)})_{b_2d_2}^{d_1a_2}\mathrm{}r^{(1)}(\lambda \lambda _{m_1^l}^{(1,l)})_{b_{m_1^l}d_{m_1^l}}^{d_{m_1^l1}a_{m_1^l}}\stackrel{~}{𝒢}_{ad_{m_1^l}}$$
(53)
where $`\stackrel{~}{𝒢}_{ab}`$ denotes the elements of the boundary matrix $`\stackrel{~}{𝒢}`$ given by

$$\stackrel{~}{𝒢}=\left(\begin{array}{ccccc}\stackrel{~}{g}_2\frac{f_2^1}{f_2^{L+1}}& \mathrm{\#}& \mathrm{\#}& \mathrm{}& \mathrm{\#}\\ 0& \stackrel{~}{g}_3\frac{f_3^1}{f_3^{L+1}}& \mathrm{\#}& \mathrm{}& \mathrm{\#}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& \stackrel{~}{g}_N\frac{f_N^1}{f_N^{L+1}}\end{array}\right)_{N1\times N1}$$
(54)
The unwanted terms are originated when the variables $`\lambda _i^{(1,l)}`$ of the multi-particle state (51) are exchanged with the spectral parameter $`\lambda `$ due to the second part of the commutation rules (47-48). It is possible to collect these terms in closed forms, thanks to the commutation rule (49) which makes possible to relate different ordered multi-particle states. It turns out that all the unwanted terms are canceled out provided that the rapidities $`\lambda _i^{(1,l)}`$ satisfy the following restriction,

$`\stackrel{~}{g}_1{\displaystyle \frac{f_1^1}{f_1^{L+1}}}\left[{\displaystyle \frac{a(\lambda _i^{(1,l)})}{b(\lambda _i^{(1,l)})}}\right]^L{\displaystyle \underset{\stackrel{j=1}{ji}}{\overset{m_1^l}{}}}b(\lambda _i^{(1,l)}\lambda _j^{(1,l)}){\displaystyle \frac{a(\lambda _j^{(1,l)}\lambda _i^{(1,l)})}{b(\lambda _j^{(1,l)}\lambda _i^{(1,l)})}}^{a_{m_1^l}\mathrm{}a_1}=`$
$`T^{(1)}(\lambda =\lambda _i^{(1,l)},\{\lambda _j^{(1,l)}\})_{a_1\mathrm{}a_{m_1^l}}^{b_1\mathrm{}b_{m_1^l}}^{b_{m_1^l}\mathrm{}b_1},i=1,\mathrm{},m_1^l`$ (55)
Now we reached a point which is fundamental to diagonalize $`T^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})`$ in order to compute the eigenvalues of $`T(\lambda )`$ and at the same time to solve Eq.(55). This becomes possible if we require that $`^{a_{m_1^l}\mathrm{}a_1}`$ is an eigenvector of the auxiliary transfer matrix with eigenvalue $`\mathrm{\Lambda }^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})`$, namely

$$T^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})_{a_1\mathrm{}a_{m_1^l}}^{b_1\mathrm{}b_{m_1^l}}^{b_{m_1^l}\mathrm{}b_1}=\mathrm{\Lambda }^{(1)}(\lambda ,\{\lambda _i^{(1)}\})^{a_{m_1^l}\mathrm{}a_1}$$
(56)
Inspection of Eq.(52) and Eq.(55) together with Eq.(56) shows that the eigenvalue of $`T(\lambda )`$ is

$$\mathrm{\Lambda }(\lambda )^{(l)}=\stackrel{~}{g}_1\frac{f_1^1}{f_1^{L+1}}[a(\lambda )]^L\underset{i=1}{\overset{m_1^l}{}}\frac{a(\lambda _i^{(1,l)}\lambda )}{b(\lambda _i^{(1,l)}\lambda )}+[b(\lambda )]^L\underset{i=1}{\overset{m_1^l}{}}\frac{1}{b(\lambda \lambda _i^{(1,l)})}\mathrm{\Lambda }^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})$$
(57)
and the nested Bethe ansatz equations (55) become

$$\stackrel{~}{g}_1\frac{f_1^1}{f_1^{L+1}}\left[\frac{a(\lambda _i^{(1,l)})}{b(\lambda _i^{(1,l)})}\right]^L\underset{\stackrel{j=1}{ji}}{\overset{m_1^l}{}}b(\lambda _i^{(1,l)}\lambda _j^{(1,l)})\frac{a(\lambda _j^{(1,l)}\lambda _i^{(1,l)})}{b(\lambda _j^{(1,l)}\lambda _i^{(1,l)})}=\mathrm{\Lambda }^{(1)}(\lambda =\lambda _i^{(1,l)},\{\lambda _j^{(1,l)}\}),i=1,\mathrm{},m_1^l$$

In order to solve the eigenvalue problem (56) it is necessary to introduce a second algebraic Bethe ansatz for the eigenvectors $`^{a_{m_1^l}\mathrm{}a_1}`$. Because the boundary matrix $`\stackrel{~}{𝒢}`$ is triangular there is no need to perform gauge transformations to find an appropriate reference state for $`T^{(1)}(\lambda ,\{\lambda _i^{(1,l)}\})`$. We can use, for instance, the usual ferromagnetic pseudovacuum build up by tensor product of elementary $`(N1)`$-dimensional $`\left(\begin{array}{c}1\\ 0\\ \mathrm{}\\ 0\end{array}\right)_{N1}`$ vectors. As a result the solution (56) becomes very similar to that of the periodic $`SU(N1)`$ vertex model in the presence of inhomogeneities. Since this problem has been extensively discussed in the literature we will only present our final results for the main eigenvalue problem (50). It turns out that the eigenvalues of the transfer matrix $`\mathrm{\Lambda }^{(l)}(\lambda )`$ is given by

$`\mathrm{\Lambda }^{(l)}(\lambda ;\{\lambda _i^{(1,l)}\},\mathrm{},\{\lambda _i^{(N1,l)}\})`$ $`=`$ $`\stackrel{~}{g}_1{\displaystyle \frac{f_1^1}{f_1^{L+1}}}[a(\lambda )]^L{\displaystyle \underset{j=1}{\overset{m_1^l}{}}}{\displaystyle \frac{a(\lambda _j^{(1,l)}\lambda )}{b(\lambda _j^{(1,l)}\lambda )}}`$ (58)
$`+[b(\lambda )]^L{\displaystyle \underset{k=1}{\overset{N2}{}}}\stackrel{~}{g}_{k+1}{\displaystyle \frac{f_{k+1}^1}{f_{k+1}^{L+1}}}{\displaystyle \underset{j=1}{\overset{m_k^l}{}}}{\displaystyle \frac{a(\lambda \lambda _j^{(k,l)})}{b(\lambda \lambda _j^{(k,l)})}}{\displaystyle \underset{j=1}{\overset{m_{k+1}^l}{}}}{\displaystyle \frac{a(\lambda _j^{(k+1,l)}\lambda )}{b(\lambda _j^{(k+1,l)}\lambda )}}`$
$`+[b(\lambda )]^L\stackrel{~}{g}_N{\displaystyle \frac{f_N^1}{f_N^{L+1}}}{\displaystyle \underset{j=1}{\overset{m_{N1}^l}{}}}{\displaystyle \frac{a(\lambda \lambda _j^{(N1,l)})}{b(\lambda \lambda _j^{(N1,l)})}}`$
The rapidities $`\{\lambda _i^{(k,l)}\}`$ $`k=1,\mathrm{},N`$ parameterize the multi-particle states of the nesting problem at step $`k`$ and are required to satisfy the following nested Bethe ansatz equations

$$\frac{\stackrel{~}{g}_1}{\stackrel{~}{g}_2}\frac{f_1^1f_2^{L+1}}{f_1^{L+1}f_2^1}\left[\frac{a(\lambda _i^{(1,l)})}{b(\lambda _i^{(1,l)})}\right]^L=\underset{\stackrel{j=1}{ji}}{\overset{m_1^l}{}}\frac{a(\lambda _i^{(1,l)}\lambda _j^{(1,l)})}{a(\lambda _j^{(1,l)}\lambda _i^{(1,l)})}\underset{\stackrel{j=1}{ji}}{\overset{m_2^l}{}}\frac{a(\lambda _j^{(2,l)}\lambda _i^{(1,l)})}{b(\lambda _j^{(2,l)}\lambda _i^{(1,l)})}$$
(59)
$`{\displaystyle \frac{\stackrel{~}{g}_k}{\stackrel{~}{g}_{k+1}}}{\displaystyle \frac{f_k^1f_{k+1}^{L+1}}{f_k^{L+1}f_{k+1}^1}}{\displaystyle \underset{j=1}{\overset{m_{k1}^l}{}}}{\displaystyle \frac{a(\lambda _i^{(k,l)}\lambda _j^{(k1,l)})}{b(\lambda _i^{(k,l)}\lambda _j^{(k1,l)})}}={\displaystyle \underset{\stackrel{j=1}{ji}}{\overset{m_k^l}{}}}{\displaystyle \frac{a(\lambda _i^{(k,l)}\lambda _j^{(k,l)})}{a(\lambda _j^{(k,l)}\lambda _i^{(k,l)})}}{\displaystyle \underset{\stackrel{j=1}{ji}}{\overset{m_{k+1}^l}{}}}{\displaystyle \frac{a(\lambda _j^{(k+1,l)}\lambda _i^{(k,l)})}{b(\lambda _j^{(k+1,l)}\lambda _i^{(k,l)})}},`$ (60)
$`k=2,\mathrm{},N2`$
$$\frac{\stackrel{~}{g}_{N1}}{\stackrel{~}{g}_N}\frac{f_{N1}^1f_N^{L+1}}{f_{N1}^{L+1}f_N^1}\underset{\stackrel{j=1}{ji}}{\overset{m_{N2}^l}{}}\frac{a(\lambda _i^{(N1,l)}\lambda _j^{(N2,l)})}{b(\lambda _i^{(N1,l)}\lambda _j^{(N2,l)})}=\underset{\stackrel{j=1}{ji}}{\overset{m_{N1}^l}{}}\frac{a(\lambda _i^{(N1,l)}\lambda _j^{(N1,l)})}{a(\lambda _j^{(N1,l)}\lambda _i^{(N1,l)})}$$
(61)
The final step is to carry out simplifications on the phase factors $`\stackrel{~}{g}_i\frac{f_i^1}{f_i^{L+1}}`$ with the help of the constraints (38) and (42). After a cumbersome algebra it is possible to show that such factors are just the eigenvalues $`g^{(i)}`$ of the boundary matrix $`𝒢_𝒜`$. To make sure that the different possibilities we have at hand for the ratios $`p_{\alpha ,\beta }^{(l)}`$ do not lead us to singular gauge matrices we choose to order them for each $`l`$-th choice of pseudovacuum by $`g^{(l)}=g^{(l+N)}`$ for $`l=1,\mathrm{},N`$. Taking into account this ordering as well as performing the shifts $`\{\lambda _j^{(k,l)}\}\{\lambda _j^{(k,l)}\}k/2`$ our result (58) for the eigenvalue becomes

$`\mathrm{\Lambda }^{(l)}(\lambda ;\{\lambda _i^{(1,l)}\},\mathrm{},\{\lambda _i^{(N1,l)}\})`$ $`=`$ $`g^{(l)}[a(\lambda )]^L{\displaystyle \underset{j=1}{\overset{m_1^l}{}}}{\displaystyle \frac{\lambda _j^{(1,l)}\lambda +\frac{1}{2}}{\lambda _j^{(1,l)}\lambda \frac{1}{2}}}`$ (62)
$`+[b(\lambda )]^L{\displaystyle \underset{k=1}{\overset{N2}{}}}g^{(l+k)}{\displaystyle \underset{j=1}{\overset{m_k^l}{}}}{\displaystyle \frac{\lambda \lambda _j^{(k,l)}+\frac{k+2}{2}}{\lambda \lambda _j^{(k,l)}+\frac{k}{2}}}{\displaystyle \underset{j=1}{\overset{m_{k+1}^l}{}}}{\displaystyle \frac{\lambda _j^{(k+1,l)}\lambda +\frac{1k}{2}}{\lambda _j^{(k+1,l)}\lambda \frac{k+1}{2}}}`$
$`+[b(\lambda )]^Lg^{(l+N1)}{\displaystyle \underset{j=1}{\overset{m_{N1}^l}{}}}{\displaystyle \frac{\lambda \lambda _j^{(N1,l)}+\frac{N+1}{2}}{\lambda \lambda _j^{(N1,l)}+\frac{N1}{2}}}`$
and the nested Bethe ansatz equations can be compactly written as

$$\left[\frac{\lambda _i^{(a,l)}+\frac{\delta _{a,1}}{2}}{\lambda _i^{(a,l)}\frac{\delta _{a,1}}{2}}\right]^L=\frac{g^{(l+a)}}{g^{(l+a1)}}\underset{b=1}{\overset{N1}{}}\underset{\stackrel{k=1}{ki}}{\overset{m_b^l}{}}\frac{\lambda _i^{(a,l)}\lambda _k^{(b,l)}+\frac{C_{a,b}}{2}}{\lambda _i^{(a,l)}\lambda _k^{(b,l)}\frac{C_{a,b}}{2}},i=1,\mathrm{},m_a^l;a=1,\mathrm{},N1$$
(63)
where $`C_{ab}`$ is the Cartan matrix elements of the $`SU(N)`$ Lie algebra.

We see that the results for the eigenvalues and the Bethe ansatz equations is similar to that expected from the $`SU(N)`$ vertex model with diagonal twists giving by the eigenvalues of the boundary matrix $`𝒢_𝒜`$. It remains to be investigated whether this interesting feature is particular of the $`SU(N)`$ symmetry or it also works in other isotropic vertex models such as those invariant by the $`O(N)`$ and $`Sp(2N)`$ Lie algebras.

## 4 Concluding remarks

In this paper we have been able to apply the quantum inverse scattering program to solve exactly the isotropic $`SU(N)`$ vertex model with non-diagonal twisted boundary conditions. We find that the eigenvectors can be constructed in terms of multi-particle states over $`N`$ possible pseudovacuums. The Bethe ansatz results for the eigenvalues are similar to that of the $`SU(N)`$ model with diagonal boundaries in which the eigenvalues of the boundary matrix $`𝒢_𝒜`$ play the role of the diagonal twists.

We expect that our results can be generalized without further difficulties to accommodate the solution of vertex models based on the $`SL(N|M)`$ super Lie algebra  with general non-diagonal twists. These will include interesting systems of correlated electrons on a lattice such as the one-dimensional supersymmetric t-J model  and the so-called Essler, Korepin and Schoutens superconducting model  with arbitrary symmetry breaking boundary conditions. With more effort we hope that our approach can be further generalized to include the trigonometric deformation of those vertex models based on the $`U_q[SL(N|M)]`$ symmetry. In these cases, however, we recall that the possible $`𝒢_𝒜`$ matrices compatible with integrability belong to a smaller group formed by one-dimensional dilatations and the discrete $`Z_{N+M}`$ symmetry.

Other interesting issue that deserves investigation is the situation when the boundary matrix $`𝒢_𝒜`$ becomes singular. For example, one would like to ask it is still possible to exhibit eigenvectors of the transfer matrix that are given by direct tensor product of $`N`$-dimensional vectors such as the reference states of sections 2 and 3. We have studied this problem in the simplest case $`N=2`$ and surprisingly we found a family of such states $`|\varphi ^{(n)}`$ which are

$$|\varphi ^{(n)}=\underset{i=1}{\overset{n}{}}\left(\begin{array}{c}\frac{g_{22}}{g_{21}}\\ 1\end{array}\right)_i\underset{i=n+1}{\overset{L}{}}\left(\begin{array}{c}\frac{g_{11}}{g_{21}}\\ 1\end{array}\right)_i,n=0,1,\mathrm{},L$$
(64)
whose corresponding eigenvalues $`\mathrm{\Lambda }^{(n)}(\lambda )`$ have also the following simple factorized form

$$\mathrm{\Lambda }^{(n)}(\lambda )=(g_{11}+g_{22})[a(\lambda )]^{Ln}[b(\lambda )]^n$$
(65)
This result prompted us to study further properties of the transfer matrix (7) when $`𝒢_𝒜`$ is a singular matrix. Our study for finite $`L`$ up to six sites reveals that the roots of the characteristic polynomial of $`T(\lambda )`$ are exactly the eigenvalues (65) whose degeneracy is the binomial coefficient $`d_n=\frac{L!}{(Ln)!n!}`$. In the case of singular boundary matrix $`T(\lambda )`$ becomes defective since it has fewer than $`2^L`$ independent eigenvectors. To each distinct eigenvalue $`\mathrm{\Lambda }^{(n)}(\lambda )`$ we find only one eigenvector which is precisely the state (64) and therefore the total number of independent states is $`L+1`$. These results are strong evidences that $`T(\lambda )`$ behaves as a non-derogatory matrix and we conjecture that its Jordan decomposition for arbitrary $`L`$ should be

$$T(\lambda )=diag(J_0,J_1,\mathrm{},J_L)$$
(66)
where $`J_n`$ is a $`d_n\times d_n`$ Jordan matrix is given by

$$J_n=\left(\begin{array}{ccccc}\mathrm{\Lambda }^{(n)}(\lambda )& 1& 0& \mathrm{}& 0\\ 0& \mathrm{\Lambda }^{(n)}(\lambda )& 1& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& \mathrm{\Lambda }^{(n)}(\lambda )\end{array}\right)_{d_n\times d_n}$$
(67)
This turns out to be a remarkable example how boundary conditions can change in a drastic way the Hilbert space of integrable models. At this point it is natural to ask what happens to the Bethe ansatz states (28) when one gradually varies the boundary matrix towards the singular manifold. In particular, if we can figure out the kind of Bethe states in each sector $`n_\pm `$ that should collapse to the eigenvectors (64). A precise answer to this question as well as possible generalizations of these results for arbitrary $`N`$ has eluded us so far.

## Acknowledgements

We would like to thank the organizers of the Workshop Flux, Charge, Topology and Statistics 2003, Amsterdam, where part of this work was carried out. The authors G.A.P Ribeiro and W. Galleas thank Fapesp (Fundação de Amparo à Pesquisa do Estado de São Paulo) for financial support. The work of M.J. Martins has been partially support by the Brazilian Research Council-CNPq and Fapesp.

## Appendix A : Completeness for $`L=2`$

This appendix is concerned with the study of the completeness of the Bethe ansatz solution (31) for $`L=2`$, i.e. that all four eigenvalues of the transfer matrix are obtained either by starting with $`|0^{(+)}`$ or with $`|0^{()}`$. Let us first begin with $`|0^{(+)}`$ whose corresponding eigenvalue $`\mathrm{\Lambda }_0^{(+)}(\lambda )`$ is clearly

$$\mathrm{\Lambda }_0^{(+)}(\lambda )=g^{(+)}[a(\lambda )]^2+g^{()}[b(\lambda )]^2$$
(A.1)
The next step is to solve the Bethe ansatz equations for the one-particle state $`\stackrel{~}{B}(\lambda _1)|0^{(+)}`$. As a result we find two possible rapidities given by

$$\lambda _1^\pm =\frac{1}{2}\frac{\left(\sqrt{g^{(+)}}\pm \sqrt{g^{()}}\right)}{\left(\sqrt{g^{(+)}}\sqrt{g^{()}}\right)}$$
(A.2)
giving us the following one-particle $`\mathrm{\Lambda }_1^\pm `$ eigenvalues

$$\mathrm{\Lambda }_1^\pm =a(\lambda )b(\lambda )\left(g^{(+)}+g^{()}\right)\pm \sqrt{g^{(+)}g^{()}}$$
(A.3)
Repeating similar exercise for the two-particle state $`\stackrel{~}{B}(\lambda _1)\stackrel{~}{B}(\lambda _2)|0^{(+)}`$ we have

$$\lambda _{1,2}=\frac{g^{(+)}+g^{()}\pm 2I\sqrt{g^{(+)}g^{()}}}{2\left(g^{()}g^{(+)}\right)}$$
(A.4)
and the corresponding eigenvalue is

$$\mathrm{\Lambda }_2(\lambda )=g^{()}[a(\lambda )]^2+g^{(+)}[b(\lambda )]^2$$
(A.5)
An exact diagonalization of the transfer matrix (7) corroborates these four possible eigenvalues for $`L=2`$. Note also that (A.5) is exactly the eigenvalue associated to reference state $`|0^{()}`$. The others eigenvalues (A.1) and (A.3) are easily obtained from $`|0^{}`$ by noticing that to each solution $`\lambda _i^{(+)}`$ one can find the corresponding $`\lambda _i^{()}`$ through the reflection $`\lambda _i^{()}=\lambda _i^{(+)}`$ symmetry. We have also investigated numerically this problem for $`L=3,4`$ and found that both references states can lead to the complete spectrum of $`T(\lambda )`$.