text:
# Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxies

## 1. Introduction

The broad band X-ray spectra (2-300 keV) of Seyfert 1 galaxies are generally well fitted by a cut-off power law continuum with superimposed secondary components like a neutral iron line and a reflection hump. The latter components have been detected in the end of the eighties with GINGA (Piro et al. pir90 (1990); Matsuoka et al. mat90 (1990); Pounds et al. pou90 (1990)) and are considered as signatures of reprocessing of the primary X-ray emission in surrounding cold matter. A high energy cut-off in the X-ray/$`\gamma `$-ray part of the spectrum was first unambiguously detected by SIGMA and OSSE in NGC 4151 (Jourdain et al., jou92 (1992); Maisack et al. mai93 (1993)). A high energy cut-off was also requested to fit the average X-ray/$`\gamma `$-ray spectra (obtained from data of different satellites) of a sample of Seyfert (Zdziarski et al. zdz95 (1995); Gondek et al. gon96 (1996)). The broad band capabilities of the Italian-Dutch satellite BeppoSAX allowed to detect high energy cut-offs in other Seyfert 1 galaxies, about half of the fifteen observed with this satellite. Values of the high energy cut-offs range from 70 keV to $``$ 300 keV (Matt mat00 (2000), hereafter M00).

A cut-off power law plus reflection model (the so-called PEXRAV model in XSPEC, Magdziarz & Zdziarski mag95 (1995)) depends on three parameters: the spectral index $`\mathrm{\Gamma }`$, the high energy cut-off $`E_c`$ and the reflection normalization $`R`$, the formers characterizing the shape of the primary continuum. The physical interpretation of these parameters is generally done in the framework of thermal Comptonization mechanism, which is commonly believed to be the origin of the X-ray emission of Seyfert galaxies. Using approximate relations (cf. section 4), it is then possible to derive, from values of $`\mathrm{\Gamma }`$ and $`E_c`$, values of the temperature $`kT_e`$ and optical depth $`\tau `$ of the “comptonizing” hot plasma (the so-called corona).

If such approximations are sufficient in isotropic geometries (e.g. a central soft photon source surrounded by a spherical hot cloud), where Comptonization process produces roughly power law spectra with a high energy cut-off, strong discrepancies may appear in anisotropic ones, especially for small optical depth and large temperature. Indeed, the (possible) anisotropy of the source of seed photons (as seen by the corona), may introduce anisotropies and modifications of the outgoing spectrum. The largest effect occurs when soft photons are emitted by a plane (disk) below the corona. In this case, photons backscattered towards the disk in the first scattering are necessarily produced in ”head–on collisions” and therefore have an energy gain larger than average, while photons scattered towards the corona (i.e. in the forward direction) have an energy gain smaller than average. As a consequence, the contribution of the first scattering order to the outgoing flux is significantly reduced producing a spectral break (the so-called anisotropy break) in the spectrum. Below the high energy cut–off the spectra is then better described by broken power laws than by a simple power law contrary to what is generally believed (Petrucci et al., pet00 (2000), hereafter P00). These anisotropic effects become however less important for larger optical depth (i.e. $`\tau >1`$) since they are “smoothed” by the larger number of Compton scatterings suffered by the soft photons.

It is also worth noting that the energy of the anisotropy break which is, in practice, close to the peak energy of the third scattering order (Haardt haa93 (1993)), depends on $`kT_e`$ but also on the temperature of the seed soft photon plasma $`kT_{bb}`$. Consequently, if in the case of isotropic geometry the spectral shape is relatively independent on $`kT_{bb}`$, which fixes only the low energy boundary of the emitted spectrum, in anisotropic geometry the shape of the Comptonization spectra depends on it.

In a recent paper (P00), we have applied accurate Comptonization models to fit high quality data of the BeppoSAX long look at NGC 5548. We have underlined the importance of the anisotropic effects and, in particular, we have shown that due to the presence of the anisotropy break, above which the intrinsic continuum is steeper, the temperature of the hot electrons estimated from Comptonization models in anisotropic geometry can be much larger that the one derived fitting the data with a power-law + cut off model.

Besides, we have shown that these different models may provide different (and even opposite) relationships between physical parameters. For instance, still in the case of NGC 5548, the spectral softening of the spectrum, observed during a small flare in the central part of the observation, goes with an increase of the high energy cut–off (i.e. of the corona temperature), in the case of PEXRAV, but with a decrease of the corona temperature when using an anisotropic Comptonization model.

In the same spirit of P00, we extend our application of accurate Comptonization models to a larger sample of objects. We report in this paper the results of this study. The paper is organized as follows. In section 2 we briefly present the main results of the BeppoSAX observations of Seyfert galaxies. We also present the sample of objects we use for the purpose of this paper. We describe the models we used and the fitting procedure in section 3. We also detail in this section the fit of the soft excess/warm absorber (hereafter WA) features observed in each object. We present the results of the fits in section 4 and discuss their physical interpretations in section 5. We then conclude in the last section.

## 2. The data

### 2.1. Overview of the BeppoSAX observations of Seyfert I galaxies

Two major BeppoSAX Core programs were devoted to classical Seyfert 1 galaxies aiming to study the broad band spectra (PI: G.C. Perola) and spectral variability (PI: L. Piro). Thirteen sources have been observed on aggregate so far, since the launch of the satellite in 1996.

Reviews summarizing the main results on these sources have been presented by Piro et al., (pir99 (1999)) and M00. For all objects, the spectral models used by the authors included a power law with exponential cut–off (i.e. $`f_EE^{1\mathrm{\Gamma }}\mathrm{e}^{E/E_\mathrm{c}}`$) characterized by a photon index $`\mathrm{\Gamma }`$ and a e–folding energy $`E_\mathrm{c}`$, a gaussian to fit the neutral Iron line near 6.4 keV and a Compton reflection continuum characterized by a reflection normalization $`R`$. Soft excess and/or WA features were added if required by the data.

The major findings are the followings:

* A high energy cut–off is generally observed, with an $`e`$-folding energy typically in the range 70–300 keV.
* All sources show the presence of a reflection component, so confirming the GINGA and ASCA results (Nandra & Pounds, nan94 (1994)).
* Only two objects (Mkn 509 and Mkn 841) show the presence of a soft excess. This sharply contrasts with the results obtained with EXOSAT (Turner & Pounds tur89 (1989)) or ROSAT (Walter & Fink, wal93 (1993)) where most of the Seyfert 1 observed seemed to show such excess (for a discussion of this point see Piro et al., pir97 (1997)). On the other hand, WA features (mostly the O VII and O VIII edges) are generally observed in the LECS range.
* A correlation between the power law index $`\mathrm{\Gamma }`$ and the amount of reflection component $`R`$ is indicated, in the same sense as found by Zdziarski et al. (zdz99 (1999), hereafter Z99). However, the issue of the reality of this correlation is still open since it may be due to the fact that the two parameters are strongly related in the fitting procedure (see for example Vaughan & Edelson vau00 (2000)).

### 2.2. Our sample

From the sample of Seyfert 1 galaxies observed by BeppoSAX , we have selected observations with good signal-to-noise ratios at high energies as measured with the PDS instrument. We therefore favored objects with hard X-ray spectra. Such observations are best suited for an analysis of the high energy spectra.

The subsample of objects studied in this paper is reported in Table 1, with the name of the different objects, the date of the BeppoSAX observations (some of these objects have been observed several times), the exposure time and the statistical significance in the PDS instrument.
We also present results from the observation of ESO 141-G55, even if its detection in the PDS instrument is relatively weak, since this source was recently observed for the first time by BeppoSAX .

NGC 4151 was observed three times by BeppoSAX , two times in 1996 during the SV-phase (in July and in December), and the last one in January 1999 during A02. However during the observation of July 96, the source experienced a strong increase of the 2-10 keV flux by a factor two on time scales of one day, with significant spectral variability as shown by Piro et al. (2000a ). According to these authors, we choose to study separately the low and high state occuring during this period (and indicated as 96L and 96H in the following). Unfortunately, during the high state, the LECS was off. Only the MECS and PDS data are therefore available for this state. We have not used, in this paper, the observation of December 96 since it does not differ substantially from the low state of July 96.

Concerning Mkn 509, its observation has been split in two parts during the AO2. Since no significant spectral variability was observed between the two data sets (Perola et al. per00 (2000)), we will only refer here to the combined spectrum.

In this paper, we will be concerned with data from three of the four instruments on board BeppoSAX : the Low Energy Concentrator Spectrometer, LECS (Parmar et al. par97 (1997)) covering the 0.15–10 keV range, the Medium Energy Concentrator Spectrometer, MECS (Boella et al. boe97 (1997)) covering the 2–10 keV range and the Phoswich detector system, PDS (Frontera et al. fro97 (1997)) covering the range 13 – 200 keV. Due to some problems that remain in the spectral analysis of LECS data above 4 keV, we shall use only data in the range 0.1–4 keV.

LECS and MECS event files and PDS pha files were downloaded from the BeppoSAX SDC archives. The spectral counts were extracted from a circular region of 4’ and 8’ radius around the source centroid in the MECS and LECS images respectively. We used the data of the three (or two, for observations done after May 1997) MECS units merged together to increase the signal-to-noise ratio. In XSPEC (Arnaud arn96 (1996)), we use the last updated responses (Sep. 1997) and background matrices of each instrument.

## 3. Model fitting

### 3.1. The primary continuum

We fit the data using two different models for the primary continuum: 1) an exponentially cut–off power law plus a reflection component from neutral material (PEXRAV model of XSPEC, Magdziarzk & Zdziardski mag95 (1995)) and 2) a thermal Comptonization spectrum from a disk+corona configuration in slab geometry. The latter was obtained using the code of Haardt (haa94 (1994), hereafter we will refer to this Anisotropic Comptonization Code as AC2). This code derived the angle–dependent spectra of the disk–corona system using an iterative scattering method, where the scattering anisotropy was taken into account only in the first scattering order. It includes also a reflection component described following White, Lightman & Zdziarski (whi88 (1988)) and Lightman & White (lig88 (1988)) and assuming neutral matter. The spectral shape of the reflected photons is averaged over angles. It is multiplied by a first normalization factor which depends on the inclination angle (see Ghisellini, Haardt & Matt ghi94 (1994) for details). In addition, the usual $`R`$ normalization is left free to vary in the fit procedure, so that, for the given inclination angle, $`R=1`$ corresponds to a solid angle, subtended by the reflector, of 2$`\pi `$ . We do not apply Comptonization models for other geometries (e.g. hemisphere) since, as shown for NGC 5548 (P00), the results follow closely those obtained with the slab geometry.

The fit parameters of AC2 are the temperature of the corona $`kT_e`$, its optical depth $`\tau `$, the temperature of the disk $`kT_{bb}`$ (assuming a black body soft emission) and the reflection normalization $`R`$. On the other hand, the PEXRAV continuum depends only on 3 parameters: the e–folding energy of the cut–off power law $`E_\mathrm{c}`$, the photon index $`\mathrm{\Gamma }`$ and the reflection normalization $`R`$.

### 3.2. Soft X-ray excess and Warm Absorber

Most of the objects present WA features and/or a Soft excess in the low part (i.e. in the LECS energy band) of their spectra. Good and realistic fits require to take these features into account. Detailed analyses of these different components are however beyond the scope of this paper, which is instead focused on the high energy continuum of the sources. Then, as a first approximation, we have added simple components (edges, gaussian, blackbody…) to the primary continuum to reproduce the main features present in the soft part of the spectra. We discuss briefly the case of each source in the following.

#### 3.2.1 NGC 4151

