text:
# Spin Fluctuation Induced Superconductivity Controlled by Orbital Fluctuation

\[

## Abstract

A microscopic Hamiltonian reflecting the correct symmetry of $`f`$-orbitals is proposed to discuss superconductivity in heavy fermion systems. In the orbitally degenerate region in which not only spin fluctuations but also orbital fluctuations develop considerably, cancellation between spin and orbital fluctuations destabilizes $`d_{x^2y^2}`$-wave superconductivity. Entering the non-degenerate region by increasing the crystalline electric field, $`d_{x^2y^2}`$-wave superconductivity mediated by antiferromagnetic spin fluctuations emerges out of the suppression of orbital fluctuations. We argue that the present scenario can be applied to recently discovered superconductors CeTIn5 (T=Ir, Rh, and Co).

\]

Unconventional superconductivity has been one of central issues in the research field of strongly correlated electron systems. Especially since the discovery of high temperature superconductivity in cuprates, much effort has been focussed on elucidating the mechanism of unconventional superconductivity, clarifying that a crucial role is played by antiferromagnetic (AF) “spin fluctuations” (SF) . The importance of AFSF is widely recognized, for instance, in $`4f`$\- and $`5f`$-electron superconducting materials such as $`\mathrm{CeCu}_2\mathrm{Si}_2`$  and $`\mathrm{UPd}_2\mathrm{Al}_3`$  as well as in organic superconductors such as $`\kappa `$-(BEDT-TTF) . Thus, it is widely believed that a broad class of unconventional superconductors originates from SF.

In $`d`$\- and $`f`$-electron systems, however, the potential importance of orbital degrees of freedom has recently been discussed intensively. In fact, “orbital ordering” is found to be a key issue for understanding microscopic aspects of the charge-ordered phase in colossal magnetoresistive manganites . This orbital ordering is primarily relevant to the insulating phase, while “orbital fluctuations” (OF) should be significant in the metallic phase. Recently, the effects of OF have attracted attention, since it is hoped to provide a new scenario for superconductivity . Especially from a conceptual viewpoint, it is important to clarify how superconductivity emerges when both SF and OF play active roles.

As a typical material for investigating superconductivity in a system with both SF and OF, let us introduce the recently discovered heavy-fermion superconductors $`\mathrm{CeTIn}_5`$ (T=Ir, Rh, and Co)  with the $`\mathrm{HoCoGa}_5`$-type tetragonal crystal structure. Due to this structure and to strong correlation effects, the similarity with cuprates has been emphasized. In particular, AFSF also plays an essential role in the Ce-115 system, since it exhibits quasi two-dimensional Fermi surfaces  and the AF phase exists next to the superconducting state in the phase diagram of $`\mathrm{CeRh}_{1\mathrm{x}}\mathrm{Ir}_\mathrm{x}\mathrm{In}_5`$ . In fact, a line node in the gap function has been observed in CeTIn5 by various experimental techniques . It may be true that superconductivity itself is due to AFSF, but an important role for the OF has been overlooked in spite of the fact that both SF and OF are originally included in the ground-state multiplet of Ce3+ ion. In actual materials, superconductivity occurs in a situation where OF is suppressed, since orbital degeneracy is lifted by the effect of the crystalline electric field (CEF). Thus, we envision a scenario where OF controls the stability of the superconductivity even though it originates from AFSF.

In this Letter, we investigate superconductivity based on the orbitally degenerate Hubbard model constructed by the tight-binding method . Solving the gap equation with the pairing interaction evaluated using the random phase approximation (RPA), we obtain several superconducting phases around the spin and orbital ordered phases. In the orbitally degenerate region where both SF and OF are developed, it is found that singlet superconductivity is suppressed due to the competition between them, while triplet superconductivity is favored since they are cooperative in this case. When a level splitting is included to lift the orbital degeneracy, $`d_{x^2y^2}`$-wave superconductivity due to AFSF is stabilized in the vicinity of the AF phase. Thus, we claim that AFSF-induced supercondutivity in Ce-115 systems is substantiated in consequence of the suppression of OF.

In order to construct the microscopic model for $`f`$-electron systems, let us start our discussion from a local basis of the Ce3+ ion. Among the 14-fold degenerate $`4f`$-electron states, due to the effect of strong spin-orbit coupling, only the $`j`$=5/2 sextuplet effectively contributes to the low-energy excitations ($`j`$ is total angular momentum). This sextuplet is further split into a $`\mathrm{\Gamma }_7`$ doublet and a $`\mathrm{\Gamma }_8`$ quadruplet due to the effect of cubic CEF, where the eigen states are given by $`|\mathrm{\Gamma }_{7\pm }`$= $`\sqrt{1/6}|\pm 5/2`$$``$$`\sqrt{5/6}|3/2`$, $`|\mathrm{\Gamma }_{8\pm }^{(1)}`$= $`\sqrt{5/6}|\pm 5/2`$+$`\sqrt{1/6}|3/2`$, and $`|\mathrm{\Gamma }_{8\pm }^{(2)}`$=$`|\pm 1/2`$. Here $`+`$ and $``$ in the subscripts denote “pseudo-spin” up and down, respectively, within each Kramers doublet.

Now we discuss the relative positions of the energy levels of $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_8^{(\tau )}`$ while taking account of certain features of CeTIn5. Since CeTIn5 has a tetragonal crystal structure and a quasi two-dimensional Fermi surface , it is natural to consider a two-dimensional square lattice composed of Ce ions. Due to the effect of anions surrounding the Ce ion, it is deduced that the energy level of $`\mathrm{\Gamma }_7`$ becomes higher than those of $`\mathrm{\Gamma }_8`$’s. Thus, in the following, the $`\mathrm{\Gamma }_7`$ orbital is neglected for simplicity. Note also that we need to include the effect of the tetragonal CEF, which lifts the degeneracy of the $`\mathrm{\Gamma }_8^{(\tau )}`$. Although a mixing between $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_8^{(1)}`$ generally occurs under a tetragonal CEF, such a mixing is expected to be small since the $`\mathrm{\Gamma }_7`$ orbital has higher energy. Thus, in this situation, the effect of a tetragonal CEF can be included in terms of a level splitting $`\epsilon `$ between the $`\mathrm{\Gamma }_8^{(1)}`$ and $`\mathrm{\Gamma }_8^{(2)}`$ orbitals. In order to explain the magnetic anisotropy observed in experiments, $`\mathrm{\Gamma }_8^{(1)}`$ must be lower than $`\mathrm{\Gamma }_8^{(2)}`$, i.e., $`\epsilon `$$`>`$0 in CeTIn5. The magnitude of $`\epsilon `$ for each CeTIn5 compound will be discussed later. It should be noted that the above level scheme is consistent with the following two facts: (i) Experimental results for CeTIn5 exhibit a larger uniform susceptibility for magnetic field perpendicular to the CeIn3 plane than for the parallel case . This significant anisotropy is well explained under the assumption that $`\mathrm{\Gamma }_7`$ is $`not`$ the lowest-energy state and $`\epsilon `$ is $`positive`$. (ii) The band-structure calculation results suggest that the almost flat band corresponding to $`\mathrm{\Gamma }_7`$ appears above the Fermi level .

In order to include the itinerant features of the $`4f`$-electrons, a simple way is to include nearest-neighbor hopping for the $`f`$-electrons by the tight-binding method . Although the hybridization with the In $`5p`$ electronic states may be important, here such an effect is considered as renormalization of the effective hopping amplitude of $`f`$-quasiparticles. Further, by adding the on-site Coulomb interaction terms among the $`f`$-electons, the Hamiltonian becomes

$`H`$ $`=`$ $`{\displaystyle \underset{\mathrm{𝐢𝐚}\tau \tau ^{}\sigma }{}}t_{\tau \tau ^{}}^𝐚f_{𝐢\tau \sigma }^{}f_{𝐢+𝐚\tau ^{}\sigma }\epsilon {\displaystyle \underset{𝐢}{}}(n_{𝐢1\sigma }n_{𝐢2\sigma })/2`$ (1)
$`+`$ $`U{\displaystyle \underset{𝐢\tau }{}}n_{𝐢\tau }n_{𝐢\tau }+U^{}{\displaystyle \underset{𝐢\sigma \sigma ^{}}{}}n_{𝐢1\sigma }n_{𝐢2\sigma ^{}},`$ (2)
where $`f_{𝐢\tau \sigma }`$ is the annihilation operator for an $`f`$-electron with pseudo-spin $`\sigma `$ in the orbital $`\mathrm{\Gamma }_8^{(\tau )}`$ at site $`𝐢`$, $`𝐚`$ is the vector connecting nearest-neighbor sites, and $`n_{𝐢\tau \sigma }`$=$`f_{𝐢\tau \sigma }^{}f_{𝐢\tau \sigma }`$. The first term represents the nearest-neighbor hopping of $`f`$-electrons with the amplitude $`t_{\tau \tau ^{}}^𝐚`$ between $`\mathrm{\Gamma }_8^{(\tau )}`$ and $`\mathrm{\Gamma }_8^{(\tau ^{})}`$ along the $`𝐚`$-direction, given by $`t_{11}^𝐱`$=$`\sqrt{3}t_{12}^𝐱`$= $`\sqrt{3}t_{21}^𝐱`$=$`3t_{22}^𝐱`$=1 for $`𝐚`$=$`𝐱`$ and $`t_{11}^𝐲`$=$`\sqrt{3}t_{12}^𝐲`$= $`\sqrt{3}t_{21}^𝐲`$=$`3t_{22}^𝐲`$=1 for $`𝐚`$=$`𝐲`$, respectively, in energy units where $`t_{11}^𝐱`$=1. Note the positive sign of the first term in $`H`$, since the M-point, not the $`\mathrm{\Gamma }`$-point, is at the bottom of the bands forming the Fermi surfaces . Note also that the present $`t_{\tau \tau ^{}}^𝐚`$ is just the same as that of the $`e_\mathrm{g}`$ electrons , but this point will be discussed elsewhere. The second term denotes the tetragonal CEF. In the third and fourth terms, $`U`$ and $`U^{}`$ are the intra- and inter-orbital Coulomb interations, respectively. In reality, $`U`$=$`U^{}`$, since they originate from the same Coulomb interactions among $`f`$-orbitals in the $`j`$=5/2 multiplet, but in this paper, we also treat the case where $`U`$$``$$`U^{}`$ in order to analyze the roles of SF and OF. Since we consider quarter-filling (one $`f`$-elecron per site), the model in the limit of $`\epsilon `$=$`\mathrm{}`$ reduce to the half-filled, single-orbital Hubbard model.

Now, in order to investigate superconductivity around the spin and/or orbital ordered phases, we calculate spin and orbital susceptibilities, $`\widehat{\chi }^\mathrm{s}(𝐪)`$ and $`\widehat{\chi }^\mathrm{o}(𝐪)`$, respectively. Within the RPA, these are given in a matrix form as

$`\widehat{\chi }^\mathrm{s}(𝐪)=[\widehat{1}\widehat{U}^\mathrm{s}\widehat{\chi }(𝐪)]^1\widehat{\chi }(𝐪),`$ (3)
$`\widehat{\chi }^\mathrm{o}(𝐪)=[\widehat{1}+\widehat{U}^\mathrm{o}\widehat{\chi }(𝐪)]^1\widehat{\chi }(𝐪),`$ (4)
where labels of row and column in the matrix appear in the order 11, 22, 12, and 21, these being pairs of orbital indices 1 and 2. Note that $`\widehat{1}`$ is the 4$`\times `$4 unit matrix. $`\widehat{U}^\mathrm{s}`$ is given by $`U_{11,11}^\mathrm{s}`$=$`U_{22,22}^\mathrm{s}`$=$`U`$, $`U_{12,12}^\mathrm{s}`$=$`U_{21,21}^\mathrm{s}`$=$`U^{}`$, otherwise zero, while $`\widehat{U}^\mathrm{o}`$ is expressed as $`U_{11,11}^\mathrm{o}`$=$`U_{22,22}^\mathrm{o}`$=$`U`$, $`U_{11,22}^\mathrm{o}`$=$`U_{22,11}^\mathrm{o}`$=$`2U^{}`$, $`U_{12,12}^\mathrm{o}`$=$`U_{21,21}^\mathrm{o}`$=$`U^{}`$, otherwise zero. The matrix elements of $`\widehat{\chi }(𝐪)`$ are defined by $`\chi _{\mu \nu ,\alpha \beta }(𝐪)`$= $`T_{𝐤,n}G_{\alpha \mu }^{(0)}(𝐤+𝐪,i\omega _n)G_{\nu \beta }^{(0)}(𝐤,i\omega _n)`$, where $`T`$ is temperature, $`G_{\mu \nu }^{(0)}(𝐤,i\omega _n)`$ is the non-interacting Green’s function for $`f`$-electrons with momentum $`𝐤`$ propagating between $`\mu `$\- and $`\nu `$-orbitals, and $`\omega _n`$=$`\pi T(2n`$$`+`$$`1)`$ with integer $`n`$. The instabilities for the spin- and orbital-ordered phases are determined by the conditions of $`\mathrm{det}[\widehat{1}\widehat{U}^\mathrm{s}\widehat{\chi }(𝐪)]`$=0 and $`\mathrm{det}[\widehat{1}+\widehat{U}^\mathrm{o}\widehat{\chi }(𝐪)]`$=0, respectively.

