text:
# Evolution, Explosion and Nucleosynthesis of Core Collapse Supernovae

## 1 Introduction

Some years ago we started a big project whose main objective was the study of the evolutionary properties of the massive stars and of their explosive yields. In the first paper of this series (Chieffi, Limongi & Straniero, 1998, hereinafter CLS98) we mainly discussed the FRANEC code and the presupernova evolution of a $`25\mathrm{M}_{}`$ star of solar metallicity. The nuclear network adopted in that paper extended up to $`{}_{}{}^{67}\mathrm{Zn}`$ and included 12 nuclei in the H burning, 25 in He burning and 149 in the more advanced phases. In the second paper (Limongi, Straniero & Chieffi, 2000, hereinafter LSC00) we presented the evolutionary properties of a set of solar metallicity stars (namely $`13,15,20,25\mathrm{M}_{}`$) with a network still extending up to Zn but with 41 nuclei in the H burning, 88 in the He burning and 179 in the more advanced burning phases. In the same paper we also presented the explosive yields (in the frame of the radiation dominated shock approximation) produced by stellar models of different initial metallicities, namely Z=0.02, Z=0.001 and Z=0. In the following papers (Limongi & Chieffi 2002; Chieffi & Limongi 2002a, hereinafter CL02a) we discussed in detail the yields produced by the first stellar generation and used them to interpret the observed surface abundances of the available ultra metal poor stars. In Chieffi & Limongi (2002b) we addressed the production of the radioactive nuclei by solar metallicity massive stars.

Here we present a new generation of stellar models and associated explosions. The present pre explosive models differ from the previous ones in several respects: first of all the input physics has been updated whenever possible, second the nuclear network has been extended up to $`{}_{}{}^{98}\mathrm{Mo}`$, so that we now follow 44 nuclei in H burning, 149 in He burning and 267 in the more advanced phases. Further, distinctive differences between the present and the previous sets of models concern a) the adoption in the current models of the Ledoux criterium to determine the stability of the H convective shell that forms towards the end of the central H burning phase and b) the adopted rate for the $`{}_{}{}^{12}\mathrm{C}(\alpha ,\gamma )^{16}\mathrm{O}`$ process. Though the effective rate for this process still shows an embarrassing large uncertainty, at variance with our previous works where we chose to adopt the rate quoted by Caughlan et al (1985), (but see Imbriani et al., 2001 and CL02a), in the present one we decided to adopt the best current estimate (Kunz et al., 2002) at its face value without any attempt to rescale it in any way.

As for the computation of the explosive burnings, a big effort has been done to build up a hydro code capable of following the passage of a shock wave through the mantle of a star. Hence we can now compute the explosive nucleosynthesis by means of either the hydro code and the radiation dominated shock approximation: a comparison between the yields obtained with the two different approximations is discussed in detail. We have also explored at some extent the dependence of the explosive yields on the final kinetic energy of the ejected mantle.

A further improvement has been the extension of the network up to the N=50 closure shell, which means the possibility to explore the production of the so called ”weak component” of the s-processes. Let us briefly remind that the solar system distribution above A=60 is the result of the combined effects of s-, r\- and p\- processes, i.e. of slow and rapid neutron captures (with respect to the decay times of the unstable nuclei near stability) plus a contribution from the p- captures. The fit to the s-only (and quasi s-only) nuclei, i.e. those produced only by slow neutron captures, showed (Ward, Newman and Clayton, 1976; Ward and Newman, 1978; Käppeler et al., 1982) that at least two different neutron fluxes are necessary to reproduce their observed abundances. A ”strong” pulsed neutron exposure (called ”main component”), usually associated to the thermal pulses occurring in the intermediate mass stars, would account for the s-nuclei above A=90 and a weaker single exposure (called ”weak component”), associated with the massive stars, would be responsible for the production of the s\- nuclei between A=60 and A=90. The nucleosynthesis of the s-nuclei in the massive stars has been addressed up to now in a number of papers, namely Couch, Schmiedekamp, and Arnett (1974); Lamb et al. (1977); Arnett and Thielemann (1985); Busso and Gallino (1985); Prantzos, Arnould and Arcoragi (1987); Langer, Arnould and Arcoragi (1989); Prantzos, Hashimoto and Nomoto (1990); Raiteri et al. (1991a, b); Baraffe, El Eid and Prantzos (1992); Raiteri et al. (1993); The, El Eid and Meyer (2000); Hoffman, Woosley and Weaver (2001); Rausher et al. (2002). In spite of this not negligible number, almost all of these papers considered only the contribution of the central He burning to the final yields, neglecting completely the possible contribution (positive and/or negative) of the further evolutionary phases including the explosive nucleosynthesis. Only two papers (Raiteri et al., 1991b, 1993) added the contribution of the C-burning shell up to the end of the 90’s. Just in two very recent papers the yields of the nuclei between A=60 and A=90 have been computed by considering the whole evolutionary history of the massive stars, including the explosive burnings, namely Hoffman, Woosley and Weaver (2001) and Rausher et al. (2002). Our work proceeds on the same guideline, in the sense that all the nuclei of interest have been fully included in both the evolutionary and explosive codes, but with the difference that our computations extend to a wider mass grid, i.e. between 13 and 35 $`\mathrm{M}_{}`$.

The paper is organized as follows: the latest version of the FRANEC code is described in section 2 while the hydrostatic evolution of the models is presented in section 3. The explosive nucleosynthesis computed by means of either the hydro code and the radiation dominated shock approximation is discussed in section 4. Section 5 presents a wide discussion of the present results. A final conclusion follows.

## 2 Stellar evolution code and input physics

The presupernova evolutions presented in this paper have been computed by means of the latest version (4.97) of the FRANEC package. The most important changes with respect to the previous release (CLS98 and LSC00) are the following.

The nuclear network has been extended up to Molybdenum and now it includes 40 isotopes (from neutrons to $`{}_{}{}^{30}\mathrm{Si}`$) in hydrogen burning, 149 isotopes (from neutrons to $`{}_{}{}^{98}\mathrm{Mo}`$) in helium burning and 267 isotopes (from neutrons to $`{}_{}{}^{98}\mathrm{Mo}`$) in all the more advanced burning phases. In total 282 isotopes (Table 1) and about 3000 reaction rates were explicitly included in the various nuclear burning stages. As in the previous versions, the nuclear network is fully coupled to the equations describing the physical structure of the star so that both the physical and chemical evolution due to the nuclear reactions are solved simultaneously. This means that the full network is used to compute the nuclear energy generation rate for each meshpoint and at each timestep. The system of differential equations is now solved by means of the Yale sparse matrix package that provides a speed up factor of $`2`$ compared to the LU decomposition (and subsequent backward and forward substitution) previously adopted. We have also removed the NSE approximation for temperatures larger than $`410^9\mathrm{K}`$ so that the full network is now used until the beginning of the core collapse of the star.

All the nuclear reaction rates are taken from the latest version of the REACLIB (Rausher & Thielemann, 2000) available also at T. Rauscher’s home page 111$`\mathrm{ℎ𝑡𝑡𝑝}://\mathrm{𝑞𝑢𝑎𝑠𝑎𝑟}.\mathrm{𝑝ℎ𝑦𝑠𝑖𝑘}.\mathrm{𝑢𝑛𝑖𝑏𝑎𝑠}.\mathrm{𝑐ℎ}/\mathrm{𝑡𝑜𝑚𝑚𝑦}`$ with some exceptions described below. REACLIB is a nuclear reaction rate library that contains fits to both experimental and theoretical rates. The main difference between this and the previous version of REACLIB (adopted in CLS98 and LSC00) consists in an upgrade of the theoretical fits, as presented by Rausher & Thielemann (2000). The experimental rates were not updated. Whenever possible we always chose the experimental rates.

The nuclear cross sections for the (n,$`\gamma `$) reactions are derived from the database of experimental values provided by Bao et al. (2000). These data cover the energy (temperature) range between 5 and 100 keV ($`5.8010^7`$ \- $`1.1610^9`$ K). Above 100 keV we use the theoretical rates computed by Rausher & Thielemann (2000), but rescaled to match the experimental values at 100 keV. The reverse rates, i.e. the ($`\gamma `$,n) ones, are those provided by Rausher, rescaled by the same multiplier factors used for the forward reactions.

The $`{}_{}{}^{12}\mathrm{C}(\alpha ,\gamma )^{16}\mathrm{O}`$ and $`{}_{}{}^{22}\mathrm{Ne}(\alpha ,\mathrm{n})^{25}\mathrm{Mg}`$ rates adopted here are the ones provided by Kunz et al. (2002) (the adopted rate) and by Jaeger et al. (2001) (the recommended rate) respectively.

The $`{}_{}{}^{18}\mathrm{O}(\alpha ,\gamma )^{22}\mathrm{Ne}`$ and $`{}_{}{}^{22}\mathrm{Ne}(\alpha ,\gamma )^{26}\mathrm{Mg}`$ rates are the recommended values given by Käppeler et al. (1994) for $`\mathrm{T9}0.6`$. Above this temperature we match the Caughlan & Fowler (1988) analytical relations to the experimental data.

The weak interaction rates are taken by Oda et al. (1994) for A=1 and $`17\mathrm{A}39`$; Fuller, Fowler & Newman (1982, 1985) for $`40\mathrm{A}45`$ (i.e. the same adopted in CLS98 and in LSC00); Langanke & Martínez-Pinedo (2000) for $`46\mathrm{A}62`$; Takahashi & Yokoi (1987) for $`\mathrm{A}63`$; REACLIB for all the nuclei not included in the above mentioned collections.

The initial distribution of the heavy elements has been assumed to be scaled solar and the relative abundances of the nuclear species above He are derived from Anders & Grevesse (1989). As usual the initial He abundance has been determined by fitting the present properties of the Sun (Straniero, Chieffi and Limongi, 1997).