The BeppoSAX data of this galaxy have already been analyzed by Piro et al. (2000a , 2000b ), assuming a simple cut–off power law for the primary underlying continuum. A prominent soft excess, known to exist in this object (Holt et al. hol80 (1980)) below the $``$2 keV cutoff due to a high absorbing column, is clearly detected in the observations of 1996 (low state) and 1999. It is generally modeled by allowing some fraction of the central source to be completely uncovered or viewed through a lower column (the so-called “leaky absorber”). In addition, there is evidence of a possibly separate soft excess component that does not vary with the 2–10 keV continuum (Pounds et al., 1986; Perola et al., 1986). Part of this emission comes from an extended X–ray region known to exist in this object (Elvis et al. elv83 (1983); Morse et al. mor95 (1995)). Recent Chandra observations with the High Energy Transmission Grating Spectrometer have revealed detail spectra of this source (Ogle et al., ogl00 (2000)) and allow for the first time precise measurements of the ionization, temperature and kinetics of the extended soft X–ray emission region.

For the scope of the present paper, we model the complex soft X–ray emission in a (relatively) simpler manner. According to Piro et al. (2000a  2000b ), we have used a dual absorber that is the source is covered completely by a medium with a column density $`N_{H1}`$ and with a covering fraction $`f_{cov}`$ by a second medium with a column density $`N_{Hcov}`$. On the other hand, the soft excess below 2 keV was modeled by a scattering component plus an ultrasoft component described by a thermal bremsstrahlung, both absorbed by the galactic column density.

Due to the lack of LECS data for the high state of 1996, we cannot constrain the different components of the soft X-ray range. We thus analyze the data only above 3.5 keV adding a simple photo-electric absorption characterized by a column density $`N_{H1}`$.

We have reported for the three observations our best fit parameters $`N_{H1}`$, $`N_{Hcov}`$ and $`f_{cov}`$, obtained with PEXRAV to model the continuum, in Table 3. We note a significant increase of $`N_{H1}`$ between 96L and 99. As already noted by Piro et al (2000a , 2000b ), such variations suggest that the structure of the dual absorber has slightly changed between the two observations.

#### 3.2.2 NGC 3783

Evidences for deep absorption features in the 0.1-1.5 keV band of NGC 3783 have already been detected in ROSAT (Turner et al. tur93 (1993)) and ASCA observations (George et al. geo95 (1995)), the main ones being the O VII and O VIII edges at 0.74 and 0.87 keV respectively. More recently, in January 2000, NGC 3783 was observed using the High Energy Transmission Grating Spectrometer on the Chandra X-ray observatory (Kaspi et al.,kas00 (2000)). These authors detected a large number of emission and absorption lines, strongly confirming the presence of warm absorbing/emitting media embedding the central X-ray source.

An emission feature near 0.5 keV has also been claimed in ASCA spectra of this source (George et al. geo98 (1998)) and has been interpreted as an O VII emission line (expected at 568 eV). Kaspi et al. (kas00 (2000)) clearly see this emission feature in the Chandra spectra, confirming for the first time the detection of this component.

According to these results, we have included in our models two absorption edges, near 0.74 and 0.87 keV, and a gaussian peaking near 0.5 keV to fit the O VII line. The best fit results are reported in Table 2. As suggested by the recent observation of Chandra, the edge near 1.3–1.4 likely corresponds to contributions by Fe L and Ne K edges whereas the second one, centered near 0.8–0.9 keV, may be due to the O VIII edge.

Concerning the line, we have checked that its addition significantly improves the fit with $`\mathrm{\Delta }\chi ^2`$=16 (according to the F-test, this difference implies a probability less than 0.01 that the line is due to random fluctuations in the counts). This feature is thus detected by BeppoSAX with a high confidence level. The corresponding best fit central energy of the line was $`E_{line}`$0.55 keV with an EW of $``$100 eV. Such EW is large, and it is unlikely that this emission features was only produced by O VII emission lines. In fact, in the more extreme case where we see the blends of the three O VII lines, we expect an EW of only 30 eV (De Rosa, private communication). The large EW and small energy of the line we obtained can be due in part to a bad modeling of the continuum near this component.

#### 3.2.3 Mkn 509

Mkn 509 was observed by BeppoSAX in 1998. The observation was split in two parts during the AO2, the first one being done in May 1998 and the second one in October 1998. As already said, a hardness ratio analysis of the counts shows marginal evidence of spectral variations between the two observations. They have thus been combined and we applied the spectral analysis to the integrated counts.

These data have already been studied in detail by Perola et al. (per00 (2000), hereafter PER00). Assuming a cut–off power law model for the intrinsic continuum of this source, the data show at low energy the presence of additional, superimposed features due to O VII and/or O VIII edges as well as a strong soft excess. PER00 described the latter by a “soft” power law component ($`A_sE^{\mathrm{\Gamma }_s}`$). They then obtained a significant improvement of the goodness of the fit ($`\mathrm{\Delta }\chi ^2`$=20), the best fit values for $`A_s`$ and $`\mathrm{\Gamma }_s`$ being 1.11$`\times 10^2`$ cm-2 sec-1 kev-1 and 2.5 respectively.

PER00 underlined the fact that the addition of this “soft” component had also the interesting effect of providing a more astrophysically sound global solution, with a reflection normalization $`R0.6`$ and Iron line width $`\sigma _{Fe}`$ 0.4 keV (compared to $`R`$2 and $`\sigma _{Fe}`$3 keV without this soft component), giving some confidence on the correctness of the model.

Following these authors, we have added this “soft” power law to our fit, fixing $`A_s`$ and $`\mathrm{\Gamma }_s`$ to the best fit values given above. We also add an edge with the best fit energy and optical depth found by PER00, i.e. $`E_{edge}=0.74`$ keV, $`\tau _{edge}=0.06`$.

#### 3.2.4 ESO 141-G55

The presence of a soft excess was already suggested for this source from ROSAT observations (Turner et al tur93a (1993)). However the BeppoSAX data do not show any evidence of such component nor WA features. It could be partly due to the low S/N of this observation. Consequently, our fits were done without any additional components in the soft band other than the absorption by the galactic column density (cf. Tables 3 and 4).

#### 3.2.5 NGC 5548

The WA observed in this source was already discussed in P00 and Nicastro et al. (nic00 (2000)) and the reader can refer to these papers for more details. For the purpose of this work, we adopt the results of P00 who have added two edges, at $`E_1`$= 0.74 keV and $`E_2`$=0.87 keV, to fit the Oxygen absorption features observed in the LECS band. The corresponding optical depths were fixed to the best fit values obtained with PEXRAV, i.e. $`\tau _1`$=0.51 and $`\tau _2`$=0.12 respectively.

#### 3.2.6 IC 4329A

The X-ray spectrum of IC 4329A is known to display strong neutral absorption, below 3 keV, due to a gas column about ten times greater than the galactic value $`N_{h_{gal}}=4.5\times 10^{20}`$ cm-2 (Elvis et al. elv89 (1989)). This gas is supposed to be associated with the disk of the host galaxy, which is oriented nearly edge-on (Petre et al. pet84 (1984)).

At low energies (below 1-2 keV), ROSAT and ASCA observations have also shown evidence of two rather strong edges, one consistent with O VII, the other with O VIII (Madejski et al. mad95 (1995); Cappi et al., cap96 (1996); Reynold rey97 (1997); George et al. geo98a (1998)).

The BeppoSAX observation of IC 4329A has been already analyzed by Perola et al. (per99 (1999)) using PEXRAV as primary continuum. To model the soft range of the spectrum, they added a neutral column of gas in addition to the galactic one and included two absorption edges. This model gives a good fit to the data and we adopt the same description in our fits. We fix the energy and optical depth of the two edges and the neutral column density to the best fit values obtained by Perola et al. (per99 (1999)), i.e. $`E_1=`$0.73 keV, $`\tau _1`$=0.52, $`E_2=`$1.03 keV, $`\tau _1`$=0.19 and $`N_h`$=0.38$`\times 10^{22}`$ cm-2.

### 3.3. Fitting procedure

Since the sample is only composed of Seyfert 1 galaxies, we expect, in the unification model framework, relatively small inclination angles. In the case of NGC 4151, there are claims of a large inclination angle ($``$60, Evans et al. eva93 (1993)). However this value is rather uncertain. To reduce the length of the fitting procedure, we thus decided to fix, for all the objects of our sample, the inclination to 30 (Nandra et al. nan97 (1997)).

The geometrical normalization of the reflection component $`R`$ was let free to vary. We recall that a value $`R=1`$ corresponds to a covering factor of the cold matter to the X–ray source of $`\mathrm{\Omega }=2\pi `$ whatever the degree of anisotropy of the X–ray emission.

We added a gaussian to reproduce the iron line complex near 6.4 kev. The central energy and width of the gaussian were let free to vary. An Iron edge in the range 7–9 keV was also added if required by the data, that is for NGC 4151, NGC 3783.

We allowed the relative MECS to PDS normalization to vary by 5% around the value of 0.86 to account for the estimated systematic uncertainty (Fiore et al. fio99 (1999)). Analogously, we let the LECS to MECS normalization free to vary over the range of acceptable values 0.7–1 (Fiore et al. fio99 (1999)).

Finally, in the case of PEXRAV, we fixed the parameters of the soft excess and/or WA components (i.e. edges, line, soft power law) to the values presented in section 3.2. For fits with AC2, we also fixed, in a first step, the column density to the best fit values obtained with PEXRAV. A second series of fits were obtained with a column density free to vary. We checked that there were no significant differences between the two sets of results (we generally obtain slightly larger column densities and changes of the soft temperature $`kT_{bb}`$ when $`N_h`$ is let free to vary). We will thus only consider the best fit results with free column density in the following.

For computing errors on the parameters of the primary continuum (that is $`kT_e`$, $`\tau `$, $`kT_{bb}`$ and $`R`$ for AC2, and $`E_c`$, $`\mathrm{\Gamma }`$ and $`R`$ for PEXRAV), we let the column density free to vary but we fixed the LECS/MECS and MECS/PDS normalization, and the Iron line and/or edges to their best fit values. Throughout the paper, errors on single parameter are quoted at a confidence level of 90 % (i.e. $`\mathrm{\Delta }\chi ^2`$=2.7) unless otherwise specified.

## 4. Results

### 4.1.  The PEXRAV model

The best fit parameter values obtained with the fitting procedure described above using the PEXRAV model are reported in Table 3. We checked, as far as possible, our PEXRAV results with those already published on these data. For IC 4329A, Mkn 509, NGC 4151, NGC 3783 and NGC 5548 our PEXRAV fits are in good agreement with those obtained by Perola et al. (per99 (1999)), Perola et al. (per00 (2000)), Piro et al. (2000a , 2000b ), De Rosa et al. (der00 (2000)) and Nicastro et al. (nic00 (2000)) respectively. In the case of NGC 5548, we also checked that a better modelization of the WA (with the CLOUDY code used by Nicastro et al. nic00 (2000)) does not affect the results on the continuum. We note also that, in the case of NGC 4151, we found smaller errors for the spectral index $`\mathrm{\Gamma }`$ ($`\mathrm{\Delta }\mathrm{\Gamma }\pm 0.01`$) in comparison to Piro et al. (2000a , where they found $`\mathrm{\Delta }\mathrm{\Gamma }\pm 0.1`$). This is simply due to the fact that the latters have computed the errors using the MECS and PDS data (above 3.5 keV) whereas we have used the LECS, MECS and PDS ones together (thus adding the soft excess parameters in the fitting procedure).

For the PEXRAV model the parameters of the primary continuum determined by the fits are $`\mathrm{\Gamma }`$ and $`E_\mathrm{c}`$. The actual physical parameters of the Comptonizing region (also reported in Table 3) have been obtained from the spectral parameters in the following way. The temperature $`kT_e`$ is simply estimated as $`kT_eE_\mathrm{c}/2`$, keeping in mind that such approximation roughly holds for $`\tau <1`$. For $`\tau 1`$, $`kT_eE_\mathrm{c}/3`$ would be more appropriate. Thus the reported values can be considered as upper limits to the temperature. Knowing the temperature, the optical depth can be computed from the spectral index derived from the PEXRAV fit using the following relation (Shapiro, Lightman & Eardley sha76 (1976); Sunyaev & Titarchuk sun80 (1980); Lightman & Zdziarski lig87 (1987)):

