text:
# HD-THEP-00-49 QCD sum rule analysis of the charmonium system: The charm quark mass

## 1 Introduction

An essential task within modern particle physics lies in the determination of quark masses which are key parameters of the standard model. In the past, QCD moment sum rule analyses have been successfully applied for extracting the charm and bottom quark masses from experimental data on the charmonium and bottonium systems respectively . The fundamental quantity in these investigations is the vacuum polarisation function $`\mathrm{\Pi }(q^2)`$:

$`\mathrm{\Pi }_{\mu \nu }(q^2)`$ $`=`$ $`i{\displaystyle d^4xe^{iqx}T\{j_\mu (x)j_\nu (0)\}}`$ (1)
$`=`$ $`(q_\mu q_\nu g_{\mu \nu }q^2)\mathrm{\Pi }(q^2),`$
where the current is represented by $`j_\mu (x)=(\overline{c}\gamma _\mu c)(x)`$. Via the optical theorem, the experimental cross section $`\sigma (e^+e^{}c\overline{c})`$ is related to the imaginary part of $`\mathrm{\Pi }(s)`$:

$$R(s)=\frac{1}{Q_c^2}\frac{\sigma (e^+e^{}c\overline{c})}{\sigma (e^+e^{}\mu ^+\mu ^{})}=12\pi \text{Im}\mathrm{\Pi }(s+iϵ).$$

Using a dispersion relation, we can express the moments by an integral over the velocity $`v=\sqrt{14m^2/s}`$:

$`_n={\displaystyle \frac{12\pi ^2}{n!}}\left(4m^2{\displaystyle \frac{d}{ds}}\right)^n\mathrm{\Pi }(s)|_{s=0}`$ (2)
$`=\left(4m^2\right)^n{\displaystyle _{s_{min}}^{\mathrm{}}}𝑑s{\displaystyle \frac{R(s)}{s^{n+1}}}=2{\displaystyle _0^1}𝑑vv(1v^2)^{n1}R(v).`$
The moments can either be calculated theoretically, including perturbation theory, Coulomb resummation and non-perturbative contributions, or can be obtained from experiment. In this way, we can relate the charm quark mass to the hadronic properties of the charmonium system.

A natural choice for the mass appearing in eq. (2) is the pole mass $`M`$. In the first part of our numerical analysis, we will use the pole mass scheme to extract the charm pole mass. However, as the pole mass suffers from renormalon ambiguities , it can only be determined up to corrections of order $`\mathrm{\Lambda }_{QCD}`$. In the second part of our analysis we shall therefore use the recently introduced potential-subtracted mass $`m_{PS}`$ . From this mass definition we can obtain the $`\overline{\mathrm{MS}}`$-mass more accurately than from the pole mass scheme.

In the next section, we shall present the perturbative expansion of the correlator which is known up to next-to-next-to-leading order (NNLO). As was found in recent analyses of the Upsilon-system , the dominant theoretical contributions arise from the threshold expansion in the framework of non-relativistic QCD (NRQCD). These contributions will be derived in section 3. Afterwards, we shall shortly discuss the non-perturbative contributions and the phenomenological spectral function. In the numerical analysis, we shall first explain the method of analysis. We will then obtain the pole mass and the $`\overline{\mathrm{MS}}`$-mass from an analysis in the pole mass and the PS-mass scheme respectively. The origin of different contributions to the error will be carefully investigated. We shall conclude with a summary and an outlook.

## 2 Perturbative expansion

The perturbative spectral function $`R(v)`$ can be expanded in powers of the strong coupling constant $`\alpha _s`$,

$$R^{Pt}(s)=R^{(0)}(s)+\frac{\alpha _s}{\pi }R^{(1)}(s)+\frac{\alpha _s^2}{\pi ^2}R^{(2)}(s)+\mathrm{}$$

From this expression the corresponding moments $`_n`$ can be calculated via the integral of eq. (2). The first two terms are known analytically and can for example be found in ref. . $`\mathrm{\Pi }^{(2)}(s)`$ is still not fully known analytically. However, the method of Padé-approximants has been exploited to calculate $`\mathrm{\Pi }^{(2)}`$ numerically, using available results at high energies, analytical results for the first eight moments and the known threshold behaviour . This information is sufficient to obtain a numerical approximation of $`\mathrm{\Pi }^{(2)}(s)`$ in the full energy range.

The numerical stability of the results can be checked in different ways. By choosing different Padé-approximants or by selecting a smaller set of input data the results for the moments remain almost unchanged. Furthermore, some contributions to the spectral density like those from internal quark loops are known analytically and are in very good agreement with the numerical spectral density. The above are strong indicators that the numerically obtained spectral density comes very close to the exact spectral density.

## 3 Coulomb resummation

The perturbative expansion which contains terms up to the order $`O(\alpha _s^2)`$ is a good approximation for high velocities. However, as one approaches the threshold region, terms of the order $`v(\alpha _s/v)^k`$ become increasingly important. These terms can be obtained in the framework of NRQCD. The correlator is then expressed in terms of a Greens function :

$$\mathrm{\Pi }(s)=\frac{N_c}{2M^2}\left(C_h(\alpha _s)G(k)+\frac{4k^2}{3M^2}G_C(k)\right),$$

where $`k=\sqrt{M^2s/4}`$ and $`M`$ represents the pole mass. The constant $`C_h(\alpha _s)`$ is a perturbative coefficient which is needed for the matching between the full and the non-relativistic theory. It naturally depends on the hard scale. The Greens function is analytically known up to NNLO  and sums up terms of order $`\alpha _s^n/v^{nk}`$ for $`n0`$ and $`k=1,2,3`$. It is crucial for the analysis that the result depends on three scales. While the hard scale $`\mu _{hard}`$ is responsible for the hard perturbative processes, the soft scale $`\mu _{soft}`$ governs the expansion of the Greens function. Furthermore, the factorisation scale $`\mu _{fac}`$ separates the contributions of large and small momenta and plays the role of an infrared cutoff.

The Greens function contains two parts: the continuum and the poles above and below threshold respectively. We are interested in both contributions separately: first, the individual corrections can be analysed and their error estimated. Second, in our numerical analysis we will reconstruct the spectral density above threshold and we thus need the corresponding spectral density at low velocities. In principle, the expressions for the energies and decay widths of the poles have been calculated. However, in the actual case of the charm quark this expansion does not converge well. We will therefore choose a different method of evaluation. Since the results for the Greens function are known analytically, we can evaluate their contribution to the moments numerically by performing the derivatives. On the other hand, by using a dispersion relation, we can obtain the continuum from the imaginary part of the correlator. From the difference we can obtain the pole contributions:

$`_n^{Poles}`$ $`=`$ $`{\displaystyle \frac{12\pi ^2}{n!}}\left(4M^2{\displaystyle \frac{d}{ds}}\right)^n\mathrm{\Pi }(s)|_{s=0}`$
$`12\pi \left(4M^2\right)^n{\displaystyle _{4M^2}^{\mathrm{}}}𝑑s{\displaystyle \frac{\text{Im}\mathrm{\Pi }(s)}{s^{n+1}}}.`$
Since we will not evaluate the poles near threshold, but rather calculate their contributions to the moments in a region where perturbation theory is expected to be valid, the convergence of the pole contributions is improved. Nevertheless, the poles will give the largest contribution to the theoretical moments and thus the dependence on the scales will remain relatively strong. The large corrections are partly due to the definition of the pole mass. These contributions can be reduced by using an intermediate mass definition. In this analysis we will use the potential-subtracted (PS) mass  where the potential below a separation scale $`\mu _{sep}`$ is subtracted:

$`\delta m(\mu _{sep})`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{|𝐪|<\mu _{sep}}{}}{\displaystyle \frac{d^3q}{(2\pi )^3}}V(q),`$
$`m_{PS}(\mu _{sep})`$ $`=`$ $`M\delta m(\mu _{sep}).`$
As will be seen in the numerical analysis, this mass definition leads to an improved scale dependence and a more precise determination of the $`\overline{\mathrm{MS}}`$-mass.

## 4 Condensate contributions

The non-perturbative effects on the vacuum correlator are parametrised by the condensates. The leading correction is the gluon condensate contribution which is known up to next-to-leading order . Furthermore, the dimension 6 and 8 contributions have been calculated and will be included in our analysis . It will turn out that the condensate contributions are suppressed when compared to former charmonium sum rule analyses and only have little influence on the mass. Besides an increase of the theoretical moments from the Coulomb contributions we will restrict the moments to $`n7`$ where the non-perturbative contributions are relatively small. Since we obtain a larger pole mass than the former analyses, the condensates, starting with a power of $`1/M^4`$, are suppressed further.

## 5 Phenomenological spectral function

Experimentally, the first six $`\psi `$-resonances have been observed. Since the widths of the poles are very small compared to the masses, the narrow-width approximation provides an excellent description of these states. To model the contributions above the 6th resonance, we use the assumption of quark-hadron-duality and integrate over the perturbative spectral density:

$$\frac{_n}{(4M^2)^n}=\frac{9\pi }{\alpha _{e.m.}^2Q_c^2}\underset{k=1}{\overset{6}{}}\frac{\mathrm{\Gamma }_k}{E_k^{2n+1}}+_{s_0}^{\mathrm{}}𝑑s\frac{R^{Pt}(s)}{s^{n+1}}.$$

To estimate the continuum contribution we will use a threshold in the range of $`3.6`$ GeV $`\sqrt{s_0}4.2`$ GeV with a central value of $`\sqrt{s_0}=3.8`$ GeV. However, the most dominant phenomenological contributions come from the first two $`\psi `$-resonances resulting in a small influence of the continuum even for low values of $`n`$.

## 6 Numerical analysis

### 6.1 Analysing method

Besides the contributions from the poles of the Greens function and the condensates, the theoretical part of the correlator contains the spectral density above threshold. Now we will discuss the different parts of the spectral density. For high velocities the spectral density is well described by the perturbative expansion. However, as one approaches smaller values of $`v`$, the perturbative expansion breaks down. The resummed spectral density, on the other hand, gives a good description for low values of $`v`$, but becomes unreliable for high velocities. For these reasons, we will introduce a separation velocity $`v_{sep}0.3`$. Above $`v_{sep}`$ we will use the perturbative spectral density. Below $`v_{sep}`$, we take the resummed spectral density adding the terms which are included in perturbation theory but not in resummation.

In fig. 1 we have displayed the different contributions. The dotted line represents the perturbative expansion in LO and NLO. The thin solid line also includes the NNLO. Whereas for high velocities the perturbative expansion is well convergent, the importance of the higher corrections increases for smaller $`v`$. The dashed line is the resummed spectral density and the thick solid line the reconstructed spectral density. For the charmonium system, there exists a range of intermediate values of $`v`$ where neither the perturbative expansion nor the resummation can be trusted. Indeed, it can be clearly seen that the reconstructed spectral density shows a gap at the separation velocity. To estimate the error we have varied $`v_{sep}`$ between 0.2 and 0.4. The analysis shows that though the introduction of the separation velocity stabilises the sum rules, the variation only has a minor influence on the mass.

### 6.2 Pole mass scheme

Since the perturbative corrections grow for large values of $`n`$, we will restrict our analysis to moments with $`n7`$. As our analysing method needs moments of $`n3`$ to reconstruct the spectral density, we will use $`3n7`$. As the central values for our scales we will choose

$`\mu _{soft}`$ $`=`$ $`1.2\text{GeV},\mu _{fac}=1.45\text{GeV},`$
$`\mu _{hard}`$ $`=`$ $`1.7\text{GeV}.`$ (3)
We have set the hard scale equal to the central value for the pole mass. For the soft scale, we would have liked to choose a somewhat smaller value, but the NNLO corrections get out of control for $`\mu _{soft}1`$ GeV. The factorisation scale lies between the two other scales. In this scheme the theoretical moments are dominated by the pole contributions. To estimate the error on the mass, we have varied the scales within reasonable ranges. The result is depicted in fig. 2.

The thick solid line gives the central value for the charm mass of $`M_c=1.70`$ GeV with the values of eq. (6.2). The error is dominated by the variation of the scales, we obtain

$`1.05\text{GeV}\mu _{soft}1.5\text{GeV}:`$ $`\mathrm{\Delta }M_c=100\text{MeV}`$
$`1.2\text{GeV}\mu _{fac}1.7\text{GeV}:`$ $`\mathrm{\Delta }M_c=50\text{MeV}`$
$`1.4\text{GeV}\mu _{hard}2.5\text{GeV}:`$ $`\mathrm{\Delta }M_c=40\text{MeV}.`$
A significant uncertainty also comes from $`\mathrm{\Lambda }_{QCD}`$.

By choosing $`\mathrm{\Lambda }_{QCD}=330\pm 30`$ MeV we get an error of $`\mathrm{\Delta }M_c=50`$ MeV. The results are summarised in table 1. Adding the errors in quadrature we obtain the charm pole mass:

$$M_c=1.70\pm 0.13\text{GeV}.$$
(4)
Instead of deriving the $`\overline{\mathrm{MS}}`$-mass from the pole mass, in the following section we will use the PS-mass to obtain a more precise value for the $`\overline{\mathrm{MS}}`$-mass.

### 6.3 Potential-subtracted mass scheme

First, we have to choose a value for the separation scale $`\mu _{sep}`$. This scale should be taken large enough to guarantee a perturbative relation to the $`\overline{\mathrm{MS}}`$-mass. On the other hand, it should be smaller than $`Mv`$. Both conditions cannot be well fulfilled at the same time. As a compromise value we will choose $`\mu _{sep}=1.0\pm 0.2`$ GeV. In this scheme the pole contributions from the Greens function turn out to be smaller than in the pole mass scheme. The contributions from the condensates get more important here and we shall restrict our analysis to $`n6`$ where these corrections are under good control.

Using the same central values for the scales (6.2) we obtain $`m_{PS}(\mu _{sep}=1.0)=1.35`$ GeV and from this value a $`\overline{\mathrm{MS}}`$-mass of $`m_c(m_c)=1.23`$ GeV. The introduction of the intermediate mass definition leads to a reduced scale dependence:

$`1.05\text{GeV}\mu _{soft}1.5\text{GeV}:`$ $`\mathrm{\Delta }m_{PS}=60\text{MeV}`$
$`1.2\text{GeV}\mu _{fac}1.7\text{GeV}:`$ $`\mathrm{\Delta }m_{PS}=40\text{MeV}`$
$`1.4\text{GeV}\mu _{hard}2.5\text{GeV}:`$ $`\mathrm{\Delta }m_{PS}=30\text{MeV}.`$
In table 2 we have listed the individual contributions to the error of $`m_{PS}`$ and $`m_c`$. Finally, we obtain for the masses

$`m_{PS}(\mu _{sep}=1.0)`$ $`=`$ $`1.35\pm 0.09\text{GeV}`$
$`m_c(m_c)`$ $`=`$ $`1.23\pm 0.09\text{GeV}.`$ (5)
When we use the pole mass from eq. (4) and calculate the $`\overline{\mathrm{MS}}`$-mass using the three-loop relation between these masses , we obtain $`m_c(m_c)=1.20\pm 0.17`$ GeV. This value is in good agreement with eq. (6.3). This is not self-evident as the dominating pole contributions are reduced in the PS-scheme and the relative influence of the individual contributions is shifted.

## 7 Conclusions

The obtained value for the pole mass lies somewhat higher than in former sum rule analyses . In  the authors used perturbation theory to NLO resulting in a value of $`M_c=1.46\pm 0.07`$ GeV. In the second investigation  the analysis has been performed in the $`\overline{\mathrm{MS}}`$-scheme with perturbation theory to NLO. Using the NLO relation to the pole mass the author obtains $`m_c(m_c)=1.26\pm 0.05`$ GeV and $`M_c=1.42\pm 0.03`$ GeV. The author has also performed an analysis using resummation in LO with a value of $`M_c=1.45\pm 0.07`$ GeV. In our analysis the increased value of the pole mass is essentially due to large Coulomb contributions which have not been included in former analyses. As a consequence, the error becomes larger as well. In a recent analysis, the pole mass has been estimated from the charmonium ground state at NNLO . Here the authors obtained a pole mass of $`M_c=1.88_{0.13}^{+0.22}`$ GeV. This value resulted from large corrections of the Coulomb potential to the ground state. During the last years, several lattice analyses obtained the following values for the $`\overline{\mathrm{MS}}`$-mass:

$`m_c(m_c)`$ $`=`$ $`1.59\pm 0.28\text{GeV}\text{[22]},`$
$`m_c(m_c)`$ $`=`$ $`1.33\pm 0.08\text{GeV}\text{[23]},`$
$`m_c(m_c)`$ $`=`$ $`1.73\pm 0.26\text{GeV}\text{[24]},`$
$`m_c(m_c)`$ $`=`$ $`1.22\pm 0.05\text{GeV}\text{[25]}.`$
Whereas the results from  are in good agreement with this analysis, the investigations from  and  obtain higher masses. For the time-being, the results are not conclusive and future lattice calculations for the charm mass might be of interest. Since other methods reveal significant uncertainties in the determination of the quark masses, the sum rules remain one of the most precise methods to extract these fundamental quantities.

Acknowledgements M.Eidemüller would like to thank S. Narison for the invitation to this very “charm”ing and interesting conference.


format: text

text:
# Temperature of the Central Stars of Planetary Nebulae and the Effect of the Nebular Optical Depth

## 1 Introduction

The temperature of the central stars of planetary nebulae (PNs) is an essential parameter for evolutionary studies as well as for an analysis of the nebulae themselves. The different methods for the determination of the stellar temperature and the problems involved were extensively discussed by Pottasch (1984) and Kaler (1985a, 1989). The most common method was suggested by Zanstra (1931) and further developed by Harman & Seaton (1966). The temperature of the ionizing star of a planetary nebula is calculated from the ratio between the flux of a recombination line and the stellar continuum flux at a given frequency. The Zanstra method assumes that all photons above the H (or the He+) Lyman limit are absorbed within the nebula and that each recombination eventually results in a Balmer photon. The total ionizing flux can then be related to the total flux of a recombination line. The Zanstra method yields then two values for the stellar temperature: the Zanstra temperature obtained from the intensity of a hydrogen recombination line, TZ(H), and the Zanstra temperature obtained from a He II recombination line, TZ(He II). When applied to observed nebulae, TZ(H) is generally lower than TZ(He II); the difference can reach values of the order of 60,000 K (Kaler 1983b), and the temperature ratio can reach a factor higher than 3 (Kaler 1983b, 1985a). This is the well known Zanstra discrepancy, which is stronger for PNs with low stellar temperatures. In fact, many PNs with TZ(H) $`<`$ 100,000 K have TZ(He II) $`>`$ TZ(H) while, for higher stellar temperatures, both Zanstra temperatures are similar (Pottasch, 1984; Gathier & Pottasch 1988, 1989). Another important point is that the calculated Zanstra temperatures do not reproduce the high temperatures predicted by evolutionary models (Kaler 1985a, 1989; Stasińska & Tylenda 1986).

Another interesting feature that could be related to the Zanstra temperature issues is the distribution of objects in a log L - log T plane. As noted by Shaw & Kaler (1989), when this distribution is based on Zanstra temperatures there is a dense crowd of planetary nebula nuclei with temperatures of $``$ 100,000 K and a strong decrease towards higher temperatures, what these authors call the “Zanstra wall”.

The causes of the Zanstra discrepancy have been discussed by several authors (Pottasch 1984; Kaler 1985a, 1989; Henry & Shipman 1986;Stasińska & Tylenda 1986; Kudritzki & Méndez 1989; Gabler, Kudritzki & Méndez 1991; Méndez, Kudritzki & Herrero 1992) and can be related to the following effects: (1) optical effects, i.e., nebulae exhibiting a Zanstra discrepancy would be optically thin to photons ionizing H yet optically thick to those ionizing He+; (2) differential dust absorption in the nebula; (3) the stellar continuum differing from the usually assumed blackbody spectrum; for example, an excess of photons with energies beyond the He+ ionization potential would result in a high He II Zanstra temperature.

Objects with a large discrepancy show fainter low-ionization lines, suggesting that the effect 1 is the correct interpretation (Kaler 1983b). Simple evolutionary models of planetary nebulae predict TZ(H) $``$ TZ(He II) during the nebula lifetime because of the variation of their optical depth (Tylenda et al. 1994; see also Schönberner & Tylenda 1990). The effect of dust on the calculated Zanstra temperature is discussed by Helfer et al. (1981). For an assumed absorption law (with a high opacity around 50 eV), Zanstra temperatures underestimate the actual stellar temperature. The effect is more intense for TZ(H) and leads to the Zanstra discrepancy. However, dust appears to be only important in some specific nebulae (Kaler 1985a). Departures from the blackbody, in particular an excess of photons beyond the He+ ionization potential, are suggested by observational studies of many central stars of planetary nebulae (Kaler 1985a; see also references in Henry & Shipman 1986). An ionizing spectrum with an excess of high energy photons (relative to a blackbody), produced by a star with a less than solar atmospheric He abundance, leads to a He II Zanstra temperature higher than TZ(H) (Henry & Shipman 1986).

Detailed theoretical analyses, using photoionization models and applied to optically thick nebulae, are presented by Stasińska & Tylenda (1986) and Henry & Shipman (1986). Stasińska & Tylenda (1986) show that for low stellar temperatures (T $``$ 100,000 K), both Zanstra temperatures are similar. For higher stellar temperatures, TZ(HI) is larger and TZ(He II) is lower than T. However, this result is opposite to what is obtained from observations. Analyzing only models with T $``$ 150,000 K, Henry and Shipman (1986) conclude that TZ(H) is a good measure of the stellar temperature. In fact, for these values of T, TZ(H) and TZ(He II) are similar (Stasińska & Tylenda 1986) and provide a good stellar temperature estimation if the nebula is completely optically thick.

Other methods for determining the stellar temperature include modelling of stellar absorption line profiles, ionic ratios, fitting of model atmospheres, the energy balance or Stoy method, stellar UV energy distribution, etc. The results from these different methods were compared with the values given by the Zanstra method in order to explain the Zanstra discrepancy; however, these methods give discordant results and have many uncertainties (see, for instance, Kaler 1985a, Stasińska & Tylenda 1986, and Kaler 1989).

Many authors have adopted TZ(He II) as representative of the stellar temperature assuming that the Zanstra discrepancy is due to an optical depth effect (Kaler 1983b; Gleizes, Acker & Stenholm 1989; Kaler, Shaw & Kwitter 1990; Kaler & Jacoby 1991; Stanghellini, Corradi & Schwarz 1993). However, TZ(He II) may not be a good indicator of the stellar temperature since it never reaches values as high as predicted by theoretical stellar evolutionary studies.

Usually the discussion of the Zanstra discrepancy in terms of an optical depth effect is based on optically thin or thick nebulae at 13.6 eV or 54.4 eV. It is also implicitly assumed that the nebula is optically thick at 54.4 eV at a distance to the central star smaller than that corresponding to 13.6 eV. Harman & Seaton (1966) suggest the following criteria for the complete absorption of H0, He0, and He+ ionizing-photons: presence of \[OI\] lines, He I images smaller than H I images, and He++ fractional abundance $``$ 0.75, respectively. For Pottasch (1984) the Zanstra method must work for $`\tau _{13.6}`$ $`>`$ 1. Some authors define a criterion to distinguish between optically thin and thick objects using the ratio of the Zanstra temperatures, ZR = TZ(He II)/TZ(H). For example, for Shaw & Kaler (1985) the nebula is optically thick to the H Lyman continuum when ZR   $`\stackrel{<}{}`$  1.2, while for ZR   $`\stackrel{>}{}`$  2.5 and He II $`\lambda `$4686/H$`\beta `$   $`\stackrel{>}{}`$  0.9 it is thin for He+ Lyman continuum photons, TZ(He II) being a lower limit for the stellar temperature. The criterion used by Kaler & Jacoby (1989), based on line intensities of low-ionization lines, states that a nebula is thick when \[O II\] $`\lambda `$3727/H$`\beta `$ $``$ 1 and \[N II\] $`\lambda `$6584/H$`\alpha `$ $``$ 1. However, planetary nebulae can present a large range of optical depths, depending on the quantity of matter. Furthermore, even considering the central stellar radiation as a blackbody, different stellar temperatures correspond to different ratios between the number of ionizing photons with energy higher than 54.4 eV and those higher than 13.6 eV. Thus, the radial ionic distribution for H and He varies with T, changing the relative sizes of the H+ and He++ zones with the nebula optical depth at 13.6 eV and 54.4 eV. As remarked by Stasińska & Tylenda (1986), the radiation transfer is much more complicated than assumed by the Zanstra method.

In brief, the Zanstra temperature is commonly used in the literature for planetary nebula modelizations as well as evolutionary analysis, and low optical depth must be at least a partial explanation for the Zanstra discrepancy. However, a more detailed analysis is required in order to explain the issues listed above. In this paper, a careful analysis of the effect of the nebula optical depth on the determination of the Zanstra temperatures is intended. The effects on the Zanstra temperatures due to deviations of a blackbody spectrum and due to an overabundance of He in the nebula are also discussed. The theoretical models used in our analysis are described in §2. The results for the Zanstra temperature ratio (ZR) and its behavior with the stellar temperature and with the nebula optical depth appear in §3, which also includes a comparison with values derived from PNs observations. An alternative method to estimate the temperature of central stars of planetary nebulae is suggested and discussed in §4. The conclusions are outlined in §5.

## 2 Planetary nebulae models and theoretical values for the Zanstra temperatures

Models for typical planetary nebulae are generated with the photoionization code AANGABA (Gruenwald & Viegas 1992). The physical conditions of the gas are determined by solving the coupled equations of ionization and thermal balance for a spherical symmetric cloud. Several processes of ionization and recombination, as well as of gas heating and cooling, are taken into account. The transfer of the primary and diffuse radiation fields is treated in the “outward-only” approximation. For the radiation-bounded models (equivalent to completely optically thick nebulae), the calculations stop when the fractional abundance H+/H reaches 10-4, defining the maximum radius for the ionized nebula, Rmax. Matter-bounded models (with a nebula radius less than Rmax) will also be discussed. The input parameters are the ionizing radiation spectrum, the gas density, and the chemical abundance for the elements included in the calculations (H, He, C, N, O, Ne, Mg, Si, S, Cl, Ar, and Fe). A range of input parameters, typical of planetary nebulae (Pottasch 1984), is assumed: T $``$ 50,000K, L = 30 - 20,000 L and nH = 102 \- 106 cm-3. In order to discuss the Zanstra temperature for very hot stars, which are predicted by evolutionary models, a maximum stellar temperature of 500,000 K is adopted. A blackbody spectrum is assumed for the ionizing radiation, but the effects due to departures from this kind of spectrum will also be discussed. Concerning the chemical abundances, average values for planetary nebulae, as given by Kingsburgh & Barlow (1994), are assumed. For elements not given by these authors, the solar value is adopted (Grevesse & Anders 1989).

The He II $`\lambda `$4686 and H$`\beta `$ fluxes obtained for the theoretical nebulae are used to derive TZ(H) and TZ(He II) by the standard Zanstra method. For each set of input parameters, TZ(H) and TZ(He II) are calculated for different values of the nebula optical depth at the H Lyman limit.

## 3 Theoretical versus “observed” Zanstra temperatures

In the following section, the assumed energy distribution of the ionizing radiation is fixed (blackbody). Thus, any discrepancy between the temperature adopted for the central star, T, and the derived Zanstra temperatures is not due to the assumed spectrum but inherent to the method.

Our results show that TZ(H) reproduces fairly well the stellar temperature for optically thick nebulae ionized by a star with T $`<`$ 150,000 K, in agreement with Henry & Shipman (1986) and Stasińska & Tylenda (1986). For higher stellar temperatures, TZ(H) is greater than T, and the difference between these two values increases with T. These results agree with those of Stasińska & Tylenda (1986). We find, however, that the deviation of TZ(H) relative to T is slightly smaller. As already discussed by Stasińska & Tylenda (1986) the deviation of Zanstra temperatures from the stellar temperature is due to the fact that each He++ recombination gives more than one photon that ionizes H, and the proportion of photons ionizing He+ increases with the stellar temperature. Furthermore, a fraction of the high energy photons are in fact absorbed by H and not by He+. A detailed discussion on the generation of H ionizing photons following the He+ and He++ recombination can be seen in Osterbrock (1989).

### 3.1 Emitting volumes of H+, He+, and He++

Before discussing the influence of the nebula optical depth on the derived Zanstra temperatures, it is useful to illustrate how the relative sizes of the H+, He+, He++ Strömgren spheres change as a function of the stellar temperature. The variation of the fractional abundances of H and He ions with the position in the nebula is shown in terms of r/Rmax in Figure 1a (left panels) where r is the distance from the center of the nebula and Rmax is the maximum dimension of the ionized region (see §2). The results given in Figure 1 correspond to models with L = 3000 $`L_{}`$ and nH = 104 cm-3. For low stellar temperatures the He++ Strömgren radius, R$`_{He^{++}}`$, is much smaller than R$`_{He^+}`$ or R$`_{H^+}`$, as expected. However, as T increases, R$`_{He^{++}}`$ approaches R$`_{H^+}`$.

For a matter-bounded nebula (with a total extent less than Rmax) the emitting zones of H+, He+, and He++ can be smaller than their corresponding Strömgren spheres. In this case, the emitted line intensities will be lower than those emitted by a radiation-bounded nebula. For a given reduction of the nebula extent, the size of the emitting zones of each of these ions will be differently affected, depending on the temperature of the central star. For nebulae with low stellar temperatures, a reduction of the nebula size affects mainly the H0, H+, and He0 zones. Thus, the smaller the nebula radius, the lower TZ(H), while TZ(He II) may still be a good indicator of the stellar temperature. For increasing stellar temperatures, the volumes of the H+ and He++ zones tend to be equal (Fig. 1a). In this case a reduction of the nebula size can result in a matter-bounded nebula where the H+ and He++ zones are almost equally affected. In this case, both TZ(H) and TZ(He II) underestimate the stellar temperature.

### 3.2 The effect of the nebula optical depth

The possible underestimate of the stellar temperature, due to the fact that a nebula may not be radiation-bounded, can be discussed as an optical depth effect. If the nebula has not enough material to be radiation-bounded, its radius is smaller than Rmax, and the nebula optical depth at a given frequency will be lower than the optical depth of a radiation-bounded nebula. The reduction of the nebula radius (creating a matter-bounded nebula) will differently affect the nebula optical depth of the H0, He0, and He+ continua. The fractional ionic distribution of H and He ions with the optical depth at the H Lyman limit ($`\tau _{13.6}`$) is shown in Figure 1b. As seen in §3.1, for a radiation-bounded nebula with low T the He++ Strömgren radius is much smaller than the H+ Strömgren radius. Thus matter-bounded nebulae can have an optical depth at the H Lyman continuum, $`\tau _{13.6}`$, close to unity, while the optical depth at the He++ Lyman limit, $`\tau _{54.4}`$, is much higher. The object is then optically thin to the H-ionizing photons and optically thick to the He+-ionizing photons, leading to ZR higher than 1. In this case, TZ(He II) provides a better estimate of the stellar temperature. As the stellar temperature increases, the optical depth at the H and He+ Lyman limits tend to have similar values; both will be reduced if the radius of the nebula is smaller than that of a radiation-bounded nebula. In this case, neither TZ(H) nor TZ is a good indicator of the stellar temperature.

The effect of the nebula optical depth on the derived Zanstra temperatures is shown in a ZR versus TZ(He II) plot (Figs. 2a and 2b) for the same models as in Figure 1. Each solid line corresponds to models with a given stellar temperature; the nebula optical depth decreases with increasing ZR. The curves are labeled by the stellar temperature in units of 1000 K. The dashed curves connect the results of completely optically thick models (radiation-bounded nebulae) with different stellar temperatures; these results correspond to the minimum ZR value for a given stellar temperature.

Recalling the ionic distribution shown in Figures 1a and 1b, the behavior of ZR shown by the curves in Figure 2 can be easily understood: (1) For T $`<`$ 150,000 K, the He++ zone is inside the H+ zone and much smaller. Matter-bounded models with decreasing $`\tau _{13.6}`$ would result in weaker H$`\beta `$ emission line, while the He II $`\lambda `$4686 line is unchanged. Thus, starting at the minimum value, corresponding to the optically thick model, ZR increases while TZ(He II) is practically constant. When the optical depth is low enough to affect the He++ zone, ZR still increases but TZ(He II) decreases and the curves turn to the left; (2) For higher stellar temperatures, the decrease of TZ(He II) with $`\tau _{13.6}`$ happens closer to the optically thick value (the volumes of the H+ and He++ zones are similar) and ZR increases slowly. In each of the solid lines in Figure 2a, the points corresponding to $`\tau _{13.6}`$ = 1, $`\tau _{13.6}`$ = 10, and $`\tau _{54.4}`$ = 1 are indicated, respectively, by crosses, triangles, and dots. Notice that for T $`>`$ 200,000 K, a nebula can be optically thin for He+-ionizing photons, even for ZR $``$ 1, though thick for photons above the H Lyman limit.

The theoretical results also show that the Zanstra method tends to underestimate the stellar temperature. The effect is larger for higher stellar temperatures, even for high optical depths. This may explain the ”Zanstra wall” in the log L - log T plot, since high-temperature stars, predicted by stellar evolutionary models, are penalized by the Zanstra method.

### 3.3 Confronting the theoretical results with the observations

Values for ZR and TZ(He II) derived from observations for a large sample of PNs are plotted in Figure 2b in order to be compared to the theoretical results. For each object several values of the Zanstra temperatures can be found in the literature. The criteria used to select the objects and the values of the Zanstra temperatures plotted in Figure 2b are the following: (1) if the Zanstra temperatures coming from different authors are similar (difference less than 20$`\%`$ from the average value), their average value is taken; (2) if the same author presents discordant data for the same object, the more recent value is taken; (3) if all the data for a given object are discordant, the object is not included. The values for both Zanstra temperatures were taken from: Martin (1981); Kaler (1983b); Pottasch (1984); Reay et al. (1984); Shaw & Kaler 1985; Viadana & de Freitas Pacheco (1985); de Freitas Pacheco, Codina & Viadana (1986); Gathier & Pottasch (1988, 1989); Gleizes at al. (1989); Jacoby & Kaler (1989); Shaw & Kaler (1989); Kaler et al. (1990); Kaler & Jacoby (1991); Méndez et al. (1992). Since we are discussing the standard Zanstra method, Zanstra temperatures corrected by the Stasińska-Tylenda effect (1986) were not included.

Most of the observational points for ZR and TZ(He II) are inside the region defined by the theoretical curves that correspond to ionizing stars with a blackbody spectrum. Our results naturally explain the trend shown by the observational values: for lower H Zanstra temperatures ($``$ 100,000 K) many planetaries may have TZ(He II) $`>`$ TZ(H), i.e., ZR $`>`$ 1, while for higher temperatures the difference between these temperatures is smaller. Such a behavior, referred to as “strange” by Pottasch (1984), induced Gathier & Pottasch (1988) to discard the optical depth explanation, since nebulae with higher stellar temperature should be older and optically thinner. The distribution of ZR versus TZ(H) presents a similar trend. Notice that a decreasing ZR ratio with increasing TZ(H) was obtained by Gathier & Pottasch (1988) with a sample including fewer objects.

The variation of the stellar temperatures with the optical depth in Figure 2 can solve some problems raised in the literature. One such problem is the temperature of the ionizing star of NGC 1360. The Zanstra temperatures for this object (34,900 K and 79,300 K, from the references given above) are much smaller than the temperature obtained from UV measurements (100,000 K; Pottasch et al. 1978). In Figure 2b the position of this object is in the region where the lines are crowded; the nebula is thus optically thin and the star can have a higher temperature than that given by TZ(He II), as suggested by the UV data. Also, the discussion (Kaler & Hartkopf 1981) regarding A43 (a thin and high-excitation nebula with a low He II Zanstra temperature star) and A50 (medium excitation, thick and high TZ) must be reviewed, since, following our results (Figure 2), the central star of A43 can have a temperature much higher than TZ(He II).

### 3.4 Other effects

In the previous section, using the results from photoionization models corresponding to a given value of the stellar luminosity and gas density and assuming a blackbody spectrum for the ionizing radiation, it was shown that the main issues concerning the Zanstra temperatures can be explained by an optical depth effect. In the following discussion the results for different values for the stellar luminosity and/or the gas density, as well as for an ionizing radiation spectrum deviating from a blackbody shape, are presented.

First, still adopting a blackbody spectrum, we discuss the results corresponding to the whole range of adopted values for the stellar luminosity and gas density. We verified that, as long as T $`<`$ 200,000 K, the behavior of ZR with TZ(He II) is the same as discussed in §§3.1 and 3.2, For T $``$ 200,000 K and $`\tau _{13.6}`$ $`>`$ 1, models with low stellar luminosities ($``$ 100 L) and low gas density ($`<`$ 103 cm-3) can give Zanstra temperatures lower than those obtained with the standard models discussed above. The differences between the Zanstra and stellar temperatures increase with increasing T and decreasing values for L and nH. For example, for T = 300,000 K, a maximum difference occurs for $`\tau _{13.6}`$ $``$ 10, L = 10 L, and nH= 100 cm-3, when both Zanstra temperatures decrease by $``$ 20 $`\%`$, increasing the difference between Zanstra and effective stellar temperatures. However, only a few PNs would have such low stellar luminosities and gas densities.

A number of authors explain the Zanstra discrepancy by an excess of photons with energy above 54.4 eV in the ionizing spectrum. This could explain the high values of TZ(He II) compared to TZ(H). As discussed by Henry $`\&`$ Shipman (1986), observations and models imply an excess of photons beyond the He+ threshold in numerous planetary nebula nuclei. Such an excess could be produced by a stellar atmosphere with subsolar He abundances and would lead to TZ(He II) higher than TZ(H) when compared with models where a blackbody is assumed.

To show the effect of a spectrum presenting an excess of high-energy photons above 54.4 eV, we discuss the results of photoionization models with an ionizing radiation spectrum of a pure H atmosphere (Wesemael et al. 1980). For example, for a completely optically thick nebula around a 150,000 K star, TZ(He II) is 6$`\%`$ higher and TZ(H) is 13$`\%`$ lower compared to the corresponding blackbody results. For decreasing optical depths, the curves tend rapidly to those corresponding to blackbody models. In brief, only for completely optically thick nebulae an excess of high energy photons will affect (in a small amount) the calculated Zanstra temperatures.

### 3.5 Outsiders

Some objects shown in Figure 2b are outside the area covered by the models. These nebulae can be either above the area limited by the curves or below it. Those above the curves limiting the high values of ZR could be explained by an error in TZ(He II) of the order of 5 - 10 $`\%`$. However, checking more carefully, it can be verified that most of these nebulae are Abell nebulae, including NGC 246, the prototype of the class (Abell 1966). Observations of the central star of some of these nebulae show characteristics of high stellar temperatures; furthermore, the nebulae have low surface brightness and large angular diameter (Abell 1966). These objects are probably in an advanced evolutionary stage. Some of them are known to have very high nebular He abundance in their inner regions (Jacoby & Ford 1983; Guerrero & Manchado 1996). Calculations for He-rich nebulae show that abundances up to He/H = 0.25 can explain the positions of the nebulae lying above the limiting curves in Figure 2b. Results for T = 150,000 K and He/H = 0.20 are shown in Figure 3 by the dot-dashed line. Note that for high optical depths the curves for the same stellar temperature but different He abundance are superposed. For T = 150,000 K, the curves separate from each other when $`\tau _{13.6}`$ $`<`$ 2.5. The reason is that increasing the He abundance, the He++ zone decreases relative to the H+ zone by 36$`\%`$ in volume for this value of T. Thus, a decrease of the nebula optical depth will only affect the He++ zone \[and consequently TZ(He II)\] when TZ(H) is very low. Thus, the area covered by the models stretches toward higher ZR, including the high ZR objects.

Large nebulae have been studied in detail by Kaler and collaborators (Kaler 1981, 1983b; Kaler & Feibelman 1985; Kaler et al. 1990). For many of these nebulae the color temperature obtained from UV observations are well above their Zanstra temperatures. This is consistent with our results (Fig. 2b) that show that the stellar temperature can be higher than TZ(He II). Because of their large diameters, the density in the Abell or other large nebulae may be smaller than the one assumed for the standard models. However, as mentioned above, results for the Zanstra temperatures with different densities are similar.

The above results do not necessarily mean that all the nebulae lying above the limiting curves are He-rich, or, inversely, that all He-rich nebulae have positions above the curves plotted in Figure 2. The same questions can be asked regarding the nebulae size. From the 19 nebulae above the curves, 14 have calculated or limiting values for the abundance; from these, only one have He abundance definitively below the average value for planetary nebulae as given by Kingsburgh & Barlow (1994). However, He-rich nebulae are also found in regions of higher optical depth in the diagrams. Regarding the size, all nebulae above the curves, except two (Hu 1-2, with 0.018 pc, and Cn 1-2, with no value calculated for the radius) have radius larger than 0.15 pc (Cahn, Kaler, & Stanghellini 1992). But large nebulae are spread everywhere in the diagrams. In brief, large and/or He-rich nebulae can be found in any location on the diagram, but most of those above the limiting curves are large and He-rich.

As discussed below, a He-rich atmosphere may provide the explanation for the outsiders with very low ZR. A high He abundance in the inner parts of a nebula can indicate a high He in the upper layers of the star, including their atmospheres. The energy spectrum of a He-rich atmosphere will show a deficit of high energy (E $`>`$ 54.4 eV) photons because of the contribution of He to the stellar continuum opacity. Henry & Shipman (1986) discarded a high He abundance in the stellar atmosphere, since they do not explain the Zanstra discrepancy shown by most PNs (ZR $`>`$ 1), concluding that the atmospheres have subsolar He abundance (and an excess of high-energy photons). Our results show that the main effect originating the Zanstra discrepancy is the optical effect. However, photoionization models assuming an ionizing spectrum with a deficit of photons with energy higher than 54.4 eV indicate that ZR is less than unity in the optically thick case, lowering the limit defined by the models. Results for models with various stellar temperatures and an ionizing spectrum showing a deficit of a factor of five in the flux of high-energy photons, relative to a blackbody, are shown by the dotted lines in Figure 3. The corresponding blackbody results are shown by the solid lines. Thus, such ionizing spectrum, presenting a deficit of high-energy photons, may explain the outsiders with low ZR.

## 4 Stellar temperatures and the He II/He I line intensity ratio

The results presented in Figures 2 and 3 can be used to obtain the nebula optical depth as well as a value for the stellar temperature more accurately than that given by the Zanstra temperatures. However, for optically thin nebulae, the curves are crowded up and the temperature is not well defined. The same occurs for high stellar temperatures, even at high optical depths. Thus, another method for obtaining the stellar temperature is required.

Since the optical depth is a major factor of the stellar temperature determination, line intensity ratios produced by ions in different ionization stages can be used to distinguish nebulae with different optical depths. When plotted against TZ(He II), many of these line ratios also show a crowding of the curves corresponding to different models. The best ratio discriminating the results for different models is the ratio between He II and He I line intensities. The results for He II $`\lambda `$4686/He I $`\lambda `$5786 and He II $`\lambda `$4686/He I $`\lambda `$4471 versus TZ(He II) are presented in Figures 4a and 4b, respectively, for the same models of Figures 1 and 2. Notice that, for a given range of T and $`\tau _{13.6}`$, particularly for low optical depths, the curves in Figure 4 are more widely spaced and provide a better determination of these parameters than the curves shown in Figure 2.

Regarding the nebulae in the crowded region of Figure. 2b, for which there is an uncertainty in the determination of T, the stellar temperature may be obtained from Figure 4. Besides the 19 nebulae above the curves in Figure 2b, 28 nebulae are in the region where the results for T are just lower limits. For all these nebulae, only 20 have measured intensities for He II and He I lines. With intensities taken from the literature (Torres-Peimbert & Peimbert 1977; Aller & Czyzak 1979, 1983; Jacoby & Ford 1983; Kaler 1983a; Kaler 1985b; Manchado, Mampaso & Pottasch 1987; Peimbert & Torres-Peimbert 1987; Kaler et al. 1990; Acker et al. 1991; de Freitas Pacheco, Maciel & Costa 1992; Stanghellini, Kaler & Shaw 1994; Kingsburgh & Barlow 1994), the stellar temperatures obtained from Figure 4 are higher than TZ(He II) by about 10$`\%`$ in general but can reach 32$`\%`$. For the selected nebulae, the maximum observed value for He II $`\lambda `$4686/He I $`\lambda `$5786 is 1.85dex for NGC 4361 (there is also a measured lower limit for its central region of 2.25dex). For He II $`\lambda `$4686/He I $`\lambda `$4471 the maximum measured ratio is 2.37dex, for NGC 2022. Higher ratios, corresponding to lower optical depths, are not measured since He I is too faint to be detected. Notice that for high T ($`>`$ 150,000 K), He I can be very faint even for $`\tau _{13.6}`$ higher than unity, since, as shown in Figure 1, the higher T the smaller the He+ region. So, for nebulae with low ZR and no detected He I line, the stellar temperature can be much higher than TZ(He II). As long as the He II Zanstra temperature is known, this alternative method can be used and may be considered as a second-order approximation to the stellar temperature, providing values closer to the real stellar temperature, even for nebulae showing a high He II/He I line ratio.

## 5 Conclusions

One of the problems concerning the understanding of PNs and their evolution is the determination of their stellar temperature. The Zanstra method is generally used, although the H and He II Zanstra temperatures may be discrepant and may underestimate the stellar temperature. Many authors suggested that the dominant mechanism explaining the Zanstra discrepancy is the optical depth. Because TZ(He II) is less affected by the nebula optical depth, these authors suggest adopting TZ(He II) as the true stellar temperature. Besides the importance of a good determination of the stellar temperature of planetary nebulae for the analysis of the emission-line spectrum and the chemical abundance determination of PNs, let us recall that the position of the central star of PN on the H-R diagram is a crucial test of evolutionary models from the AGB to the white dwarf stages.

Here the nebula optical depth effect is analyzed in detail using photoionization models. From the theoretical H$`\beta `$ and He II $`\lambda `$4686 line intensities, the Zanstra temperatures are calculated and compared to the adopted stellar temperatures for a variety of models with different optical depths. Our results show that the nebula optical depth is the main factor explaining the behavior of the Zanstra temperatures with the stellar temperature. The Zanstra discrepancy mainly occurring for low stellar temperatures is clearly explained by the changes induced by the optical depth in the relative ionic distribution in the nebula (Figs. 1 and 2). Another consequence is that even the He II Zanstra temperature underestimates T, mainly for nebulae with high-temperature stars. The results showing that the stellar temperature can be higher than TZ(He II) are consistent with UV data for the central stars of PNs (Pottasch et al. 1978; Kaler & Feibelman 1985). The relation between the Zanstra temperature ratio ZR and TZ(He II) (Fig. 2) can be used to obtain a more accurate estimate of the stellar temperature. However, for nebulae with very high stellar temperature and/or small optical depths the theoretical results for different stellar temperatures and optical depths are crowded and the method is uncertain. For these nebulae, a better determination of these parameters can be obtained from a plot of He II/He I versus TZ(He II) (Figs. 4).

An important source of uncertainty is related to the kind of observations used to determine the stellar temperature. Observed line intensities do not always refer to the whole nebula since line ratios are usually obtained from observations with a narrow slit crossing the nebula. Furthermore, the nebula may be inhomogeneous and the optical depth anisotropic. That is, the nebula can be optically thick in some direction but optically thin in others. This is a very important point not addressed in this paper. Line intensity ratios obtained with a narrow slit, or in a given position of a nebula, may not correspond to the ratio for the entire nebula (Gruenwald, Viegas, & Broguière 1997; Gruenwald & Viegas 1998).

Cases of very high or very low ZR can be explained by the coupled effect of optical depth and over- or under-abundance of He in the stellar atmosphere, which affects the ionizing spectrum. A nebular overabundance of He relative to solar values and a low optical thickness explain the very high discrepancy shown by some Abell planetaries. On the other hand, ZRs less than unity are characteristic of optically thick PNs with undersolar He abundance.

Stanghellini et al. (1993) state that since bipolar nebulae have TZ(He II) $``$ TZ(H), they are thicker than other nebulae. However, central stars of bipolar nebulae are known to have high temperatures (Corradi & Schwarz 1995). From our results (Figure 2) these nebulae are expected to have TZ(He II) $``$ TZ(H), even if they are not completely optically thick.

Finally, the “Zanstra wall” in the log L - log T diagram (Shaw and Kaler 1989) is related to the fact that the Zanstra method underestimates the stellar temperature and this effect is larger for high-temperature stars. Consequently, the lack of high-temperature stars (predicted by the evolutionary models) in the log L - log T diagram can be easily understood.

###### Acknowledgements.

We are thankful to M. Peimbert for his valuable suggestions. This work was partially supported by FAPESP, CNPq, and PRONEX/FINEP.


text:
# 1 Introduction

## 1 Introduction

Recent interest in Turbo Codes  and Low Density Parity Check Codes  has fuelled development of Factor Graphs and associated Sum-Product Algorithm  (SPA), with applications to error-correction, signal processing, statistics, neural networks, and system theory. Meanwhile the possibility of Quantum Computing has sparked much interest , and Quantum Bayesian Nets have been proposed to help analyse and design Quantum Computers . This paper links these areas of research, showing that quantum resources can achieve maximum-likelihood (ML) decoding of classical information. The natural Hilbert Space of a quantum particle encodes a probability vector, and the joint-state of quantum particles realises the ’products’ associated with SPA. SPA summary is omitted as quantum bits (qubits) naturally encode the total joint-probability state. Dependencies between vector indices become ’entanglement’ in quantum space, with the Factor Graph defining dependency (entanglement) between qubits. Graph function nodes are implemented as unitary matrix 222 ’Unitary’ means that $`𝐔`$ satisfies $`\mathrm{𝐔𝐔}^{}=𝐈`$, where $``$ means ’conjugate transpose’. -vector products followed by quantum measurement. This is the Quantum Product Algorithm (QPA). As QPA avoids summary it avoids problems encountered by SPA on graphs with short cycles. Moreover, whereas SPA is iterative, using message-passing and activating each node more than once, QPA does not iterate but must successfully activate each node only once. However the (severe) drawbacks with QPA are as follows: 1) Each function node must be repeatedly activated until it successfully ’prepares’ it’s local variable nodes (qubits) in the correct entangled state - any activation failure destroys evolution in all variable nodes already entangled with local variables. 2) Once a complete Factor Graph has successfully evolved, final quantum measurement only delivers the ML codeword with a certain (largest) probability. Repeated successful evolutions then determine the ML codeword to within any degree of confidence. This second drawback can be overcome by suitable ”ML Amplification” of QPA output, prior to measurement.

Section 2 presents QPA, highlighting its ability to deliver the optimal output joint-state, unlike SPA. Quantum systems describe the exact joint-state by appropriate ’entanglement’ with and measurement of ancillary qubits. Section 3 considers a simple example of QPA on Quantum Factor Graphs, showing that iteration on graphs with cycles is unnecessary because QPA avoids premature summary. Section 4 shows how to amplify the likelihood of measuring the ML codeword from QPA output. Unfortunately QPA must be repeated many times and/or executed in parallel to have a hope of successful completion. Suitable distributed QPA scheduling is discussed in Section 5, and it is argued that successful QPA completion is conceivable using asynchronous distributed processing on many-node Factor Graphs. This paper does not deal with phase properties of quantum computers. It is expected that the inclusion of phase and non-diagonal unitary matrices will greatly increase functionality of the Quantum Factor Graph.

The aim of this paper is not to propose an immediately realisable implementation of a quantum computer. Rather, it is to highlight similarities between graphs for classical message-passing, and graphs that ’factor’ quantum computation. The paper also highlights the differences between the two graphs: whereas classical graphs can only ever compute over a tensor product space, the quantum graph can compute over the complete entangled (tensor-irreducible) space.

## 2 The Quantum Product Algorithm (QPA)

### 2.1 Preliminaries

Consider the Factor Graph of Fig 1.

Let $`𝐔_𝐟=\left(\begin{array}{cccc}f_0& 0& 0& 0\\ 0& f_1& 0& 0\\ 0& 0& f_2& 0\\ 0& 0& 0& f_3\end{array}\right)`$, and $`𝐔_𝐠=\left(\begin{array}{cccc}g_0& 0& 0& 0\\ 0& g_1& 0& 0\\ 0& 0& g_2& 0\\ 0& 0& 0& g_3\end{array}\right)\begin{array}{c}\text{, where }|g_k|^2=1|f_k|^2\text{,}\\ \text{and }f_k^{}g_k+f_kg_k^{}=0\text{}k\text{.}\\ \text{}\text{’ means complex conjugate.}\end{array}`$ Let $`𝐔_{\mathrm{𝐟𝐠}}=\left(\begin{array}{cc}𝐔_𝐟& 𝐔_𝐠\\ 𝐔_𝐠& 𝐔_𝐟\end{array}\right)`$. $`𝐔_{\mathrm{𝐟𝐠}}`$ is unitary, and the $`𝐔_𝐟`$ of Fig 1 and subsequent figures always implies the action of $`𝐔_{\mathrm{𝐟𝐠}}`$ together with the measurement of an ancillary qubit, $`z`$, as described below. A qubit, $`x_i`$, can be in states $`0`$ or $`1`$ or in a statistical superposition of $`0`$ and $`1`$. Let qubits $`x_0,x_1`$ be initialised (by the black boxes) to states $`x_0=(\alpha _0,\beta _0)^T`$ and $`x_1=(\alpha _1,\beta _1)^T`$, where $`\alpha _i,\beta _i`$ are complex probabilities such that $`|\alpha _i|^2+|\beta _i|^2=1`$. For instance, $`x_0`$ is in states $`0`$ and $`1`$ with probabilities $`|\alpha _0|^2`$ and $`|\beta _0|^2`$, respectively. Let an ancillary qubit, $`z`$, be initialised to state $`0`$, i.e. $`z=(1,0)`$. Then the initial joint probability product-state of qubits $`x_0,x_1,z`$ is $`𝐀=(\alpha _0,\beta _0)^T(\alpha _1,\beta _1)^T(1,0)^T=(\alpha _0\alpha _1,\beta _0\alpha _1,\alpha _0\beta _1,\beta _0\beta _1,0,0,0,0)^T=(s_0,s_1,s_2,s_3,0,0,0,0)^T`$, where $`|s_0|^2+|s_1|^2+|s_2|^2+|s_3|^2=1`$, and ’$``$’ is the tensor product. The element at vector index $`v`$ is the probability that the qubits are in state $`v`$. For instance, qubits $`x_0x_1z`$ are in joint-state $`010`$ with probability $`|s_2|^2`$. Subsequent measurement of a subset of the qubits projects the measured qubits to a fixed substate with a certain probability, and ’summarises’ the vector for the remaining non-measured qubits. Thus QPA is as follows,

* Compute $`𝐒=𝐔_{\mathrm{𝐟𝐠}}𝐀`$.
* Measure qubit $`z`$. With probability $`p_f=|s_0f_0|^2+|s_1f_1|^2+|s_2f_2|^2+|s_3f_3|^2`$ we collapse $`z`$ to $`0`$, and $`x_0,x_1`$ to joint-state $`𝐒_𝐟=\mu _0(s_0f_0,s_1f_1,s_2f_2,s_3f_3)^T`$. With probability $`p_g=|s_0g_0|^2+|s_1g_1|^2+|s_2g_2|^2+|s_3g_3|^2`$ we collapse $`z`$ to $`1`$, and $`x_0,x_1`$ to joint-state $`𝐒_𝐠=\mu _1(s_0g_0,s_1g_1,s_2g_2,s_3g_3)^T`$. $`\mu _0`$ and $`\mu _1`$ are normalisation constants. $`p_f+p_g=1`$. $`𝐒_𝐟`$ is our desired QPA result. Successful QPA completion is self-verified when we measure $`z=0`$.

In contrast, classical SPA computes $`𝐒_𝐟=𝐔_𝐟𝐀`$ (with probability 1) and must then perform a subsequent ’summary’ step on $`𝐒_𝐟`$ before returning a result for each variable separately. This result is,
$`x_0=|\mu _0|^2(|s_0f_0|^2+|s_2f_2|^2,|s_1f_1|^2+|s_3f_3|^2)^T`$, $`x_1=|\mu _1|^2(|s_0f_0|^2+|s_1f_1|^2,|s_2f_2|^2+|s_3f_3|^2)^T`$.
For instance, for $`x_0=0`$ we sum the two classical 333Classical SPA probabilities in this paper are always represented as the magnitude-squared of their quantum counterparts probabilities of $`𝐒_𝐟`$ where $`x_0=0`$ to get $`|s_0f_0|^2+|s_2f_2|^2`$. Similarly, for $`x_0=1`$ we summarise to $`|s_1f_1|^2+|s_3f_3|^2`$. It is in this sense that SPA is a ’tensor-approximation’ of QPA.

We identify the following successively accurate computational scenarios (decoding modes) for a space of $`N`$ binary-state variables:

* Hard-Decision operates on a probability space,
$`(\alpha _0,\beta _0)(\alpha _1,\beta _1)\mathrm{}(\alpha _{N1},\beta _{N1}),\text{ }\alpha ,\beta \{0,1\}`$
* Soft-Decision operates on a probability space,
$`(\alpha _0,\beta _0)(\alpha _1,\beta _1)\mathrm{}(\alpha _{N1},\beta _{N1}),\text{ }\alpha ,\beta \{\text{Real Numbers }01\}`$
* Quantum Soft-Decision operates on a probability space,
$`(\alpha _0,\beta _0)(\alpha _1,\beta _1)\mathrm{}(\alpha _{N1},\beta _{N1}),\text{ }\alpha ,\beta \{\text{Complex Numbers}\}`$
* Entangled-Decision operates on a probability space,
$`(s_0,s_1,s_2,\mathrm{},s_{2^N1})`$, $`s\{\text{Complex Numbers}\}`$

All four of the above Decision modes satisfies the probability restriction that the sum of the magnitude-squareds of the vector elements is 1. Both Quantum Soft-Decision and Entangled-Decision make use of the natural quantum statistical properties of matter, including the property of Superposition. Moreover, Entangled-Decision operates over exponentially larger space. Classical SPA operates in Soft-Decision mode. QPA operates in Entangled-Decision mode. In the previous discussion it was assumed that the QPA was operating on input of the form, $`(\alpha _0,\beta _0)^T(\alpha _1,\beta _1)^T(1,0)^T`$. More generally, QPA can operate on input and deliver output in Entangled-Decision mode. This is in strong contrast to SPA which must summarise both input and output down to Soft-Decision mode. It is this approximation that forces SPA to iterate and to sometimes fail on graphs with cycles.

Consider the following example. If the diagonal of $`𝐔_𝐟`$ is $`(1,0,0,1)`$, then $`𝐔_𝐟`$ represents XOR, and Fig 1 decodes to codeset $`𝐂=\{00,11\}`$ (i.e. $`x_0+x_1=0,\text{ mod }2`$). $`𝐂`$ has distance 2, which is optimal for length 2 binary codes: in general if $`𝐔_𝐟`$ cannot be tensor-decomposed then it represents a code $`𝐂`$ with good distance properties. Initially, let $`x_0=(\sqrt{0.4},\sqrt{0.6})^T`$, $`x_1=(\sqrt{0.6},\sqrt{0.4})^T`$. Then $`𝐀=(\sqrt{.24},0.6,0.4,\sqrt{0.24},0,0,0,0)^T`$, and $`𝐒_𝐟=\frac{1}{\sqrt{2}}(1,0,0,1)^T`$. $`p_f=0.48`$, so, on average, $`48`$ $`𝐒_𝐟`$ outputs are computed for every $`100`$ QPA attempts. The ML codeword is both $`00`$ and $`11`$, and when $`𝐒_𝐟`$ is measured, $`00`$ and $`11`$ are equally likely to be returned. In contrast, classical SPA for the same input returns $`x_0=x_1=(\frac{1}{2},\frac{1}{2})`$, implying (wrongly) an equally likely decode to any of the words $`00,01,10,11`$. So even in this simplest example the advantage of QPA over SPA is evident.

### 2.2 Product Space for Classical SPA

Because $`x_0`$ and $`x_1`$ are separated in Fig 1, their classical joint-state only represents tensor product states (Soft-Decision mode). An equivalent Factor Graph to that of Fig 1 could combine $`x_0`$ and $`x_1`$ into one quaternary variable which would reach all non-product quaternary states. But this requires ’thickening’ of graph communication lines and exponential increase in SPA computational complexity. Consequently only limited variable ’clustering’ is desirable, although too little clustering ’thins out’ the solution space to an insufficient highly-factored product space. This is the fundamental Factor Graph trade-off - good Factor Graphs achieve efficient SPA by careful variable ’separation’, ensuring the joint product space is close enough to the exact (non-summarised) non-product space.

### 2.3 Entangled Space for QPA

In contrast, although $`x_0`$ and $`x_1`$ are physically separated in Fig 1, quantum non-locality must take into account correlations between $`x_0`$ and $`x_1`$. Their joint-state now occurs over the union of product and (much larger) non-product (entangled) space (Entangled-Decision mode). An entangled joint-state vector cannot be tensor-factorised over constituent qubits. QPA does not usually output to product space because the joint-state of output qubits is usually entangled. In fact QPA is algorithmically simpler than SPA, as SPA is a subsequent tensor approximation of QPA output at each local function.

### 2.4 Example

Let the diagonal of $`𝐔_𝐟`$ be $`(1,0,0,1)`$. Initialise $`x_0`$ and $`x_1`$ to joint-product-state, $`x_0=\frac{1}{\sqrt{3}}(1,\sqrt{2})^T`$, $`x_1=\frac{1}{\sqrt{2}}(1,1)^T`$. With probability $`p_f=0.5`$ QPA measures $`z=0`$ and computes the joint-state of $`x_0,x_1`$ as $`𝐒_𝐟=\frac{1}{\sqrt{3}}(1,0,0,\sqrt{2})^T`$. A final measurement of qubits $`x_0`$ and $`x_1`$ yields codewords $`11`$ and $`00`$ with probability $`\frac{2}{3}`$, and $`\frac{1}{3}`$, respectively. In contrast SPA summarises $`𝐒_𝐟`$ to $`x_0=x_1=\frac{1}{3}(1,2)`$. Although a final ’hard-decision’ on $`x_0`$ and $`x_1`$ chooses, correctly, the ML codeword $`x_0=x_1=1`$, the joint-product-state output, $`\frac{1}{3}(1,2)^T\frac{1}{3}(1,2)^T=\frac{1}{9}(1,2,2,4)^T`$ assigns, incorrectly, a non-zero probability to words $`01`$ and $`10`$.

### 2.5 A Priori Initialisation

To initialise $`x_0`$ to $`(\alpha _0,\beta _0)^T`$, we again use QPA. Let the diagonal of $`𝐔_𝐟`$ (for the left-hand black box of Fig 1) be $`(\alpha _0,\beta _0)`$. Then the diagonal of $`𝐔_𝐠`$ is $`\pm i(\frac{\alpha _0\sqrt{1|\alpha _0|^2}}{|\alpha _0|},\frac{\beta _0\sqrt{1|\beta _0|^2}}{|\beta _0|})^T`$. Measurement of $`z=0`$ initialises $`x_0`$ to $`(\alpha _0,\beta _0)^T`$, and this occurs with probability $`p_f=0.5`$. $`x_1`$ is initialised likewise.

### 2.6 Comments

The major drawback of QPA is the significant probability of QPA failure, occurring when $`z`$ is measured as $`1`$. This problem is amplified for larger Quantum Factor Graphs where a different $`z`$ is measured at each local function; QPA evolution failure at a function node not only destroys the states of variables connected with that function, but also destroys all states of variables entangled with those variables. QPA is more likely to succeed when input variable probabilities are already skewed somewhat towards a valid codeword. Section 3 shows how QPA can operate successfully even when SPA fails.

## 3 Quantum Product Algorithm on Factor Graphs with Cycles

This section shows that graph cycles do not compromise QPA performance. Consider the Factor Graph of Fig 2.

Functions $`𝐔_{\mathrm{𝐟𝟎}}`$ and $`𝐔_{\mathrm{𝐟𝟏}}`$ are both $`8\times 8`$ XOR diagonal matrices with diagonal elements $`(10010110)`$. Acting on the combined four-qubit space, $`x_0x_1x_2x_3`$, they are the functions $`𝐔_{\mathrm{𝐟𝟎}}𝐈_\mathrm{𝟐}`$ and $`𝐈_\mathrm{𝟐}𝐔_{\mathrm{𝐟𝟏}}`$, respectively, with diagonal elements $`(1001011010010110)`$ and $`(1100001100111100)`$, respectively, where $`𝐈_\mathrm{𝟐}`$ is the $`2\times 2`$ identity matrix. QPA on Fig 2 performs the global function $`𝐔_𝐅=\left(𝐔_{\mathrm{𝐟𝟎}}𝐈_\mathrm{𝟐}\right)\left(𝐈_\mathrm{𝟐}𝐔_{\mathrm{𝐟𝟏}}\right)`$ on four-qubit space, with diagonal elements $`(1000001000010100)`$, forcing output into codeset $`𝐂=\{0000,0110,1011,1101\}`$. Functions $`𝐔_{\mathrm{𝐟𝟎}}`$, $`𝐔_{\mathrm{𝐟𝟏}}`$, and $`𝐔_𝐅`$ ’sieve’ the input joint-state, where $`𝐔_𝐅`$ is the combination of two ’sub-sieves’, $`𝐔_{\mathrm{𝐟𝟎}}`$ and $`𝐔_{\mathrm{𝐟𝟏}}`$. QPA iteration (i.e. successfully completing a sub-function more than once on the same qubits) has no purpose, as only one needs apply a particular sieve once. So graph cycles have no bearing on QPA. (However iteration may be useful to maintain the entangled result in the presence of quantum decoherence and noise). To underline cycle-independence, consider the action of SPA, then QPA on Fig 2.

Initialise as follows (using classical probabilities),

$$x_0x_1x_2x_3=(0.1,0.9)^T(0.6,0.4)^T(0.6,0.4)^T(0.6,0.4)^T$$

Hard-decision gives $`x_0x_1x_2x_3=1000`$, which can then be decoded algebraically to codeword $`0000`$. However optimal soft-decision would decode to either $`x_3x_2x_1x_0=1011`$ or $`1101`$, with equal probability. Because of the small graph cycle SPA fails to decode correctly, and settles to the joint-product-state,
$`x_0x_1x_2x_3=(0.108,0.892)^T(0.521,0.479)^T(0.521,0.479)^T(0.601,0.399)^T`$. A final hard-decision on this output gives non-codeword $`x_0x_1x_2x_3=1000`$ which can then be decoded algebraically, again to codeword $`0000`$. In contrast, successful QPA outputs the optimal entangled joint-state,
$`S_F=\frac{1}{\sqrt{2040}}(\sqrt{216},0,0,0,0,0,\sqrt{96},0,0,0,0,\sqrt{864},0,\sqrt{864},0,0)^T`$. Final measurement of $`S_F`$ always outputs a codeword from $`𝐂`$, and with probability $`\frac{2864}{2040}`$ outputs either $`1011`$ or $`1101`$. QPA evolves on Fig 2 correctly with probability $`0.204`$. Therefore $`1000`$ attempts produce around $`204`$ correctly entangled joint-states.

To underline QPA advantage, consider the single variable extension of Fig 2 in Fig 3, where $`x_4`$ is initialised to $`(\sqrt{0.5},\sqrt{0.5})^T`$.

As $`x_4=x_0x_3`$, and our original code, $`𝐂=\{0000,0110,1011,1101\}`$, always had $`x_0=x_3`$, then $`x_4`$ should always be $`0`$. But SPA on Fig 3 computes $`x_4=(0.421,0.579)^T`$ and subsequent hard-decision gives $`x_4=1`$. In contrast, successful QPA computes the optimal non-product joint-state,

$$S_F^{}=\frac{1}{\sqrt{2040}}(\sqrt{216},0,0,0,0,0,\sqrt{96},0,0,0,0,\sqrt{864},0,\sqrt{864},0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)^T$$

Final measurement of $`S_F^{}`$ always outputs $`x_4=0`$. QPA evolves on Fig 3 correctly with probability $`0.2040.5=0.114`$.

## 4 Maximum-Likelihood (ML) Amplification

### 4.1 Preliminaries

The ML codeword is the one most likely to be measured from QPA output, with probability, $`p_M`$, say. For instance, if QPA output of Fig 1 is $`S_f=\frac{1}{\sqrt{3}}(1,0,0,\sqrt{2})^T`$, say, then $`11`$ is the ML codeword, and it is measured with probability $`p_M=\frac{2}{3}`$. Numerous executions of QPA on the same input will verify that $`11`$ is, indeed, the ML codeword. However these numerous executions must output to a length $`2^N`$ final averaging probability vector (for $`N`$ qubits). We do not want to store such an exponential vector. Instead, therefore, we ’amplify’ the statistical advantage of $`11`$ over $`00`$ prior to measurement, thereby making $`11`$ significantly more likely to be read. This is achieved by computing the square of each quantum vector element as follows. Consider two independent QPA executions on the same input, both outputting $`S_f`$. Associate these outputs with qubits $`x_{0,0},x_{1,0}`$, and $`x_{0,1},x_{1,1}`$. The joint-state of qubits $`x_{0,0},x_{1,0},x_{0,1},x_{1,1}`$ is,

$$V_0=S_fS_f=\frac{1}{3}(1,0,0,\sqrt{2},0,0,0,0,0,0,0,0,\sqrt{2},0,0,2)^T$$

Consider the unitary permutation matrix

$$𝐏=\left(\begin{array}{c}1000000000000000\\ 0000010000000000\\ 0000000000100000\\ 0000000000000001\\ 0100000000000000\\ 0010000000000000\\ 0001000000000000\\ 0000100000000000\\ 0000001000000000\\ 0000000100000000\\ 0000000010000000\\ 0000000001000000\\ 0000000000010000\\ 0000000000001000\\ 0000000000000100\\ 0000000000000010\end{array}\right)$$

Only the ’1’ positions in the first four rows are important. Performing $`𝐏`$ on $`x_{0,0},x_{1,0},x_{0,1},x_{1,1}`$, gives,

$$𝐏V_0=\frac{1}{3}(1,0,0,2,0,0,\sqrt{2},0,0,0,0,0,0,\sqrt{2},0,0)^T$$

We then measure qubits $`x_{0,1},x_{1,1}`$. With probability $`p_{a_0}=\frac{5}{9}`$ we read $`x_{0,1}=x_{1,1}=0`$, in which case $`x_{0,0}`$ and $`x_{1,0}`$ are forced into joint state $`S_{f,1}=\frac{1}{\sqrt{5}}(1,0,0,2)`$, which is the element-square of $`S_f`$. A measurement of $`S_{f,1}`$ returns $`11`$ with probability $`p_M=\frac{4}{5}`$, which is a significant improvement over $`p_M=\frac{2}{3}`$. Likewise we compute the element fourth-powers of $`S_f`$ by preparing two independent qubit pairs in $`S_{f,1}`$ and permuting the (umeasured) joint state vector $`V_1=S_{f,1}S_{f,1}`$ to give $`𝐏V_1`$, and then measuring the second pair of qubits. With probability $`p_{a_1}=\frac{17}{25}`$ we read this pair as $`00`$, in which case the first two qubits are forced into the joint-state $`S_{f,2}=\frac{1}{\sqrt{1}7}(1,0,0,4)`$, which is the element fourth-power of $`S_f`$. A measurement of $`S_{f,2}`$ returns $`11`$ with probability $`p_M=\frac{16}{17}`$, which is a further improvement over $`p_M=\frac{2}{3}`$. In this way we amplify the likelihood of measuring the ML codeword. To compute the element $`2^k`$$`^{\text{th}}`$-power, $`S_{f,k}`$, we require, on average, $`\frac{2}{p_{a_k}}`$ independent preparations, $`S_{f_{k1}}`$, each of which requires, on average, $`\frac{2}{p_{a_{k1}}}`$ independent preparations, $`S_{f_{k2}}`$, and so on.

We can perform QPA on large Factor Graphs, then amplify the result $`k`$ times to ensure a high likelihood of measuring the ML codeword, as described above. However the above amplification acts on the complete graph with one operation, $`𝐏`$. It would be preferable to decompose $`𝐏`$ into $`4\times 4`$ unitary matrices which only act on independent qubit pairs $`x_{i,0}`$ and $`x_{i,1}`$, thereby localising amplification. Consider, once again, Fig 1. From the point of view of $`x_{0,1}`$, $`x_{0,0}`$ appears to be in summarised state 444$`x_{0,0}`$ is generally not in this summarised state, due to phase considerations, but the viewpoint is valid for our purposes as long as subsequent unitary matrix operations on $`x0`$ only have one non-zero entry per row. , $`s_f=\frac{1}{\sqrt{3}}(1,\sqrt{2})^T`$. Similarly, from the point of view of $`x_{0,0}`$, $`x_{0,1}`$ appears to be in state $`s_f`$. Thus $`x_{0,0},x_{0,1}`$ appear to be in joint product state $`v_0=\frac{1}{3}(1,\sqrt{2},\sqrt{2},2)^T`$. Consider unitary permutation matrix,

$$𝐐=\left(\begin{array}{c}1000\\ 0001\\ 0100\\ 0010\end{array}\right)$$

We compute $`𝐐v_0=\frac{1}{3}(1,2,\sqrt{2},\sqrt{2})^T`$ on qubits $`x_{0,0},x_{0,1}`$ and measure qubit $`x_{0,1}`$. With probability $`p_{a_0}=\frac{5}{9}`$ we read $`x_{0,1}=0`$, in which case $`x_{0,0}`$ is forced into joint state $`s_{f,1}=\frac{1}{\sqrt{5}}(1,2)`$, which is the element-square of $`s_f`$. Due to the exact form of our joint-state vector, $`S_f`$, this single measurement is enough to also force $`x_{0,0}x_{1,0}`$ into joint state $`S_{f,1}`$. However, for a general function $`S_f`$, we should perform $`𝐐`$ on every qubit pair, $`x_{i,0}x_{i,1}`$, then measure $`x_{i,1}`$ $`i`$. This is equivalent to performing $`𝐏^{}=𝐐𝐐`$ on (re-ordered) joint-state vector $`x_{0,0}x_{0,1}x_{1,0}x_{1,1}`$, and this is identical to performing $`𝐏`$ on $`x_{0,0}x_{1,0}x_{0,1}x_{1,1}`$. The probability of measuring $`x_{1,0}=x_{1,1}=0`$ is the same whether $`𝐏`$ or $`𝐐`$ is used. The same process is followed to achieve element $`2^k`$th powers.

### 4.2 The Price of Amplification

There is a statistical cost to qubit amplification. Let $`s=(\alpha ,\beta )^T`$ be the initial state of a qubit $`x`$, where, for notational convenience, we assume that $`\alpha `$ and $`\beta `$ are both real. Then $`\alpha ^2+\beta ^2=1`$ and, given $`2^k`$ qubits all identically prepared in state $`s`$, the likelihood of preparing one qubit in (unnormalised) state $`s_k=(\alpha ^{2^k},\beta ^{2^k})^T`$ is $`\gamma _k`$, where,

$$\gamma _k=\gamma _{k1}^2\frac{r_{k+1}}{r_k^2},\text{ }\gamma _0=1$$

and $`r_k=\alpha ^{2^k}+\beta ^{2^k}`$. For a qubit in state $`s_k`$, the probability of selecting the ML codebit is,

$$P_{Mk}=\frac{\alpha ^{2^{k+1}}}{\alpha ^{2^{k+1}}+\beta ^{2^{k+1}}}$$

(assuming $`\alpha \beta `$). We can plot $`\gamma _k`$ against $`P_{Mk}`$ for various $`\alpha ^2`$ as $`k`$ varies, as shown in Fig 4.

Each of the 25 lines in Fig 4 refers to a different value of $`\alpha ^2`$, for $`\alpha ^2`$ from $`0.5`$ up to $`0.98`$ in steps of $`0.2`$. The initial state, $`s`$, when $`k=0`$, occurs with probability $`\gamma _k=1`$, and is marked on the right-hand side of Fig 4 for each of the 25 lines. After one amplification step, $`k=1`$, and another 25 points are marked on the graph to the left of the points for $`k=0`$, indicating that a successful amplification step has occurred with probability $`\gamma _k1`$. Similarly points for $`k=2`$, $`k=3`$,..etc are marked successively to the left on Fig 4. The $`y`$-axis shows the ML advantage, $`P_{Mk}`$, which can be achieved with probability $`\gamma _k`$ after $`k`$ steps for each value of $`\alpha ^2`$. For instance, when $`s=(\alpha ,\beta )^T=(\sqrt{0.62},\sqrt{0.38})^T`$, then an ML advantage of $`P_{Mk}=0.9805`$ can be ensured after $`k=3`$ steps, and this can be achieved with probability $`\gamma _k=0.0223`$ given $`2^3=8`$ independently prepared qubits, all in state $`s`$. Amplification is more rapid if $`s`$ already has significant ML advantage (i.e. when $`\alpha `$ is high). In contrast if $`\alpha ^2=0.5`$ then no amplification of that qubit is possible. This is quite reasonable as, in this case, both states $`0`$ and $`1`$ are equally likely, so there is no ML state. Successive measurement of zero of all second qubits of each qubit pair self-verifies that we have obtained successful amplification. If, at any step, $`k`$, the second qubit of the qubit pair is measured as one then amplification fails and the graph local to this qubit which has been successfully entangled up until now is completely destroyed.

## 5 Distributed QPA on Many-Node Factor Graphs

### 5.1 Preliminaries

In classical systems it is desirable to implement SPA on Factor Graphs which ’tensor-approximate’ the variable space using many small-state variables (e.g. bits), linked by small-dimensional constituent functions, thereby minimising computational complexity. In quantum systems it is similarly desirable to implement QPA on Factor Graphs using many small-state variables (e.g. qubits), linked by small-dimensional constituent unitary functions. Any Quantum Computation can be decomposed into a sequence of one or two-bit ’universal’ gate unitary operations 555This also implies that any classical Factor Graph can be similarly decomposed. . Computational complexity is minimised by using small-dimensional unitary matrices for constituent functions. Moreover, fine granularity of the Factor Graph allows distributed node processing. This appears to be essential for large Quantum Factor Graphs to have acceptable probability of successful global evolution, as we will show. Distributed QPA allows variable nodes to evolve entanglement only with neighbouring variable nodes so that, if a local function measurement or amplification is unsuccessful, only local evolution is destroyed. Remember that local evolution is OFTEN unsuccessful, as failure occurs when a local ancillary qubit, $`z`$, is measured as $`1`$, or when a local amplifying qubit is measured as 1. Therefore node localities with high likelihood of successful evolution (i.e. with positively skewed input probabilities) are likely to evolve first. These will then encourage other self-contradictory node localities to evolve successfully. In contrast, non-distributed QPA on large Factor Graphs using one large global function is very unlikely to ever succeed, especially for graphs encoding low-rate codes.

To illustrate the advangtage of distributed QPA, consider the low rate code of Fig 5, where $`𝐔_{\mathrm{𝐟𝐢𝐣}}=\text{diag}(1,0,0,1)`$. Both top and bottom graphs represent the code $`𝐂=\{0000,1111\}`$, where $`𝐔`$ is a combination of XOR sub-matrices, $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$,$`𝐔_{\mathrm{𝐟𝟏𝟐}}`$, and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$. The top graph distributes processing. We allow $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$ to operate independently and in parallel. Moreover, if $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ fails to establish, then it does not destroy any successful evolution of $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$, as the two localities are not currently entangled. Once both $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$ have completed successfully, the subsequent probability of successful completion of $`𝐔_{\mathrm{𝐟𝟏𝟐}}`$ is, in general, likely to increase. So distributing QPA increases likelihood of successful evolution of the complete Factor Graph. We now demonstrate this graphically. Let qubits $`x_0,x_1,x_2,x_3`$ of Fig 5 initially be in states $`x_0=(\alpha _0,\beta _0)^T`$, $`x_1=(\alpha _1,\beta _1)^T`$, $`x_2=(\alpha _2,\beta _2)^T`$, $`x_3=(\alpha _3,\beta _3)^T`$, where, for notational convenience, we assume all values are real. Then $`\alpha _i^2+\beta _i^2=1`$, $`i`$. The probability of successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ is $`p_{f01}=(\alpha _0\alpha _1)^2+(\beta _0\beta _1)^2`$, and probability of successful completion of $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$ is $`p_{f23}=(\alpha _2\alpha _3)^2+(\beta _2\beta _3)^2`$. Therefore the probability of successful completion of both $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$ after exactly $`q`$ parallel attempts (no less) is,

$$p_{03}\left(q\right)=\left(1p_{f01}\right)^{q1}\left(1p_{f23}\right)^{q1}p_{f01}p_{f23}+\left(1p_{f01}\right)^{q1}\left(1\left(1p_{f23}\right)^{q1}\right)p_{f01}+\left(1p_{f23}\right)^{q1}\left(1\left(1p_{f01}\right)^{q1}\right)p_{f23}$$

Given successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$, the probability of subsequent successful completion of $`𝐔_{\mathrm{𝐟𝟏𝟐}}`$ is,

$$p_{f12}^{}=\frac{(\alpha _0\alpha _1\alpha _2\alpha _3)^2+(\beta _0\beta _1\beta _2\beta _3)^2}{p_{f01}p_{f23}}$$

Therefore the probability of successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$, immediately followed by successful completion of $`𝐔_{\mathrm{𝐟𝟏𝟐}}`$ is, $`p_{03}(q)=p_{03}(q1)p_{f12}^{}`$, and the probability of successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$, immediately followed by completion failure of $`𝐔_{\mathrm{𝐟𝟏𝟐}}`$ is, $`\overline{p_{03}}(q)=p_{03}(q1)(1p_{f12}^{})`$. Therefore the probability of successful completion after exactly $`t`$ steps of $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$ in parallel, followed by $`𝐔_{\mathrm{𝐟𝟏𝟐}}`$, is,

$$p_e(t)=\underset{q=2}{\overset{t}{}}p_{03}(q)\underset{𝐯𝐃(𝐭𝐪)}{}\underset{u𝐯}{}\overline{p_{03}}(u)$$

where $`𝐃(𝐤)`$ is the set of unordered partitions of $`k`$. Therefore the probability of successful completion after at most $`t`$ steps of $`𝐔_{\mathrm{𝐟𝟎𝟏}}`$ and $`𝐔_{\mathrm{𝐟𝟐𝟑}}`$ in parallel, followed by $`𝐔_{\mathrm{𝐟𝟏𝟐}}`$, is,

$$p_m(t)=\underset{i=2}{\overset{t}{}}p_e(i)$$

In contrast, for non-distributed QPA, the probability of successful completion, after at most $`t`$ steps, of $`𝐔`$, (the bottom graph of Fig 5) is $`P(t)=1(1(\alpha _0\alpha _1\alpha _2\alpha _3)^2(\beta _0\beta _1\beta _2\beta _3)^2)^t`$. Figs 6 and 7 show plots of $`p_m(t)`$ and $`P(t)`$ versus $`t`$ for $`\alpha _0=\alpha _1=\alpha _2=\alpha _3=w`$ as $`w`$ varies, and $`\alpha _0=u,\alpha _1=\alpha _2=\alpha _3=w=0.9`$ as $`u`$ varies, respectively. For Fig 7, low values of $`u`$ indicate a contradiction between $`x_0`$ and the other three variables. In particular the contradiction is so pronounced when $`\alpha _0=0.0`$ that successful QPA completion is highly unlikely. More generally, this indicates that severe internal Factor Graph contradictions are fatal to QPA (as they are for SPA). Both Fig 6 and 7 indicate that, due to initial latency of distributed processing, non-distributed QPA appears marginally faster for the first few steps. However, after a few steps distributed QPA in general becomes marginally faster. In fact results are unfairly biased towards the non-distributed case, as it is assumed that attempts to complete $`𝐔`$ and $`𝐔_{\mathrm{𝐟𝐢𝐣}}`$ have the same space-time-complexity cost, whereas $`𝐔`$ is far more costly. Hence, even for this smallest example, Distributed QPA outperforms non-Distributed QPA.

The example of Fig 5 only achieves marginal advantage using Distributed QPA because the example has so few nodes. The advantage is more pronounced in Fig 8.

Fig 8 represents the code $`𝐂=\{000000000,111111111\}`$ 666This code is trivial but demonstrates a ’worst-case’ low-rate scenario. In general, codes of higher rate, with or without cycles, decode more quickly. , where $`𝐔_{\mathrm{𝐢𝐣𝐤}}=\text{diag}(1,0,0,1,0,1,1,0)`$. We allow $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$ to operate independently and in parallel. If $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$ fails to establish, then it does not destroy any successful evolution of $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$ or $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$, as the three localities are not currently entangled. Once $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$ have completed successfully, the probability of successful subsequent completion of $`𝐔_{\mathrm{𝐟𝟐𝟓𝟖}}`$ is, in general, amplified. Let qubits $`x_i,0i<9`$ of Fig 8 initially be in states $`x_i=(\alpha _i,\beta _i)^T`$, where, for notational convenience, we assume all values are real. Then $`\alpha _i^2+\beta _i^2=1`$, $`i`$. Let the probability of successful completion after at most $`t`$ steps of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$ in parallel, followed by $`𝐔_{\mathrm{𝐟𝟐𝟓𝟖}}`$, be $`p_m(t)`$, and the probability of successful completion, after at most $`t`$ steps, of a non-distributed version of Fig 8 be $`P(t)`$. Appendix A derives $`p_m(t)`$ and $`P_t`$ for this case. Figs 9 and 10 show plots of $`p_m(t)`$ and $`P(t)`$ versus $`t`$ for $`\alpha _i=w`$, $`i`$, as $`w`$ varies, and $`\alpha _0=u,\alpha _i=w=0.9`$, $`i`$, $`i0`$, as $`u`$ varies, respectively. For Fig 10 low values of $`u`$ indicate contradiction between $`x_0`$ and the other eight qubits. The contradiction is so pronounced when $`\alpha _0=0.0`$ that successful QPA completion is highly unlikely.

Figs 11 and 12 show plots of $`p_m(t)`$ and $`P(t)`$ versus $`t`$ for $`\alpha _0=\alpha _1=u,\alpha _i=w=0.9`$, $`i`$, $`i0,1`$, and $`\alpha _0=\alpha _8=u,\alpha _i=w=0.9`$, $`i`$, $`i0,8`$, respectively, as $`u`$ varies. Both figures indicate contradictions between two qubits and the rest, but the scattered nature of contradictions for Fig 12 ($`x_0`$ and $`x_8`$ are connected to different local functions) enhances Distributed QPA performance compared to Fig 11.

Figs 6-12 indicate that distributed QPA completes significantly faster than non-distributed QPA, in particular for cases requiring many steps, $`t`$. Even more so as the presented results are unfairly biased towards the non-distributed case, as it is assumed that attempts to complete non-distributed $`𝐔`$ or each constituent $`𝐔_{\mathrm{𝐟𝐢𝐣𝐤}}`$ have the same space-time-complexity cost, whereas $`𝐔`$ is far more costly. We conclude that Distributed QPA is essential for large Quantum Factor Graphs.

### 5.2 Free-Running Distributed QPA

Consider the notional Factor Graph of Fig 13. Each (square) function node activates time-independently on its local (circular) variable nodes. Functions successfully completed are marked with an ’X’. After a certain time, say, three ’areas of success’ evolve, due to general agreement between input variable states at these localities. This means that variables on the perimeter of each region of success are ’encouraged’ to agree with the ’general view’ of the associated region of success. Unfortunately, in the bottom left of the graph is a variable (dark circle) which strongly contradicts with the rest of the graph. No area of success evolves around it, and it is difficult for other areas of success to ’swallow’ it. Assuming the contradiction is not too strong then, eventually, after numerous attempts, the complete graph is marked with ’X’s and the Graph evolves successfully. At this point the contents of each qubit variable can be amplified, and final measurement of all qubits provides the ML codeword with high probability. The advantage of a free-running strategy, where each function node is free to activate asynchronously, is that regions of general agreement develop first and influence other areas of the graph to ’follow their opinion’. Fig 13 also shows that one ’bad’ (contradictory) qubit can be a fatal stumbling block to successful evolution of the whole graph (as it can for SPA on classical graphs). Thus Distributed QPA requires Fault-Tolerance, where only an arbitrary subset of entangled nodes are required as a final result (node redundancy). The free-running schedule of Fig 13 naturally avoids the ’bad’ qubits, but sufficient evolution occurs when enough function nodes complete. Alternatively, bad qubits could be set to $`(\sqrt{0.5},\sqrt{0.5})`$ after a time-out. A more detailed proposal of Fault-Tolerant QPA is left for future work.

Fig 13 also serves to illustrate the ’template’ for a Reconfigurable Quantum Graph Array. One can envisage initialising an array of quantum variables so that two local variables can be strongly or weakly entangled by identifying the mutual square function nodes with strongly or weakly-entangling matrices, respectively. In particular, two neighbouring nodes may be ’locally disconnected’ by setting the function node joining them to a tensor-decomposable matrix, (i.e. zero-entangling). The quantum computer is then measurement-driven. The concept of measurement-driven quantum computation has also recently been pursued in , where a uniform entanglement is set-up throughout the array 777It is interesting that this entanglement is strongly related to Rudin-Shapiro and quadratic constructions  prior to computation via measurement.

Fig 14 shows the system view of QPA. A continual stream of pure qubits needs to be initialised and then entangled, and then amplified, so as to ensure at least one successful entangled and amplified output from the whole apparatus.

## 6 Phase QPA

The above discussions have ignored the capacity of Quantum Systems to carry phase information. In fact QPA, as presented so far, is immune to phase modification, as classical probabilities have no phase component. However QPA should be generalised to cope with phase shift in order to decode quantum information. This is the subject of ongoing research.

## 7 Conclusion and Discussion

The Quantum Product Algorithm (QPA) on a Factor Graph has been presented for Maximum-Likelihood (ML) Decoding of Classical ’soft’ information using quantum resources. The relationship of QPA to the Sum-Product Algorithm (SPA) has been indicated, where avoidance of summary allows QPA to overcome small graph cycles. Quantum Factor Graphs use small unitary matrices which each act on only a few qubits. QPA is measurement-driven and is only statistically likely to succeed after many attempts. The ML codeword is obtained with maximum likelihood by measuring the entangled vector resulting from successful QPA. To ensure a high probability of measuring the ML codeword QPA output can be amplified prior to measurement. The complete ML decoder is only successful after many attempts. Finally, free-running Distributed QPA is proposed to improve the likelihood of successful QPA completion. The free-running distributed structure suggests further benefit will be obtained by introducing Fault-Tolerance in the form of redundant function and variable nodes. Phase aspects of QPA have yet to be explored. This paper has been written to demonstrate the exponential capacity of quantum systems, and their natural suitability for graph decompositions such as the Factor Graph. The paper has not tried to deal with quantum noise and quantum decoherence, but one can expect the Factor Graph form to ’gracefully’ expand to cope with the extra redundancy necessary to protect qubits from decoherence and noise. When viewed in the context of entangled space, it is surprising how successful classical message-passing algorithms are, even though they are restricted to operate in tensor product space. This suggests that methods to improve the likelihood of successful QPA completion may include the possibility of hybrid QPA/SPA graphs, where SPA operates on non-cyclic and resolvable parts of the graph, leaving QPA to cope with small cycles or unresolved areas of the graph.

## 8 Appendix A - Deriving $`p_m(t)`$ and $`P(t)`$ for Fig 8

The probability of successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, is $`p_{f012}=(\alpha _0\alpha _1\alpha _2)^2+(\beta _0\beta _1\beta _2)^2`$, and similarly for $`p_{f345}`$ and $`p_{f678}`$. Let $`h_{012}=(1p_{f012})^{q1}`$, $`h_{345}=(1p_{f345})^{q1}`$, $`h_{678}=(1p_{f678})^{q1}`$. Then the probability of successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$ after exactly $`q`$ parallel attempts is,

$$\begin{array}{c}p_{036}\left(q\right)=h_{012}h_{345}h_{678}p_{f012}p_{f345}p_{f678}+\left(1h_{012}\right)h_{345}h_{678}p_{f345}p_{f678}\\ +h_{012}\left(1h_{345}\right)h_{678}p_{f012}p_{f678}+h_{012}h_{345}\left(1h_{678}\right)p_{f012}p_{f345}\\ +\left(1h_{012}\right)\left(1h_{345}\right)h_{678}p_{f678}+\left(1h_{012}\right)h_{345}\left(1h_{678}\right)p_{f345}\\ +h_{012}\left(1h_{345}\right)\left(1h_{678}\right)p_{f012}\end{array}$$

Given successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$, the probability of subsequent successful completion of $`𝐔_{\mathrm{𝐟𝟐𝟓𝟖}}`$ is,

$$p_{f258}^{}=\frac{(\alpha _0\alpha _1\alpha _2\alpha _3\alpha _4\alpha _5\alpha _6\alpha _7\alpha _8)^2+(\beta _0\beta _1\beta _2\beta _3\beta _4\beta _5\beta _6\beta _7\beta _8)^2}{p_{f012}p_{f345}p_{678}}$$

Therefore the probability of successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$, immediately followed by successful completion of $`𝐔_{\mathrm{𝐟𝟐𝟓𝟖}}`$ is, $`p_{08}(q)=p_{036}(q1)p_{f258}^{}`$, and the probability of successful completion of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$, immediately followed by completion failure of $`𝐔_{\mathrm{𝐟𝟐𝟓𝟖}}`$ is, $`\overline{p_{08}}(q)=p_{036}(q1)(1p_{f258}^{})`$. The probability of successful completion after exactly $`t`$ steps of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$ in parallel, followed by $`𝐔_{\mathrm{𝐟𝟐𝟓𝟖}}`$, is,

$$p_e(t)=\underset{q=2}{\overset{t}{}}p_{08}(q)\underset{𝐯𝐃(𝐭𝐪)}{}\underset{u𝐯}{}\overline{p_{08}}(u)$$

where $`𝐃(𝐤)`$ is the set of unordered partitions of $`k`$. Therefore the probability of successful completion after at most $`t`$ steps of $`𝐔_{\mathrm{𝐟𝟎𝟏𝟐}}`$, $`𝐔_{\mathrm{𝐟𝟑𝟒𝟓}}`$, and $`𝐔_{\mathrm{𝐟𝟔𝟕𝟖}}`$ in parallel, followed by $`𝐔_{\mathrm{𝐟𝟐𝟓𝟖}}`$, is,

$$p_m(t)=\underset{i=2}{\overset{t}{}}p_e(i)$$

In contrast, the probability of successful completion, after at most $`t`$ steps, of a non-distributed version of Fig 8 is $`P(t)=1(1(\alpha _0\alpha _1\alpha _2\alpha _3\alpha _4\alpha _5\alpha _6\alpha _7\alpha _8)^2(\beta _0\beta _1\beta _2\beta _3\beta _4\beta _5\beta _6\beta _7\beta _8)^2)^t`$.