By using $`\widehat{\chi }^\mathrm{s}(𝐪)`$ and $`\widehat{\chi }^\mathrm{o}(𝐪)`$, the superconducting gap equation is given as

$$𝚫^\xi (𝐤)=\underset{𝐤^{}}{}\widehat{V}^\xi (𝐤𝐤^{})\widehat{\varphi }(𝐤^{})𝚫^\xi (𝐤^{}),$$
(5)
where $`𝚫^\xi (𝐤)`$ =\[$`\mathrm{\Delta }_{11}^\xi (𝐤)`$, $`\mathrm{\Delta }_{22}^\xi (𝐤)`$, $`\mathrm{\Delta }_{12}^\xi (𝐤)`$, $`\mathrm{\Delta }_{21}^\xi (𝐤`$)\]t is the gap function in the vector representation for a singlet ($`\xi =\mathrm{S}`$) or triplet ($`\xi =\mathrm{T}`$) pairing state, and the matrix elements of the singlet- and triplet-pairing potentials, respectively, are given by

$`V_{\alpha \beta ,\mu \nu }^\mathrm{S}(𝐪)=[(3/2)\widehat{W}^\mathrm{s}(𝐪)+(1/2)\widehat{W}^\mathrm{o}(𝐪)+\widehat{U}^\mathrm{s}]_{\alpha \mu ,\nu \beta },`$ (6)
$`V_{\alpha \beta ,\mu \nu }^\mathrm{T}(𝐪)=[(1/2)\widehat{W}^\mathrm{s}(𝐪)+(1/2)\widehat{W}^\mathrm{o}(𝐪)\widehat{U}^\mathrm{s}]_{\alpha \mu ,\nu \beta }.`$ (7)
Here, the spin and orbital susceptibilities are included as $`\widehat{W}^\mathrm{s}(𝐪)`$=$`\widehat{U}^\mathrm{s}+\widehat{U}^\mathrm{s}\widehat{\chi }^\mathrm{s}(𝐪)\widehat{U}^\mathrm{s}`$ and $`\widehat{W}^\mathrm{o}(𝐪)`$=$`\widehat{U}^\mathrm{o}+\widehat{U}^\mathrm{o}\widehat{\chi }^\mathrm{o}(𝐪)\widehat{U}^\mathrm{o}`$. The elements of the pair correlation function $`\widehat{\varphi }(𝐤)`$ are given by $`\varphi _{\alpha \beta ,\mu \nu }(𝐤)`$=$`T_nG_{\alpha \mu }^{(0)}(𝐤,i\omega _n)G_{\nu \beta }^{(0)}(𝐤,i\omega _n)`$. The superconducting transition is obtained for each respective irreducible representation by solving Eq. (5), where the maximum eigenvalue becomes unity.

Here we mention the essential symmetries of the Cooper pairs in multi-orbital systems. The pairing states are classified into four types owing to spin $`SU(2)`$ symmetry and space inversion symmetry; (1) spin-singlet and orbital-symmetric with even-parity, (2) spin-triplet and orbital-symmetric with odd-parity, (3) spin-singlet and orbital-antisymmetric with odd-parity, and (4) spin-triplet and orbital-antisymmetric with even-parity. Note that $`SU(2)`$ symmetry does not exist in orbital space, since the lattice and the local wavefunctions rotate simultaneously. For (1) and (2), $`f`$-electrons in the same band form the Cooper pair, while for (3) and (4), such a pair is formed only between different bands. Except for the special case in which two Fermi surfaces connect with each other, orbital-antisymmetric pairing is unstable due to depairing effects destroy inter-band pairs. In fact, even after careful calculations, we do not find any region for (3) and (4) in the parameter space considered here. Thus, in the following, we discuss only the orbital-symmetric pair.

Let us discuss first the effect of $`\epsilon `$ on the correlation between SF and OF. Since we are considering a realistic situation corresponding to CeTIn5, we restrict ourselves to the case of $`U`$=$`U^{}`$ for a while. In Fig. 1, the principal components of $`\widehat{\chi }^\mathrm{s}(𝐪)`$ and $`\widehat{\chi }^\mathrm{o}(𝐪)`$ are shown in $`𝐪`$ space for $`U`$=$`U^{}`$=0.9$`U_\mathrm{m}`$, where $`U_\mathrm{m}`$ denotes spin instability determined from $`\mathrm{det}[\widehat{1}\widehat{U}^\mathrm{s}\widehat{\chi }(𝐪)]`$=0. In the orbital-degenerate case with $`\epsilon `$=0, as shown in Figs. 1(a) and (b), the overall magnitude of OF is comparable with that of SF. As for the $`𝐪`$ dependence, common structures are found in $`\chi ^\mathrm{s}(𝐪)`$ and $`\chi ^\mathrm{o}(𝐪)`$. Namely, there are two peaks around $`(\pi /2,\pi /2)`$ and $`(\pi /2,0)`$ owing to the nesting properties of the Fermi surface (see the inset). Thus, the orbital-degenerate region is characterzied by competition between SF and OF. With increasing $`\epsilon `$, the Fermi surface approaches a shape having the nesting vector $`(\pi ,\pi )`$ (see the inset). For $`\epsilon `$=4, as shown in Fig. 1(c) and (d), OF is almost completely suppressed, and SF around $`(\pi ,\pi )`$ become dominant. Note here that the increase of $`\epsilon `$ makes the lower energy state favorable and suppresses excitations to the upper energy state, indicating the suppression of OF. Thus, even at this stage, it is understood that the suppression of OF leads to the development of SF for the paramagnetic system which is wavering between spin- and orbital-ordering.

Next, we consider how a superconducting phase emerges when $`\epsilon `$ is increased. In Fig. 2(a), the maximum eigenvalues for several irreducible representations are depicted as a function of $`\epsilon `$ for $`U`$=$`U^{}`$=2.5. The calculations are carried out for a fixed $`T`$=0.02, and the first Brillouin zone is divided into $`128\times 128`$ meshes. In the orbital-degenerate region, several eigenvalues are very close to each other owing to multi-peak structures in $`\widehat{\chi }^\mathrm{s}(𝐪)`$ and $`\widehat{\chi }^\mathrm{o}(𝐪)`$. With increasing $`\epsilon `$, as is easily understood from the growth of AFSF mentioned above, the eigenvalue for $`B_{1\mathrm{g}}`$ symmetry becomes dominant and finally at $`\epsilon `$$``$3, the $`B_{1\mathrm{g}}`$ superconducting phase is stabilized. On further increasing $`\epsilon `$, the AF instability eventually occurs, since the system asymptotically approaches the half-filled single-orbital Hubbard model. We emphasize that the increase of $`\epsilon `$ brings transitions successively in the order of PM, $`B_{1\mathrm{g}}`$ superconducting, and AF phases.

Based on the above calculated results, let us try to explain the differences among three Ce-115 compounds, CeRhIn5 (Néel temperature $`T_\mathrm{N}`$=3.8K), CeIrIn5 (superconducting transition temperature $`T_\mathrm{c}`$=0.4K), and CeCoIn5 ($`T_\mathrm{c}`$=2.3K) . One property which distinguishes these compounds is two-dimensionality, as expressed by the ratio $`c/a`$ between lattice constants. It is naturally expected that two-dimensionality becomes stronger in the order of CeIrIn5, CeRhIn5, and CeCoIn5 . Since the magnetic state should be stabilized with increasing three-dimensionality, CeIrIn5 should be most favorable for the occurence of antiferromagnetism among three compounds, but this is obviously inconsistent with experimental results. This inconsistency is resolved by introducing another important ingredient $`\epsilon `$. One can show that the increase of magnetic anisotropy just corresponds to the increase of $`\epsilon `$. Analysis of experimental results for the anisotropy of magnetic susceptibilities  leads to the conclusion that $`\epsilon `$ becomes larger in the order of CeCoIn5, CeIrIn5, and CeRhIn5. Taking into account the effect of $`\epsilon `$ in addition to dimensionality, we arrive at the picture schematically shown in Fig. 2(b). Due to the enhancement of AFSF induced by increasing $`\epsilon `$, as shown in Fig. 1, it may be understood that CeRhIn5 rather than CeIrIn5 is more favorable to antiferromagnetism. Namely, it is considered that CeRhIn5 with the largest $`\epsilon `$ is antiferromagnet, while CeCoIn5 and CeIrIn5 with smaller $`\epsilon `$ than CeRhIn5 exhibit superconductivity. The difference in $`T_\mathrm{c}`$ between CeCoIn5 and CeIrIn5 may be attributed to the extent of two-dimensionality. Concerning the superconductivity and antiferromagnetism of CeTIn5 compounds, the combination of two dimensionality and the crystalline field splitting is necessary to reach a consistent picture.

The orbitally degenerate model, Eq. (1), shows new, rich superconducting properties both for singlet and triplet pairings. In Fig. 3(a) and (b), we show the phase diagrams in the $`U`$-$`U^{}`$ plane for $`\epsilon `$=0 and $`U`$-$`\epsilon `$ plane for $`U`$=$`U^{}`$, respectively. In these figures, four characteristic superconducting phases are observed: (1) a SF-mediated spin-singlet superconducting phase with $`B_{1\mathrm{g}}`$-symmetry in the vicinity of a spin-ordered phase; (2) a $`g`$-wave superconducting phase with $`A_{2\mathrm{g}}`$-symmetry; (3) a spin-triplet superconducting phase with $`E_\mathrm{u}`$-symmetry due to the cooperation between SF and OF around $`U`$=$`U^{}`$ and $`\epsilon `$=0 and (4) an OF-mediated spin-singlet superconducting phase with $`A_{1\mathrm{g}}`$-symmetry for $`U`$$``$$`U^{}`$ around an orbital-ordered phase. With respect to (1), it is thought that $`B_{1\mathrm{g}}`$-superconductivity is induced by SF around $`𝐪`$=$`𝐐`$=$`(\pi ,\pi )`$. Concerning the triplet superconductivity (3), note that for this region all superconducting instabilities are quite close to each other because of the multi-peak structures in $`\widehat{\chi }^\mathrm{s}(𝐪)`$ and $`\widehat{\chi }^\mathrm{o}(𝐪)`$ as shown in Fig. 1. We also note that the factors in front of the SF term in $`\widehat{W}^\mathrm{s}(𝐪)`$ are $`3/2`$ for singlet and $`+1/2`$ for triplet pairing, while the factor for the OF term in $`\widehat{W}^\mathrm{o}(𝐪)`$ is $`+1/2`$ for both pairings. Namely, OF is cooperative with SF for spin-triplet pairing, while SF and OF compete with each other for spin-singlet pairing . Thus, the spin-triplet pairing phase appears in the region for $`\epsilon `$$``$0 and $`U`$$``$$`U^{}`$. The mechanism of the OF-mediated singlet supercondutivity (4) is interesting. In this state, two quasi-particles with different pseudo-spin form an on-site intra-orbital singlet pair since $`U^{}U`$, leading to $`A_{1\mathrm{g}}`$ superconductivity. In short, singlet superconductivity with $`B_{1\mathrm{g}}`$\- or $`A_{1\mathrm{g}}`$-symmetry is stabilized when either SF or OF is dominant, while triplet superconductivity is favored when there is cooperation between SF and OF.

In summary, we have studied superconductivity in the Ce-115 systems based on the orbitally degenerate Hubbard model, and found that with increasing the tetragonal CEF splitting $`d_{x^2y^2}`$-wave superconductivity due to AFSF emerges out of the suppression of OF in the vicinity of the AF phase. The concept of AFSF-induced superconductivity controlled by OF qualitatively explains the difference among Ce-115 materials.

The authors thank H. Harima, A. Hasegawa, T. Moriya, Y. Ōnuki, R. E. Walstedt, and K. Yamada for discussions. K.U. is financially supported by Grant-in-Aid for Scientific Research Areas (B) from the Ministry of Education, Science, Sports, Culture, and Technology.


format: text