$$\mathrm{\Gamma }1\left[\frac{9}{4}+\frac{m_ec^2}{kT_e\tau (1+\tau /3)}\right]^{1/2}\frac{3}{2}$$
(1)
This equation is valid for $`\tau >1`$, and we have checked a posteriori that such condition is roughly matched in all cases.

### 4.2.  The slab model

The fits with the AC2 yield directly and consistently the temperature $`kT_e`$ and optical depth $`\tau `$, which was our first motivation for starting this approach. The best fit parameters are reported in Table 4. It is worth noting that another parameter enters the fitting procedure in an important way, that is the temperature of the soft photons, $`kT_{bb}`$. Its values (see Table 4) range, in most cases, between 10 and 40 eV and are relatively well constrained with errors of the order of 5-50% whereas we have only data above 100 eV. For large values of $`kT_{bb}`$ (i.e. $`kT_{bb}>`$30 eV), the high energy tail of the black body shape may be detected in the low energy part of the BeppoSAX data and $`kT_{bb}`$ can then be roughly constrained. We recall that the intensity of the seed photon component with respect to the comptonized one is also fixed by the fitted parameters (mainly $`\tau `$). When no tail is observed, the soft-X-ray data impose, at least, an upper limit of $`kT_{bb}`$. For anisotropic geometries, the high energy part of the Comptonization spectrum also depends on $`kT_{bb}`$, through the location of the anisotropy break. The fitting procedure has thus to adjust $`kT_{bb}`$, together with $`\tau `$, $`kT_e`$ and $`R`$, to reproduce the soft and hard X-ray data resulting in a complex interdependence in the fitting procedure, between the soft and hard part of the spectra.

It is worth noting that, in the case of NGC 4151 96L and 96H, the best fit values of $`kT_{bb}`$ are relatively large (in fact we have only lower limit, the AC2 working only for $`kT_{bb}<100`$ eV). It is true however that in the case of the high state 96H, we don’t have any LECS data. Consequently, as explained above, $`kT_{bb}`$ is not constrained by the soft part of the spectra and large values may be reached. But it is surprising that in the low state 96L, we still find large $`kT_{bb}`$, at least 3 times larger than in 1999. A possible explanation is that, in the observation of 1996, the presence of a larger reflection component (also found with the PEXRAV fits) may hide the anisotropy break, thus diminishing the constraints on $`kT_{bb}`$, and allowing larger values.

### 4.3. Comparison between the slab model and PEXRAV

A purely statistical comparison between the two models is not possible since both give acceptable values of $`\chi ^2`$. The slab model has one more parameter and gives generally lower values of $`\chi ^2`$, except for NGC 4151. However, we have checked, for this source, that the best fit $`\chi ^2`$ probability are similar for the two models. Besides, the two best fits have residuals without significant features and differences between each other. The two models give thus an equivalent representation of the data.

#### 4.3.1 The temperature and optical depth of the corona

For all objects, the estimated corona temperatures (respectively optical depths) from PEXRAV are substantially smaller (respectively larger) than those inferred with an anisotropic Comptonization model. This trend is similar to the the result obtained for NGC 5548 by P00. In Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesa we have plotted the corona temperatures obtained with AC2 versus the temperatures obtained with PEXRAV. Large differences (up to a factor 8) are found between the two estimates. Note however that for IC 4329A, for which the high energy data are comparatively better than for the rest of the sample, the two models give results which are quite close to each other.

The reason of these differences is relatively simple. Within PEXRAV type models, the slope of the power law, determined with small errors by the LECS and MECS data, cannot change, by hypothesis, at higher energies. A cut–off around 100 keV is then required to fit the PDS data. In Comptonization models, the LECS and MECS data determine the slope below the anisotropy break. Above this break the intrinsic spectrum is steeper and thus can fit the PDS data without an additional steepening beyond 100 keV, allowing for a larger temperature and (consequently) a smaller optical depth to keep, roughly, the same power law slope.

#### 4.3.2 The reflection component

The reflection normalizations obtained with AC2 cluster around 1, except for two extreme states of NGC 4151. They are, in all cases but one (ES0 141-G55), larger than those found using the simple cut–off power-law model. We have reported in Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesb the amplitude of the reflection component derived from both models. In some cases we obtain differences of factors 4–5. Note that with the exception of NGC 4151 for all other sources the slab model yields values of $`R`$ consistent with 1.