The equation of state adopted in CLS98 and LSC00 for temperatures below $`10^6\mathrm{K}`$ has been replaced by the latest EOS and EOSPLUS tables given by Rogers, Swenson & Iglesias (1996) and Rogers (2001). Above $`\mathrm{T}=10^6\mathrm{K}`$ we are still using the same EOS described in LSC00.

The opacity tables are the same used in CLS98 and in LSC00: the radiative opacity coefficients are derived from Kurucz (1991) for $`T10^4\mathrm{K}`$, from Iglesias, Rogers and Wilson (1992) (OPAL) for $`10^4\mathrm{K}<T10^8\mathrm{K}`$, and from the Los Alamos Opacity Library (LAOL) (Huebner et al., 1977) for $`10^8\mathrm{K}<T10^{10}\mathrm{K}`$. The opacity coefficients due to the thermal conductivity are derived from Itoh et al. (1983).

As in CLS98 and LSC00, the borders of the convective zones are defined by means of the Schwarzschild criterion except for the H convective shell that forms towards the end of the central H-burning, where the Ledoux criterion has been adopted (Imbriani et al., 2001). According to this stability criterion, the H ”convective” shell remains essentially stable against the growth of convective motions so that no mixing occurs in this region. This choice is based on some observational constraints mainly related to the relative number of red and blue supergiants (Langer & Maeder, 1995). No mass loss has been taken into account.

As a last comment let us emphasize that the present computations are much more numerically accurate than the CLS98 and LSC00 ones since they have been obtained with a higher degree of both spatial and temporal resolution: the typical number of meshpoints is increased from 600 to 1500 while the number of timesteps in the advanced burning phases is roughly doubled (from 10000 to 20000).

## 3 The Hydrostatic Evolution

We computed the evolution of 6 models, namely 13, 15, 20, 25, 30 and $`35\mathrm{M}_{}`$, covering most of the mass interval that is limited at the lower end by the stars that ignite C in a partially degenerate environment and at the high one by the dominion of the Wolf Rayet stars. The initial chemical composition for all the models is Y=0.285 and Z=0.02. We have already largely discussed the evolutionary properties of the various models from the pre main sequence phase up to the final core collapse in earlier papers of this series (i.e. CLS98 and LSC00) and hence we will not repeat them here. For sake of completeness, let us mention that no breathing pulses develop during the last phase of central He burning because of the very fine spatial and temporal resolution (about 3000 models and 2000 meshpoints just to follow the central He burning phase). The key properties of these evolutions are reported in Table 2 and may be compared to the older ones reported in Table 2 of LSC00. Such a comparison shows that the new models do not differ significantly with respect to the older ones: the largest differences occur for the $`20\mathrm{M}_{}`$, but this is due to the quite unusual behavior of the old $`20\mathrm{M}_{}`$. The two key quantities that play a fundamental role in the determination of the final explosive yields, i.e. the electron mole number $`Y_\mathrm{e}`$ and the final mass-radius relation within the region where these explosive burnings occur, are shown in Figure 1 as a function of the radial coordinate: the old and new relations are shown as dotted and solid lines respectively while the thin and thick lines refer, respectively, to $`Y_\mathrm{e}`$ and to the mass coordinate M. This figure shows that $`Y_\mathrm{e}`$ is more or less systematically lower by $`10^4`$ in the new models while the M-R relation is now slightly steeper for the two smaller masses but somewhat shallower for the two more massive ones. However, the total mass exposed to the various explosive burnings (Chieffi, Limongi & Straniero, 2000) did not change by more than $`30\%`$ at most.

## 4 The Explosion

The present modeling of the core collapse supernovae in spherical symmetry did not yield to successful explosions yet (at least with the presently available microphysics). Even improvements in the neutrino transport (Mezzacappa et al., 2001) and the inclusion of general relativistic corrections (Liebendörfer et al., 2001) in the model calculations did not change this situation. For this reason the explosive nucleosynthesis calculations for core collapse supernovae are still based on explosions induced by injecting an arbitrary amount of energy in a (also arbitrary) mass location and then following the development of the blast wave with an hydro code. The amount of energy injected is tuned in order to obtain a prefixed amount of kinetic energy of the ejecta at the infinity (usually $`10^{51}\mathrm{erg}`$). Two different techniques are adopted at present: the thermal bomb and the piston. In the first approach the extra energy is deposited instantaneously in the form of internal energy at the base of the mantle of the exploding star. In the second one, conversely, the extra energy is given in the form of kinetic energy by assuming the inner edge of the mantle to move like a piston in the gravitational field of the compact remnant and to follow a ballistic trajectory. A third technique that may be adopted to compute the final explosive yields is the radiation dominated shock approximation (Weaver & Woosley, 1980; Arnett 1996; CL02a) that we adopted in our previous works. In the following we will firstly present the yields obtained by adopting this last technique and compare them to the analogous ones presented in LSC00. We then present our hydro code together to the properties of the explosions. A comparison between the explosive yields obtained with the hydro code and the ones computed in the radiation dominated shock approximation will close this section.

### 4.1 The radiation dominated shock approximation (RDA)

A full set of explosive yields have been computed in the radiation dominated shock approximation by (arbitrarily) imposing a final kinetic energy of 1.2 foe ($`1\mathrm{foe}=1\times 10^{51}\mathrm{erg}`$), the ejection of $`0.05\mathrm{M}_{}`$ of $`{}_{}{}^{56}\mathrm{Ni}`$ and a time delay of 0.5 s between the beginning of the core collapse and the rejuvenation of the shock wave. Figure 2 shows the comparison, for the four masses in common, between the LSC00 and the present yields. Since the numerical technique adopted in both the older and newer simulations of the explosion is identical, the differences are mainly due to the different presupernova structures and to the different nuclear reaction rates. This figure clearly shows that the differences remain confined within a factor of two for most of the isotopes, raising up to one order of magnitude for just few isotopes. Note that the largest differences between the new and the old computations do not show any specific clear trend with the mass and since the improvements to either the code and the input physics were applied all together, it is not easy to disentangle which is/are the specific changes that produced such differences. However, it is clear that the differences concerning the nuclei that are produced by the explosive burnings mainly reflect the differences in both the final Mass-Radius relation and the Ye profile. Viceversa the differences concerning the nuclei produced by the hydrostatic burnings will mainly reflect the changing in the reaction rates (the physical evolution of the stellar interiors did not change appreciably in the central H and He burning phases).

### 4.2 The hydrodynamical explosion

#### 4.2.1 The hydro code

The propagation of a shock front through the mantle of a star is followed by solving the hydrodynamical equations in spherical symmetry and in lagrangean form:

$$\frac{v}{t}=\frac{Gm}{r^2}4\pi r^2\frac{P}{m}4\pi \frac{(r^2Q)}{m}$$
(1)
$$\frac{r}{m}=\frac{1}{4\pi r^2\rho }$$
(2)
$$\frac{e}{t}=\frac{P}{\rho ^2}\frac{\rho }{t}4\pi r^2Q\frac{v}{m}$$
(3)
supplemented by the boundary condition:

$$\frac{v_1}{t}=\frac{Gm}{r^2}$$
(4)
that imposes the inner edge of the exploding mantle to move on a ballistic trajectory under the gravitational field of the compact remnant.

The technique adopted to solve the system of equations 1, 2 and 3, is the same described by Richtmeyer & Morton (1967) and Mezzacappa & Bruenn (1993). Few additional zones are added to the computational domain in order to allow both the forward and the reverse shock waves to leave the structure. Once the shock enters one of these zones, they are removed from the computational domain, simulating in this way the escape of the shock wave out of the expanding mantle. The chemical evolution of the matter is computed by coupling the same nuclear network adopted in the hydrostatic calculations to the system of equations 1, 2 and 3. The blast wave is started by imparting an initial velocity $`v_0`$ at a mass coordinate of $`1\mathrm{M}_{}`$, i.e. well inside the iron core, and $`v_0`$ is properly tuned in order to obtain a given final kinetic energy of the ejecta.

#### 4.2.2 Propagation of the shock, fall back and explosive nucleosynthesis

The main properties of the explosion of the presupernova models presented in this paper can be summarized by discussing one specific model taken as a representative case. Figure 3 shows the time history of the shock propagation in a $`25\mathrm{M}_{}`$ model characterized by an initial velocity $`v_0=1.55510^9\mathrm{cm}/\mathrm{s}`$ and by a final kinetic energy at the infinity of $`1.14410^{51}\mathrm{erg}`$. Once the shock forms, it propagates outward in mass increasing locally both the temperature and the density, triggering therefore the explosive nucleosynthesis. Behind the shock front both the pressure and the density (and hence the temperature) are fairly flat (occurrence largely exploited in the radiation dominated shock approximation) and progressively lower as the matter kicked by the shock wave expands and cools down. In $`4\mathrm{s}`$ the shock reaches the outer edge of the CO core; at this time the temperature in the shocked region has dropped down to $`10^9\mathrm{K}`$ and the explosive nucleosynthesis almost vanishes. At $`\mathrm{t}100\mathrm{s}`$ the inner zones of the mantle begin to fall back onto the compact remnant while the shock is still within the He core. The fall back lasts $`110^5\mathrm{s}`$ though most of the mass that will eventually return to the remnant actually falls back in just $`200\mathrm{s}`$. The shock reaches the He/H interface at $`\mathrm{t}370\mathrm{s}`$ and since this region corresponds to a very steep increase of $`\rho r^3`$, a reverse shock forms (Bethe, 1990). This reverse shock propagates inward in mass and decelerates somewhat the expanding matter that encounters in its way back: however, at least in our computations, this reverse shock does not affect significantly the amount of fall back. The main outgoing shock reaches the surface of the star at $`\mathrm{t}210^5\mathrm{s}`$, while the reverse shock escapes from the inner edge of the mantle at $`\mathrm{t}610^6\mathrm{s}`$. The further evolution of the mantle is characterized by an homologous expansion with velocity ranging between 1000 (inner zones) and 3000 Km/s (outer zones).