text:
# The relation between the mid-infrared emission and Black Hole Mass in Active Galactic Nuclei: A direct way to probe black hole growth?

## 1. Introduction

The cosmological evolution of the black hole (BH) population is a fundamental property to constrain models of galaxy formation and evolution. In the last few years we have seen significant progress, especially for local galaxies where BH masses are now being measured. More importantly, the discovery of two fundamental correlations between the BH mass and the bulge luminosity (Kormendy & Richstone 1995; Magorrian et al. 1998) and the BH mass and the stellar velocity dispersion (Gebhardt et al. 2000; Ferrarese & Merritt 2000) suggests a tight link in the evolution of the BH and its galactic host.

To produce reliable statistics of BH masses at high redshift, one would need unrealistic amounts of observing time if the methods employed to determine BH masses in the Local Universe (see Section 2) were to be used. Thus, understanding the cosmological evolution of BHs requires finding a good indicator of the active galactic nuclei (AGN) luminosity, independent from the effects of obscuration and the viewing angle, and directly related to the BH mass.

Among the suggested AGN power indicators are the infrared luminosity, the narrow line region \[O iii\]$`\lambda 5007`$ line emission, the hard X-ray ($`210`$keV) emission, and the radio emission (e.g., Mulchaey et al. 1994). The $`210`$keV hard X-ray emission is only a good indicator of the intrinsic luminosity of the AGN for those cases where it is transmitted through the torus, that is, in Compton thin galaxies.

QSOs emit a significant fraction of their bolometric luminosity in the infrared (e.g., Sanders et al. 1989). The mid- and far-infrared emission in Sys and radio quiet quasars is predominantly thermal in origin, namely dust emission (e.g., Rieke 1978; Barvainis 1987; McAlary & Rieke 1988; Sanders et al. 1989; and recently Haas et al. 2000; Polletta et al. 2000). Moreover, Spinoglio, & Malkan (1989) found that the $`12\mu `$m flux is approximately a constant fraction of the bolometric flux in Seyfert (Sy) galaxies and quasars. These considerations suggest that the main parameter governing the mid-infrared behavior is the power of the central engine, almost independently of the geometry, dust content and dust properties of the central region.

For type 2 objects, the infrared emission at $`\lambda 5\mu `$m may still be significantly affected by viewing angle effects (obscuration) and/or stellar contribution (Alonso-Herrero et al. 2001 and references therein). At $`\lambda 20\mu `$m there may be significant contributions from star formation and/or the underlying galaxy. These contributions are not easily disentangled from the AGN emission unless high spatial resolution imaging is used. Thus the nuclear $`10\mu `$m luminosity appears as a good choice to represent the AGN bolometric luminosity.

In this letter we explore for the first time the relation between the BH mass and the thermal $`10\mu `$m nuclear emission for a heterogeneous sample of low-redshift AGN: Sy 1s, Sy 2s, and PG quasars, and explore the possibility of using the mid-infrared luminosity of AGN as a representation of their BH mass. We also discuss implications for the BH growth and demographics over cosmological times.

## 2. The Data and Correlations

We have compiled BH masses of local AGNs – Sy galaxies and PG quasars – derived using three methods and for which nuclear $`10\mu `$m fluxes are available (see Table 1 for the sample and references for the data). The most direct way is from spatially resolved stellar kinematics (see Ho 1999 for a review), and measurements of masers in disks around massive BHs (see Moran, Greenhill, & Herrnstein 1999 for a review).

In the case of bright active galaxies where the stellar emission is swamped by the AGN emission, the best technique to derive BH masses is reverberation mapping (see a review by Peterson 1993). For sample of AGNs analyzed in Kaspi et al. (2000) we have used the BH masses (and errors) derived from the mean of the FWHM velocity measurements.

The third and most indirect method uses the empirical relation between the BH mass and the stellar velocity dispersion (the $`M_{\mathrm{BH}}`$$`\sigma _{}`$ relation) found by Ferrarese & Merritt (2000): $`\mathrm{log}M_{\mathrm{BH}}=4.80(\pm 0.54)\mathrm{log}\sigma _{}2.9(\pm 1.3)`$. The $`M_{\mathrm{BH}}`$ errors here are dominated by the empirical correlation uncertainty rather than by the observational errors in $`\sigma _{}`$. In Figure 1 we only plot the errors associated with the slope of the relation, that is, $`\mathrm{\Delta }\mathrm{log}M_{\mathrm{BH}}=0.54\times \mathrm{log}\sigma _{}`$.

We also explore the $`210`$ keV hard X-ray luminosity vs. BH mass relation for Compton thin AGNs based upon the good correlation found between the hard X-ray and the $`10\mu `$m luminosities in Sys and PG quasars (Alonso-Herrero et al. 2001, Krabbe et al. 2001). The $`210`$keV hard X-ray luminosities (corrected for absorption) are from Nandra & Pounds (1994), Lawson & Turner (1997), Bassani et al. (1999), George et al. (2000), Iwasawa et al. (2000), and the compilation in Alonso-Herrero et al. (1997). We have included a few low luminosity AGNs with available hard X-ray fluxes corrected for absorption (although mid-infrared fluxes were not available). These are NGC 3031, NGC 4594, NGC 6251 ($`M_{\mathrm{BH}}`$ from the compilation of Ho 2002), NGC 3079 ($`M_{\mathrm{BH}}`$ from Moran et al. 1999) and NGC 5194 ($`M_{\mathrm{BH}}`$ derived from the $`M_{\mathrm{BH}}`$$`\sigma _{}`$ relation, where $`\sigma _{}`$ is from Nelson & Whittle 1995).