The case of NGC 4151 is peculiar, showing a large variation in $`R`$ between the 96 and 99 observations and an indication of a variation of $`R`$ also between the low and high states of the July 96 observation. These changes are present irrespective of the model used in the spectral fitting (see also Piro et al. 2000a , 2000b ). However the slab model yields higher absolute values, allowing for some reflection also in the 99 observation for which PEXRAV yields $`R=0.01\pm 0.1`$. According to the slab model, the corona temperature is very high in the 99 observation, which may offer a physical connection with the low value of $`R`$ as dicussed in section 5.1.
We can also remark that the relatively small reflection component and the hardness of the spectrum of this object may be difficult to reconcile with simple corona geometries (like a plane or a patchy corona for example). More complex ones (as dynamics coronae for example, Malzac et al. mal00 (2000)) may be needed. Another possibility may be that the reflecting material in NGC 4151 is highly ionized. In this case, fitting with models which do not take into account the complex ionization pattern of the reflector can severely underestimate the reflection normalization (Ballantyne, Ross & Fabian bal00 (2000); Done & Nayakshin don01 (2001); .

As an illustration of these results, we have plotted in Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesa and Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesb the unfolded best fit spectra of IC4329A and NGC 4151 (observation of 1999) obtained with PEXRAV (dashed lines) and AC2 (solid lines). We have also plotted in each case the reflection hump component to better estimate its contribution to the total spectrum. In the case of IC 4329A, the two models give roughly the same result at high energy. For NGC 4151 however, large differences are expected between PEXRAV and AC2, with a factor of 2–3 at 300 keV and larger than 10 above 500 keV.

## 5. Discussion

A Comptonizing corona above a passive optically thick layer (approximating an accretion disk with negligible internal heating) is a complex system, in which the soft seed photons for the Comptonization process are produced by reprocessing of the Comptonized emission in the underlying passive layer. Such a corona is ”self-cooled” and, for a given geometrical configuration of the system, can only be in equilibrium if the temperature and optical depth satisfy a definite relation (Haardt & Maraschi haa91 (1991)). These relations which essentially correspond to roughly constant Compton parameters $`y4\left({\displaystyle \frac{kT_e}{m_ec^2}}\right)\left[1+4\left({\displaystyle \frac{kT_e}{m_ec^2}}\right)\right]\tau (1+\tau )`$, with $`y`$ 0.5 and 2 for the slab and hemisphere respectively, have been accurately computed by Stern et al. (ste95 (1995); see also for a review Svensson sve96 (1996)).

In Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesa and Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesb the values of $`\tau `$ vs. $`kT_e`$, obtained for the 6 sources considered here using the spectral models AC2 and PEXRAV, are compared with the theoretical relations expected for a Comptonizing region in energy balance for a plane (filled rectangles, solid line) and hemispherical (filled hemispheres, dashed line) configuration respectively.

The best fit results obtained with AC2 (Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesa) are relatively close to the theoretical expectations for the slab case even if they tend to fall preferentially above the solid line. The data therefore indicate a Comptonization parameter larger than for a pure slab geometry, that is a more “photon-starved” configuration. A similar result was already derived by P00 for NGC 5548. For the latter source, best fits with a hemispherical Comptonisation model were also perfomed. The derived coronal parameters were in fact in good agreement with the theoretical predictions for the hemisphere. However an implausibly large value of the reflection component (of the order of 2) was required, suggesting that the real configuration is more complex than these two ideal cases. For the present sample we did not perform fits with hemispherical models.

The values of $`kT_e`$ and $`\tau `$ obtained with PEXRAV (cf. Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesb) show the same trend as found with AC2 that is larger values of $`\tau `$ for smaller values of $`kT_e`$. They also fall above the theoretical expectation for a slab and even above the theoretical expectation for a hemisphere (the theoretical curves of Stern et al. ste95 (1995) have been extended to large $`\tau `$ assuming a constant value of the Compton parameter for each curve i.e. y $``$0.5 and 2 for the slab and hemispherical configuration respectively). Therefore this set of parameters suggests, for each source of our sample, a configuration more “photon-starved” than a hemisphere.

An important difference with respect to the slab model is that the optical depths derived from PEXRAV fits are generally larger than 1. The corona should then be optically thick reducing or cancelling the effects of anisotropy. Therefore also the PEXRAV model has an internal consistency. However a corona with large optical depth may wash out discrete features from the underlying disk (e.g. lines and reflection itself, Petrucci et al. pet01b (2001)) more than would be desirable. This problem could be alleviated if the corona was ”patchy” (and possibly dynamic cf. section 4.3.2).

We conclude that both a hot, optically thin corona with significant anisotropic effects and a less hot, optically thick, patchy corona with negligible anisotropy are consistent with the available data.

### 5.1. Correlations between physical parameters

A correlation between the reflection normalization $`R`$ and the photon index $`\mathrm{\Gamma }`$ has been claimed by Z99 from the study, using PEXRAV to model the primary continuum, of a large number of GINGA observations of Seyfert and galactic black hole objects. Yet, the validity of this correlation is under debate since $`R`$ and $`\mathrm{\Gamma }`$ are strongly correlated in the fitting procedure (Vaughan & Edelson vau00 (2000)). It seems however more significant and less dispute in the case of the galactic black hole objects where the errors on both parameters are smaller.

We have plotted in Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxies the reflection normalization $`R`$ versus the photon index $`\mathrm{\Gamma }`$ we obtained with PEXRAV. In our case, we do not find a clear correlation between the two parameters. A Spearman rank-order correlation test gives inconclusive results with $`r_s`$0.3. However our sample is biased in favor of objects with hard spectra, to ensure a good detection in the PDS instrument. When the 13 objects actually observed by BeppoSAX are taken into account a stronger correlation is observed (M00), in agreement with Z99.

In the case of the AC2 model we have no simple way of characterizing the spectral shape. Also results are available for only 6 sources and we cannot include the other Seyfert 1s observed by BeppoSAX. Therefore a direct comparison with the Z99 correlation is not possible. Below, however, we examine the correlation of $`R`$ with the corona temperature which can be computed for the two models.

#### 5.1.1  Relation between $`R`$ or $`\mathrm{\Gamma }`$ and the temperature of the corona

Z99 interpret the correlation between $`R`$ and $`\mathrm{\Gamma }`$ in the framework of thermal reprocessing models where $`R`$ is directly proportional to the solid angle subtended by the cold matter surrounding the corona. The latter solid angle is not necessarily proportional to the soft photon flux reentering the corona, but it is likely that the larger is $`R`$, the larger is the ratio of soft luminosity to hard luminosity in the corona, i.e. the smaller is the Comptonization parameter $`y`$, resulting in softer spectra (larger $`\mathrm{\Gamma }`$).

It is interesting that the BeppoSAX data can provide the temperatures associated with sources with different values of $`R`$ or $`\mathrm{\Gamma }`$, which were not available from the GINGA data used by Z99. The result is shown in Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxies where we plot $`R`$ vs. the corona temperature $`E_c/2`$ obtained with PEXRAV. Here again, in order to have a larger sample, we added the best fit values obtained by M00 for the other Seyfert 1s observed with BeppoSAX. This plot suggests a positive correlation between the two parameters, that is, the temperature of the corona is larger for larger values of the reflection normalization. In fact a plot of $`E_C/2`$ vs. $`\mathrm{\Gamma }`$ shows that steeper spectra correspond to higher temperatures and (as a corollary) to lower optical depths (cf. Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesa and Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesb).

The slab model analysis yields different trends. A plot of $`R`$ vs. $`kT_e`$ for the latter model shows that sources with larger $`R`$ tend to have lower $`kT_e`$ (cf. Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxies), a behaviour opposite to what is found with PEXRAV. The suggested correlation depends however strongly on the different states of NGC 4151.

To summarize, an analysis in terms of PEXRAV indicates the following correlations: larger reflection component (steeper $`\mathrm{\Gamma }`$) – larger temperature (cf. Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxies) while the analysis in terms of anisotropic Comptonization yields: larger reflection – lower temperature (cf. Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxies).

The latter behaviour is naturally expected in pair free Compton cooled coronae (see Fig Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxiesa) since the transition to lower Compton parameter (i.e. higher $`l_s`$/$`l_h`$ as suggested by a larger $`R`$) implies a decrease in the temperature for constant optical depth, as (relatively well) verified by our sources.

The behaviour indicated by PEXRAV could instead be understood better in the case of a pair dominated corona in pair equilibrium. In such a case, for constant hard compactness, the larger the ratio $`l_s/l_h`$, the larger the temperature of the corona. This behavior is simply imposed by the pair equilibrium. Indeed, an increase of the cooling corresponds to a decrease of the number of particles in the hard tail of the thermal particle distribution. To reach a new equilibrium the temperature must increase. (Zdziarski zdz85 (1985); Svensson sve82 (1982); Ghisellini & Haardt ghi94b (1994) ; Coppi cop99 (1999)).

A possible limitation of the “pair dominated” interpretation comes from the high values of the compactness $`l_h`$ required by the physical parameters determined with PEXRAV. We have reported on Fig. Testing Comptonization models using BeppoSAX observations of Seyfert 1 galaxies the $`\tau `$$`T_e`$ relations expected, for a pair dominated corona, for different values of $`l_h`$, but varying the ratio $`l_h/l_s`$. We have also included the observations of M00. The data require values of $`l_h`$ in the range 100–1000. Given the (unabsorbed) luminosity of the sources of our sample ($`L_{0.1200keV}`$ of the order of $`10^{44}`$ erg.cm-2.s-1), this implies black hole mass upper limits of the order of $`10^7`$ solar masses (for $`l_h`$=100). This is relatively smaller than (but roughly consistent with) black hole mass estimations obtained from reverberation mapping or resolved kinematics methods for most of the objects of our sample (see Gebhardt et al. geb00 (2000) and Nelson nel00 (2000) and references therein). The worse case is that of NGC 4151, where the PEXRAV best fit values of $`kT_e`$ and $`\tau `$ require very large hard compactness (larger than 1000), implying a black hole mass of $`10^5`$ M (reverberation mapping methods estimate a black hole mass of $`10^7M_{\mathrm{}}`$, Clavel et al. cla87 (1987); Ulrich & Horne ulr96 (1996)). In this case, however, the complexity of the soft part of the spectrum may lead to an uncorrect modelisation of the continuum. We therefore conclude that the pair dominated coronae interpretation is also still valid for this object. Other observations (with better constraints on the compactness and the physical parameter of the corona) are clearly needed to confirm or disprove these results.

## 6. Conclusion

The aim of this paper was to test anisotropic Comptonization models over the high signal to noise BeppoSAX observations of a sample of six Seyfert 1 (seven observations). We use two types of model: a detailed Comptonization code in slab geometry, which treats carefully the anisotropy effects in Compton processes, and a simple cut–off power law plus reflection model (PEXRAV model of XSPEC). The latter is generally used as a zero order approximation to Comptonization spectra. If this is a relatively good approximation for isotropic geometry, it fails to reproduce the real shape of Comptonization spectra in anisotropic ones (like a plane, hemispherical or spherical hot region above a flat disk). In this case, due to the deficiency of the first order scattering component towards the observer, it better resembles a broken power law with convex curvature. This leads to strong differences between the best fit parameter values obtained with these two models. The main results of this work are the following:

* The data are well fitted by both models and there is no statistical evidence for a model to be better than the other. Both models give results in agreement with a X-ray source geometry more “photon starved” than the slab case.
* Our best fits with Comptonization models in slab geometry give a temperature generally much larger and an optical depth much smaller than derived from the power law + cut–off fits, using standard Comptonization formulae. The estimate of the reflection normalization is also larger with the slab corona model.
* The two models also lead to different relationships between physical parameters. For instance, PEXRAV indicates a correlation between the reflection normalization and the corona temperature whereas the slab corona model suggests an anticorrelation. This has major consequences for the physical interpretation of the data. Thus, the PEXRAV results suggest that the hot corona may be pair dominated whereas the slab corona model is in better agreement with a low pair density solution.

These results should have observational consequences. First, the large differences in the fitted temperatures lead to widely different predictions as to the fluxes emitted at higher energies, above the PDS range. This underlines the need of better soft $`\gamma `$–ray observations, as those possibly provided by INTEGRAL, in order to directly confirm or disprove the temperatures inferred with anisotropic Comptonization models.

On the other hand, the large values of $`\tau `$ predicted by PEXRAV would spread out the reflection hump in the entire X-ray range, avoiding any detection of this component. The fact that we obtain not negligible values of $`R`$ with PEXRAV may however be reconciled with the large values of $`\tau `$ if the corona is “patchy”. Indeed, in this case, reflection could come from the uncovered part of the cold, reflecting material. Part of the reflection could also come from the torus. The variability of the reflection component may allow to differentiate these two possibilities. Large values of $`\tau `$ would also modify the profile of the Fe line component produced in the disk, leading to an attenuation and broadening of the line. Forthcoming observations with CHANDRA and XMM-Newton are expected to bring substantial progress on this issue.

Acknowledgements: POP acknowledges a grant of the European Commission under contract number ERBFMRX-CT98-0195 (TMR network ”Accretion onto black holes, compact stars and protostars”). This work was partially supported by the Italian MURST trough the grant COFIN98-02-15-41 (JM) and by the Agenzia Spaziale Italiana (ASI) trough the grant ASI-ARS-99-74 (LM,FH).

Figure Captions

Figures


format: text

text:
# IMSc/2001/04/13 UWThPh-2001-14 IFIC/01-18 The solar LMA neutrino oscillation solution in the Zee model

## 1 Introduction

The results of the atmospheric  and solar neutrino experiments (for a recent presentation of the solar results see, e.g., Ref. ) indicate that neutrinos have non-zero masses. A suitable environment for obtaining Majorana neutrino masses is to extend the Higgs sector of the Standard Model . Extending the Standard Model with a singly charged gauge singlet scalar and adding a second Higgs doublet allows to write down an explicitly lepton number-violating interaction in the Higgs potential and leads to 1-loop neutrino masses . In the following we will only consider the lepton sector of the Zee model.

The Zee model is traded in two versions in the literature: the original model , which we will call general Zee model (GZM) in this paper, and a simpler version of the Zee model where only one of the two Higgs doublets couples in the lepton sector . The latter version, which we will call restricted Zee model (RZM), can naturally be achieved with a discrete symmetry ; it has the advantage that the family lepton number-violating interactions mediated by the couplings of the three physical neutral scalars are absent. The interesting point is that the RZM leads to a symmetric neutrino mass matrix with zeros in the diagonal444We always work in a basis where the charged lepton mass matrix is diagonal.; this mass matrix is called *Zee mass matrix* in the literature. One can easily check that all phases in the Zee mass matrix can be absorbed into the left-handed neutrino fields and, therefore, in neutrino oscillations no CP violation is observable if the Zee mass matrix is the correct neutrino mass matrix.

The Zee mass matrix has a special feature: it has been shown  that it allows to accommodate only bimaximal mixing, i.e., solar and atmospheric mixing angles are both very close to $`45^{}`$. While this is in excellent agreement with the atmospheric neutrino data, results from the solar neutrino experiments are not that well compatible with a solar mixing angle $`\theta 45^{}`$, though this value is also not excluded . The purpose of the present paper is to show that in the GZM it is no problem to accommodate the large mixing angle (LMA) MSW solution of the solar neutrino deficit, while at the same time all constraints on the additional couplings in the lepton sector stemming from the second Higgs doublet are respected. In the GZM also the diagonal elements of the neutrino mass matrix are non-zero in general. Since the GZM is a quite rich and intricate model, we restrict ourselves rather to an “existence proof” of the LMA MSW solution within the GZM instead of discussing the GZM in full generality.

Let us outline our procedure:

* We assume that $`U_{e3}=0`$ ($`U`$ is the neutrino mixing matrix) and the atmospheric mixing angle is exactly $`45^{}`$ . Then in the mixing matrix the only free parameter is the solar mixing angle $`\theta `$. This gives a certain form of the mass matrix with 4 complex parameters (Section 2).
* The general problem is still quite intricate, so we set one of these parameters of the mass matrix equal to zero and assume that the remaining ones are real. In this scenario we can relate in a simple way the three real parameters with the physical quantities $`\theta `$, the solar mixing angle, $`\mathrm{\Delta }m_{}^2=|m_1^2m_2^2|`$, the solar mass-squared difference, $`\mathrm{\Delta }m_{\mathrm{atm}}^2=|m_3^2(m_1^2+m_2^2)/2|`$, the atmospheric mass-squared difference and $`m_3`$ ($`m_j`$ with $`j=1,2,3`$ are the neutrino masses) and it is possible to have the LMA angle solution of the solar neutrino problem555This means that $`\theta `$ is large but safely smaller than $`45^{}`$. (Section 3).
* Now the GZM is brought into play. After a discussion of the neutrino mass matrix in this model (Section 4) we assume that all quantities in the model are real and we set to zero all but two of the additional Yukawa couplings present in the GZM; this has the purpose of avoiding as many family lepton number-violating neutral scalar interactions a possible. We show that with these two additional coupling constants the restricted mass matrix mentioned in Point 2, which allows for the LMA MSW solution, can be accommodated (Section 4).
* We complete our procedure by a numerical discussion of the parameters of our scenario and by estimates of the rates of (family) lepton number-violating processes (Section 6).

We also review the features of the RZM as a limit of the GZM (Section 7) and summarize the results (Section 8). In the appendix we present the general formulas for the 1-loop Majorana neutrino mass matrix induced by charged scalar loops.

## 2 Neutrino mixing and the mass matrix

The Majorana neutrino mass matrix $`_\nu `$ is diagonalized with a unitary matrix $`V`$ by

$$V^T_\nu V=\widehat{m}=\mathrm{diag}(m_1,m_2,m_3).$$
(1)
With the assumptions mentioned in the introduction in Point 1, we can write the matrix $`V`$ as

$$V=e^{i\widehat{\alpha }}Ue^{i\widehat{\beta }}\text{with}U=\left(\begin{array}{ccc}\hfill \mathrm{cos}\theta & \hfill \mathrm{sin}\theta & 0\\ \hfill \mathrm{sin}\theta /\sqrt{2}& \hfill \mathrm{cos}\theta /\sqrt{2}& 1/\sqrt{2}\\ \hfill \mathrm{sin}\theta /\sqrt{2}& \hfill \mathrm{cos}\theta /\sqrt{2}& 1/\sqrt{2}\end{array}\right).$$
(2)
The phases in $`V`$ suggest the definitions

$$_\nu ^{}=U\widehat{\mu }U^T\text{with}_\nu ^{}e^{i\widehat{\alpha }}_\nu e^{i\widehat{\alpha }}\text{and}\widehat{\mu }=e^{2i\widehat{\beta }}\widehat{m}.$$
(3)
Then $`_\nu ^{}`$ can be expressed by $`\widehat{\mu }`$ and the parameters of $`U`$ in the following way:

$$_\nu ^{}=\left(\begin{array}{ccc}c^2\mu _1+s^2\mu _2& cs(\mu _1\mu _2)/\sqrt{2}& cs(\mu _1\mu _2)/\sqrt{2}\\ cs(\mu _1\mu _2)/\sqrt{2}& (s^2\mu _1+c^2\mu _2+\mu _3)/2& (s^2\mu _1c^2\mu _2+\mu _3)/2\\ cs(\mu _1\mu _2)/\sqrt{2}& (s^2\mu _1c^2\mu _2+\mu _3)/2& (s^2\mu _1+c^2\mu _2+\mu _3)/2\end{array}\right)$$
(4)
with $`c\mathrm{cos}\theta `$ and $`s\mathrm{sin}\theta `$ (see, e.g., Refs. ). Thus, the mass matrix has the structure

$$_\nu ^{}=\left(\begin{array}{ccc}\hfill a& \hfill b& \hfill b\\ \hfill b& \hfill c& \hfill d\\ \hfill b& \hfill d& \hfill c\end{array}\right).$$
(5)
Consequently, for having $`U_{e3}=0`$ and an atmospheric mixing angle of exactly $`45^{}`$ the following conditions on $`_\nu ^{}`$ are necessary:

Condition 1: $`_{\nu e\mu }^{}+_{\nu e\tau }^{}=0,`$ (6a)
Condition 2: $`_{\nu \mu \mu }^{}=_{\nu \tau \tau }^{},`$ (6b)
Condition 3: $`\mathrm{tan}2\theta =2\sqrt{2}{\displaystyle \frac{b}{ac+d}}𝐑.`$ (6c)
Finally, with the parameterization (5) the complex masses $`\mu _j`$ are found as

$`\mu _1`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(a+cd\pm \left[(ac+d)^2+8b^2\right]^{1/2}\right),`$ (7a)
$`\mu _2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(a+cd\left[(ac+d)^2+8b^2\right]^{1/2}\right),`$ (7b)
$`\mu _3`$ $`=`$ $`c+d.`$ (7c)
## 3 A simplified mass matrix

Having discussed the general form of the mass matrix which leads to the mixing matrix (2), we now investigate the consequences of the following simplifying assumptions:

$$a,b,d𝐑\text{and}c=0.$$
(8)
In the next section, this scenario will be reproduced in the framework of the Zee model. With the reality assumptions, the quantities $`\mu _j`$ (3) are identical with the neutrinos masses apart from possible signs. The experimentally accessible quantities are expressed as

$`\mathrm{tan}2\theta `$ $`=`$ $`2\sqrt{2}{\displaystyle \frac{b}{a+d}},`$ (9a)
$`\mathrm{\Delta }m_{}^2`$ $`=`$ $`|ad|\left[(a+d)^2+8b^2\right]^{1/2},`$ (9b)
$`\mathrm{\Delta }m_{\mathrm{atm}}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}(a^2d^2)+2b^2,`$ (9c)
$`m_3`$ $`=`$ $`|d|`$ (9d)
by the parameters $`a`$, $`b`$, $`d`$. We have chosen $`m_3`$ as representative of the absolute neutrino mass values, since it is simply given by Eq. (9d). Without loss of generality we will adopt henceforth the following conventions: $`0^{}\theta 90^{}`$, $`m_1<m_2`$, $`\mu _3=m_3`$. It follows from the last relation and from Eq. (7c) that $`d`$ is positive. Note that in Eq. (9c) no absolute value of the right-hand side of the equation is necessary, because it must be positive. The argument goes as follows. Suppose that $`(a^2d^2)/2+2b^2=\mathrm{\Delta }m_{\mathrm{atm}}^2`$. Then it follows that $`d^2a^22\mathrm{\Delta }m_{\mathrm{atm}}^2`$. Therefore, $`d^2a^2`$ is positive, which allows to derive the inequality $`\mathrm{\Delta }m_{}^2d^2a^22\mathrm{\Delta }m_{\mathrm{atm}}^2`$ from Eq. (9b). This is a contradiction to the values of the mass-squared differences, fitted from the data .

In Eqs. (9a), (9b), (9c) and (9d), four physical quantities are expressed by three parameters. Therefore, a consistency condition exists, which is given by

$$\mathrm{\Delta }m_{\mathrm{atm}}^2=\frac{1}{2}\eta \mathrm{\Delta }m_{}^2|\mathrm{cos}2\theta |+\frac{1}{4}\left(\sqrt{m_3^2+\eta \mathrm{\Delta }m_{}^2|\mathrm{cos}2\theta |}+\eta ^{}m_3\right)^2\mathrm{tan}^22\theta .$$
(10)
The signs $`\eta `$ and $`\eta ^{}`$ occurring in this equation are $`\eta =\text{sign}(a^2d^2)`$ and $`\eta ^{}=\text{sign}a`$. In the context of the Zee model we will finally need the relations

$`a^2`$ $`=`$ $`m_3^2+\eta \mathrm{\Delta }m_{}^2|\mathrm{cos}2\theta |,`$ (11a)
$`b^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }m_{\mathrm{atm}}^2{\displaystyle \frac{1}{4}}\eta \mathrm{\Delta }m_{}^2|\mathrm{cos}2\theta |.`$ (11b)
Looking at the consistency condition (10) and assuming that $`m_3^2`$ is of the order of $`\mathrm{\Delta }m_{}^2|\mathrm{cos}2\theta |`$ or smaller, we obtain $`\mathrm{sin}^22\theta /|\mathrm{cos}2\theta |\mathrm{\Delta }m_{\mathrm{atm}}^2/\mathrm{\Delta }m_{}^2`$, which amounts to bimaximal mixing. Since we want to show that the Zee model allows to avoid bimaximal mixing we concentrate on

$$m_3^2\mathrm{\Delta }m_{}^2|\mathrm{cos}2\theta |.$$
(12)
With this assumption it is easy to obtain an approximate expression for $`m_3^2`$. One can check that for $`\eta ^{}=1`$ one arrives again at bimaximal mixing. Using $`\eta ^{}=1`$, we have $`a>0`$ and we easily calculate

$$m_3^2=\mathrm{\Delta }m_{\mathrm{atm}}^2\left\{\frac{1}{\mathrm{tan}^22\theta }\eta \frac{1}{2}\frac{\mathrm{\Delta }m_{}^2}{\mathrm{\Delta }m_{\mathrm{atm}}^2}\frac{|\mathrm{cos}2\theta |}{\mathrm{sin}^22\theta }+\frac{1}{16}\left(\frac{\mathrm{\Delta }m_{}^2}{\mathrm{\Delta }m_{\mathrm{atm}}^2}\right)^2\mathrm{sin}^22\theta +\mathrm{}\right\}.$$
(13)
Using this equation the condition (12) implies

$$|\mathrm{cos}2\theta |\frac{\mathrm{\Delta }m_{}^2}{\mathrm{\Delta }m_{\mathrm{atm}}^2}.$$
(14)
Equation (13) together with $`a>0`$ and $`d>0`$, allows to estimate from Eq. (11a) that

$$ad\frac{1}{2}\eta \frac{\mathrm{\Delta }m_{}^2}{\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}}\mathrm{sin}2\theta .$$
(15)
From Eqs. (9a) and (9b) and the convention $`m_1<m_2`$, we can express the masses $`m_1`$ and $`m_2`$ as

$`m_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(|ad|+{\displaystyle \frac{\mathrm{\Delta }m_{}^2}{|ad|}}\right),`$ (16a)
$`m_2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(|ad|+{\displaystyle \frac{\mathrm{\Delta }m_{}^2}{|ad|}}\right).`$ (16b)
Inspection of Eqs. (7a) and (7b) reveals that our conventions fix the signs of $`\mu _{1,2}`$: $`\text{sign}\mu _1=\text{sign}\mu _2=\eta `$. Then, with Eq. (15), an estimate of the neutrino masses which neglects the solar mass-squared difference is given by

$$m_1m_2\frac{\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}}{\mathrm{sin}2\theta }\text{and}m_3\frac{\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}}{|\mathrm{tan}2\theta |}.$$
(17)
This equation tells us that $`m_3<m_{1,2}`$, at least in the regime of large solar mixing.

## 4 Neutrino masses in the general Zee model

In the previous section we have discussed the mass matrix determined by Eqs. (5) and (8) without reference to any specific model of neutrino masses. Now we introduce the Zee model  and discuss the neutrino mass matrix in the case that both scalar doublets of the Zee model couple in the lepton section. The Yukawa Lagrangian is given by

$$_Y=\underset{a=1}{\overset{2}{}}\overline{L}\mathrm{\Gamma }_a\varphi _a\mathrm{}_R+L^TC^1i\tau _2fLh^++\mathrm{h}.\mathrm{c}.,$$
(18)
where $`f`$ is an antisymmetric $`3\times 3`$ matrix . The mass matrix of the charged leptons arises at tree level through

$$\varphi _a_0=\frac{v_a}{\sqrt{2}}\left(\begin{array}{c}0\\ 1\end{array}\right)\text{and}M_{\mathrm{}}=\frac{1}{\sqrt{2}}\underset{a=1}{\overset{2}{}}v_a\mathrm{\Gamma }_a$$
(19)
with $`v\sqrt{|v_1|^2+|v_2|^2}246`$ GeV. The physical charged scalar fields $`H_1^+`$, $`H_2^+`$ with masses $`M_1`$, $`M_2`$, respectively, and the would-be Goldstone boson $`\varphi _w^+`$ are obtained by the unitary transformation

$$\left(\begin{array}{c}\varphi _1^+\\ \varphi _2^+\\ h^+\end{array}\right)=\left(\begin{array}{ccc}\frac{v_1}{v}& \frac{v_2^{}}{v}W_{11}& \frac{v_2^{}}{v}W_{12}\\ \frac{v_2}{v}& \frac{v_1^{}}{v}W_{11}& \frac{v_1^{}}{v}W_{12}\\ 0& W_{21}& W_{22}\end{array}\right)\left(\begin{array}{c}\varphi _w^+\\ H_1^+\\ H_2^+\end{array}\right).$$
(20)
As anticipated in the introduction, we assume to be in a basis where the charged lepton mass matrix is diagonal, i.e., $`M_{\mathrm{}}=\widehat{M}_{\mathrm{}}`$, where the hat symbolizes that this mass matrix is diagonal. In this basis we have

$$\mathrm{\Gamma }_1=\frac{1}{v_1}(\sqrt{2}\widehat{M}_{\mathrm{}}v_2\mathrm{\Gamma }_2).$$
(21)
In our parameterization the flavour-changing Higgs couplings are given by the off-diagonal elements of $`\mathrm{\Gamma }_2`$. Furthermore, off-diagonal elements of the 2$`\times `$2 unitary matrix $`W`$ are present because of the vacuum expectation value of the lepton number-violating term

$$\lambda h^+\varphi _1^{}\stackrel{~}{\varphi }_2+\text{h.c.}$$
(22)
in the Higgs potential.

The physical charged Higgses couple to the leptons in the following way:

$`_Y(H^+)=`$ (23)
$`\overline{\nu }_L\left(\sqrt{2}{\displaystyle \frac{v_2^{}}{v_1v}}\widehat{M}_{\mathrm{}}+{\displaystyle \frac{v}{v_1}}\mathrm{\Gamma }_2\right)\mathrm{}_R(W_{11}H_1^++W_{12}H_2^+)`$
$`+\overline{(\nu _L)^c}(2f)\mathrm{}_L(W_{21}H_1^++W_{22}H_2^+)+\text{h.c.}`$
With the formulas in the appendix we obtain

$$_\nu =\underset{j=1}{\overset{2}{}}A_L^j\widehat{M}_{\mathrm{}}I(M_j^2,\widehat{M}_{\mathrm{}}^2)A_{R}^{j}{}_{}{}^{}+\text{transp.}$$
(24)
with

$$A_R^j=\left(\sqrt{2}\frac{v_2^{}}{v_1v}\widehat{M}_{\mathrm{}}+\frac{v}{v_1}\mathrm{\Gamma }_2\right)W_{1j}\text{and}A_L^j=2fW_{2j}.$$
(25)
The infinity in Eq. (24) cancels  because $`_{j=1}^2W_{2j}W_{1j}^{}=0`$. Defining

$$J(M_1^2,M_2^2,m^2)=\frac{M_1^2}{M_1^2m^2}\mathrm{ln}\frac{M_1^2}{m^2}\frac{M_2^2}{M_2^2m^2}\mathrm{ln}\frac{M_2^2}{m^2},$$
(26)
from Eq. (24) we obtain the final result

$`_\nu =2W_{21}W_{11}^{}\times {\displaystyle \frac{1}{(4\pi )^2}}`$ (27)
$`\times \left\{\sqrt{2}{\displaystyle \frac{v_2}{v_1^{}v}}\left(\widehat{M}_{\mathrm{}}^2\widehat{J}ff\widehat{J}\widehat{M}_{\mathrm{}}^2\right){\displaystyle \frac{v}{v_1^{}}}\left(\gamma \widehat{M}_{\mathrm{}}\widehat{J}ff\widehat{J}\widehat{M}_{\mathrm{}}\gamma ^T\right)\right\}`$
with $`\widehat{J}\text{diag}(J(M_1^2,M_2^2,m_\alpha ^2))`$ ($`\alpha =e,\mu ,\tau `$) and $`\gamma \mathrm{\Gamma }_2^{}`$. Note that for $`m^2M_{1,2}^2`$ the function $`J`$ simplifies to

$$J2\mathrm{ln}(M_1/M_2).$$
(28)
For the product $`W_{21}W_{11}^{}`$ of elements of the charged-scalar mixing matrix $`W`$, we obtain the relation

$$W_{21}W_{11}^{}=\frac{\lambda ^{}v}{\sqrt{2}(M_1^2M_2^2)}.$$
(29)
It shows explicitly that the Majorana neutrino masses are proportional to the coupling $`\lambda `$ in the Higgs potential (22).

In the GZM considered here, there are family lepton number-violating processes induced by the charged and the neutral scalar interactions. Experimental bounds constrain the coupling matrices $`f`$ and $`\mathrm{\Gamma }_2`$.

## 5 The simplified mass matrix within the Zee model

In this section, our aim is to reproduce the neutrino mass matrix defined by Eqs. (5) and (8) within the GZM. In order to save the amount of writing we introduce the notation

$$_\nu =A\left(\widehat{r}^2ff\widehat{r}^2\right)B\left(\gamma \widehat{r}ff\widehat{r}\gamma ^T\right)\text{with}\widehat{r}=\text{diag}(m_e,m_\mu ,m_\tau )/v.$$
(30)
Both constants $`A`$ and $`B`$ are of the order of 1 GeV, resulting from dividing the electroweak scale by $`16\pi ^2`$:

$$B=2W_{21}W_{11}^{}\times \frac{1}{(4\pi )^2}\mathrm{ln}\frac{M_1^2}{M_2^2}\times \frac{v^2}{v_1^{}}\text{and}A=\sqrt{2}\frac{v_2}{v}B.$$
(31)
Since the off-diagonal elements of $`\gamma `$ introduce flavour-changing neutral interactions, we adopt the philosophy to set to zero as many of them as possible. As we will see, it turns out that

$$\gamma _{e\tau }0\text{and}\gamma _{\tau \tau }0,$$
(32)
and all other elements of $`\gamma `$ being equal to zero, is sufficiently general to avoid bimaximal mixing, which necessarily happens for $`\gamma =0`$ . It can easily be checked that with this assumption Condition 2 (6b) is fulfilled by having $`_{\nu \mu \mu }^{}=_{\nu \tau \tau }^{}=0`$. Furthermore, we assume that all quantities we deal with are real: $`A`$, $`B`$, the elements of the matrix $`f`$, and $`\gamma _{e\tau }`$ and $`\gamma _{\tau \tau }`$. Thus we identify $`_\nu `$ with $`_\nu ^{}`$.

We have to reproduce with the mass matrix (30) the parameters $`a`$ (11a), $`b`$ (11b), $`d`$ (9d) and Condition 1 (6a). Thus we have five coupling constants, three in $`f`$ and $`\gamma _{e\tau }`$, $`\gamma _{\tau \tau }`$ and four relations. It is convenient to express the other four coupling constants as a function of $`\gamma _{\tau \tau }`$:

$`f_{e\mu }`$ $`=`$ $`{\displaystyle \frac{b}{A(r_\mu ^2r_e^2)}}\left\{1+{\displaystyle \frac{ad}{2b^2}}{\displaystyle \frac{A(r_\tau ^2r_e^2)Br_\tau \gamma _{\tau \tau }}{A(r_\tau ^2r_\mu ^2)Br_\tau \gamma _{\tau \tau }}}\right\},`$ (33a)
$`f_{e\tau }`$ $`=`$ $`{\displaystyle \frac{b}{A(r_\tau ^2r_e^2)Br_\tau \gamma _{\tau \tau }}},`$ (33b)
$`f_{\mu \tau }`$ $`=`$ $`{\displaystyle \frac{d}{A(r_\tau ^2r_\mu ^2)Br_\tau \gamma _{\tau \tau }}},`$ (33c)
$`\gamma _{e\tau }`$ $`=`$ $`{\displaystyle \frac{a}{2Br_\tau b}}\left\{A(r_\tau ^2r_e^2)Br_\tau \gamma _{\tau \tau }\right\}.`$ (33d)
## 6 Numerical estimates

Let us now estimate the values of the coupling constants. For definiteness we take $`\mathrm{\Delta }m_{\mathrm{atm}}^2=3\times 10^3`$ eV2 and $`A=2`$ GeV. Furthermore, we need the values of $`r_\tau ^25.22\times 10^5`$ and $`m_\mu ^2/m_\tau ^23.535\times 10^3`$. Defining $`x=B\gamma _{\tau \tau }/Ar_\tau `$ and assuming that666This avoids some finetuning for $`\gamma _{\tau \tau }`$. Note, however, that $`x=0`$ is also possible. $`x1`$ and $`(1m_e^2/m_\tau ^2x)/(1m_\mu ^2/m_\tau ^2x)1`$, we obtain

$`f_{e\mu }`$ $``$ $`1.0\times 10^4\times {\displaystyle \frac{\text{sign}b}{\mathrm{sin}^22\theta }},`$ (34a)
$`f_{e\tau }`$ $``$ $`3.7\times 10^7\times {\displaystyle \frac{\text{sign}b}{1x}},`$ (34b)
$`f_{\mu \tau }`$ $``$ $`5.2\times 10^7\times {\displaystyle \frac{1}{|\mathrm{tan}2\theta |(1x)}},`$ (34c)
$`\gamma _{e\tau }`$ $``$ $`7.2\times 10^3\times {\displaystyle \frac{v_2}{v}}{\displaystyle \frac{1x}{\mathrm{tan}2\theta }}.`$ (34d)
These equations serve to see the orders of magnitude and any effects of $`\mathrm{\Delta }m_{}^2`$ are neglected. As can be seen from Eq. (15), a considerable amount of finetuning is involved in order to reproduce the solar mass-squared difference.

Now we concentrate on the LMA MSW solution of the solar neutrino problem, where $`\theta `$ is in the first octant. With Eq. (9a) it follows that $`b>0`$. In this case a representative value of the mixing angle $`\theta `$ is given by the best fit value $`\mathrm{sin}^22\theta 0.75`$ ($`\mathrm{tan}2\theta 1.71`$) of Ref. , with the corresponding mass-squared difference $`\mathrm{\Delta }m_{}^23.2\times 10^5`$ eV2. We note that in this case we have $`|f_{e\tau }/f_{e\mu }|m_\mu ^2/m_\tau ^2`$, which is similar to the case of the Zee mass matrix . On the other hand, in our scenario we have $`f_{e\tau }f_{\mu \tau }`$, whereas for the Zee mass matrix the relation $`|f_{e\tau }||f_{\mu \tau }|`$ holds . As far as $`\gamma `$ is concerned, with $`x1`$ we can have $`\gamma _{e\tau }\gamma _{\tau \tau }10^3÷10^2`$.

Due to our assumption (32), flavour-changing neutral scalar interactions at the tree level are very constrained. Among the charged lepton decays we only have $`\tau ^{}e^{}\mathrm{}^+\mathrm{}^{}`$ with $`\mathrm{}=e,\mu `$. A generous estimate of the branching ratio of this decay for $`\mathrm{}=\mu `$ is obtained by

$$\text{Br}(\tau ^{}e^{}\mu ^+\mu ^{})\frac{1}{G_F^2M_0^4}\gamma _{e\tau }^2\left(\frac{m_\mu }{v}\right)^2<10^8,$$
(35)
where we have taken $`G_F^2M_0^410^2`$ ($`M_0100`$ GeV) and $`M_0`$ is a generic neutral Higgs mass. At the 1-loop level, neutral Higgs exchange also induces the decay $`\tau ^{}e^{}\gamma `$. Making again an estimate, we obtain

$$\text{Br}(\tau ^{}e^{}\gamma )\frac{\alpha }{48\pi }\frac{1}{G_F^2M_0^4}(\gamma _{e\tau }\gamma _{\tau \tau })^2<10^{10}.$$
(36)
Both estimates are well compatible with the experimental upper bounds on these branching ratios of the order of $`10^6`$ . More detailed discussions of these decays are found in Refs. . We have used $`|\gamma _{e\tau }|\stackrel{<}{}\mathrm{\hspace{0.25em}10}^2`$ and $`|\gamma _{\tau \tau }|\stackrel{<}{}\mathrm{\hspace{0.25em}10}^2`$ in Eqs. (35) and (36).

Also the charged scalars participate in various charged lepton decays as intermediate particles. Numerous decays of the type $`\mathrm{}_1^{}\mathrm{}_2^{}+2\text{neutrinos}`$ proceed at tree level via the Lagrangian (23). According to the couplings in this Lagrangian we can distinguish between $`f`$, $`\gamma `$ and $`\widehat{r}`$ vertices, and we can have decay amplitudes with all possible combinations of these vertices, except $`\gamma `$-$`\gamma `$, which is forbidden due to the restricted form of $`\gamma `$, Eq. (32). E.g., to the Standard Model amplitude of ordinary muon decay there is an amplitude with couplings $`r_\mu r_e`$ and another one with $`f_{e\tau }f_{\mu \tau }`$; another example is $`\tau ^{}\overline{\nu }_e\overline{\nu }_\mu \mu ^{}`$ with couplings $`f_{e\tau }r_\mu `$, which is lepton number violating. The branching ratios of all these decays are negligible because of the smallness of the coupling constants $`f_{\alpha \beta }`$ and $`\gamma _{e\tau }`$ and the ratios $`r_\alpha `$. Also negligible are radiative decays $`\mathrm{}_1^{}\mathrm{}_2^{}+\gamma `$ induced by charged Higgs loops . In this case one always has two $`f`$-vertices in the loop graph, except in the case of the amplitude for $`\tau ^{}e^{}+\gamma `$, where there are two contributions, proportional to $`f_{e\mu }f_{\mu \tau }`$ and $`\gamma _{e\tau }r_e`$. Recent reviews of the restrictions on the coupling constants $`f_{\alpha \beta }`$ are found in Refs. .

Scalar contributions to the anomalous magnetic moments of the electron and muon involving $`f`$  and $`\widehat{r}`$-couplings are totally negligible because these constants are too small. A contribution from the $`\gamma `$-couplings to the electron magnetic moment coming from $`\tau `$ exchange is proportional to $`\gamma _{e\tau }\gamma _{\tau e}`$  and is thus zero in view of $`\gamma _{\tau e}=0`$.

In our scenario we have $`_{\nu ee}=a0`$. The matrix element $`_{\nu ee}`$ is identical with the effective neutrino mass $`m_\nu `$ probed in neutrinoless double-beta decay. Therefore, this decay is allowed and the effective neutrino mass is given by

$$|m_\nu |m_3\frac{0.05\text{eV}}{\mathrm{tan}2\theta }.$$
(37)
This represents an order of magnitude which is accessible in future experiments (for recent reviews see, e.g., Ref. ).

Having seen that the LMA MSW solution of the solar neutrino problem can be accommodated in our scenario, we now proceed to the small mixing angle (SMA) MSW solution. In this case we take as illustration the best fit value of Ref. , $`\mathrm{sin}^22\theta 2.3\times 10^3`$ ($`\mathrm{tan}2\theta 4.8\times 10^2`$), which has a corresponding $`\mathrm{\Delta }m_{}^20.5\times 10^5`$ eV2. From Eq. (34a) we see that now $`f_{e\mu }0.05`$ becomes relatively large and barely compatible with the requirement that the Zee boson does not have an effect on the muon decay rate so that it does not destroy the agreement in electroweak precision tests . Thus in our simple scenario we cannot incorporate safely the SMA solution.

## 7 The limit $`\mathrm{\Gamma }_20`$

To make contact with Refs. , we explore the effect of $`\mathrm{\Gamma }_2=0`$. This means that all diagonal elements of the neutrino mass matrix are zero  and, in our notation, we have $`a=0`$, i.e., $`m_3^2=\mathrm{\Delta }m_{}^2|\mathrm{cos}2\theta |`$ and $`\eta =1`$ (see (Eq. (11a)). These relations and inspection of the consistency condition (10) leads to

$$m_3\frac{1}{2}\frac{\mathrm{\Delta }m_{}^2}{\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}},m_1m_2\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}\text{and}|\mathrm{cos}2\theta |\frac{1}{4}\frac{\mathrm{\Delta }m_{}^2}{\mathrm{\Delta }m_{\mathrm{atm}}^2}.$$
(38)
From the last relation we read off that $`\theta `$ is 45 for all practical purposes. The consequences for the coupling matrix $`f`$ are

$$\left|\frac{f_{e\tau }}{f_{e\mu }}\right|\frac{m_\mu ^2}{m_\tau ^2}\text{and}\left|\frac{f_{e\tau }}{f_{\mu \tau }}\right|\sqrt{2}\frac{\mathrm{\Delta }m_{\mathrm{atm}}^2}{\mathrm{\Delta }m_{}^2}.$$
(39)
The first of these two relations is obtained by taking the ratio of Eq. (33b) and Eq. (33a), where $`a=0`$ is taken into account. In the second relation, Eq. (11b) has been used. The results (38) and (39) agree with those of Refs. .

## 8 Summary

In this paper we have discussed neutrino masses and mixing in the general Zee model, where both Higgs doublets couple in the lepton sector. In this endeavour, we were motivated by the result that for $`\mathrm{\Gamma }_2=0`$, where $`\mathrm{\Gamma }_2`$ denotes the Yukawa coupling matrix of the second Higgs doublet, the Zee mass matrix leads to bimaximal mixing. It is true that a mixing angle of $`45^{}`$ is perfect for the description of the atmospheric neutrino data, but it does not represent a very good fit for the solar neutrino data. The general Zee model is a rather rich and intricate model. Therefore, in order to simplify the analysis, we have assumed that all quantities appearing in the neutrino mass matrix are real and that in the second Yukawa coupling matrix $`\mathrm{\Gamma }_2`$, which is set to zero usually, only two elements, the $`e\tau `$ and $`\tau \tau `$ elements, are non-zero. We have shown that this is sufficient to accommodate solutions of the solar neutrino problem with a large mixing angle instead of maximal mixing; on the other hand, in the atmospheric sector maximal mixing remains. At the same time, all dangerous processes induced by the new couplings are sufficiently suppressed. Actually, we could have even set the $`\tau \tau `$ element of $`\mathrm{\Gamma }_2`$ equal to zero and used a single non-zero element in this coupling matrix. However, the $`\tau \tau `$ element represents a possibility to achieve a further suppression of the potentially more dangerous $`e\tau `$ element. As in the case $`\mathrm{\Gamma }_2=0`$, the neutrino mass $`m_3`$ is the smallest neutrino mass, and the resulting mass spectrum is rather of the type which is called “inverted hierarchy”. However, whereas for $`\mathrm{\Gamma }_2=0`$ one has $`m_3m_1m_2`$, now $`m_3\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}/\mathrm{tan}2\theta `$ is of the same order of magnitude as the other two masses.

Numerous finetunings are involved in our scenario. Implementing the condition (6a) and setting all but two elements of $`\mathrm{\Gamma }_2`$ equal to zero represent rather severe finetunings. These procedures were useful for exactly having $`U_{e3}=0`$ and an atmospheric mixing angle of $`45^{}`$, and for avoiding a dangerous class of lepton family number-violating interactions. In order to reproduce the neutrino masses within our scenario, we need $`|f_{e\tau }||f_{\mu \tau }||f_{e\mu }|(m_\mu /m_\tau )^2`$ and a quite drastic finetuning between $`a`$ and $`d`$ in the mass matrix (5), in order to accommodate the solar mass-squared difference. Of course, the mass-squared difference $`\mathrm{\Delta }m_{}^210^8`$ eV2 of the quasi-vacuum oscillation solution  would require a stronger finetuning than $`\mathrm{\Delta }m_{}^210^5`$ eV2 of the large mixing angle solution. Examination of the Zee model with a general $`\mathrm{\Gamma }_2`$ and with small deviations of the neutrino mass matrix from the form (5) should reveal how much our finetunings could be relaxed. In the present work we have confined ourselves to prove that it is possible to reproduce the LMA MSW solution within the Zee model.

In conclusion, even in the case that bimaximal mixing is ruled out, the Zee model will remain a viable and interesting scenario in order to accommodate the neutrino oscillation solutions of the solar and atmospheric neutrino problems.

## Acknowledgements

This work was supported by Spanish DGICYT under grant PB98-0693, by the European Commission TMR network HPRN-CT-2000-00148 and the ESF network 86. T.S. is supported by the Marie Curie Training Site Program for Particle Physics Beyond the Standard Model, contract number HPMT-2000-00124.

## Appendix A Charged scalar exchange and the 1-loop neutrino mass

We proceed from the general Lagrangian

$$_Y(\nu _L,\mathrm{},\phi ^+)=(\overline{\nu }_LA_R\mathrm{}_R+\overline{(\nu _L)^c}A_L\mathrm{}_L)\phi ^++\text{h.c.},$$
(A1)
where the $`A_{L,R}`$ are 3$`\times `$3 coupling matrices in the case of 3 families and $`\mathrm{}_{L,R}`$ denote chiral charged fermion fields. The neutrino Majorana mass Lagrangian is defined by

$$_{\nu \mathrm{mass}}=\frac{1}{2}\nu _L^TC^1_\nu \nu _L+\text{h.c.}$$
(A2)
With the scalar interactions (A1) we obtain

$$_\nu ^{1\mathrm{loop}}=A_L\widehat{M}_fI(M_\phi ^2,\widehat{M}_f^2)A_R^{}+\text{transp.},$$
(A3)
where we have defined

$`I(M^2,m^2)={\displaystyle \frac{1}{(4\pi )^2}}`$
$`\times \left\{{\displaystyle \frac{1}{ϵ}}+\mathrm{ln}(4\pi )+\gamma _E1+{\displaystyle \frac{1}{M^2m^2}}\left(M^2\mathrm{ln}{\displaystyle \frac{M^2}{\mu ^2}}m^2\mathrm{ln}{\displaystyle \frac{m^2}{\mu ^2}}\right)\right\}.`$
We have used dimensional regularization; thus, $`ϵ=(4n)/2`$ with $`n`$ being the number of space-time dimensions, $`\gamma _E`$ is Euler’s constant and $`\mu `$ is an arbitrary mass scale.


text:
# Vertical pairing of identical particles suspended in the plasma sheath

## Abstract

It is shown experimentally that vertical pairing of two identical microspheres suspended in the sheath of a radio-frequency (rf) discharge at low gas pressures (a few Pa), appears at a well defined instability threshold of the rf power. The transition is reversible, but with significant hysteresis on the second stage. A simple model, which uses measured microsphere resonance frequencies and takes into account besides Coulomb interaction between negatively charged microspheres also their interaction with positive ion wake charges, seems to explain the instability threshold quite well.

Recent interest in the properties of complex plasmas – charged monodisperse microparticles suspended in an electron-ion environment, is partially due to the possibility to model condensed matter phenomena on an “atomic” level . The particles are charged negatively in a radio-frequency (rf) discharge (up to $`10^310^4`$ elementary charges) and levitate in the lower plasma sheath, where gravity can be balanced by an inhomogeneous vertical electric field. The particles repel each other via a screened Coulomb interaction. Inside a confining horizontal potential they arrange themselves either into an ordered structure – the “plasma crystal”, or form disordered, liquid- or gas-like, states .

At low gas pressures one has typically only a very few layers (as few as one) for these structures. Ion streaming motion in the plasma sheath produces a nonequilibrium environment, which is reflected in the properties of the complex plasma structures and makes them rather different from the known equilibrium ones . One of the striking observations at sufficiently low gas densities (such that the ion-neutral mean free path is not small compared with a typical distance between microspheres) is an unusual “stacking” of the particles such that adjacent horizontal layers are located on top of each other . This vertical “polarization” of the dusty plasma was ascribed to the attractive forces which result from the focusing of the ion flux under each particle (ion wake effect, see, e.g., Ref. ). Thus, the attractive interaction between the horizontal layers in a plasma crystal is asymmetric, such that the attractive force is communicated only downstream in the direction of the ion flow . The same arguments were applied to predict a “binding” force for a possible “dust molecule” formation . This theoretical prediction was verified recently by experimental studies of the competition between repulsive and attractive forces acting on two different dust particles levitated on different levels in the plasma sheath.

In this Letter we present experiments on a new instability which is observed for two identical microspheres suspended initially on the same level in the plasma sheath. In this work we will report on the quantitative investigations involving (for simplicity) only two microspheres, which were undertaken to clarify the basic behavior for many identical particles. The instability appears first as a continuous (forward) bifurcation to a state of the vertically separated microspheres, and then to a discontinuous vertical pairing. Depending on gas pressure, the last stage – the pairing – can be either strongly hysteretic (for lower pressure) or weakly hysteretic (for higher pressure).

The experiments were performed in a standard Gaseous Electronics Conference (GEC) rf Reference Cell  with the lower electrode powered at 13.56 MHz and the upper electrode grounded (see Fig. 1). Argon gas at various pressures between 1 and 7 Pa was used, and the discharge power (or a rf peak-to-peak voltage, $`U_{\mathrm{pp}}`$) was the control parameter of the instability described. The electron temperature and density were measured at the center of the discharge with a rf-compensated passive Langmuir probe. At these conditions the electron temperature $`T_e`$ was found to be within 2-5 eV, while the electron density $`n_e`$ ranged from $`10^7`$ to $`8\times 10^8`$ cm-3, so that the electron Debye length varied from $`\lambda _{\mathrm{D}e}0.5`$ mm (high $`U_{\mathrm{pp}}`$) to $`\lambda _{\mathrm{D}e}5`$ mm (low $`U_{\mathrm{pp}}`$). These values of the Debye length are comparable to, or larger than the separation distances measured, and thus Debye shielding will play only a minor role in this rather low-density plasma. The microspheres, suspended in the plasma, were illuminated by a laser sheet of about $`100\mu `$m thickness, and their imaging was performed by external CCD cameras from the top and from the side (in both cases via 300 mm objectives). We used polystyrene (density $`1.05`$ g/cm3) microspheres with diameter $`7.6\pm 0.1\mu `$m and mass $`M2.4\times 10^{10}`$ g. In order to confine the microspheres horizontally we placed a quartz glass cylinder of 20 mm height and 50 mm diameter on the lower electrode (see Fig. 1). In contrast to a copper ring normally used to confine particles, the glass cylinder produces extremely flat, almost “square well” confining potential. The latter conclusion is based on our observation of unusually flat shape of 2D plasma crystal, particularly up to its edges, contrary to the trap formed by a standard metal ring. The “square well” confining potential provides us the opportunity to study two identical microspheres suspended in the plasma sheath on the same initial height.

In order to quantitatively investigate the pairing instability and to verify that it is indeed symmetry breaking, we conducted experiments with two identical microspheres at different pressures. Figure 2 shows, how, at a pressure of 2 Pa, the relative position of one microsphere, which eventually moves below and then under the other (i.e., downstream) with respect to the eventual upstream microsphere, varies as the power (or $`U_{\mathrm{pp}}`$) decreases. As long as $`U_{\mathrm{pp}}`$ exceeds the threshold value $`U_{\mathrm{pp}}^{\mathrm{th}}55`$ V, both particles are located on the same level (initial position 1). When $`U_{\mathrm{pp}}`$ decreases below the threshold, vertical separation starts (position 2) and grows continuously up to position 3. A further small $`U_{\mathrm{pp}}`$ decrease causes the lower particle to “jump” to a position directly beneath the upper one and to create a vertical pair (final position 4). By reversing the process, i.e., increasing $`U_{\mathrm{pp}}`$, the particles revert back to their initial horizontal configuration. The transition is strongly hysteretic with respect to $`U_{\mathrm{pp}}`$. Further measurements showed that the threshold, $`U_{\mathrm{pp}}^{\mathrm{th}}`$, rapidly increases with pressure, while the relative width of the hysteresis decreases with pressure.

Figure 3 represents the forward (positions 1-5) and the reverse (positions 6,7) transitions at $`p=7`$ Pa. The instability (vertical separation) starts at $`U_{\mathrm{pp}}^{\mathrm{th}}200`$ V and develops continuously until position 4. Then a small $`U_{\mathrm{pp}}`$ decrease is accompanied by the discontinuous vertical pairing (position 5). When the voltage is increased back, the particles remain vertically paired until position 6 is reached, and then the lower particle “jumps” to position 7. In Fig. 4 the vertical separation distance is plotted against the control parameter, $`U_{\mathrm{pp}}^{\mathrm{th}}U_{\mathrm{pp}}`$, at $`p=7`$ Pa. As in the low pressure case, the transition is divided into two stages: (i) continuous transition from the horizontal configuration, and (ii) discontinuous, hysteretic transition to the final vertical pairing. It is worth noting that at the onset of the continuous transition strong symmetrical fluctuations of the vertical particle separation are observed in the experiments. This means that the pairing instability may be initiated with equal probability by either microsphere in the pair. Therefore, this transition is related to the symmetry breaking.

In the experiments both $`U_{\mathrm{pp}}`$ (at the fixed $`p`$) and $`p`$ (at the fixed $`U_{\mathrm{pp}}`$) were used as the control parameters of the instability. But in a subsequent analysis, presented below, it is much more convenient to use the resonance frequencies of vertical, $`\omega _z`$, and horizontal, $`\omega _r`$, oscillations of a single particle in the sheath as the control parameters. The resonance frequencies are certain functions of $`U_{\mathrm{pp}}`$ and $`p`$, and can be determined experimentally.

We first examine the simplest model without wake effects. The model uses the measured values of $`\omega _z`$ and $`\omega _r`$ as the control parameters and provides an order of magnitude estimate for the instability threshold. We consider a pair of two identical particles of mass $`M`$ and (negative) charge $`Q`$, electrostatically confined in the harmonic potential well. If the particles are separated horizontally by $`R`$ and vertically by $`\delta `$ (i.e., each particle is displaced at $`R/2`$ and $`\delta /2`$ from the center), then their energy in the confinement is $`\frac{1}{4}M(\omega _z^2\delta ^2+\omega _r^2R^2)`$. Since the Debye length is relatively large, the bare coupling energy of two unshielded particles, $`Q^2(R^2+\delta ^2)^{1/2}`$, should be added to the confinement energy. This system has two stable configurations – horizontally (when $`\omega _z>\omega _r`$) or vertically (when $`\omega _z<\omega _r`$) aligned at distances of $`R=(2Q^2/M\omega _r^2)^{1/3}`$ or $`\delta =(2Q^2/M\omega _z^2)^{1/3}`$, respectively.

Our model is characterized in terms of the measured resonance frequencies. Thus, it is necessary to trace these frequencies as the control parameter, $`U_{\mathrm{pp}}^{\mathrm{th}}U_{\mathrm{pp}}`$, is varied. The typical dependence of $`\omega _r`$ and $`\omega _z`$ on the value of $`U_{\mathrm{pp}}`$ at $`p=`$ 1 Pa is shown in Fig. 5. One can see that the frequencies converge rapidly as $`U_{\mathrm{pp}}`$ decreases. Comparison of Figs. 4 and 5 shows (albeit for different pressures) that the transition starts continuously when $`\omega _z`$ still considerably exceeds $`\omega _r`$. This implies that the harmonic approximation with the isotropic particle field is too simple. Thus, we have to include the wake effect in its simplest form. For consistency in the nonlinearity analysis, we must also take into consideration the anharmonicity of the potential well in the vertical direction . The corresponding expression for the total potential energy of a particle pair in a 2D confinement is

$$𝒲_{\mathrm{pair}}\frac{M}{4}\left(\omega _r^2R^2+\omega _z^2\delta ^2+\frac{\beta ^{}}{8}\delta ^4\right)+\frac{Q^2}{\sqrt{R^2+\delta ^2}}\frac{Qq}{2\sqrt{R^2+(\mathrm{\Delta }+\delta )^2}}\frac{Qq}{2\sqrt{R^2+(\mathrm{\Delta }\delta )^2}}.$$
(1)
In writing Eq. (1) the simplest two-microsphere model of the wake is used (see Fig. 6). We treat the excessive positive charge of the wake, $`q`$, as point-like, located under the particle at distance $`\mathrm{\Delta }`$. We assume also that since both the particle and the wake charges vary only slightly with variation in $`R`$ and $`\delta `$, so that $`Q`$ and $`q`$ are approximately constant. The last two terms in Eq. (1) represent the particle-wake interaction, and cross-terms of the particle-wake and wake-wake ($`qq`$) interactions are neglected . Since the coefficients multiplying the terms $`\delta ^3`$ cancel out for a pair of identical particles, the lowest order vertical anharmonicity in Eq. (1) is the fourth order term with coefficient $`\beta ^{}>0`$ but whose value is here unknown, since it was not feasible to measure it by technique used elsewhere .

The dependence $`R(\delta )`$ is determined from the equilibrium condition in the radial direction, $`𝒲_{\mathrm{pair}}/R=0`$:

$$\frac{2Q^2}{(R^2+\delta ^2)^{3/2}}Qq\left(\frac{1}{[R^2+(\mathrm{\Delta }+\delta )^2]^{3/2}}+\frac{1}{[R^2+(\mathrm{\Delta }\delta )^2]^{3/2}}\right)=M\omega _r^2.$$
(2)
The equilibrium value of $`\delta `$ is given by the condition $`𝒲_{\mathrm{pair}}/\delta =0`$. A linear combination of this condition and Eq. (2) gives

$$M(\omega _{z}^{}{}_{}{}^{2}\omega _{r}^{}{}_{}{}^{2})\delta +\frac{M\beta ^{}}{4}\delta ^3+Qq\mathrm{\Delta }\left(\frac{1}{[R^2+(\mathrm{\Delta }+\delta )^2]^{3/2}}\frac{1}{[R^2+(\mathrm{\Delta }\delta )^2]^{3/2}}\right)=0.$$
(3)
Expanding Eq.(3) to third power in $`\delta `$ and substituting $`R(\delta )`$ from Eq. (2) provide a typical stationary equation for the vertical displacement as the order parameter:

$$\left(\frac{\omega _z^2}{\omega _r^2}\mathrm{\Omega }^2\right)\stackrel{~}{\delta }+\frac{\beta ^{}R_0^2}{4\omega _r^2}\stackrel{~}{\delta }^3\frac{35(\mathrm{\Omega }^21)\stackrel{~}{\mathrm{\Delta }}^2}{6(1+\stackrel{~}{\mathrm{\Delta }}^2)^2}\left(\frac{1+\frac{8}{7}\frac{q/Q}{(1+\stackrel{~}{\mathrm{\Delta }}^2)^{5/2}}}{1\frac{q/Q}{(1+\stackrel{~}{\mathrm{\Delta }}^2)^{5/2}}}\right)\stackrel{~}{\delta }^3=0.$$
(4)
Here we have introduced the dimensionless parameters, $`\stackrel{~}{\delta }=\delta /R_0`$ and $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }/R_0`$ \[normalized by $`R_0=R|_{\delta =0}`$ from Eq. (2)\], and the critical frequency ratio,