In order to explore the dependence of the explosion on the initial velocity $`v_0`$ we have computed several hydro simulations that are summarized in Table 3. The legend is as follows: the initial mass (column 1), a label to easily identify the model in the subsequent tables (column 2), the initial velocity $`v_0`$ (column 2), the times at which the shock reaches the CO core (column 3), the fall back starts (column 4) and stops (column 5), the shock reaches the He/H interface (or equivalently the time at which the reverse shock forms) (column 6), the forward shock escapes from the surface (column 7) and the reverse shock escapes from the inner edge of the mantle (column 8), the final kinetic energy at the infinity (column 9), the mass of the compact remnant (column 10) and the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ synthesized (column 11). The final explosive yields produced by most of these hydro simulations are reported in Tables 4-9 once the unstable isotopes have been decayed into their stable closest neighbor. The yields of all the $`\gamma `$-ray emitters (having a yield larger than $`10^{10}`$ solar masses) included in the network are collected in Table 10.

An interesting result that comes out from these hydro simulations is that the correlation between the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ ejected and the final kinetic energy of the ejecta weakens as the mass of the progenitor increases. Let us consider, for each mass, two limiting final kinetic energies $`\mathrm{E}_{\mathrm{kin}}`$: the first one, $`\mathrm{E}_{\mathrm{kin}}^{\mathrm{max}}`$, is the maximum kinetic energy that prevents the ejection of part of the iron core while the second one, $`\mathrm{E}_{\mathrm{kin}}^{\mathrm{min}}`$, is the minimum kinetic energy necessary to eject at least $`10^3\mathrm{M}_{}`$ of $`{}_{}{}^{56}\mathrm{Ni}`$. The range of kinetic energies defined in this way is much larger for the lower mass models than for the more massive ones. For example the $`\mathrm{\Delta }(\mathrm{E}_{\mathrm{kin}}^{\mathrm{max}}\mathrm{E}_{\mathrm{kin}}^{\mathrm{min}})0.63\mathrm{foe}`$ for the $`13\mathrm{M}_{}`$ model while it reduces to $`\mathrm{\Delta }(\mathrm{E}_{\mathrm{kin}}^{\mathrm{max}}\mathrm{E}_{\mathrm{kin}}^{\mathrm{min}})0.2\mathrm{foe}`$ for the $`35\mathrm{M}_{}`$ model. Viceversa the maximum amount of the $`{}_{}{}^{56}\mathrm{Ni}`$ ejected remains roughly constant over the mass range 13-30 $`\mathrm{M}_{}`$ and increases a while in the 35 $`\mathrm{M}_{}`$ case. This means that explosions leading to the ejection of significantly different amounts of $`{}_{}{}^{56}\mathrm{Ni}`$ imply also a significative scatter in the final kinetic energy of the mantle if the mass is not very large, while they would imply a substantially constant final kinetic energy if the star is massive enough.

To eject a prefixed amount of $`{}_{}{}^{56}\mathrm{Ni}`$ it is usually necessary to compute several hydro runs iterating on the initial velocity $`v_0`$. It would be much easier to run just an hydro model and obtain the given amount of $`{}_{}{}^{56}\mathrm{Ni}`$ (plus all the other yields) by imposing the mass cut a posteriori by hand. By means of the present runs we can check if this procedure is reliable or not. Let us check if, e.g., the yields provided by case 15B may be obtained from case 15F. Hence let us impose in the 15F case a mass cut such that the yield of $`{}_{}{}^{56}\mathrm{Ni}`$ is the same obtained in case 15B; let us call this hand made model 15FB. The comparison between the ”real” case 15B and the ”fictitious” case 15FB is shown in Figure 4 where the (percentage) differences between the two cases are shown in the upper panel. It is remarkable that the differences remain confined within 20% for all the isotopes except that $`{}_{}{}^{20,21}\mathrm{Ne}`$, $`{}_{}{}^{25}\mathrm{Mg}`$ and $`{}_{}{}^{46}\mathrm{Ca}`$. Also for these nuclei, however, the maximum difference never exceeds 50%. To span a range of different mass cuts, the central panel in Figure 4 shows the comparison between case 15C and a ”fictitious” case 15FC while the lower panel shows the comparison between case 15D and 15FD. As expected, the discrepancies progressively reduce as the imposed mass cut gets closer to the real one. Figure 5 shows a similar comparison for the $`35\mathrm{M}_{}`$ model. In this case the maximum differences never exceeds 5%.

### 4.3 Comparison between RDA and Hydro calculations

It is worth comparing now the explosive yields obtained with the RDA technique and the ones obtained by means of the hydro code. Since, at variance with the RDA calculations, the hydro code does not take into account the time delay ($`\tau _{delay}`$) between the core collapse and the explosion, we have recalculated the RDA explosions for all the models by imposing $`\tau _{delay}=0`$. Figure 6 shows, for the $`13\mathrm{M}_{}`$, the logarithmic ratio of the elemental yields obtained with the hydro code and the RDA technique for the same amount of ejected $`{}_{}{}^{56}\mathrm{Ni}`$. By the way, the elemental abundances are obtained by firstly decaying all the unstable nuclei to their parental stable neighbor and then summing up all the stable nuclei of each given element. Although the two techniques can be considered very different, it is surprising that the difference in the final yields are always confined within 0.1 dex for all the elements except for K, Ti, Co, Ni, Cu and Zn, i.e. the elements produced only (or mainly) in the zone that undergoes explosive Si burning with complete Si exhaustion. Note however that, even in the worst case, for these elements the differences between the two techniques remain confined within 0.5 dex, i.e. less than a factor of 3. A closer inspection to the temporal evolution of one meshpoint exposed to the complete explosive Si burning shows that the RDA technique overestimates both the peak temperature and the peak density and leads to a slower expansion compared to the hydro calculation. This is a well known feature of the radiation dominated shock approximation and it is also readily visible in, e.g., WW95 (their Figure 9). Although this occurrence can be generalized also to all the other zones, it mostly influences the region exposed to temperatures larger than 5 billion degrees, i.e. the ones characterized by the alpha rich freeze out. Such differences become less important as the mass of the progenitor star increases, hence the yields obtained with the two techniques become smaller as well. More specifically,the maximum difference is confined within 0.15 and 0.10 dex for the $`25\mathrm{M}_{}`$ and the $`35\mathrm{M}_{}`$ respectively. Figure 7 shows that even the isotopic differences between hydro and RDA calculations are nearly always confined within 0.1 dex with some notable exceptions, namely, $`{}_{}{}^{20}\mathrm{Ne}`$, $`{}_{}{}^{25}\mathrm{Mg}`$, $`{}_{}{}^{39}\mathrm{K}`$, $`{}_{}{}^{43,44,46}\mathrm{Ca}`$, $`{}_{}{}^{47,48}\mathrm{Ti}`$, $`{}_{}{}^{50}\mathrm{V}`$, $`{}_{}{}^{59}\mathrm{Co}`$, $`{}_{}{}^{58,60,61,62}\mathrm{Ni}`$, $`{}_{}{}^{63,65}\mathrm{Cu}`$, $`{}_{}{}^{64,66,67,70}\mathrm{Zn}`$ and $`{}_{}{}^{84}\mathrm{Sr}`$. $`{}_{}{}^{39}\mathrm{K}`$, $`{}_{}{}^{48}\mathrm{Ti}`$, $`{}_{}{}^{59}\mathrm{Co}`$, $`{}_{}{}^{58}\mathrm{Ni}`$, $`{}_{}{}^{63,65}\mathrm{Cu}`$ and $`{}_{}{}^{64,66,67}\mathrm{Zn}`$ are the largest (if not the only) contributors to their respective elemental abundances and hence the possible origin of the differences obtained for these isotopes has already been discussed above. $`{}_{}{}^{47}\mathrm{Ti}`$, $`{}_{}{}^{43}\mathrm{Ca}`$ and $`{}_{}{}^{60,61,62}\mathrm{Ni}`$ are also produced only, or mainly, in the region of the alpha rich freeze out hence the differences can be readily understood. $`{}_{}{}^{20}\mathrm{Ne}`$ and $`{}_{}{}^{25}\mathrm{Mg}`$ are produced during the hydrostatic evolution in the carbon convective shell and then destroyed in all the zones shocked to a temperature larger than $`210^9\mathrm{K}`$. As a consequence their final abundance will depend on the radial coordinate at which the peak temperature of $`210^9\mathrm{K}`$ is reached. Such a critical radius is more external in the RDA compared to the hydro calculations. $`{}_{}{}^{46}\mathrm{Ca}`$, $`{}_{}{}^{50}\mathrm{V}`$, $`{}_{}{}^{70}\mathrm{Zn}`$ and $`{}_{}{}^{84}\mathrm{Sr}`$ are produced by the explosion; they show a gaussian-like profile with a maximum corresponding to a mass coordinate characterized by a given peak temperature that is $`210^9\mathrm{K}`$ for $`{}_{}{}^{46}\mathrm{Ca}`$, $`{}_{}{}^{50}\mathrm{V}`$ and $`{}_{}{}^{70}\mathrm{Zn}`$, and $`310^9\mathrm{K}`$ for $`{}_{}{}^{84}\mathrm{Sr}`$. Both the mass location of this maximum and its width depend on the relation between the peak temperature and the interior mass. In general the RDA calculations tend to underproduce all of these three isotopes with respect to the hydro explosions, this effect being larger for lower explosion energies in the hydro models. $`{}_{}{}^{44}\mathrm{Ca}`$ is made by itself and by $`{}_{}{}^{44}\mathrm{Ti}`$ in proportions depending on the explosion energy - $`{}_{}{}^{44}\mathrm{Ti}`$ is synthesized in the alpha rich freeze out while $`{}_{}{}^{44}\mathrm{Ca}`$ is produced by the explosion in more external zones (like, e.g., $`{}_{}{}^{46}\mathrm{Ca}`$). For highest energies most of the zone undergoing alpha rich freeze out is ejected and $`{}_{}{}^{44}\mathrm{Ca}`$ is completely made by $`{}_{}{}^{44}\mathrm{Ti}`$. For lower explosion energies a large amount of fall back occurs and $`{}_{}{}^{44}\mathrm{Ca}`$ is made only by itself. Both techniques give similar yields for $`{}_{}{}^{44}\mathrm{Ca}`$ for the lowest explosion energies while the hydro technique overproduces it $`(^{44}\mathrm{Ti})`$ because of a higher degree of alpha rich freeze out. All the differences discussed above progressively and significantly reduce as the mass of the progenitor increases.