When available we have used distances derived from surface brightness fluctuations (Tonry et al. 2001). For the rest of the galaxies distances were computed using $`H_0=75\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

## 3. Discussion and Conclusions

The relations between the BH mass and the $`10\mu `$m nuclear luminosity and $`210`$keV luminosity for our AGN sample are presented in Figure 1 (left and right panel respectively). The latter relation has been previously analyzed by Wandel & Mushotzky (1986) and Awaki et al. (2001) for a limited number of AGNs. Ho (2002) and Wu & Han (2001) have investigated the relationship between the radio nuclear and total power and the BH mass and found a good correlation for bright AGNs. Ho (2002) however finds that weakly active galaxies do not follow the same relation.

In a series of papers, McLeod & Rieke (1995), McLeod, Rieke, & Storrie-Lombardi (1999) and McLeod & McLeod (2001) have demonstrated that there is a relation between the nuclear $`B`$ magnitude of quasars and the near-infrared luminosity of the host galaxy (that is, the stellar mass of the galaxy). This relation is also present in luminous infrared galaxies when the far-infrared luminosity is used to represent the bolometric luminosity (McLeod et al. 1999). McLeod & McLeod (2001) interpreted these relations as the maximum nuclear luminosity possible in a host galaxy of a given near-infrared luminosity (mass). In view of this relation and the good correlation between the BH mass and the bulge luminosity of the host galaxy it would not be unreasonable to expect a relation between the BH mass and the AGN luminosity represented by the $`10\mu `$m or the hard X-ray luminosities, the latter only for Compton thin objects.

In the hard X-ray vs. $`M_{\mathrm{BH}}`$ relation the lines ($`\eta _{\mathrm{Edd}}=0.01`$ and $`\eta _{\mathrm{Edd}}=1`$) are the Eddington ratios computed in Awaki et al. (2001, their figure 3) assuming a bolometric correction of $`\frac{L_{\mathrm{bol}}}{L_\mathrm{X}}=27.2`$ (note that these authors pointed out that such bolometric correction may be lower for low luminosity AGNs). All the objects located to the right of the $`\eta _{\mathrm{Edd}}=0.01`$ line are low-luminosity AGNs. Although this is a good relation, it is only useful for Compton thin galaxies (column densities $`N_\mathrm{H}10^{24}\mathrm{cm}^2`$). Unfortunately for a significant fraction of type 2 objects the $`210`$keV X-ray emission is not transmitted through the torus (e.g., Bassani et al. 1999), and the $`210`$keV luminosities cannot be used to trace the obscured BH growth.

The scatter in the relations may be accounted for by efficiency variations in the accretion process in different kinds of AGN and/or evolutionary effects. The evolution will be relevant if the accretion rate is higher in the early growing stage of an AGN, as proposed by Mathur, Kuraszkiewicz, & Czerny (2001) to explain the properties of narrow line Sy 1 galaxies (NLS1s). This class of galaxies would represent an early stage of the AGN evolution based upon their lower BH to bulge mass ratios. The four NLS1s (shown as open triangles in Figure 1) in our sample show a similar behavior to the other Sy galaxies in both panels of Figure 1, if only with slightly higher accretion rates than the average. However the small number of NLSy1s does not allow us to reach a firm conclusion.

Another factor that may introduce some scatter in the mid-infrared vs. BH mass relation could be the presence of nuclear star formation in Sy 2 galaxies. If star formation were present (especially for those galaxies observed with the relatively large ISO aperture), we would expect an excess of $`10\mu `$m emission for a given BH mass. We do not find such excesses in our sample of galaxies. On the other hand, if the optical depths were high for type 2 objects even at $`10\mu `$m, then for a given BH mass they should lie below type 1 objects. Type 1 objects do not generally lie above type 2 objects in the $`10\mu `$m luminosity vs. BH mass diagram, suggesting that the reprocessing of the $`10\mu `$m nuclear emission is not severely affected by geometric and optical depth effects (see also Alonso-Herrero et al. 2001).

Hosokawa et al. (2001) and Kawaguchi, Shimura, & Mineshige (2001) constructed a disk-corona model to account for the X-ray through infrared spectral energy distributions (SEDs) of quasars and Sy galaxies. These models demonstrate that the AGN SEDs are sensitive to both the BH mass ($`M_{\mathrm{BH}}`$) and accretion rate ($`\dot{M}`$). In the optical, UV and soft X-ray spectral ranges decreasing $`\dot{M}`$ and/or increasing $`M_{\mathrm{BH}}`$ will move the peak of the SED toward longer wavelengths (i.e., figure 1 in Hosokawa et al. 2001). On the other hand, the hard X-ray and mid-infrared spectral shapes remain approximately constant for varying $`M_{\mathrm{BH}}`$ and $`\dot{M}`$. This is expected, as in this model the infrared bump is produced by dust heated by the AGN, and modelled with constant spectral indices, regardless of $`M_{\mathrm{BH}}`$ or $`\dot{M}`$. The infrared power varies so that the luminosity ratio of the big blue bump to the infrared bump remains constant.

The predicted $`10\mu `$m luminosity as a function of $`M_{\mathrm{BH}}`$ from Hosokawa et al. (2001) models is shown in Figure 1 (left panel) for two accretion rates555For reference in this model a BH with an accretion rate of $`\dot{M}=12\times L_{\mathrm{Edd}}/c^2`$ will represent a disk with Eddington luminosity, $`L_{\mathrm{Edd}}=1.5\times 10^{38}\frac{M_{\mathrm{BH}}}{\mathrm{M}_{}}`$ (in erg s$`{}_{}{}^{1})`$ $`L_{\mathrm{Edd}}`$. $`\dot{M}=50.5\times L_{\mathrm{Edd}}/c^2`$. The third accretion rate $`\dot{M}=0.05\times L_{\mathrm{Edd}}/c^2`$ curve has been extrapolated from the above mass accretion rates. The $`10\mu `$m luminosities of PG quasars are well reproduced with accretion rates of $`\dot{M}=0.55\times L_{\mathrm{Edd}}/c^2`$ (as found also by Kawaguchi et al. 2001 from the optical to X-ray properties), whereas Sy galaxies in general require lower accretion rates (see also e.g., Su & Malkan 1989).

Figure 2 shows the ratio of the optical $`B`$-band, near-infrared $`J`$-band, $`10\mu `$m and hard X-ray to the AGN bolometric luminosity as a function of $`M_{\mathrm{BH}}`$ at two different redshifts $`z=0`$ (left panel) and $`z=3`$ (in the observer’s rest-frame, right panel) from Hosokawa et al. (2001) model predictions. The accretion rate are $`\dot{M}=5.0\times L_{\mathrm{Edd}}/c^2`$ and $`\dot{M}=0.5\times L_{\mathrm{Edd}}/c^2`$. From this figure it is clear that the ratios of the $`10\mu `$m and hard X-ray luminosities to the AGN bolometric luminosity remain approximately constant for increasing $`M_{\mathrm{BH}}`$ at a given $`\dot{M}`$ over redshifts of $`z=03`$. Moreover, the ratio of the $`10\mu `$m luminosity to the AGN bolometric luminosity is also almost constant for the two accretion rates. This shows that the $`10\mu `$m luminosites can be used to trace the BH growth by comparing the mid-infrared luminosity functions at different epochs.

The optical to bolometric luminosity ratio varies significantly as a function of the BH mass for redshifts $`z=0`$ and $`z=3`$, and the two accretion rates. Although the $`J`$-band to bolometric luminosity ratio does only vary slightly with $`M_{\mathrm{BH}}`$ at $`z=0`$ and $`\dot{M}=5.0\times L_{\mathrm{Edd}}/c^2`$, it does significantly at $`z=3`$ for both accretion rates. Leaving aside the unavoidable obscuration (viewing angle) effects of type 2 AGNs, the complex dependence of the UV, optical and near-infrared SEDs with $`M_{\mathrm{BH}}`$ and $`\dot{M}`$ would make it more intricate the interpretation of the time evolution of luminosity functions to infer the BH growth.

There is now growing evidence for a missing population of obscured type 2 AGNs, necessary to reproduce the hard X-ray background (e.g., Fabian & Iwasawa 1999), as well as to reconcile the observed local comoving density of BHs with that predicted from the QSO luminosity function at $`z=3`$ (Haehnelt, Natarajan, & Rees 1998). This indicates that a significant fraction of the accretion by BHs may be obscured by dust. The two relations discussed here suggest that the hard X-ray luminosity and more importantly, the mid-infrared luminosity can be effective tools to probe the processes taking place in obscured nuclei of galaxies. These relations can be used to conduct a detailed census of the BH masses over a range of redshifts, and measure directly their growth rate as long as an estimate of the accretion rate is available. The present X-ray space missions are already producing results (e.g., Brandt et al. 2002 and references therein). The launch of SIRTF will provide further capabilities to obtain $`M_{\mathrm{BH}}`$ for obscured AGNs at redshifts of up to $`z4`$.

Summarizing, for a heterogeneous sample of local AGNs we find: (i) the $`M_{\mathrm{BH}}`$ correlates with the $`10\mu `$m nuclear luminosity in local Sy galaxies and PG quasars; (ii) the available data suggest no significant differences between type 1 and type 2 objects, implying that the reprocessing of the $`10\mu `$m nuclear emission is not severely affected by geometric and optical depth effects; (iii) there is a good relation between the $`M_{\mathrm{BH}}`$ and the $`210`$keV hard X-ray luminosity, but only for Compton thin galaxies; (iv) the $`M_{\mathrm{BH}}`$ vs. $`L(10\mu \mathrm{m})`$ relation will allow to estimate the obscured BH mass growth over cosmological times using future observations with SIRTF.

We thank Dave Alexander, Roberto Maiolino and Laura Ferrarese for helpful discussions, and Paulina Lira for providing us with data prior to publication. We are also grateful to an anonymous referee for useful comments. AAH acknowledges support from NASA Contract 960785 through the Jet Propulsion Laboratory.


text:
# References

1. The possibility that nature becomes supersymmetric at the weak scale offers a well motivated and exciting scenario for physics beyond the standard model. It is well motivated because it allows a dynamical generation of the weak scale, and an understanding of why this scale is much less than the Planck scale. Furthermore, it gives rise to a highly successful numerical prediction for gauge coupling unification. It is exciting because weak scale supersymmetry will be thoroughly tested at hadron colliders over the next decade.

Nevertheless, the underlying structure of the fields and interactions of the weak scale supersymmetric theory contains three puzzles:

* In non-supersymmetric field theories there are three distinct types of fields, corresponding to particles with spin $`0`$, $`1/2`$ and $`1`$. In this case, there is no doubt as to what distinguishes the Higgs field from the lepton doublet field. In contrast, in supersymmetric field theories there are just two types of fields: vector multiplets and chiral multiplets. Hence it is now no longer clear what distinguishes a Higgs field, $`H`$, from a matter field, $`M`$. In particular, the down type Higgs and lepton doublets have identical gauge and spacetime properties; what distinguishes the Higgs boson from the sneutrino?
* Phenomenologically, interactions cannot be the most general allowed by known gauge and spacetime symmetries. The superpotential must contain interactions of the form $`MMH`$ for quark and lepton masses, and most probably $`MMHH`$ at the weak scale for neutrino masses. Yet other forms of superpotential interactions, such as $`MH`$, $`M^3`$, and $`M^4`$ are highly constrained by neutrino masses and proton decay. Some of these interactions are either highly suppressed or forbidden.
* A superpotential interaction of the form $`HH`$ is a special and intriguing case. On one hand it must be highly suppressed since a coefficient of order the Planck scale or unified mass scale would remove the Higgs doublets from the low energy theory. On the other hand it cannot vanish since otherwise there is a massless charged fermion coming from the Higgs/vector multiplets. Indeed the theory is only realistic if the coefficient, $`\mu `$, is of order the weak scale, leading to the well known $`\mu `$ problem. Why should this supersymmetric parameter be of order the supersymmetry breaking scale?

The other great mystery of low energy supersymmetry is the origin of supersymmetry breaking. Like the supersymmetric interactions, a great deal about the structure of the supersymmetry breaking interactions is governed by the requirement of consistency with experiment. However, nothing determines the “messenger” scale, $`M_m`$: the highest scale at which the supersymmetry breaking interactions of squarks, sleptons, Higgs and gauginos are local. In supergravity theories this locality is maintained up to the Planck scale , but other methods of transporting supersymmetry breaking to the standard model superpartners, such as gauge mediation , gaugino mediation , and boundary condition supersymmetry breaking , have messenger scales that can be anywhere between the weak and Planck scales.

In this paper we study the consequences of imposing an $`R`$ symmetry on the supersymmetric interactions of effective theories with weak scale supersymmetry. This symmetry will allow us to address all three of the puzzles listed above in the the context of supersymmetry breaking with an arbitrary messenger scale. We will describe a continuous symmetry, $`U(1)_R`$, although a discrete subgroup is sufficient for our purposes. The soft supersymmetry breaking operators break the $`R`$ symmetry, since they include Majorana gaugino masses, and therefore holomorphic, $`A`$ term, scalar interactions. Using this $`R`$ symmetry, we find a new mechanism for solving the $`\mu `$ problem for arbitrary messenger scales, and this requires a unique choice for the $`R`$ quantum numbers of matter and Higgs fields. Furthermore, we find that this is also the unique choice which accounts for the absence of $`MH`$, $`M^3`$ and $`M^4`$ superpotential interactions, while allowing consistency with quark-lepton unification. Finally, this $`R`$ symmetry forces a distinction between Higgs and lepton doublets.

There are two well-known classes of solutions to the $`\mu `$ problem. One class corresponds to modifying physics at the Planck scale by adding non-renormalizable operators ; the other changes the physics at the weak scale, as is the case in the Next-to-Minimal Supersymmetric Standard Model (NMSSM) . The first approach requires supersymmetry breaking to be mediated at the Planck scale, while the second requires a departure in the weak scale theory from the Minimal Supersymmetric Standard Model (MSSM). Neither of these requirements apply to our mechanism.

In our mechanism, we suppose that all supersymmetric interactions respect some global symmetry, $`G`$, which forbids $`\mu `$ and commutes with both gauge and flavor symmetries. To allow for the the possibility of unification of the matter, we require the $`G`$ charges of all matter multiplets to be equal. We consider that some field acquires a vacuum expectation value (vev) at scale $`\mathrm{\Lambda }`$ between the weak and Planck scales. Supersymmetry breaking interactions, mediated at any scale above $`\mathrm{\Lambda }`$, break $`G`$ and cause a deformation in the vacuum resulting in the generation of $`\mu `$ of order the supersymmetry breaking scale. Such theories have been constructed for the case of mediation at the Planck scale . However, for such high mediation scales $`\mu `$ is not a problem since the Giudice-Masiero mechanism is available. We study the case of an arbitrary mediation scale. Related solutions for the $`\mu `$ problem in the context of gauge mediation have been proposed , but differ from our mechanism in the origin of the vacuum deformation.

2. If we require that our solution to the $`\mu `$ problem apply to arbitrary mediation scales, we find that the operator giving rise to $`\mu `$ is uniquely determined. We assume this operator must be present at tree level — if $`\mu `$ arises radiatively then the soft $`\mu BH_uH_d`$ term is typically generated at too high a level. Furthermore, $`\mu `$ must arise from a renormalizable operator. If the operator had a coefficient suppressed by powers of the Planck scale, it would not be possible to get a $`\mu `$ parameter of the desired size for arbitrary values of the messenger scale. Thus the operator generating $`\mu `$ is unique:

$$W_\mu =\lambda XH_uH_d,$$
(1)
where $`X`$ is a standard model gauge singlet chiral superfield. Our mechanism requires that the supersymmetric interactions break some symmetry at a scale $`\mathrm{\Lambda }`$, giving a mass to $`X`$ and forcing $`A_X=F_X=0`$. Here $`A_X`$ and $`F_X`$ represent the lowest and highest components of the chiral superfield $`X`$, respectively. Non-zero values for $`A_X`$ and $`F_X`$ are generated by supersymmetry breaking.

We discover that our global symmetry, $`G`$, must be an $`R`$ symmetry from the following argument. The order of magnitude of $`A_X`$ and $`F_X`$ generated after supersymmetry breaking can be understood from the $`G`$ charges of both $`X`$ and the soft supersymmetry breaking operators. If $`G`$ is a non-$`R`$ symmetry then $`A_X`$ and $`F_X`$ have the same $`G`$ transformation, so that both $`\mu `$ and $`\mu B`$ are generated at the same order in supersymmetry breaking. Hence $`B`$ is of order $`\mathrm{\Lambda }`$ and much too large: our mechanism requires $`G`$ to be an $`R`$ symmetry.

Now, under the $`R`$ symmetry, there are only two types of supersymmetry breaking terms appearing in the scalar potential. There are holomorphic terms, which we denote by $`A`$ and assign $`R`$ charge $`2`$, and there are non-holomorphic terms, denoted $`m^2`$, which have $`R`$ charge zero. A successful solution to the $`\mu `$ problem requires $`\mu `$, and therefore $`A_X`$, to be linear in supersymmetry breaking. This requires that $`X`$ has $`R`$ charge $`+2`$ or $`2`$, so that $`A_X`$ can be generated proportional to $`A^{}`$ or $`A`$. Interestingly, this automatically leads to $`F_X`$ of exactly the right order, since now $`F_X`$ has $`R`$ charge $`0`$ or $`4`$ and is naturally generated at second order in supersymmetry breaking proportional to $`|A|^2,m^2`$ or $`A^2`$.

Let us consider our two possible $`R`$-charge assignments separately. If $`X`$ has $`R`$ charge $`2`$, then $`H_uH_d`$ has total $`R`$ charge $`+4`$. Due to the presumed existence of unification, all matter fields carry the same $`R`$ charge. The Yukawa couplings then force equal $`R`$ quantum numbers for the two Higgs doublets, $`R(H_u)=R(H_d)`$, so that $`R(H)=2`$ and $`R(M)=0`$. However, in this case the superpotential interaction $`MH`$ is allowed, which will push some MSSM matter fields to have masses of order the Planck scale. Therefore, we prefer to consider the charge assignment $`R(X)=2`$. Again, the Yukawa couplings require $`R(H_u)`$=$`R(H_d)`$, so we have

$$R(H)=0R(M)=1.$$
(2)
This $`R`$ symmetry is extremely powerful: as well as distinguishing between the lepton and Higgs doublet and forbidding $`MH`$, it also forbids $`M^3`$ and $`M^4`$, leading to baryon and lepton number conservation from operators of dimension four and five. Therefore, our solution to the $`\mu `$ problem has forced us to forbid dangerous dimension four and five operators that might lead to too rapid proton decay. Moreover, if the soft supersymmetry breaking operators provide the only source of $`R`$ breaking, then $`R`$ parity remains unbroken, leading to stability of the lightest superpartner.111 $`R`$ parity forbids the generation of dimension four baryon and lepton number violating operators, even after supersymmetry is broken. While dimension five proton decay operators could be generated with small coefficients proportional to supersymmetry breaking, they are phenomenologically irrelevant.

With a mild assumption about the origin of neutrino masses, a completely independent argument will lead us to an identical conclusion for the global symmetry $`G`$. The Yukawa interactions, $`M^2H`$, possess a non-$`R`$ Peccei-Quinn symmetry (PQ: $`M(1),H(2)`$) as well as the $`R`$ symmetry of Eq. (2). In fact, $`G`$ must be a linear combination of $`R`$ and PQ, or one of its subgroups. The existence of small neutrino masses strongly suggests that the superpotential also contains $`M^2H^2`$. Provided that this interaction is not generated by supersymmetry breaking , this immediately implies that $`H`$ is neutral under $`G`$, and hence $`G`$ must be the $`R`$ symmetry given in Eq. (2).

3. The simplest model which realizes the above general mechanism for generating the $`\mu `$ term is given by the superpotential

$$W_0=fX(Y^2\mathrm{\Lambda }^2),$$
(3)
where $`R(X)=2`$ and $`R(Y)=0`$. Here we imagine that the scale $`\mathrm{\Lambda }`$ is much larger than the weak scale, although this is not necessary for our mechanism to work. We also impose a discrete symmetry $`YY`$, so that the gauge hierarchy is not destabilized by the generation of a large tadpole operator for a singlet field.222 We assume that violations of global symmetries by nonperturbative gravitational effects  are sufficiently suppressed. Without supersymmetry breaking, the minimum of the potential lies at $`A_X=0`$ and $`A_Y=\mathrm{\Lambda }`$, satisfying $`F_X=F_Y=0`$. There is no flat direction at this level, and all the fields have masses of order $`\mathrm{\Lambda }`$.

When we add supersymmetry breaking terms, with a scale $`\stackrel{~}{m}`$ of order the weak scale, the vevs will shift. The most general soft supersymmetry breaking terms are given by

$$_{\mathrm{soft},0}=m_X^2|X|^2m_Y^2|Y|^2(a_fXY^2a_\mathrm{\Lambda }\mathrm{\Lambda }^2X+\mathrm{h}.\mathrm{c}.).$$
(4)
Here, $`R(m_X^2)=R(m_Y^2)=0`$ and $`R(a_f)=R(a_\mathrm{\Lambda })=2`$, and we have used $`X`$ and $`Y`$ to denote the scalar fields of the respective chiral superfields. By minimizing the scalar potential, we obtain $`X(a_\mathrm{\Lambda }^{}a_f^{})/4|f|^2\stackrel{~}{m}`$ and $`F_X[(a_\mathrm{\Lambda }+a_f)(a_\mathrm{\Lambda }^{}a_f^{})/4|f|^2+m_Y^2]/2f\stackrel{~}{m}^2`$. The vevs of $`X`$ and $`F_X`$ are both of the order of the weak scale as indicated by the previous general analysis. Therefore, if we introduce couplings to the Higgs doublet $`W=W_0+\lambda XH_uH_d`$ and $`_{\mathrm{soft}}=_{\mathrm{soft},0}m_{H_u}^2|H_u|^2m_{H_d}^2|H_d|^2(a_\lambda XH_uH_d+\mathrm{h}.\mathrm{c}.)`$, $`\mu `$ and $`\mu B`$ terms of order $`\stackrel{~}{m}`$ and $`\stackrel{~}{m}^2`$ are generated as

$$\mu =\lambda X\frac{\lambda (a_\mathrm{\Lambda }^{}a_f^{})}{4|f|^2},$$
(5)
$$\mu B=\lambda F_X+a_\lambda X\frac{(2fa_\lambda \lambda a_\mathrm{\Lambda }\lambda a_f)(a_\mathrm{\Lambda }^{}a_f^{})}{8f|f|^2}\frac{\lambda m_Y^2}{2f},$$
(6)
where $`\mu `$ and $`\mu B`$ are defined by $`W=\mu H_uH_d`$ and $`=\mu BH_uH_d+\mathrm{h}.\mathrm{c}.`$. Although $`\mu =0`$ for $`a_\mathrm{\Lambda }=a_f`$, it is natural to expect $`a_\mathrm{\Lambda }a_f`$ since the two parameters run differently under renormalization group evolution; a simple realistic example will be given later. An important point here is that there exists a parameter region where the additional Higgs coupling in Eq. (1) does not change the gross dynamics at the high scale. For instance, if $`m_{H_u}^2`$ and $`m_{H_d}^2`$ are sufficiently large, the vev of order $`\mathrm{\Lambda }`$ is still entirely contained in the $`Y`$ field, and the vevs for the Higgs doublets stay smaller than the weak scale. Note that here the various supersymmetry breaking parameters are evaluated at the scale $`\mathrm{\Lambda }`$. Thus, both $`m_{H_u}^2`$ and $`m_{H_d}^2`$ can be positive without conflicting with electroweak symmetry breaking. Below the scale $`\mathrm{\Lambda }`$, the heavy fields $`X`$ and $`Y`$ are integrated out and the effective theory contains only the Higgs doublets with the $`\mu `$ and $`\mu B`$ parameters given by Eqs. (5, 6). The soft supersymmetry breaking parameters must be further evolved from $`\mathrm{\Lambda }`$ to $`\stackrel{~}{m}`$ using the renormalization group equations of this effective theory (MSSM), to evaluate electroweak symmetry breaking.

An important requirement for our mechanism is that $`\mathrm{\Lambda }`$ must be smaller than the messenger scale of supersymmetry breaking, $`\mathrm{\Lambda }<M_m`$. Therefore, the superpotential Eq. (3) itself is not sufficient for a complete solution of the $`\mu `$ problem (except for the supergravity mediation case), since we have introduced by hand a mass parameter, $`\mathrm{\Lambda }`$, smaller than the fundamental scale. A complete solution, however, is obtained if we generate the scale $`\mathrm{\Lambda }`$ by the dynamics of strong gauge interactions. Consider, for example, the $`SU(2)_S`$ gauge theory with four doublet chiral superfields $`Q_i`$ ($`i=1,\mathrm{},4`$) with the following superpotential:

$$W_{0,\mathrm{tree}}=fX\left(Y^2(QQ)\right)+f^{}X^a(QQ)_a.$$
(7)
This superpotential explicitly breaks a flavor $`SU(4)_F`$ symmetry of the $`Q_i`$ down to $`SP(4)_F`$; $`(QQ)`$ and $`(QQ)_a`$ ($`a=1,\mathrm{},5`$) denote singlet and five-dimensional representations of $`SP(4)_F`$ given by suitable combinations of gauge invariants $`Q_iQ_j`$. The strong dynamics of the $`SU(2)_S`$ gauge theory is described by the effective superpotential

$$W_{0,\mathrm{eff}}=W_{0,\mathrm{tree}}+S\left((QQ)^2+(QQ)_a^2\mathrm{\Lambda }^4\right),$$
(8)
where $`S`$ is an additional Lagrange multiplier chiral superfield . For a relatively large value of the coupling $`f^{}`$, the vacuum lies at $`(QQ)=\mathrm{\Lambda }^2`$ and $`(QQ)_a=0`$, so that the superpotential $`W_{0,\mathrm{eff}}`$ is effectively reduced to Eq. (3). Note that the original tree-level superpotential, Eq. (7), does not contain any mass parameters and is invariant under the $`U(1)_R`$ symmetry with $`R(X)=R(X^a)=2`$ and $`R(Q_i)=0`$. In fact, it is the most general superpotential consistent with the combined $`R`$ and $`SP(4)_F`$ symmetries. (A linear term in $`X`$ is forbidden either by requiring that the superpotential not contain any mass parameters, or by imposing an anomalous discrete $`Z_3`$ symmetry under which all the fields are transformed by $`\mathrm{exp}(2\pi i/3)`$.) It is also important that $`U(1)_R`$ does not have an anomaly for $`SU(2)_S`$ (i.e. $`\mathrm{\Lambda }`$ does not carry $`U(1)_R`$ charge), so that the previous general argument is not affected by the strong $`SU(2)`$ gauge dynamics.

We now consider an application of our mechanism to realistic theories. We find that the mechanism fits beautifully into the framework where small neutrino masses are generated by integrating out right-handed neutrino fields through the see-saw mechanism . We consider the following theory. In addition to the usual three generations of standard-model quark and lepton superfields, $`Q`$, $`U`$, $`D`$, $`L`$ and $`E`$, we introduce three right-handed neutrino superfields $`N`$. Here, we have omitted generation indices. The Yukawa couplings are given by

$$W_{\mathrm{Yukawa}}=y_uQUH_u+y_dQDH_d+y_eLEH_d+y_\nu LNH_u.$$
(9)
We also introduce the $`U(1)_X`$ gauge symmetry, contained in $`SO(10)/SU(5)`$, under which various fields transform as $`Q(1)`$, $`U(1)`$, $`D(3)`$, $`L(3)`$, $`E(1)`$ and $`N(5)`$. This gauge symmetry is broken by the vevs of the fields $`\mathrm{\Phi }(10)`$ and $`\overline{\mathrm{\Phi }}(10)`$ through the superpotential

$$W_{\mathrm{Breaking}}=fX\left(\mathrm{\Phi }\overline{\mathrm{\Phi }}(QQ)\right)+f^{}X^a(QQ)_a.$$
(10)
Here $`(QQ)`$ and $`(QQ)_a`$ are gauge invariants consisting of $`Q_i`$, the doublets under the strong $`SU(2)_S`$ gauge interaction (see discussion around Eqs. (7, 8)). Note that the above superpotentials, Eqs.(9, 10), do not contain any mass parameters and are invariant under the $`U(1)_R`$ symmetry, $`R(Q)=R(U)=R(D)=R(L)=R(E)=R(N)=1`$, $`R(H_u)=R(H_d)=R(\mathrm{\Phi })=R(\overline{\mathrm{\Phi }})=R(Q_i)=0`$ and $`R(X)=R(X^a)=2`$. For a relatively large $`f^{}`$, the dynamics of the $`SU(2)_S`$ gauge interaction cause the condensation of $`(QQ)=\mathrm{\Lambda }^2`$, which is transmitted to the vevs of the $`\mathrm{\Phi }`$ and $`\overline{\mathrm{\Phi }}`$ fields, $`\mathrm{\Phi }=\overline{\mathrm{\Phi }}=\mathrm{\Lambda }`$. Here, the equality $`\mathrm{\Phi }=\overline{\mathrm{\Phi }}`$ is forced by the $`U(1)_X`$ $`D`$-term condition. The vevs for all the other fields are zero at this stage: $`X=X^a=(QQ)_a=0`$. After introducing soft supersymmetry breaking operators, the vevs of the fields shift. In particular, non-vanishing vevs for $`X`$ and $`F_X`$ are generated as $`X\stackrel{~}{m}`$ and $`F_X\stackrel{~}{m}^2`$, as long as holomorphic soft supersymmetry breaking parameters are not subject to the special relation $`a_{X\mathrm{\Phi }\overline{\mathrm{\Phi }}}=a_{X(QQ)}`$ at the scale $`\mathrm{\Lambda }`$. In fact, it is quite natural to expect that the $`A`$ terms for $`X\mathrm{\Phi }\overline{\mathrm{\Phi }}`$ and $`X(QQ)`$ are different since they are renormalized differently above the scale $`\mathrm{\Lambda }`$; for example, they receive contributions from $`U(1)_X`$ and $`SU(2)_S`$ gauginos, respectively. Therefore, by introducing the couplings

$$W_{\mathrm{Masses}}=\kappa \overline{\mathrm{\Phi }}N^2+\lambda XH_uH_d,$$
(11)
the Majorana masses $`M_R`$ for the right-handed neutrinos of order $`M_R\overline{\mathrm{\Phi }}\mathrm{\Lambda }`$ and $`\mu `$ and $`\mu B`$ parameters of order $`\mu B\stackrel{~}{m}`$ are generated. As in the previous example, the superpotential $`W=W_{\mathrm{Yukawa}}+W_{\mathrm{Breaking}}+W_{\mathrm{Masses}}`$ is the most general renormalizable superpotential consistent with the gauge $`SU(3)_C\times SU(2)_L\times U(1)_Y\times U(1)_X\times SU(2)_S`$ and global $`U(1)_R\times SP(4)_F`$ symmetries of the theory, after removing a linear term in $`X`$ as before.

We finally discuss how our general mechanism works explicitly in various supersymmetry breaking scenarios. In gauge mediated supersymmetry breaking, our mechanism requires that the mass of the messenger fields, $`M_m`$, is larger than $`\mathrm{\Lambda }`$. Since the holomorphic supersymmetry breaking terms ($`A`$ terms) required for the $`\mu `$-term generation are small at the messenger scale $`M_m`$, they must be generated by renormalization group evolution from $`M_m`$ to $`\mathrm{\Lambda }`$. This can be accomplished, for example, by giving non-trivial $`SU(2)_S`$ or $`U(1)_X`$ quantum numbers to the messenger fields. In the case of gaugino mediation and boundary condition supersymmetry breaking, our mechanism requires that the compactification scale is larger than $`\mathrm{\Lambda }`$. In these cases, the relevant $`A`$ terms of order the weak scale may already exist at the compactification scale, so we do not necessarily have to rely on renormalization group evolution for their generation.333 In anomaly mediation , the holomorphic supersymmetry breaking parameter associated with the scale $`\mathrm{\Lambda }`$ is large — of the order of the gravitino mass $`10\mathrm{TeV}`$. Thus $`\lambda 10^2`$ is needed to generate a $`\mu `$ parameter of the correct size. Then, to avoid a too large $`\mu B`$ term, a cancellation between two contributions, such as $`fa_\lambda `$ and $`\lambda a_\mathrm{\Lambda }`$ in Eq. (6), is required at the $`1\%`$ level. In any of these mediation mechanisms, $`A`$ is real in the basis where the gaugino masses are real (except for the case of gaugino mediation with tree-level $`A`$ terms), and our origin for $`\mu `$ and $`\mu B`$ then leads to a real $`B`$ parameter: the supersymmetric $`CP`$ problem is solved.

4. In this paper we have proposed an origin for the parameters $`\mu `$ and $`\mu B`$ of the minimal supersymmetric standard model, which is applicable for any messenger scale, $`M_m`$. Although quite general, it does require specific symmetries and interactions. Both $`\mu `$ and $`\mu B`$ parameters arise from the superpotential interaction $`XH_uH_d`$. A stage of symmetry breaking occurs at some scale $`\mathrm{\Lambda }<M_m`$, giving a mass of order $`\mathrm{\Lambda }`$ to $`X`$, while determining $`A_X=F_X=0`$. Providing the form for the superpotential is guaranteed by an $`R`$ symmetry, with the quantum numbers of Eq. (2), the soft supersymmetry breaking operators, with coefficients $`A`$ and $`m^2`$, lead to a small readjustment of the vacuum, giving

$$\mu A^{},\mu B|A|^2,m^2.$$
(12)
This $`R`$ symmetry provides a distinction between Higgs and matter superfields, and forbids superpotential interactions that would otherwise lead to baryon number violation at too rapid a rate. Although $`R`$ is broken by supersymmetry breaking, the discrete $`R`$ parity survives so that the lightest superpartner is stable. In the case that the original $`R`$ symmetry is continuous, an $`R`$ axion will be produced by the underlying dynamics which breaks supersymmetry. If all the $`R`$ breaking effects are generated spontaneously (including the constant term in the superpotential needed to cancel the cosmological constant), the dominant mass contribution to the $`R`$ axion will come from the QCD anomaly of the $`R`$ symmetry. In this case, the $`R`$ axion provides a solution to the strong $`CP`$ problem .

At first sight our $`R`$ symmetry appears to be in conflict with grand unification: since $`R`$ forbids $`H_uH_d`$, it also forbids the corresponding mass term for the colored Higgs triplets of unified theories. However, this turns out to be a virtue — such mass terms need to be forbidden to avoid too large a proton decay rate mediated by triplet Higgsino exchange. The colored partners of $`H_{u,d}`$ must become heavy by acquiring mass terms coupling them to other colored states of the theory. This occurs in the missing partner  and Dimopoulos-Wilczek  mechanisms; however, although these mechanisms are consistent with an underlying $`U(1)_R`$ symmetry, in the simplest such models $`U(1)_R`$ is broken at the unification scale, so our $`\mu `$ generation mechanism may not work in these cases. In contrast, in Kaluza-Klein grand unification  the desired colored Higgs mass terms arise while preserving $`U(1)_R`$ symmetry, so that our $`\mu `$ generation mechanism works well in this case.

The $`U(1)_R`$ symmetry is so crucial in providing an understanding of the form for the interactions in the superpotential, it is important to seek its origin. Higher dimensional theories are particularly interesting since they have an enlarged set of supersymmetry transformations, which results in a global $`R`$ symmetry in the equivalent four dimensional description. In the case of a five dimensional grand unified theory, compactification breaks the unified gauge symmetry and also the $`SU(2)_R`$ symmetry to $`U(1)_R`$, so that precisely the $`R`$ charges considered here may arise .

Acknowledgements

Y.N. thanks the Miller Institute for Basic Research in Science for financial support. This work was supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, of the U.S. Department of Energy under Contract DE-AC03-76SF00098, and in part by the National Science Foundation under grant PHY-00-98840.


text:
# Quantum Jumps, EEQT and the Five Platonic Fractals

## 1 Introduction

Recently, Numerical Algorithms for quantum jumps were discovered. First they were introduced in quantum optics as convenient and effective tools for numerical simulations of the Liouville equation . A short history and more references can be found in Ref. . Ph. Blanchard and A. Jadczyk, in a series of papers on EEQT - Event Enhanced Quantum Theory (cf.  and references therein) - developed a new approach that unifies what John von Neumann called U- and R-processes  into a single piecewise deterministic process (PDP) where continuous evolution of a system is cyclically interrupted by discontinuous ”jumps”.

Normally one would expect that time evolution of a physical system is described by a differential equation. That is how laws of physics are usually expressed. Here however we have a surprise: in EEQT a history of an individual quantum system, coupled to a monitoring device, as in every real world experiment, is described by a process or an algorithm rather than by a differential equation. The PDP is similar to those studied in the science of economics, where periods of smooth fluctuations are interrupted by market crashes . A quantum jump is what corresponds to a market crash - a discontinuity, a ’catastrophe’. But discontinuities and catastrophes have their own laws, and here comes the concept of a piecewise deterministic Markov process111The Markovian property is not really important. \- a PDP.

PDPs are the simple and elegant ways to describe the world in terms of cyclic, rather than linear, time; that is the world of cyclically, though somewhat irregularly, recurring catastrophes. In the cases studied in the present paper the catastrophes come from the coupling of a quantum system to a system of two-state ”detectors”. In this case the catastrophes are not really catastrophes for the detectors - detectors just flip, which is exactly what detectors do for living. But these flips bring catastrophes for the quantum system, because with each flip of the detector, with each ”event”, as we call it, the quantum system state vector breaks its continuous evolution, and instantaneously jumps to a different state - it ”rejuvenates” and it starts another cycle of a peaceful, continuous evolution - till the next catastrophe.

Note. The term instantaneously in the last sentence may suggest that our formalism is incompatible with Einstein’s relativity. That this is not so has been demonstrated in Ref. (for a somewhat different approach cf also ). The point is that for a relativistic theory the role of time is being played by the Fock-Schwinger ”proper time” as a Floquet variable. The reader should bear in mind that the EEQT algorithm is explicitly nonlocal: to simulate a history of an individual system integrations over entire space (or space-time) are needed.

The EEQT algorithm generating quantum jumps is similar in its nature to a nonlinear iterated function system (IFS)  (see also  and references therein) and, as such, it generically produces a chaotic dynamics for the coupled system. Here the probabilities assigned to the maps are derived from quantum transition probabilities and thus depend on the actual point, but such generalizations of the IFS’s have been also studied (cf.  and references therein). In the present paper we describe the algorithm generating quantum fractals, that is self-similar patterns on the projective plane $`_1()S^2`$, when a continuous in time ”measurement” of several spin directions at once takes place.

Note: The term quantum fractals has been used before by Casati et al.  in a different context.

Note: The operators for different spin directions do not commute, but this does not contradicts Heisenberg’s uncertainty relations as these deal with statistical ensembles averages, while here we are describing an individual quantum system. In fact, realizing the chaotic behavior of a quantum state vector, when several noncommuting observables are being simultaneously monitored, can help us to understand the mechanisms of statistical uncertainties.

As stated above the algorithm of EEQT describes a piecewise deterministic random process - periods of a smooth evolution interspersed with catastrophic ”jumps.” Of course, once we have individual description, we can also get the laws for statistical ensembles. Here we get a nice, linear, Liouville evolution equation for measures - as it is usual in studying chaotic dynamics. The fact that there is a unique PDP generating the Liouville equation in the framework of EEQT has been proven in Ref. .

## 2 EEQT - Quantum Fractals

### 2.1 Geometry

We have published, as an OpenSource project , the algorithm implemented in Java that generates the Five Platonic Fractals - that is fractals generated by five most symmetric detector configurations. The algorithm generates self-similar patterns on a sphere of a unit radius. The points on the sphere represent (pure) states of the simplest quantum system - the spin $`1/2`$ rotator. This spin $`1/2`$ quantum system is coupled, continuously in time, to a finite number of symmetrically distributed spin-direction detectors. Thus the symmetry of the pattern reflects the symmetry of the detector directions distribution. Each spin direction is characterized by a vector $`𝐧`$ of unit length. Here we study the most symmetrical configurations, therefore we chose direction vectors $`𝐧_i`$ pointing from the origin to the vertices of one of the five platonic solids. We consider the following five detectors configurations:

1. tetrahedron: 4 detectors along the directions $`𝐧[i],i=1,\mathrm{},4`$
{{0, 0, 1.}, {a, 0, -a}, {-a, a, -a}, {-a, -a, -a}}
2. octahedron: 6 detectors along the directions $`𝐧[i],i=1,\mathrm{},6`$
{{0, 0, 1.}, {1., 0, 0}, {0, 1., 0},
{-1., 0, 0}, {0, -1., 0}, {0, 0, -1.}}
3. cube: 8 detectors along the directions $`𝐧[i],i=1,\mathrm{},8`$
{{0, 0, 1.}, {a, 0, a}, {-a, a, a}, {-a, -a, a},
{a, a, -a}, {a, -a, -a}, {-a, 0, -a}, {0, 0, -1.}}
4. icosahedron: 12 detectors along the directions $`𝐧[i],i=1,\mathrm{},12`$
{{0, 0, 1.}, {0.a, 0, a}, {a, a, a}, {-a, a, a},
{-a, -a, a}, {a, -a, a}, {a, a, -a},
{a, -a, -a}, {-a, a, -a}, {-a, 0, -a},
{-a, -a, -a}, {0, 0, -1.}}
5. dodecahedron: 20 detectors along the directions $`𝐧[i],i=1,\mathrm{},20`$
{{0, 0, 1.}, {a, 0, a}, {-a, a, a}, {-a, -a, a},
{a, a, a}, {a, -a, a}, {-a, a, a},
{a, a, a}, {a, -a, a}, {-a, -a, a},
{a, a, -a}, {a, -a, -a}, {-a, a, -a},
{-a, a, -a}, {-a, -a, -a}, {-a, -a, -a},
{a, a, -a}, {a, -a, -a}, {-a, 0, -a},
{0, 0, -1.}}

where the array of real numbers $`a[i],i=1,\mathrm{},17`$ is given in the following table.

| $`a[1]=\frac{3\sqrt{5}}{6}`$ | $`a[2]=\frac{5\sqrt{5}}{10}`$ | $`a[3]=\frac{1}{3}`$ | $`a[4]=\frac{\sqrt{5}1}{2\sqrt{3}}`$ |
| --- | --- | --- | --- |
| $`a[5]=\frac{1}{\sqrt{5}}`$ | $`a[6]=\frac{\sqrt{2}}{3}`$ | $`a[7]=\sqrt{\frac{5\sqrt{5}}{10}}`$ | $`a[8]=\frac{1}{\sqrt{3}}`$ |
| $`a[9]=\frac{2}{3}`$ | $`a[10]=\frac{5+\sqrt{5}}{10}`$ | $`a[11]=\frac{\sqrt{5}}{3}`$ | $`a[12]=\sqrt{\frac{2}{3}}`$ |
| $`a[13]=\sqrt{\frac{5+\sqrt{5}}{10}}`$ | $`a[14]=\frac{3+\sqrt{5}}{6}`$ | $`a[15]=\frac{2}{\sqrt{5}}`$ | $`a[16]=\sqrt{\frac{3+\sqrt{5}}{6}}`$ |
| $`a[17]=\frac{2\sqrt{2}}{3}`$ |  |  |  |

Note: There is no deep reason why we chose these configurations - simplicity and beauty are the main factors here. Notice that, to enable easy zooming onto the attractor, we have chosen the orientations in such a way that in each case the North Pole, with coordinates (0,0,1), of the sphere is occupied by one of the vertices.

### 2.2 The algorithm

Here we describe the algorithm. Comments on its meaning and on the derivation can be found in the endnotes in Section 4

#### 2.2.1 The Hilbert space

The quantum system is represented in a two-dimensional complex Hilbert space, which we realize as $`^2`$ \- the set of all column vectors

$$a=\left(\begin{array}{c}a_1\\ a_2\end{array}\right)$$
(1)
where $`a_1`$ and $`a_2`$ are complex numbers, and with scalar product $`(a,b)`$ defined by

$$(a,b)=\overline{a}_1b_1+\overline{a}_2b_2$$
(2)
where the bar $`\overline{c}`$ stands for the complex conjugation of the complex number $`c.`$

#### 2.2.2 Spin directions

We choose the Pauli matrices $`\sigma _x,\sigma _y,\sigma _z`$ to represent spin directions along $`x,y,z`$ axes respectively.

$$\sigma _1=\sigma _x=\left(\begin{array}{cc}0,& 1\\ 1,& 0\end{array}\right),\sigma _2=\sigma _y=\left(\begin{array}{cc}0,& \mathrm{i}\\ \mathrm{i},& 0\end{array}\right),\sigma _3=\sigma _z=\left(\begin{array}{cc}1,& 0\\ 0,& 1\end{array}\right)$$

Together with the identity matrix

$$\sigma _0=I=\left(\begin{array}{cc}1,& 0\\ 0,& 1\end{array}\right)$$
(3)
they span the whole $`2\times 2`$ complex matrix algebra. In computations it is important to make use of the fact that Pauli matrices (after multiplication by ”-i”) represent the quaternion algebra, that is:

$`\sigma _1^2=\sigma _2^2=\sigma _3^2=I`$
$`\sigma _1\sigma _2=\sigma _2\sigma _1=\mathrm{i}\sigma _3,\sigma _2\sigma _3=\sigma _3\sigma _2=\mathrm{i}\sigma _1,\sigma _3\sigma _1=\sigma _1\sigma _3=\mathrm{i}\sigma _2`$
To each direction $`𝐧`$ in space there is associated spin matrix

$$\sigma (𝐧)=n_1\sigma _1+n_2\sigma _2+n_3\sigma _3=\left(\begin{array}{cc}n_3,& n_1\mathrm{i}n_2\\ n_1+\mathrm{i}n_2,& n_3\end{array}\right)$$
(4)
satisfying automatically $`\sigma (𝐧)^2=I,`$ and with eigenvalues $`+1,1`$. Vectors $`\left(\begin{array}{c}1\\ 0\end{array}\right)`$ and $`\left(\begin{array}{c}0\\ 1\end{array}\right)`$ are eigenvectors of $`\sigma _3`$ to eigenvalues $`+1`$ and $`1`$ respectively and thus correspond to ”North” and ”South” spin orientations respectively. Let $`P(𝐧)`$ denote the projection operator that projects onto eigenstate of $`\sigma (𝐧)`$ to the eigenvalue $`+1.`$ Then $`P(𝐧)`$ is given by the formula:

$$P(𝐧)=\frac{1}{2}(I+\sigma (𝐧)).$$
(5)
Indeed, $`P(𝐧)`$ is Hermitian and has eigenvalues $`\frac{1}{2}(1\pm 1)=`$ $`1`$ or $`0`$ \- thus it is the orthogonal projection, and it projects onto the eigenstate of $`\sigma (𝐧)`$ with spin direction $`𝐧.`$

#### 2.2.3 Fuzzy projections

For each $`𝐧`$ let $`P(𝐧,ϵ)`$ be the fuzzy projection opearator defined by the formula:

$$P(𝐧,ϵ)=\frac{1}{2}(I+ϵ\sigma (𝐧))$$
(6)
where $`0ϵ1`$ ), or better: $`1ϵ,`$ is a parameter that measures the ”fuzziness.” The extreme cases are not the very interesting ones: for $`ϵ=0`$ we get the identity operator - maximal fuzziness and no information whatsoever, while for $`ϵ=1`$ we get the sharp projection $`P(𝐧)=P(𝐧,ϵ=1).`$

We restrict the range of the parameter $`ϵ`$ to the interval $`[0,1]`$ because only in this range $`P(𝐧,ϵ)`$ is a positive operator. It is easy to see that this is so. Indeed, a Hermitian matrix is positive when its eigenvalues are positive, and the eigenvalues of $`P(𝐧,ϵ)`$ are $`(1\pm ϵ)/2,`$ thus $`1ϵ1.`$ On the other hand negative $`ϵ`$ for $`𝐧`$ is the same as positive $`ϵ`$ for $`𝐧,`$ thus we restrict the range of $`ϵ`$ to $`[0,1].`$

It is the operators $`P(𝐧,ϵ)`$ that will act on quantum states to implement ”quantum jumps” whenever detectors ”flip.”

The overall coefficient in the definition (5), chosen to be $`\frac{1}{2}`$ here, is not important because in applications each of the operators $`P(𝐧,ϵ)`$ is multiplied by a coupling constant, and, in our case, when we are not interested in timing of the jumps, the value of the coupling constant plays no role.

#### 2.2.4 Jumps are implemented by fuzzy projections

Let us now discuss the algebraic operation that is associated with each quantum jump. Suppose before the jump the state of the quantum system is described by a projection operator $`P(𝐫)`$, $`𝐫`$ being a unit vector on the sphere. That is, suppose, before the detector flip, the spin ”has” direction $`𝐫`$. Now, suppose the detector $`P(𝐧,ϵ)`$ flips, and the spin right after the flip has some other direction, $`𝐫^{}`$. What is the relation between $`𝐫`$ and $`𝐫^{}`$? It is easy to see that the action of the operator $`P(𝐧,ϵ)`$ on a quantum state vector is given, in terms of operators, by the formula:

$$\lambda (ϵ,𝐧,𝐫)P(𝐫^{})=P(𝐧,ϵ)P(𝐫)P(𝐧,ϵ),$$
(7)
where $`\lambda (ϵ,𝐧,𝐫)`$ is a positive number. It is a simple (though somewhat lengthy) matrix computation that leads to the following result:

$$\lambda (ϵ,𝐧,𝐫)=\frac{1+ϵ^2+2ϵ(𝐧𝐫)}{4}$$
(8)
$$𝐫^{}=\frac{(1ϵ^2)𝐫+2ϵ(1+ϵ(𝐧𝐫))𝐧}{1+ϵ^2+2ϵ(𝐧𝐫)}$$
(9)
where $`(𝐧𝐫)`$ denotes the scalar product $`𝐧𝐫=n_1r_1+n_2r_2+n_3r_3.`$

#### 2.2.5 Transition probabilities

Given the actual state $`𝐫`$ of the quantum system, and the configuration of the detectors $`\{𝐧[i],i=1,2,\mathrm{},N\}`$ we compute probabilities $`p[i],_{i=1}^Np[i]=1`$ for the $`ith`$ detector to flip. Then we select randomly, with the calculated probability distribution, the flipping detector, and we implement the jump by changing $`𝐫`$ to $`𝐫^{}`$ according to the formula (7), with $`𝐧=𝐧[i]`$. The probabilities $`p[i]`$ are computed using the theory of piecewise deterministic Markov processes applied to the case of quantum measurements - as developed within EEQT.111 It is of interest that Born’s probabilistic interpretation of quantum mechanics as well as the standard formula for quantum mechanical transition probabilities, can be derived in this way and there is no need of adding it as a separate postulate. According to EEQT the probabilities $`p[i]`$ are given by the formula:

$$p[i]=const\text{Tr}\left(P(𝐫)P(𝐧[i],ϵ)^2P(𝐫)\right)$$
(10)
where $`const`$ is the normalizing constant. Using cyclic permutation under the trace, as well as the fact that $`P(𝐫)^2=P(𝐫)`$ we find, taking trace of both sides of the formula (7), that $`p[i]`$ are proportional to $`\lambda (ϵ,𝐧[i],𝐫)`$ given by (8), thus

$$p[i]=\frac{1+ϵ^2+2ϵ(𝐧[i]𝐫)}{N(1+ϵ^2)}.$$
(11)
Note that, owing to the fact that $`_{k=1}^N𝐧[k]=\mathrm{𝟎}`$ we have $`_{i=1}^Np[i]=1`$, as it should be.

## 3 The Five Platonic Fractals

As we noted above, different values of $`ϵ`$ give different fuzziness. To produce representative pictures, one for each solid, we adjusted $`ϵ`$ so that, as a rule, the more vertices, the higher value of $`ϵ`$ \- thus higher resolution of details. While rendering the pictures, for obtaining grayscale value for a given pixel we were using either the formula $`\mathrm{log}(data+1)`$ or, to get more details at the peak values, even $`log(log(data+1)+1).`$

### 3.1 Other Polyhedra

In principle our algorithm should create quantum fractals for each of the regular polyhedra. The only restriction on the array of vectors $`𝐧[i]`$ is that they are all of unit length, and their sum is a zero vector. We added, for comparison with the Platonic solids configurations, two additional simple yet regular figures: double tetrahedron and icosidodecahedron. Notice that tetrahedron is self-dual, while dodecahedron and icosahedron are dual to each other. Double tetrahedron array is obtained by combining $`𝐧[i]`$ with $`𝐧[i]`$ \- that is with the inverted configuration.

Icosidodecahedron has particularly simple and elegant expression for its 30 vertices: they are of the form: $`(\pm 1,0,0)`$ and its cyclic permutations, and $`\frac{1}{2}(\pm 1,\varphi ,\frac{\pm 1}{\varphi })`$ and its cyclic permutations, where $`\varphi =\frac{1+\sqrt{5}}{2}=1.61803\mathrm{}`$ is the golden ratio. All of its edges are of length $`\varphi `$. For its 30 vertices $`ϵ=0.85`$ was needed to resolve the atrractor’s fine structure.

## 4 Notes

#### 4.0.1 Complex projective plane serves as a canvas

Pure states of the quantum spin (we are discussing spin $`1/2`$ here) are described by unit vectors in our Hilbert space $`^2`$. But, as it is standard in quantum theory, proportional vectors describe the same state - the overall phase of the vector has no physical significance. Therefore, in geometrical terms, the set of all pure states is nothing but the projective complex space $`_1()`$ which happens to be the same as the sphere $`S^2.`$ That is why the fractal pattern, in our case, is being drawn on a spherical canvas. There is, however, another possible interpretation of the same algorithm. It is well know that there is an intrinsic relation between Minkowski space and the space of $`2\times 2`$ complex Hermitean matrices. 222In fact, by using Cayley transform, this relation identifies the space of unitary matrices with the compactified Minkowski space. In coordinates the map is given by $`p=\{p^\mu \}p\text{/}p^\mu \sigma _\mu ,`$ so that $`\mathrm{det}(p\text{/})=p^2=p^\mu p_\mu .`$ In particular our projection operators $`P(𝐫)`$, $`𝐫^2=1`$ correspond to null directions. In other words our canvas, the sphere $`S^2`$, can be also thought of as the projective light cone \- the space of light directions. Quantum jumps would then correspond to sudden changes of directions of light or light-like entities. Indeed the operators $`P(𝐧,ϵ)`$ are positive, thus proportional to Lorentz boosts. The formula (9) for jumps implemented by these operators should be compared with the formula for a Lorentz boost, with velocity $`\beta `$ in the direction $`𝐧`$:

$$y^0=x^0\mathrm{cosh}\alpha +(\mathrm{𝐱𝐧})\mathrm{sinh}\alpha ,𝐲=𝐱[(\mathrm{𝐱𝐧})\mathrm{cosh}\alpha +x^0\mathrm{sinh}\alpha ]𝐧,$$

where $`\beta =\mathrm{tanh}\alpha `$ is the velocity. It follows that the jumps described by the equation (9) can indeed be interpreted as Lorentz boosts witht velocity $`\beta =2ϵ/(1+ϵ^2).`$

The formula (7) can be easily generalized for $`P(𝐫)`$ being a generic Hermitean matrix (thus representing a four-vector $`p`$ of space- or time-like character as well. However there is no such generalization for the formula (10), so that the probabilities would have to be assigned equal, and a physical interpretation, even tentative one, is missing in such a case (cf Subsection 4.11 below).

### 4.1 Pure states as projection operators

For a spin $`1/2`$ quantum system pure quantum states are uniquely described by projection operators $`P(𝐧)`$, where $`𝐧`$ is a unit length direction vector, starting at the origin, and ending on one of the points of the unit sphere. All pure states are of this form. Indeed, every orthogonal projection, except of the two trivial ones: $`0`$ and $`I,`$ are of the form $`P(𝐧)`$ for some $`𝐧.`$

### 4.2 Fuzzy projections

Detecting the particle spin is somewhat different than detecting its position components. To measure the position of a particle we can use a photographic plate or a bubble chamber. In such a case a simple mathematical model of a detector is obtained by associating with each active center of the detector a fuzzy, bell-type function with its half-width corresponding to the active region of the center ( for instance about $`1\mu m`$ for an AgBr grain of a photographic emulsion). What would correspond to such a ”fuzzy” projection operator in the case of a spin measurement? We do not have much of a choice here. Due to symmetry reasons there is only one formula possible, namely one given by Eq. (6).

### 4.3 Importance of fuzziness

It is importance to note that in our generalization of the projection postulate, as the result of a jump, not all of the old state is forgotten. The new states depends, to some degree, on the old state. Here EEQT differs in an essential way from the naive von Neumann projection postulate of quantum theory. The parameter $`ϵ`$ becomes important. If $`ϵ=1`$ \- the case where $`P(𝐧,ϵ)=P(𝐧)`$ is a projection operator - the new state, after the jump, is always the same, it does not matter what was the state before the jump. There is no memory of the previous state, no ”learning” is possible, no ”lesson” is taken. This kind of a ”projection postulate” was rightly criticized in physical literature as being contradictory to the real world events, contradicting, for instance, the experiments when we take photographs of elementary particles tracks. But when $`ϵ`$ is just close to the value $`1`$, but smaller than $`1,`$ the contradiction disappears. This has been demonstrated in the cloud chamber model , where particles leave tracks much like in real life, and that happens because the multiplication operator by a Gaussian function does not kill the information about the momentum content of the original wave function. Notice that our fuzzy projection operators $`P(𝐧,ϵ)`$ have the properties similar to those of Gaussian functions, namely333This can also be interpreted as Lorentz formula for addition of relativistic velocities- cf. Section 4.0.1 below.

$$P(𝐧,ϵ)^2=\frac{1+ϵ^2}{2}P(𝐧,\frac{2ϵ}{1+ϵ^2}).$$
(12)
### 4.4 Geometrical meaning of the parameter $`ϵ`$

It is instructive to have a visual picture of the map $`𝐫𝐫^{}`$ of the sphere $`S^2`$ implemented by the operator $`P(𝐧,ϵ).`$ To this end let us assume the vector $`𝐧`$ is pointing North, i.e. $`𝐧=(0,0,1)`$. Then the result $`𝐫^{}`$ of applying the operator $`P(𝐧,ϵ)`$ to a point $`𝐫`$ on the sphere is on the same longitude as the original point $`𝐫`$, but its latitude $`\theta `$ changes - it moves towards the North Pole along its meridian, the new latitude being given by the formula:

$$\theta ^{}=\mathrm{arccos}\left(\frac{(1ϵ^2)\mathrm{cos}(\theta )+2ϵ(1+ϵ\mathrm{cos}(\theta ))}{1+ϵ^2+2ϵ\mathrm{cos}(\theta )}\right)$$
(13)
Remark: Here $`\theta `$ is not exactly the ”geographical latitude”. It is zero at the ”North Pole” ($`\sigma _3=+1`$), $`90`$ degrees at the equator, and $`180`$ degrees at the ”South Pole” ($`\sigma _3=1`$).

Each map $`P(𝐧,ϵ)`$ maps the sphere onto itself in an injective way. For $`𝐧=(0,0,1)`$ the map is easy to picture. All points of the sphere move towards the North Pole along their meridians, except of the two fixed points: North and South Pole. All of the Northern hemisphere, and a strip below the equator, shrinks, while the other part, near the South Pole, stretches. The amount of stretching can be found by plotting the function $`\theta \theta ^{}(\theta )`$ \- it has a maximum at $`\theta `$ corresponding to $`z=ϵ`$. Thus the parameter $`ϵ`$ gets a simple interpretation: it is the value of $`z`$ coordinate for which shrinking of meridians is replaced by stretching - an equilibrium point. This point is always on the southern hemisphere. For $`ϵ`$ close to zero, where the map is close to the identity map, the equilibrium point is close to the equator. Then, as $`ϵ`$ approaches the value of $`1.0`$, corresponding to the sharp projection operator, the equilibrium latitude gets closer and closer to the South Pole. In the limit of $`ϵ=1`$ all of the sphere shrinks to the North Pole, only the South Pole remains where it was.

### 4.5 Quantum Fractals and IFS

Our algorithm is, in fact, a version of a nonlinear iterated function system (IFS). Such algorithms are known to produce complex geometrical structures by repeated application of several non-commuting affine maps. The best known example is the Sierpinski triangle generated by random application of $`3\times 3`$ matrices $`A[i],i=1,2,3`$ to the vector: $`v_0=(x_0,y_0,1)`$ where $`A[i]`$ are given by $`A[i]=((0.5,0,ax_i),(0,0.5,ay_i),(0,0,1))`$ and $`ax_1=1.0,ay_1=1.0,ax_2=1.0,ay_2=0.5,ax_3=0.5,ay_3=1.0.`$ (Our $`3\times 3`$ matrices encode affine transformations - usually separated into a $`2\times 2`$ matrix and a translation vector.) At each step one of the three transformations $`A[i],i=1,2,3`$ is selected with probability $`p[i]=1/3`$. After each transformation the transformed vector is plotted on the $`(x,y)`$ plane. Theoretical papers on IFSs usually assume that the system is hyperbolic that is that each transformation is a contraction, i.e. the distances between points get smaller and smaller. It was shown in  that this assumption can be essentially relaxed when transformations are non-linear and act on a compact space - as is in the case of quantum fractals we are dealing with. But it is not known whether the results of  apply in our case.

### 4.6 Continuous part of the evolution

In our discussion of quantum fractals, we neglect completely the continuous part of the time evolution. Here we are not interested in ”when” jumps happen. We are only interested in the final pattern produced by a long sequence of jumps. Normally spin interacts with a magnetic field (if present), which causes our spin sphere to rotate around the direction of the magnetic field vector. Such a rotation, between jumps, would smear out our pattern. To study the pattern we neglect the rotation part. Once we discard the continuous evolution, the timing of jumps does not really matter, so we neglect this part of EEQT algorithm (timing is very important in simulations of particle detectors, arrival and tunneling times. Here, we simulate jumps as fast as our computer can crunch the numbers.

### 4.7 Nonunitarity

There is one important comment that applies here. Even if we neglect the continuous time evolution due to magnetic field, the very presence of the detectors causes a non-unitary time evolution of the pure state of the quantum system. This evolution is also called, by some physicists (Dicke, Elitzur, Vaidman) , ”interaction-free”. We do not want to enter into this subject here, except for one remark: for symmetric geometric configurations that we are considering here, the EEQT algorithm implies that this, continuous in time, non-unitary evolution can be neglected as well. In fact, it follows from the EEQT model that the ”interaction-free” or, as we call it, ”binamical part” of the evolution is determined by the generator $`\kappa \mathrm{\Lambda }`$, where $`\mathrm{\Lambda }=P(𝐧[i],ϵ)^2.`$ In our case, when our detectors are symmetrically placed, so that $`_{i=1}^N𝐧[i]=\mathrm{𝟎},`$ the formula (12) implies $`\mathrm{\Lambda }=\frac{N(1+ϵ^2)}{4},`$ thus the ”binamical” part is just decreasing the norm of the state vector, while leaving its direction unchanged. Thus it does not affect the geometric pattern of jumps (it is responsible for the mean frequency of jumps, but here timing is not important). For a recent review of EEQT, cf.

### 4.8 Quantum characteristic exponent

Averaging our nonlinear PDP over individual histories one gets a linear Liouville equation for the density matrix of the total system. Tracing over the classical subsystem is, in our case, easily performed and we then get:

$$\dot{\rho }=\kappa (\underset{i=1}{\overset{N}{}}P(𝐧[i],ϵ)\rho P(𝐧[i],ϵ)\frac{1}{2}\{\underset{i=1}{\overset{N}{}}P(𝐧[i],ϵ)^2,\rho \})$$
(14)
where $`\{,\}`$ stands for anti-commutator, and $`\kappa `$ is a coupling constant. For $`\rho `$ written as $`\rho =\frac{1}{2}(I+\sigma (𝐦))`$, $`𝐦^2=m_1^2+m_2^2+m_3^21`$ after some calculations we get a very simple time evolution: $`𝐦(t)=\mathrm{exp}\left(N\kappa ϵ^2t/3\right)𝐦(0).`$ The quantum characteristic exponent, as defined in Ref , is thus $`\frac{2N}{3}\kappa `$ \- not a very useful quantity in our case. The Hausdorff dimension of the limit set, for the tetrahedral case, has been numerically estimated in Ref.  and shown to decrease from 1.44 to 0.49 while $`ϵ`$ increases from 0.75 to 0.95. We hope that by publishing the generating algorithm we will create interest in confirming these as well as obtaining new results in this field.

### 4.9 How to measure the wave function itself

The fractal patterns are produced by a jumping point on the space of pure quantum states - thus are not directly observable. That is why we consider our model as simply a toy to play with. But the model can be developed further on so as to predict observable effects. One way to do it is by adding another layer of detectors, densely spaced, that would detect the fractal pattern. Here we come to the famous question: how to measure the wave function itself. The idea as to how to model such a measurement within the EEQT formalism has been indicated in . We hope to implement the appropriate algorithms at a later time. On the other hand, interpreting the patterns as resulting from Lorentz boosts applied to directions in real space, as discussed in Section 4.0.1, we may try to find similar fractal patterns in the distribution of Galaxies (compare for instance Figure (6) with recent paper on the topology of the universe .

### 4.10 Detectors

A detector is represented by a two-state classical system. It can be in one of the two states, denoted $`0`$ and $`1`$. The fact that it is “classical” means that its two states define two different superselection sectors that can be mixed statistically, but there are no observables that connect these sectors. We will assume that it can ”flip” from 0 to 1 or from 1 to 0 when coupled to a quantum spin. Each flip represents an event; specifically: a detection event. The interpretation is that when the detector flips, the experimental question ”is the spin oriented along the vector $`𝐧`$?” gets an affirmative answer. Note that $`𝐧`$ and $`𝐧`$ are two different experimental questions. They corresponds to two opposite spin directions.

A realistic detector should also exhibit a relaxation time, that is, after each flip it should take some time before it is ready to flip again. We could easily model this phenomenon in our model, but here we are interested in patterns that are created, not in the timing of its appearance .

### 4.11 Hyperbolicity

In a recent paper  Łozinski, Słomczynski and Zyczkowski studied iterated function systems on the space of mixed states, when probabilities that are associated with maps are given independently of the maps. In section III of their paper they give a short discussion dealing with iterated function systems on the space of pure states as well. They start with the following definition of a (pure states) quantum iterated function system(QIFS):

Definition (QIFS) Let $`_N`$ be a complex Hilbert space of dimension $`N.`$ Let $`𝒫_N`$ be the space of one-dimensional subspaces of $`_N`$. Given a unit vector $`\varphi _N`$, let $`P_\varphi `$ be the orthogonal projection onto the subspace spanned by $`\varphi .`$. Specify two sets of $`k`$ linear invertible operators:

* $`V_i:_N_N`$ ($`i=1,\mathrm{},k`$), which generates maps $`F_i:𝒫_N𝒫_N`$ ($`i=1,\mathrm{},k`$ ) by $`\varphi V_i\varphi /V_i\varphi `$ for any $`\varphi 𝒫_N`$, and
* $`W_i:_N_N`$ ($`i=1,\mathrm{},k`$), forming an operational resolution of identity, $`_{i=1}^kW_i^{}W_i=𝕀`$, which generates probabilities $`p_i:𝒫_N[0,1]`$ ($`i=1,\mathrm{},k`$) by $`p_i\left(\varphi \right):=W_i\left(\varphi \right)^2.`$

Comments

1. The authors define the system to be hyperbolic if the maps $`F_i`$ are contractions with respect to the Fubini-Study distance $`d(P_\varphi ,P_\psi )=\mathrm{arccos}\left(\sqrt{Tr(P_\varphi P_\psi )}\right),`$ i.e. there exists constants $`0<L_i<1`$ such that $`d(F_i(P_\varphi ),F_i(P_\psi ))L_id(P\varphi ,P_\psi )`$ for all $`\varphi ,\psi _N.`$ Then they state a proposition (Proposition 1 in ) that guarantees existence of an invariant measure for a hyperbolic system. It seems that the assumptions of this proposition can not be satisfied. It is well known  that a smooth injective map of a compact orientable manifold is automatically a bijection. On the other hand, iterating the map several times if necessary, the distance between any two image points can make smaller than any given positive number, and therefore, a fortiori less than the maximal distance between the points of $`𝒫_N`$, which contradicts surjectivity.444After informing the authors about this problem, they kindly replied that they saw it too, that their Proposition 1 has been deleted from the paper they submitted to Physical Review, and that they are going to replace the electronic version of the paper on the eprint server as well.
2. The authors of  assume that the probabilities $`p_i`$ are given independently of the mapping operators $`V_i.`$ This is not the case in the EEQT scheme. In EEQT $`V_i`$ need not be invertible, but the probabilities $`p_i`$ are determined by the $`V_i`$-s automatically, and in such a way that whenever there is a danger of dividing by zero, the associated probability is automatically zero.

One of us (A.J) would like to thank L.K-J for invaluable help, and to Palle Jorgensen for encouragement and useful comments.