$$\mathrm{\Omega }\sqrt{1+\frac{6Qq\stackrel{~}{\mathrm{\Delta }}^2}{M\omega _r^2R_0^3(1+\stackrel{~}{\mathrm{\Delta }}^2)^{5/2}}}.$$
(5)
The instability criterion $`(\omega _z/\omega _r)_{\mathrm{cr}}=\mathrm{\Omega }`$ is obtained from Eq. (4) by setting the coefficient of $`\stackrel{~}{\delta }`$ equal to zero. At $`\omega _z/\omega _r>\mathrm{\Omega }`$, the energy $`𝒲_{\mathrm{pair}}(\delta )`$ has a minimum at $`\delta =0`$ and the equilibrium configuration is horizontal. In the opposite case, $`\omega _z/\omega _r<\mathrm{\Omega }`$, the $`\delta =0`$ state becomes unstable, and the particles start separating vertically. From Fig. 5 we note that the separation starts at $`(\omega _z/\omega _r)_{\mathrm{cr}}4`$ (when $`U_{\mathrm{pp}}U_{\mathrm{pp}}^{\mathrm{th}}`$). Assuming $`Qq`$ and substituting Eq. (2) at $`\delta =0`$ into Eq. (5), one gets $`\mathrm{\Omega }`$ as a function of $`\stackrel{~}{\mathrm{\Delta }}`$ only. Then at $`\stackrel{~}{\mathrm{\Delta }}1`$ we obtain $`\mathrm{\Omega }3`$. Thus, the threshold of the pairing instability is adequately described by the model. A saturation of the order parameter, $`\stackrel{~}{\delta }`$, above the instability threshold in the first continuous stage of the transition (see Fig. 4) could be reached if the sum of both coefficients of $`\stackrel{~}{\delta }^3`$ in Eq. (4) is positive. By using the expression $`\mathrm{\Omega }(\stackrel{~}{\mathrm{\Delta }})`$ and the value of the ratio $`Q/q1`$ we achieve this condition at $`\beta ^{}R_0^2/\omega _r^27`$, or $`\beta ^{}/(2\pi )^220`$ Hz2/mm2. Our recent experiments  indicate that at pressures below $`10`$ Pa the anharmonic coefficients become relatively large: For example, $`\beta /(2\pi )^220`$ Hz2/mm2 at $`p1`$ Pa and $`U_{\mathrm{pp}}70`$ V. At smaller $`U_{\mathrm{pp}}`$, the whole vertical structure of the sheath changes dramatically, because the Debye length increases by the order of magnitude. Hence, the anharmonic coefficients should also increase rapidly as $`U_{\mathrm{pp}}`$ decreases, and it is plausible that $`\beta ^{}`$ will provide the instability saturation. Then, as follows from Eq.(4), $`\stackrel{~}{\delta }\sqrt{\mathrm{\Omega }\omega _z/\omega _r}`$ close to the threshold. However, failing an accurate experimental determination of the shape of the sheath potential for the actual conditions of the experiment, further testing of the model via the experimental value of $`\delta `$ is not feasible here.

The modelling of the second (discontinuous) stage of the instability is rather complicated problem: When the lower particle approaches the wake, its actual “shape” becomes crucial, and the approximation of the point-like wake charge cannot be used. This very nonlinear second stage will almost certainly involve nonlinear wake dynamics, and these can probably only be adequately treated by three-dimensional numerical simulations. We would like to point out that the effect of the wake is stronger at lower $`p`$ and $`U_{\mathrm{pp}}`$, when $`\lambda _{\mathrm{D}e}`$ is larger. At the same time, the ion mean free path should exceed $`\lambda _{\mathrm{D}e}`$ (otherwise, the wake effect weakens due to the ion scattering on neutrals). Both these requirements are well satisfied in our experiments.

V. S. is grateful for the support of the Alexander von Humboldt Foundation, Germany, during his stay at MPE. Many useful critical remarks of an anonymous Referee are gratefully acknowledged.