## 5 Discussion and Conclusions

The yields provided by a generation of solar metallicity stars should be used, in principle, to interpret the ejecta of a core collapse supernova of similar initial metallicity as well as a part of a database of yields to be included into a galactic chemical evolutionary code. In principle there is no reason to require that the ejecta of a generation of solar metallicity stars preserve solar relative proportions because the Solar System distribution is the result of the cumulative contribution of many generations of stars of very different metallicities. Nonetheless, it is generally assumed the production factors (PFs) of a generation of solar metallicity stars to be essentially flat. This is the consequence of the (reasonable) assumption that the average metallicity Z grows slowly and continuously with respect to the evolutionary timescales of the stars that contribute to the environment enrichment; if this is true, the stars that mostly contribute to the abundances of the various nuclei at a given metallicity $`\mathrm{Z}_0`$ are those whose initial Z is quite close to $`\mathrm{Z}_0`$ (just to give an idea, the stars in the metallicity range $`\mathrm{Z}_0>\mathrm{Z}\mathrm{Z}_0/2`$ contribute to roughly $`50\%`$ of the metallicity $`\mathrm{Z}_0`$, while stars in the range $`\mathrm{Z}_0>\mathrm{Z}\mathrm{Z}_0/10`$ contribute to the metallicity $`\mathrm{Z}_0`$ for more than the $`90\%`$). Since the assumption mentioned above is certainly well verified at least for the massive stars, and very probably also for the intermediate mass ones, it follows that most of the Solar System distribution is the result (as a first approximation) of the ejecta of ”quasi” solar metallicity stars. This is the reason why it is desirable that a generation of solar metallicity stars provides yields in roughly solar proportions or, in other words, that the PFs of the various nuclei remain roughly flat. Since Oxygen is produced only by the core collapse supernovae and it is also the most abundant element produced by these stars, it is convenient to use its PF as the one that better marks the overall increase of the average ”metallicity” and to verify if the other nuclei follow or not its behavior. Arbitrarily we chose a factor of two as a suitable warning threshold in the sense that we will assume that all the nuclei whose PF falls within a factor of two of the Oxygen one are compatible with a flat distribution while those outside this compatibility range deserve a closer look- up and may potentially constitute a problem.

Let us firstly concentrate on the elemental PFs (Figure 8). The symbols refer to the 6 masses while the two lines refer to a generation of massive stars having a Salpeter mass function $`(dn/dmm^{2.35})`$ but different choices for the mass cut. The dotted line (hereinafter case ”Flat”) refers to the case in which all the core collapse supernovae are assumed to eject 0.05 $`\mathrm{M}_{}`$ of $`{}_{}{}^{56}\mathrm{Ni}`$ independently on the initial mass, while the solid one (hereinafter case ”Trend”) is obtained by assuming the following dependence of the mass cut (actually the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ ejected) on the initial mass: $`13\mathrm{M}_{}`$ $`(0.15\mathrm{M}_{}\mathrm{of}{}_{}{}^{56}\mathrm{Ni})`$, $`15\mathrm{M}_{}`$ $`(0.10\mathrm{M}_{}\mathrm{of}{}_{}{}^{56}\mathrm{Ni})`$, $`20\mathrm{M}_{}`$ $`(0.08\mathrm{M}_{}\mathrm{of}{}_{}{}^{56}\mathrm{Ni})`$, $`25\mathrm{M}_{}`$ $`(0.07\mathrm{M}_{}\mathrm{of}{}_{}{}^{56}\mathrm{Ni})`$, $`30\mathrm{M}_{}`$ $`(0.05\mathrm{M}_{}\mathrm{of}{}_{}{}^{56}\mathrm{Ni})`$, $`35\mathrm{M}_{}`$ $`(0.05\mathrm{M}_{}\mathrm{of}{}_{}{}^{56}\mathrm{Ni})`$. Since there is not a subset of models in our grid of explosions ejecting the required $`{}_{}{}^{56}\mathrm{Ni})`$ abundances as a function of the mass, we have obtained the desired models by choosing by hand the mass cut in models 13B, 15C, 20C, 25C, 30C and 35E. The hydro models shown in Figure 8 are those used to build up the ”Flat” case and hence are those that eject 0.05 $`\mathrm{M}_{}`$ of $`{}_{}{}^{56}\mathrm{Ni}`$. Since we are focusing on the relative scaling of the PFs with respect to that of the Oxygen, all the PFs have been renormalized by imposing $`\mathrm{Log}(\mathrm{PF}_{\mathrm{Oxygen}})=0`$. This means that all the nuclei having a Log(PF) close to zero in Figure 8 preserve a scaled solar proportion with respect to the O. To facilitate the comparison with the yields provided by other authors we report in the first row of Table 11 our Oxygen PFs for the various masses while in rows 2 to 4 of the same Table we quote the yields given by Rausher et al. (2002, RHHW02), Thielemann, Nomoto and Hashimoto (1996, TNH96) and Woosley & Weaver (1995, WW95). The WW95 yields are all of the ”A” variety while the TNH96 ones are those computed with the original Ye. The first things worth noting in Figure 8 are the following: a) the yields produced by a generation of massive stars having a Salpeter mass function depend mainly on the yields of the masses between 20 and 25 $`\mathrm{M}_{}`$ (Weaver, Zimmerman & Woosley, 1978); b) the only elements that vary significantly between the cases ”Flat” (dotted) and ”Trend” (solid) are Fe and Ni and, at a smaller extent, also Ti, Co and Zn; c) the majority of the elements from C to Sr have PFs compatible (in the sense mentioned above) with that of the Oxygen: the exceptions are N, F, Na, K, Ti, Fe and Ni. N and F are largely underabundant (as expected) because our yields do not include neither the contributions from the intermediate mass stars (N) nor the neutrino induced reactions (F). The overproduction of Na (and partly of Al) is not totally unexpected because these two elements are ”secondary”, i.e. are elements whose production directly depends on the initial metallicity: the higher the initial metallicity the larger their final production. Hence their slight overproduction could simply indicate that the average metallicity to start with in order to get the solar system distribution should not be solar but slightly subsolar. By the way, similar results have been obtained by RHHW02. Let us also remark that the yields of these two elements strongly depend on the C abundance left by the central He burning (Imbriani et al., 2001) and in this paper we did not make any attempt to tune the $`{}_{}{}^{12}\mathrm{C}(\alpha ,\gamma )^{16}\mathrm{O}`$ reaction rate (nor the behavior of the convective core in the central He burning phase). K (produced mainly as $`{}_{}{}^{39}\mathrm{K}`$ by both the explosive O burning and the complete explosive Si burning) is certainly significantly underproduced in our models but such a disagreement could be alleviated if the neutrino induced processes were efficient enough (WW95 found, for example, that the inclusion of the neutrino irradiation could increase the final abundance of $`{}_{}{}^{39}\mathrm{K}`$ up to a factor of two, see their Table 4). Ti, Fe and Ni depend significantly on the adopted mass cut and a changing from the ”Flat” to the ”Trend” case increases the yields of all these elements pushing them towards a closer scaled solar distribution. It must be noted, however, that the ”Trend” case leaves less room for the Type Ia contribution because of the larger amount of Fe $`(^{56}\mathrm{Ni})`$ produced (see below). There are few other things worth mentioning: since C has a PF close to zero, these models would predict that most of the C in the Solar System Distribution would come from the core collapse supernovae. The Iron peak elements that are synthesized only by the incomplete explosive Si burning (i.e. V, Cr and Mn) show PF’s in good agreement with that of the O and do not show any dependence on the mass cut because the changes in the mass cut are confined (in the present test) to the region of the complete explosive Si burning. Sc does not show any dependence on the mass cut because in the range $`2025\mathrm{M}_{}`$ the yield of this element is dominated by the $`{}_{}{}^{45}\mathrm{Sc}`$ itself (that is produced in the central He burning and in the C burning shell) and not by the $`{}_{}{}^{45}\mathrm{Ca}`$ (that is synthesized by the complete explosive Si burning). Ti is produced mainly as $`{}_{}{}^{48}\mathrm{Cr}`$ and comes essentially from the complete and incomplete explosive Si burning. Co is produced directly as $`{}_{}{}^{59}\mathrm{Co}`$ (by the central He burning and the C burning shell) and as $`{}_{}{}^{59}\mathrm{Ni}`$ (by the complete explosive Si burning): in the ”Flat” case the Co yield reflects mainly the hydrostatic production of $`{}_{}{}^{59}\mathrm{Co}`$ while in the ”Trend” case there is also the contribution from the $`{}_{}{}^{59}\mathrm{Ni}`$. The Ni abundance is dominated by $`{}_{}{}^{58}\mathrm{Ni}`$ and $`{}_{}{}^{60}\mathrm{Ni}`$: the first of the two is destroyed by both the central He burning and by the convective C shell and is largely produced by both the complete and incomplete explosive Si burning; the second one is produced by the central He burning, the C and Ne explosive burnings and by the complete explosive Si burning. Cu is produced in similar amount as $`{}_{}{}^{63}\mathrm{Cu}`$, $`{}_{}{}^{63}\mathrm{Ni}`$ and $`{}_{}{}^{65}\mathrm{Cu}`$; all these three isotopes are produced by the central He burning and the convective C shell. Zn is produced in comparable amounts as $`{}_{}{}^{64}\mathrm{Zn}`$, $`{}_{}{}^{66}\mathrm{Zn}`$ and $`{}_{}{}^{68}\mathrm{Zn}`$. The first of these three nuclei has a complex production history, since it is firstly largely produced by the central He burning but then it is significantly destroyed by the convective C- shell; it is also produced in the layers experiencing an $`\alpha `$ rich freeze out in complete explosive Si burning. It is therefore clear that its final outcome will depend on a delicate balance between size of the convective core, outer border of the last convective C-shell and amount of fall back, if we limit ourself to the uncertainties connected with the physical evolution of the stars. In the present computations roughly 30% of the $`{}_{}{}^{64}\mathrm{Zn}`$ yield comes from the complete explosive Si burning (for the mass cut adopted in the ”Trend” case). The other two nuclei are, viceversa, produced by both the central He burning and the convective C shell. The PFs of the s-weak component, i.e. the elements from Ga to Sr, are compatible with those of the O and hence these models could provide the bulk of the abundances of these elements leaving a marginal role to the intermediate mass stars. The elements heavier than Sr are (correctly) largely underproduced because they come from the intermediate mass stars that are not included in the mass interval analyzed in this paper.

Let us now turn to the isotopic PFs. Figure 9 is the isotopic version of Figure 8. As one could expect the main differences between the two integrated PFs concern the nuclei synthesized by the complete explosive Si burning, i.e. $`{}_{}{}^{44}\mathrm{Ti}`$ $`(^{44}\mathrm{Ca})`$, $`{}_{}{}^{48}\mathrm{Cr}`$ $`(^{48}\mathrm{Ti})`$, $`{}_{}{}^{56}\mathrm{Ni}`$ $`(^{56}\mathrm{Fe})`$, $`{}_{}{}^{57}\mathrm{Ni}`$ $`(^{57}\mathrm{Fe})`$, $`{}_{}{}^{59}\mathrm{Ni}`$ $`(^{59}\mathrm{Co})`$, $`{}_{}{}^{58}\mathrm{Ni}`$, $`{}_{}{}^{60}\mathrm{Ni}`$, $`{}_{}{}^{61}\mathrm{Ni}`$, $`{}_{}{}^{62}\mathrm{Ni}`$, $`{}_{}{}^{64}\mathrm{Zn}`$ and $`{}_{}{}^{66}\mathrm{Zn}`$. Most of the nuclei in the range $`{}_{}{}^{12}\mathrm{C}`$ \- $`{}_{}{}^{64}\mathrm{Ni}`$ are compatible with a flat distribution of the PFs, with 17 nuclei being significantly underproduced and 4 nuclei largely overproduced. While the nuclei underproduced by the massive stars could come from other astrophysical sites (like, e.g., the $`{}_{}{}^{14}\mathrm{N}`$ that is very probably produced by the intermediate mass stars), it is clear that it is much more difficult to account for the nuclei overproduced; our models predict only 4 nuclei overabundant by more than a factor of two with respect to the O: they are $`{}_{}{}^{23}\mathrm{Na}`$, $`{}_{}{}^{40}\mathrm{K}`$, $`{}_{}{}^{46}\mathrm{Ca}`$ and $`{}_{}{}^{62}\mathrm{Ni}`$. The PF of the $`{}_{}{}^{23}\mathrm{Na}`$ has been already discussed above. $`{}_{}{}^{40}\mathrm{K}`$ is an unstable nucleus with an half life of 1.25 Gyr; since the age of the Sun is 4.6 Gyr, the present solar abundance of $`{}_{}{}^{40}\mathrm{K}`$ must be increased by almost a factor of 10 to obtain its abundance at the moment of the solar system formation. If we apply such a correction to the solar abundance of $`{}_{}{}^{40}\mathrm{K}`$, its PF lowers to $`\mathrm{PF}(^{40}\mathrm{K})0.4`$ becoming now significantly underproduced with respect to O. $`{}_{}{}^{46}\mathrm{Ca}`$ is produced by both the convective C-shell and the Ne explosive burning: in both cases the fuel that feeds its production is the neutron flux that, in turn, largely depends on the amount of $`{}_{}{}^{22}\mathrm{Ne}`$ available and hence on the initial metallicity Z; this could mean that, as already suggested for Na, the schematic assumption that the solar system distribution may be attributed to the ejecta of solar metallicity stars fails significantly in this case. $`{}_{}{}^{62}\mathrm{Ni}`$ is produced by the complete explosive Si burning and hence its PF largely depends on the adopted mass cut: it is indeed overproduced only in the ”Trend” case.

Also the PFs of the nuclei between Cu and Sr are in reasonable agreement with that of the O, confirming the largely accepted idea that these elements are significantly produced by massive stars. 7 isotopes in this range scatter by more than a factor of 2 relative to O: $`{}_{}{}^{64}\mathrm{Zn}`$, $`{}_{}{}^{76}\mathrm{Ge}`$, $`{}_{}{}^{80}\mathrm{Se}`$, $`{}_{}{}^{82}\mathrm{Se}`$, $`{}_{}{}^{78}\mathrm{Kr}`$, $`{}_{}{}^{71}\mathrm{Ga}`$ and $`{}_{}{}^{80}\mathrm{Kr}`$. The first 5 are underproduced while the last two are overproduced. $`{}_{}{}^{76}\mathrm{Ge}`$ and $`{}_{}{}^{82}\mathrm{Se}`$ have a strong r-process contribution, $`{}_{}{}^{78}\mathrm{Kr}`$ has probably a strong p-process contribution while the abundance of $`{}_{}{}^{80}\mathrm{Kr}`$ strongly depends on the lifetime of $`{}_{}{}^{79}\mathrm{Se}`$ that is still fairly uncertain for temperatures larger than $`710^8\mathrm{K}`$. It has been long recognized the block of elements between Cu and the neutron closure shell at N=50 cannot be produced mainly by the main s\- component but that they require the existence of a weak s-component (early attributed to the central He burning of massive stars). However most of the matter exposed to the central He burning phase will be also exposed to further burnings, i.e. the C-burning and the passage of the shock wave. The contribution of these three phases to the yields are reported in Table 11. A clear trend exists: the main production site of most of these nuclei progressively shifts from the explosive burnings towards the C convective shell and central He burning as the mass of the progenitor increases. In particular, in the 13 and 15 $`\mathrm{M}_{}`$ models all these isotopes are produced by the explosion itself. Viceversa, in the 20 - 35 $`\mathrm{M}_{}`$ range most of these nuclei are produced by either the C-burning shell and/or the central He burning, the relative contribution privileging the C- burning shell for the smaller masses and the central He burning for the larger ones. The negative percentages that appear in Table 11 refer to the percentage destroyed by that evolutionary phase: for example, in the $`30\mathrm{M}_{}`$, the $`\mathrm{C}_{\mathrm{sh}}`$ burning destroys half of the $`{}_{}{}^{64}\mathrm{Zn}`$ produced by the central He burning.

As for the contribution of the Type Ia supernovae to the Solar System distribution, we have tentatively adopted the yields provided by Iwamoto et al. (1999). In this paper a variety of models spanning different values of central ignition density, flame speed and deflagration to detonation critical density have been presented. By evaluating the pros and the cons of their various models, Iwamoto et al. (1999) concluded that the CS15DD2 is the one that better reproduces the overall properties of type Ia supernovae (we refer the reader to the quoted paper for a deeper insight in the properties of these models). Though we have adopted the CS15DD2 model as the typical one describing the yields provided by the type Ia supernovae, even the adoption of the other models presented in that paper would not change the main results discussed below. The percentage frequency of Type Ia supernovae relative to the core collapse ones has been fixed by imposing $`\mathrm{PF}_{\mathrm{Fe}}`$ = $`\mathrm{PF}_\mathrm{O}`$. In the ”Trend” case we find that this contribution is 12% while in the ”Flat” case it increases to 16%. Figure 10 shows the quoted PF’s with (filled dots connected by a solid line) and without (open dots connected by a dashed line) the Type Ia contribution for the ”Trend” case. From this figure it is quite clear that the Type Ia supernovae contribute only to the Solar System abundances of the nuclei in the range Ti-Ni. Moreover, the inclusion of the Type Ia supernovae brings both $`{}_{}{}^{50}\mathrm{Ti}`$ and $`{}_{}{}^{54}\mathrm{Cr}`$ within the compatibility band that now includes all the nuclei in the range Ti-Ni, with the exception of $`{}_{}{}^{62}\mathrm{Ni}`$ whose discrepancy is, viceversa, slightly increased. Note that, at variance with current believing, Type Ia supernova do not contribute appreciably to the synthesis of $`{}_{}{}^{48}\mathrm{Ca}`$ which hence remains largely unexplained. Our models already predict $`2.48\times 10^6`$ $`\mathrm{M}_{}`$ of $`{}_{}{}^{48}\mathrm{Ca}`$ per Type II supernovae while the Type Ia supernovae (Iwamoto et al., 1999) produce at most $`1.64\times 10^6`$ $`\mathrm{M}_{}`$ of $`{}_{}{}^{48}\mathrm{Ca}`$ per Type Ia explosion: since there is room just for only 12% of Type Ia supernova, it is readily understandable why these Type Ia yields do not change appreciably the total yield of this nucleus.

As a last comment, a changing of the slope of the Salpeter mass function from $`\alpha `$=2.35 to $`\alpha `$=0 would change the percentage of the Type Ia supernovae relative to the core collapse ones from 12% to 28%; such a changing would not alter significantly the present results.

We warmly thank Tony Mezzacappa for very helpful discussions, Franz Käppeler for useful discussions about the lifetime of $`{}_{}{}^{79}\mathrm{Se}`$ and Roberto Gallino for having pushed us to start these computations. We also thank Brad Gibson and John Lattanzio for their kind hospitality in Melbourne.


format: text

text:
# 1 Introduction

## 1 Introduction

Consider the open spin-$`\frac{1}{2}`$ XXZ quantum spin chain with nondiagonal boundary terms, defined by the Hamiltonian

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}\{{\displaystyle \underset{n=1}{\overset{N1}{}}}(\sigma _n^x\sigma _{n+1}^x+\sigma _n^y\sigma _{n+1}^y+\mathrm{cosh}\eta \sigma _n^z\sigma _{n+1}^z)`$ (1.1)
$`+`$ $`\mathrm{sinh}\eta [\mathrm{coth}\xi _{}\sigma _1^z+{\displaystyle \frac{2\kappa _{}}{\mathrm{sinh}\xi _{}}}(\mathrm{cosh}\theta _{}\sigma _1^x+i\mathrm{sinh}\theta _{}\sigma _1^y)`$
$``$ $`\mathrm{coth}\xi _+\sigma _N^z{\displaystyle \frac{2\kappa _+}{\mathrm{sinh}\xi _+}}(\mathrm{cosh}\theta _+\sigma _N^x+i\mathrm{sinh}\theta _+\sigma _N^y)]\},`$
where $`\sigma ^x,\sigma ^y,\sigma ^z`$ are the usual Pauli matrices, $`\eta `$ is the bulk anisotropy parameter, $`\xi _\pm ,\kappa _\pm ,\theta _\pm `$ are arbitrary boundary parameters, and $`N`$ is the number of spins. This is the prototypical integrable quantum spin chain with boundary. It is related to many other models, including the sine-Gordon field theory . Moreover, this model has applications in various branches of physics, including condensed matter and statistical mechanics.

This model has resisted solution for many years (see, e.g., ). The main difficulty is that, in contrast to the special case of diagonal boundary terms (i.e., $`\kappa _\pm =0`$) considered in , a simple pseudovacuum (reference) state does not exist. For example, the state with all spins up $`\left(\genfrac{}{}{0pt}{}{1}{0}\right)^N`$ is not an eigenstate of the Hamiltonian.

We recently formulated  a method of deriving the Bethe Ansatz solution of integrable spin chain (vertex-type) models which does not rely on the existence of a pseudovacuum state. In particular, we used this method to solve the model (1.1) for the special case

$`\kappa _+=\kappa _{},\xi _+=\xi _{},\theta _+=\theta _{}=0,N=\text{odd}.`$ (1.2)
Here we propose the solution for a more general case. Indeed, in terms of the boundary parameters $`\alpha _{},\beta _{}`$ introduced below in Eq. (3.30), we find an expression for the eigenvalues of the transfer matrix corresponding to the Hamiltonian

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}\{{\displaystyle \underset{n=1}{\overset{N1}{}}}(\sigma _n^x\sigma _{n+1}^x+\sigma _n^y\sigma _{n+1}^y+\mathrm{cosh}\eta \sigma _n^z\sigma _{n+1}^z)`$ (1.3)
$`+`$ $`\mathrm{sinh}\eta [\mathrm{coth}\alpha _{}\mathrm{tanh}\beta _{}\sigma _1^z+csch\alpha _{}sech\beta _{}(\mathrm{cosh}\theta _{}\sigma _1^x+i\mathrm{sinh}\theta _{}\sigma _1^y)`$
$``$ $`\mathrm{coth}\alpha _+\mathrm{tanh}\beta _+\sigma _N^z+csch\alpha _+sech\beta _+(\mathrm{cosh}\theta _+\sigma _N^x+i\mathrm{sinh}\theta _+\sigma _N^y)]\},`$
where the boundary parameters are subject to the linear relation

$`\alpha _{}+\beta _{}+\alpha _++\beta _+=\pm (\theta _{}\theta _+)+\eta k,`$ (1.4)
where $`k`$ is an even integer if $`N`$ is odd, and is an odd integer if $`N`$ is even. In the recent paper , similar results have been obtained by a different approach.

The outline of this article is as follows. In Section 2, we briefly review the construction of the model’s transfer matrix, and list some of its important properties. In Section 3, we find the eigenvalues of the transfer matrix by the three-step procedure formulated in . The first two steps, which lead to a functional relation for the transfer matrix, are the same as in , except for the introduction of the parameters $`\theta _{}`$. The principal new results appear at the third step, where we succeed to recast the functional relation in terms of a determinant for the more general case (1.4). We conclude with a brief discussion of our results in Section 4.

## 2 The transfer matrix

The fundamental transfer matrix $`t(u)`$ corresponding to the model (1.1) is given by

$`t(u)=tr_0K_0^+(u)T_0(u)K_0^{}(u)\widehat{T}_0(u),`$ (2.1)
where the monodromy matrices are given by

$`T_0(u)=R_{0N}(u)\mathrm{}R_{01}(u),\widehat{T}_0(u)=R_{10}(u)\mathrm{}R_{N0}(u),`$ (2.2)
and the $`R`$ matrix is the solution of the Yang-Baxter equation given by

$`R(u)=\left(\begin{array}{cccc}\mathrm{sinh}(u+\eta )& 0& 0& 0\\ 0& \mathrm{sinh}u& \mathrm{sinh}\eta & 0\\ 0& \mathrm{sinh}\eta & \mathrm{sinh}u& 0\\ 0& 0& 0& \mathrm{sinh}(u+\eta )\end{array}\right).`$ (2.7)
Moreover, the $`K^{}`$ matrix is the solution of the boundary Yang-Baxter equation  given by

$`K^{}(u)=\left(\begin{array}{cc}\mathrm{sinh}(\xi _{}+u)& \kappa _{}e^\theta _{}\mathrm{sinh}2u\\ \kappa _{}e^\theta _{}\mathrm{sinh}2u& \mathrm{sinh}(\xi _{}u)\end{array}\right),`$ (2.10)
which evidently depends on three boundary parameters $`\xi _{},\kappa _{},\theta _{}`$. It is related to the symmetric matrix $`K^{}(u)|_{\theta _{}=0}`$ used in  by a gauge transformation,

$`K^{}(u)=K^{}(u)|_{\theta _{}=0}^1,`$ (2.11)
with

$`=\left(\begin{array}{cc}e^{\frac{1}{2}\theta _{}}& 0\\ 0& e^{\frac{1}{2}\theta _{}}\end{array}\right).`$ (2.14)
The matrix $`K^+(u)`$ is equal to $`K^{}(u\eta )`$ with $`(\xi _{},\kappa _{},\theta _{})`$ replaced by $`(\xi _+,\kappa _+,\theta _+)`$. Finally, $`tr_0`$ denotes trace over the (two-dimensional) “auxiliary space” 0. Further details about the construction of this transfer matrix can be found in .

The transfer matrix constitutes a one-parameter commutative family of matrices

$`[t(u),t(v)]=0.`$ (2.15)
The Hamiltonian (1.1) is related to the first derivative of the transfer matrix,

$`=c_1{\displaystyle \frac{}{u}}t(u)|_{u=0}+c_2𝕀,`$ (2.16)
where

$`c_1={\displaystyle \frac{1}{4\mathrm{sinh}\xi _{}\mathrm{sinh}\xi _+\mathrm{sinh}^{2N1}\eta \mathrm{cosh}\eta }},c_2={\displaystyle \frac{\mathrm{sinh}^2\eta +N\mathrm{cosh}^2\eta }{2\mathrm{cosh}\eta }},`$ (2.17)
and $`𝕀`$ is the identity matrix. The two relations (2.15),(2.16) signal that the model is integrable. Moreover, it is evident that in order to determine the energy eigenvalues, it suffices to determine the eigenvalues of the transfer matrix.

The transfer matrix has the periodicity property

$`t(u+i\pi )=t(u),`$ (2.18)
as well as crossing symmetry

$`t(u\eta )=t(u),`$ (2.19)
and the asymptotic behavior (for $`\kappa _\pm 0`$)

$`t(u)\kappa _{}\kappa _+\mathrm{cosh}(\theta _{}\theta _+){\displaystyle \frac{e^{u(2N+4)+\eta (N+2)}}{2^{2N+1}}}𝕀+\mathrm{}\text{for}u\mathrm{}.`$ (2.20)
## 3 Bethe Ansatz solution

We now proceed to find an expression for the transfer matrix eigenvalues using the method formulated in . This method consists of three main steps:

### 3.1 Step 1: fusion hierarchy

The first step is to obtain the model’s so-called fusion hierarchy. The transfer matrix (2.1) is actually the first ($`j=\frac{1}{2}`$) member of an infinite hierarchy of commuting transfer matrices $`t^{(j)}(u)`$ corresponding to spin-$`j`$ (i.e., $`(2j+1)`$-dimensional) auxiliary spaces, $`j=\frac{1}{2},1,\frac{3}{2},\mathrm{}`$. Using the fusion procedure for $`R`$  and $`K`$  matrices, one finds that these higher-level transfer matrices obey the relations

$`t^{(j)}(u)`$ $`=`$ $`\stackrel{~}{\zeta }_{2j1}(2u+(2j1)\eta )[t^{(j\frac{1}{2})}(u)t^{(\frac{1}{2})}(u+(2j1)\eta )`$ (3.1)
$``$ $`{\displaystyle \frac{\mathrm{\Delta }(u+(2j2)\eta )\stackrel{~}{\zeta }_{2j2}(2u+(2j2)\eta )}{\zeta (2u+2(2j1)\eta )}}t^{(j1)}(u)],`$
with $`t^{(0)}=𝕀`$, and $`j=1,\frac{3}{2},\mathrm{}`$. The quantum determinant $`\mathrm{\Delta }(u)`$ is given by

$`\mathrm{\Delta }(u)`$ $`=`$ $`\left[\mathrm{sinh}(u+\eta +\xi _{})\mathrm{sinh}(u+\eta \xi _{})+\kappa _{}^2\mathrm{sinh}^2(2u+2\eta )\right]`$ (3.2)
$`\times `$ $`\left[\mathrm{sinh}(u+\eta +\xi _+)\mathrm{sinh}(u+\eta \xi _+)+\kappa _+^2\mathrm{sinh}^2(2u+2\eta )\right]`$
$`\times `$ $`\mathrm{sinh}2u\mathrm{sinh}(2u+4\eta )\zeta (u+\eta )^{2N},`$
and

$`\stackrel{~}{\zeta }_j(u)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{j}{}}}\zeta (u+k\eta ),\stackrel{~}{\zeta }_0(u)=1,`$ (3.3)
$`\zeta (u)`$ $`=`$ $`\mathrm{sinh}(u+\eta )\mathrm{sinh}(u\eta ).`$ (3.4)
These relations are the same as those for the case of symmetric $`K`$ matrices ($`\theta _{}=0`$) . We remark that the spin-$`j`$ matrix $`K_{1\mathrm{}2j}^{}(u)`$ is related to the corresponding matrix with $`\theta _{}=0`$ by a generalization of the gauge transformation (2.11),

$`K_{1\mathrm{}2j}^{}(u)=_1\mathrm{}_{2j}K_{1\mathrm{}2j}^{}(u)|_{\theta _{}=0}_{2j}^1\mathrm{}_1^1.`$ (3.5)
### 3.2 Step 2: truncation at roots of unity

The second step is to observe  that for anisotropy values

$`\eta ={\displaystyle \frac{i\pi }{p+1}},p=1,2,\mathrm{},`$ (3.6)
(and hence $`qe^\eta `$ is a root of unity, satisfying $`q^{p+1}=1`$), the level-$`\frac{p+1}{2}`$ transfer matrix can be expressed in terms of a transfer matrix of one level lower,

$`t^{(\frac{p+1}{2})}(u)=\alpha (u)\left[t^{(\frac{p1}{2})}(u+\eta )+\beta (u)𝕀\right].`$ (3.7)
The quantities $`\alpha (u)`$ and $`\beta (u)`$ are given by the corresponding expressions (4.31) in , except with $`\sigma _{}(u)e^{(p+1)\theta _{}}\sigma _{}(u)`$ and $`\rho _{}(u)e^{(p+1)\theta _{}}\rho _{}(u)`$, as a consequence of (3.5).

This result provides an example of McCoy’s dictum “Complicated is simple” . Indeed, the essential point of this step is to exploit the higher symmetry which occurs at roots of unity to help solve the model.

Combining the fusion hierarchy (3.1) and the truncation identity (or “closing relation”) (3.7) for the $`\eta `$ values (3.6), we arrive at a functional relation for the fundamental transfer matrix $`t(u)t^{(\frac{1}{2})}(u)`$ (and hence, for the corresponding eigenvalues $`\mathrm{\Lambda }(u)`$) of order $`p+1`$ :

$`\mathrm{\Lambda }(u)\mathrm{\Lambda }(u+\eta )\mathrm{}\mathrm{\Lambda }(u+p\eta )`$ (3.8)
$``$ $`\delta (u\eta )\mathrm{\Lambda }(u+\eta )\mathrm{\Lambda }(u+2\eta )\mathrm{}\mathrm{\Lambda }(u+(p1)\eta )`$
$``$ $`\delta (u)\mathrm{\Lambda }(u+2\eta )\mathrm{\Lambda }(u+3\eta )\mathrm{}\mathrm{\Lambda }(u+p\eta )`$
$``$ $`\delta (u+\eta )\mathrm{\Lambda }(u)\mathrm{\Lambda }(u+3\eta )\mathrm{\Lambda }(u+4\eta )\mathrm{}\mathrm{\Lambda }(u+p\eta )`$
$``$ $`\delta (u+2\eta )\mathrm{\Lambda }(u)\mathrm{\Lambda }(u+\eta )\mathrm{\Lambda }(u+4\eta )\mathrm{}\mathrm{\Lambda }(u+p\eta )\mathrm{}`$
$``$ $`\delta (u+(p1)\eta )\mathrm{\Lambda }(u)\mathrm{\Lambda }(u+\eta )\mathrm{}\mathrm{\Lambda }(u+(p2)\eta )`$
$`+`$ $`\mathrm{}=f(u),`$
where $`\delta (u)`$ is defined by

$`\delta (u)={\displaystyle \frac{\mathrm{\Delta }(u)}{\zeta (2u+2\eta )}}.`$ (3.9)
Moreover, the function $`f(u)`$ is given by

$`f(u)`$ $`=`$ $`{\displaystyle \frac{(1)^{p(N+1)}}{2^{2p(N+1)}}}\mathrm{sinh}^{2N}((p+1)u){\displaystyle \frac{\mathrm{cosh}^2((p+1)u+\frac{i\pi }{2}ϵ)}{\mathrm{cosh}^2((p+1)u)}}`$ (3.10)
$`\times `$ $`\{n(u;\xi _{},\kappa _{})n(u;\xi _+,\kappa _+)+n(u;\xi _{},\kappa _{})n(u;\xi _+,\kappa _+)`$
$`+2(1)^N(\kappa _{}\kappa _+)^{p+1}\mathrm{sinh}^2(2(p+1)u)\mathrm{cosh}((p+1)(\theta _{}\theta _+))\},`$
where $`ϵ=2\text{frac}(p/2)`$ equals 0 if $`p`$ is even, and equals 1 if $`p`$ is odd; and the function $`n(u;\xi ,\kappa )`$ is defined by

$`n(u;\xi ,\kappa )=\mathrm{sinh}\left((p+1)(\xi +u)\right)+{\displaystyle \underset{l=1}{\overset{\left[\frac{p+1}{2}\right]}{}}}c_{p,l}\kappa ^{2l}\mathrm{sinh}\left((p+1)u+(p+12l)\xi \right),`$ (3.11)
with

$`c_{p,l}={\displaystyle \frac{(p+1)}{l!}}{\displaystyle \underset{k=0}{\overset{l2}{}}}(plk).`$
For instance, for the case $`p=3`$, the functional relation is given by

$`\mathrm{\Lambda }(u)\mathrm{\Lambda }(u+\eta )\mathrm{\Lambda }(u+2\eta )\mathrm{\Lambda }(u+3\eta )\delta (u\eta )\mathrm{\Lambda }(u+\eta )\mathrm{\Lambda }(u+2\eta )`$
$`\delta (u)\mathrm{\Lambda }(u+2\eta )\mathrm{\Lambda }(u+3\eta )\delta (u+\eta )\mathrm{\Lambda }(u)\mathrm{\Lambda }(u+3\eta )\delta (u+2\eta )\mathrm{\Lambda }(u)\mathrm{\Lambda }(u+\eta )`$
$`+\delta (u)\delta (u+2\eta )+\delta (u\eta )\delta (u+\eta )=f(u).`$ (3.12)
### 3.3 Step 3: determinant representation

Following the strategy used in  to solve RSOS models, the third and final step is to rewrite the functional relation as the determinant of a $`(p+1)\times (p+1)`$ matrix. Let us assume that this matrix has the same form as the one for the diagonal case ($`\kappa _\pm =0`$) and for the case (1.2). That is, we assume the functional relation can be cast in the form

$`det\left(\begin{array}{cccccccc}\mathrm{\Lambda }_0& h_1^{}& 0& 0& \mathrm{}& 0& 0& h_0\\ h_1& \mathrm{\Lambda }_1& h_0^{}& 0& \mathrm{}& 0& 0& 0\\ 0& h_2& \mathrm{\Lambda }_2& h_1^{}& \mathrm{}& 0& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& 0& \mathrm{}& h_{p1}& \mathrm{\Lambda }_{p1}& h_{p2}^{}\\ h_{p1}^{}& 0& 0& 0& \mathrm{}& 0& h_p& \mathrm{\Lambda }_p\end{array}\right)=0,`$ (3.19)
where $`\mathrm{\Lambda }_k=\mathrm{\Lambda }(u+\eta k)`$, $`h_k=h(u+\eta k)`$, $`h_k^{}=h^{}(u+\eta k)`$,

$`h^{}(u)=h(u2\eta ),`$ (3.20)
and the function $`h(u)`$ is yet to be determined. We find that the functional relation (3.8) can indeed be recast in the form (3.19), provided that $`h(u)`$ satisfies the three conditions

$`h(u+i\pi )`$ $`=`$ $`h(u)`$ (3.21)
$`h(u+\eta )h(u\eta )`$ $`=`$ $`\delta (u),`$ (3.22)
$`{\displaystyle \underset{j=0}{\overset{p}{}}}h(u+j\eta )+{\displaystyle \underset{j=0}{\overset{p}{}}}h(uj\eta )`$ $`=`$ $`f(u).`$ (3.23)
The results  for $`h(u)`$ in the diagonal case and in the case (1.2) suggest that, in general, $`h(u)`$ has form

$`h(u)=\mathrm{sinh}^{2N}(u+\eta ){\displaystyle \frac{\mathrm{sinh}(2u+2\eta )}{\mathrm{sinh}(2u+\eta )}}g_{}(u)g_+(u),`$ (3.24)
where the functions $`g_{}(u)`$ contain all the dependence on the boundary parameters.

Then second condition (3.22) together with (3.2), (3.9) and (3.24) imply that

$`g_{}(u)g_+(u)g_{}(u)g_+(u)`$ $`=`$ $`(\mathrm{sinh}^2u\mathrm{sinh}^2\xi _{}+\kappa _{}^2\mathrm{sinh}^22u)`$ (3.25)
$`\times `$ $`(\mathrm{sinh}^2u\mathrm{sinh}^2\xi _++\kappa _+^2\mathrm{sinh}^22u).`$
This suggests that $`g_{}(u)`$ obey the functional equation

$`g_{}(u)g_{}(u)=(\mathrm{sinh}^2u\mathrm{sinh}^2\xi _{}+\kappa _{}^2\mathrm{sinh}^22u).`$ (3.26)
Assuming that the functions $`g_{}(u)`$ are given by

$`g_{}(u)=2\kappa _{}\mathrm{sinh}(u+\alpha _{})\mathrm{cosh}(u+\beta _{}),`$ (3.27)
then (3.26) implies that the parameters $`\alpha _{}`$, $`\beta _{}`$ obey

$`\mathrm{sinh}^2\alpha _{}\mathrm{cosh}^2\beta _{}={\displaystyle \frac{1}{4\kappa _{}^2}}\mathrm{sinh}^2\xi _{},\mathrm{cosh}^2\alpha _{}\mathrm{sinh}^2\beta _{}={\displaystyle \frac{1}{4\kappa _{}^2}}\mathrm{cosh}^2\xi _{}.`$ (3.28)
A similar reparametrization appears in . Below we shall argue that (3.27) is essentially the unique solution of (3.25).

The third condition (3.23) together with (3.10) and (3.24) imply that 111Note that the right-hand-side of (3.29) depends on $`N`$ only through its parity $`(1)^N`$.

$`{\displaystyle \underset{j=0}{\overset{p}{}}}g_{}(u+j\eta )g_+(u+j\eta )+{\displaystyle \underset{j=0}{\overset{p}{}}}g_{}(uj\eta )g_+(uj\eta )`$ (3.29)
$`=`$ $`{\displaystyle \frac{(1)^p}{2^{2p}}}\{n(u;\xi _{},\kappa _{})n(u;\xi _+,\kappa _+)+n(u;\xi _{},\kappa _{})n(u;\xi _+,\kappa _+)`$
$`+2(1)^N(\kappa _{}\kappa _+)^{p+1}\mathrm{sinh}^2(2(p+1)u)\mathrm{cosh}((p+1)(\theta _{}\theta _+))\},`$
where $`n(u;\xi ,\kappa )`$ is given by Eq. (3.11). We find that this requirement can be satisfied for $`p=\text{odd}`$, with $`g_{}(u)`$ given by (3.27) and 222The requirement of including the special case (1.2), which corresponds to $`\alpha _{}=\alpha _+`$, $`\beta _{}=\beta _+`$, $`k=0`$, helps to resolve the sign ambiguity in passing from (3.28) to (3.30).

$`\mathrm{sinh}\alpha _{}\mathrm{cosh}\beta _{}`$ $`=`$ $`{\displaystyle \frac{1}{2\kappa _{}}}\mathrm{sinh}\xi _{},\mathrm{cosh}\alpha _{}\mathrm{sinh}\beta _{}={\displaystyle \frac{1}{2\kappa _{}}}\mathrm{cosh}\xi _{},`$
$`\mathrm{sinh}\alpha _+\mathrm{cosh}\beta _+`$ $`=`$ $`{\displaystyle \frac{1}{2\kappa _+}}\mathrm{sinh}\xi _+,\mathrm{cosh}\alpha _+\mathrm{sinh}\beta _+={\displaystyle \frac{1}{2\kappa _+}}\mathrm{cosh}\xi _+,`$ (3.30)
provided that the various parameters obey the linear constraint

$`\alpha _{}+\beta _{}+\alpha _++\beta _+=\pm (\theta _{}\theta _+)+\eta k,`$ (3.31)
where $`k`$ is an even integer if $`N`$ is odd, and is an odd integer if $`N`$ is even.

In short, the functional relations can be cast in the determinant form (3.19) for $`p=\text{odd}`$ with

$`h(u)`$ $`=`$ $`\mathrm{sinh}^{2N}(u+\eta ){\displaystyle \frac{\mathrm{sinh}(2u+2\eta )}{\mathrm{sinh}(2u+\eta )}}`$ (3.32)
$`\times `$ $`4\kappa _{}\kappa _+\mathrm{sinh}(u+\alpha _{})\mathrm{cosh}(u+\beta _{})\mathrm{sinh}(u+\alpha _+)\mathrm{cosh}(u+\beta _+),`$
where $`\alpha _{}`$, $`\beta _{}`$ are defined by (3.30) and satisfy the constraint (3.31).

We now proceed as in , and assume that the matrix in (3.19) has the null vector $`(Q_0,Q_1,\mathrm{},Q_p)`$. That is,

$`\mathrm{\Lambda }_0Q_0h_1^{}Q_1h_0Q_p`$ $`=`$ $`0,`$
$`h_kQ_{k1}+\mathrm{\Lambda }_kQ_kh_{k1}^{}Q_{k+1}`$ $`=`$ $`0,k=1,\mathrm{},p1,`$
$`h_{p1}^{}Q_0h_pQ_{p1}+\mathrm{\Lambda }_pQ_p`$ $`=`$ $`0.`$ (3.33)
We make the Ansatz $`Q_k=Q(u+\eta k)`$, where $`Q(u)`$ is given by

$`Q(u)={\displaystyle \underset{j=1}{\overset{M}{}}}\mathrm{sinh}(uu_j)\mathrm{sinh}(u+u_j+\eta ),`$ (3.34)
which has the crossing symmetry $`Q(u)=Q(u\eta )`$. The zeros $`u_j`$ of $`Q(u)`$ are still to be determined. Eqs. (3.33) and (3.20) imply that the eigenvalues are given by

$`\mathrm{\Lambda }(u)=h(u){\displaystyle \frac{Q(u\eta )}{Q(u)}}+h(u\eta ){\displaystyle \frac{Q(u+\eta )}{Q(u)}}.`$ (3.35)
We verify that this result is consistent with both the periodicity (2.18) and crossing (2.19) properties of the transfer matrix. The requirement that $`\mathrm{\Lambda }(u)`$ be analytic at $`u=u_j`$ yields the Bethe Ansatz equations

$`{\displaystyle \frac{h(u_j)}{h(u_j\eta )}}={\displaystyle \frac{Q(u_j+\eta )}{Q(u_j\eta )}},j=1,\mathrm{},M.`$ (3.36)
The asymptotic behavior (2.20), together with the result (3.35) for the eigenvalues and the constraint (3.31), imply that the number $`M`$ of Bethe roots is given by

$`M={\displaystyle \frac{1}{2}}(N1+k),`$ (3.37)
where $`k`$ is the integer appearing in (3.31). We leave to a future investigation the interesting question of determining further restrictions on the value of $`k`$, which presumably is related to the question of completeness.

We now argue that (3.27) is essentially the unique solution of (3.25). Indeed, if $`\stackrel{~}{g}_{}(u)`$ are also solutions of (3.25), then

$`\stackrel{~}{g}_{}(u)\stackrel{~}{g}_+(u)=g_{}(u)g_+(u)\varphi (u),`$ (3.38)
where $`g_{}(u)`$ are given by (3.27), and $`\varphi (u)`$ satisfies

$`\varphi (u)\varphi (u)=1.`$ (3.39)
The periodicity condition (3.21) implies that $`\varphi (u)`$ has the same periodicity

$`\varphi (u+i\pi )=\varphi (u).`$ (3.40)
We infer from (3.39) and (3.40) that $`\varphi (u)`$ is a CDD-like factor

$`\varphi (u)={\displaystyle \underset{j}{}}{\displaystyle \frac{\mathrm{sinh}(u+v_j)}{\mathrm{sinh}(uv_j)}}.`$ (3.41)
The requirement that $`\mathrm{\Lambda }(u)`$ be analytic then restricts $`\varphi (u)`$ to the form

$`\varphi (u)={\displaystyle \frac{q(u\eta )}{q(u)}},\text{where}q(u)={\displaystyle \underset{j}{}}\mathrm{sinh}(uv_j)\mathrm{sinh}(u+\eta +v_j).`$ (3.42)
This is equivalent to having additional Bethe roots, which can be included in $`Q(u)`$ (3.34).

Although the above results for the eigenvalues (3.32), (3.35), (3.36) have been obtained under the assumption that $`\eta `$ is restricted to the values (3.6) with $`p`$ odd, we expect that these results remain valid for generic values of $`\eta `$. Indeed, we have explicitly verified that these expressions reproduce the correct eigenvalues for $`N=0`$ (with $`M=0`$) and $`N=1`$ (with $`M=1`$) for arbitrary $`\eta `$.

## 4 Conclusion

Our proposed expression for the eigenvalues $`\mathrm{\Lambda }(u)`$ of the transfer matrix (2.1) corresponding to the Hamiltonian (1.3) is given by (3.35), where $`h(u)`$ and $`Q(u)`$ are given by (3.32) and (3.34), respectively; the Bethe Ansatz equations are given by (3.36), with $`M`$ given by (3.37); and the parameters $`\alpha _{}`$, $`\beta _{}`$ (which are related to $`\xi _{}`$, $`\kappa _{}`$ by (3.30)) must satisfy the constraint (3.31).

It remains an open question whether a solution with Bethe Ansatz equations of the “conventional” form (3.36) can be found which does not require a constraint among the boundary parameters. (Although the solution proposed in  does not require any constraint among the boundary parameters, it holds only for the $`\eta `$ values (3.6), and the Bethe Ansatz equations are not of the conventional form.)

For the special case (1.2), an analysis of the thermodynamic ($`N\mathrm{}`$) limit and an extension to higher-dimensional representations has recently been given in . For the more general case discussed here, it should now be possible to address such questions, and also to find generalizations to higher rank-algebras, both for the trigonometric and elliptic cases.

## Acknowledgments

I am grateful to all the organizers of the Annecy workshop “Recent Advances in the Theory of Quantum Integrable Systems” for the opportunity to present this work. I am also grateful to many of the participants, in particular A. Belavin, A. Doikou, B. McCoy, F. Ravanini, V. Rittenberg and S. Ruijsenaars, for their questions or comments. This work was supported in part by the National Science Foundation under Grant PHY-0098088.