text:
# Structure Function Scaling in the Taurus and Perseus Molecular Cloud Complexes.

## 1 Introduction

Due to the complexity of the Navier–Stokes equations, mathematical work on turbulence is often inspired by experimental and observational measurements. Since geophysical and laboratory flows are predominantly incompressible, turbulence studies have been limited almost entirely to incompressible flows (or to infinitely compressible ones, described by the Burgers equation). Little attention has been paid to highly–compressible, or super–sonic turbulence.

The cold interstellar medium (ISM) of galaxies such as the Milky Way is both highly turbulent and highly super–sonic, within a range of scales from hundreds of parsecs to approximately a tenth of a parsec. Large observational surveys of the cold ISM of our galaxy and of dust within it have become available in the last few years and new surveys of unprecedented sensitivity and resolution will be obtained in the near future by FIR satellites such as SIRTF and Herschel.

The cold ISM provides a good laboratory for super–sonic turbulence and results from new observational surveys should motivate new mathematical work. Furthermore, the interpretation of the astronomical data requires a knowledge of basic properties of super–sonic turbulence.

Previous works have tried to investigate the properties of ISM turbulence by estimating the second order structure function or the power spectrum of the velocity field in dark clouds (Kleiner & Dickman, 1987; Hobson, 1992; Miesch & Bally, 1994; Brunt & Heyer, 2002a, b), as sampled by molecular emission lines, with a number of different methods. The centroid velocity at each map position is used as an estimate of the local radial velocity. This centroid velocity results from a convolution of density, velocity and excitation temperature along the line of sight and is not easily related to the three dimensional velocity structure. Other works have instead studied the column density distribution of dark clouds, estimated from the integrated intensity of molecular emission lines or from dust thermal emission. Fractal dimensions have been estimated (Beech, 1987; Bazell & Désert, 1988; Scalo, 1990; Dickman et al., 1990; Falgarone et al., 1991; Zimmermann et al., 1992; Henriksen, 1991; Hetem & Lepine, 1993; Vogelaar & Wakker, 1994; Elmegreen & Falgarone, 1996). Multifractal (Chappell & Scalo, 2001) and wavelet (Langer et al., 1993) analysis have also been proposed as a way of characterizing the projected density structure of molecular clouds.

The multifractal analysis applied by Chappell & Scalo (2001) to dust continuum images is related to the scaling of the moments of the projected density. In the present work we study the scaling of integrated intensity differences (structure functions) of the J=1–0 transition of 13CO from the Taurus and the Perseus molecular cloud complexes. If we had computed the moments of intensity, instead of intensity differences, then or analysis would be equivalent to the multifractal analysis by Chappell & Scalo (2001).

We have already shown in previous works that super–sonic turbulence in a roughly isothermal gas, such as the cold ISM, generates a complex density field with density contrasts of several orders of magnitude and with statistical properties consistent with observational data from molecular clouds (Padoan et al., 1998b, 1999; Padoan & Nordlund, 1999). We usually refer to this process as turbulent fragmentation. Turbulent fragmentation of star forming clouds is unavoidable, since super–sonic turbulent motions are ubiquitously observed. The results of this work provide further evidence of the importance of turbulent fragmentation in star forming clouds.

In § 2 we compute the relative scaling of the structure functions and in § 3 we discuss our results. The hierarchical structure model is briefly presented in § 4 and its validity for our data is verified with the so called $`\beta `$ and $`\gamma `$ tests. We draw our conclusions in § 5.

## 2 Structure Functions of Projected Density Field

We use two observational J=1-0 13CO spectral maps of the Perseus MC complex (Billawala et al., 1997) and of the Taurus MC complex (Mizuno et al., 1995). These MC complexes have an extension of approximately 30 pc, and rms radial velocity $`\sigma _v=2.0`$ km/s in Perseus and $`\sigma _v=1.0`$ km/s in Taurus. The rms sonic Mach number is therefore approximately 14 in Perseus and 7 in Taurus, assuming a gas kinetic temperature of approximately 10 K. The angular resolution is approximately 0.09 pc in the Perseus map, assuming a distance of 300 pc and 0.08 pc in the Taurus map, assuming a distance of 140 pc.

In each map we compute an image of the integrated intensity, $`I(𝐫)`$, defined as:

$$I(𝐫)=\underset{v}{}T(v,𝐫)dv,$$
(1)
where $`T(v,𝐫)`$ is the antenna temperature at the velocity channel $`v`$ and map position $`𝐫`$. The integrated intensity is set to zero in each map position with peak antenna temperature smaller than five times the rms noise averaged over the whole map. The integrated intensity of J=1-0 13CO is roughly proportional to the total gas column density since in most lines of sight of the two maps analyzed here this transition is not optically thick. However, the conversion from J=1-0 13CO integrated intensity and gas column density also depends on the abundance ratio between 13CO and H2, and on the precise spatial distribution of the excitation temperature (see Padoan et al. (1998a) for a detailed computation of the effect of the spatial distribution of excitation temperature on the estimation of the J=1-0 13CO column density).

The structure functions of the integrated intensity image $`I(𝐫)`$ are defined as:

$$S_p(r)=|I(𝐫^{})I(𝐫^{}+𝚫𝐫^{})|^p=|\mathrm{\Delta }I|^p$$
(2)
where $`p`$ is the order and the average is extended to all map positions $`𝐫^{}`$ and all position differences $`𝚫𝐫^{}`$ such that $`|𝚫𝐫^{}|=r`$.

It is possible to compute the moments by first deriving the probability density function (PDF) of $`\mathrm{\Delta }I`$, $`P(\mathrm{\Delta }I,r)`$, from the observational map, and then computing:

$$S_p(r)=_{\mathrm{\Delta }I}\mathrm{\Delta }I^pP(\mathrm{\Delta }I,r)d(\mathrm{\Delta }I)$$
(3)
The PDF can be obtained by smoothing appropriately the histogram obtained from the observational data, as suggested by Lévêque & She (1997). Here, however, we prefer to average the moments directly over the whole observational sample, as in (2), without first deriving the PDF.

The structure functions are plotted in Figure 1, up to the order $`p=20`$. They are well approximated by a power law, at least at low orders, between approximately 0.3 and 3 pc, for both Taurus and Perseus. We have performed a least square fit in that range of scales for all orders. The least square fits are shown in Figure 1 as solid lines. We call $`\eta (p)`$ the exponent of the power law fit of the $`p`$ order moment111The Greek letters $`\zeta `$ and $`\tau `$ are usually reserved for the exponents of the moments of velocity differences and energy dissipation respectively.:

$$S_p(r)r^{\eta (p)}$$
(4)
We then compute the scaling exponents relative to the third order, $`\eta (p)/\eta (3)`$, following the idea of extended self–similarity (Benzi et al., 1993; Dubrulle, 1994).

Values of $`\eta (p)/\eta (3)`$ for both Taurus and Perseus are plotted in Figure 2, up to the 20th order, and are compared with those predicted for the scaling of the structure functions of the velocity field ($`\zeta (p)/\zeta (3)`$) by the Kolmogorov model (Kolmogorov, 1941), the She & Lévêque model (formula (7) below) and the model by Boldyrev (2002). The latter is an extension of the She–Lévêque model to super–sonic turbulence, verified with numerical simulations of super–sonic turbulence by (Boldyrev et al., 2002a, b). The scaling exponents of the structure functions of the integrated intensity images of Taurus and Perseus follow very closely the velocity scaling found in Boldyrev (2002), within 5% accuracy up to the 20th order for Taurus and within 10% accuracy for Perseus.

For practical purposes it may be useful to know the intrinsic value of the scaling exponents, not just that relative to the third order one. For Taurus we find $`\eta (3)=1.10`$ and for Perseus $`\eta (3)=1.18`$, both close to the Kolmogorov’s value for the velocity structure function, $`\zeta (3)=1.0`$. The second order exponents are $`\eta (2)=0.77`$ for Taurus and $`\eta (2)=0.83`$ for Perseus, not far from the value of $`\eta (2)=0.7`$ predicted analytically (Boldyrev et al., 2002a) from the velocity scaling in Boldyrev (2002).

## 3 Statistical Significance of High Order Moments

The velocity scalings of laboratory, geophysical or numerical flows are hardly verified up to the 10th order (Müller & Biskamp, 2000; She et al., 2001). In the present work we have used astronomical data to constrain moments beyond the 10th order. A way to verify the statistical significance of high order moments is to plot the integrand of the expression (3), and verify that at its peak position the PDF is defined by a significant number of samples. The peak position of the integrand of (3) grows with increasing order $`p`$ and converges to the value of the cutoff of the PDF. For a finite sample size the cutoff of the PDF may depend on the sample size, in which case high order moments would be inaccurate. However, as pointed out in Lévêque & She (1997), if a definite cutoff is present, for example if it is well defined by a significant number of samples, then the moment calculation up to any order should be accurate, and the fast convergence of high order moments is not an artifact of the finite sample size.

We have plotted in Figure 3, left panel, the tail of the PDF of $`|\mathrm{\Delta }I|=|I(𝐫^{})I(𝐫^{}+𝚫𝐫^{})|`$ from the Taurus data, for $`|𝚫𝐫^{}|=2.6`$ pc. Vertical dotted segments mark the peak position of the integrand of (3) for moments of increasing order, from $`p=3`$ to $`p=15`$. The moments converge at $`p=15`$, and are therefore accurate up to that order. However, the PDF shows a very sharp cutoff. The cutoff is definitely significant since it is defined by a few hundred samples. The ‘real’ PDF describing the physical process may therefore have a similarly sharp cutoff, in which case the convergence around the 15th order found in the observational data would be significant, and orders above the 15th would be accurately estimated by the present data set. The moment convergence corresponds to the asymptote of the function $`\eta (p)/\eta (3)`$ plotted in Figure 2. A similar result regarding the PDF cutoff is obtained for the Perseus data (right panel of Figure 3). However, in this case the peak of the integrand of (3) for every order is very shallow. Although formally the peak position has converged already at the 5th order (right panel of Figure 3), the moments converge around the 8th order. The convergence may again be well defined, as for the Taurus data, since the PDF cutoff is defined by a large sample size.

## 4 The Hierarchical Structure Model for Super–Sonic Turbulence

The scaling of the velocity structure functions in incompressible turbulence are best described by the She–Lévêque (1994) formula:

$$\frac{\zeta (p)}{\zeta (3)}=\gamma p+C(1\beta ^p)$$
(5)
where

$$C=\frac{13\gamma }{1\beta ^3}$$
(6)
is interpreted as the Hausdorff codimension of the support of the most singular dissipative structures. In incompressible turbulence the most dissipative structures are organized in filaments along coherent vortex tubes with Hausdorff dimension $`D=1`$ and so $`C=2`$. Furthermore, $`\beta ^3=2/3`$ (She & Lévêque, 1994), which yields the She–Lévêque’s velocity scaling for incompressible turbulence:

$$\frac{\zeta (p)}{\zeta (3)}=p/9+2\left[1\left(\frac{2}{3}\right)^{p/3}\right]$$
(7)
Boldyrev (2002) has proposed to apply the scaling of (5) to super–sonic turbulence, with the assumption that the Hausdorff dimension of the support of the most singular dissipative structures is $`D=2`$ ($`C=1`$), since dissipation of super–sonic turbulence occurs mainly in sheet–like shocks. Using the physical interpretation of (5) by Dubrulle (1994), $`\gamma =1/9`$ and $`C=(2/3)/(1\beta ^3)`$. With $`C=1`$, therefore, $`\beta ^3=1/3`$ and one obtains the Boldyrev’s velocity scaling:

$$\frac{\zeta (p)}{\zeta (3)}=p/9+1\left(\frac{1}{3}\right)^{p/3}$$
(8)
This velocity scaling has been found to provide a very accurate prediction for numerical simulations of super–sonic and super–Alfvénic turbulence (Boldyrev et al., 2002a, b). The computation of the corresponding structure functions for the gas density distribution is discussed in Boldyrev et al. (2002a), where an exponent $`\eta (2)=0.7`$ is predicted for the second order structure function of the projected density. In this paper we have found the ‘provoking’ result that the projected density field of super–sonic turbulence in the ISM has a structure function that is almost indistinguishable from Boldyrev’s velocity scaling.

This result should inspire the mathematical work. As shown by Dubrulle (1994) and by She & Waymire (1995), a hierarchy of structures producing the scaling relation (5) can be obtained by a random multiplicative process with Log–Poisson statistics. The result of this paper suggests that the density field of super–sonic turbulence is fragmented by a multiplicative process with Log–Poisson statistics.

### 4.1 The $`\beta `$ and $`\gamma `$ Tests for the Hierarchical Model

In the previous section the scaling of the structure functions of integrated intensity of two maps of the Taurus and Perseus molecular cloud complexes has been compared with the velocity scaling of super–sonic turbulence. A different approach would be to use observational, experimental or numerical data to determine the parameters $`\beta `$ and $`\gamma `$ from a fit to (5). There is, however, an independent method to first estimate the parameter $`\beta `$ and then obtain $`\gamma `$ from (5), which also tests the validity of two major assumptions of the hierarchical structure model at the same time. This method is presented in She et al. (2001); we apply it in the following.

The basic assumption in the derivation of (5) by She & Lévêque (1994) is the existence of the universal scaling behavior:

$$F_{p+1}(r)=A_pF_p(r)^\beta F_{(\mathrm{})}(r)^{1\beta }$$
(9)
where

$$F_p(r)=S_{p+1}(r)/S_p(r)=|\mathrm{\Delta }I|^{p+1}/|\mathrm{\Delta }I|^p$$
(10)
is usually referred to as the $`p`$th order intensity of fluctuations and $`A_p`$ are constants independent of $`r`$ and found to be also independent of $`p`$ in a number of cases.

The $`\beta `$ test verifies the validity of this basic assumption with a log-log plot of $`F_{p+1}/F_2`$ versus $`F_p/F_1`$. If the plot is a straight line the assumed hierarchical symmetry (9) is satisfied and the data pass the $`\beta `$ test. The plot is shown in Figure 5, top panel, for both Perseus (diamond symbols) and Taurus (asterisk symbols). We have combined together values of $`F_p(r)`$ for all values of $`r`$ used to compute the moment scaling exponents $`\eta (p)`$ (approximately the range of scales between 0.3 and 3 pc). Both Taurus and Perseus data apparently pass the $`\beta `$ test.

The constants $`A_p`$ are independent of $`p`$, as for the hierarchy of intensity of fluctuations of velocity in the shell model analyzed by Lévêque & She (1997), in the turbulent Couette–Taylor flow studied by She et al. (2001) and in laboratory data obtained by Chavarria et al. (1995). The slope of the intensity of fluctuation function plotted in Figure 5 provides an estimate of the value of $`\beta =0.79`$ for Taurus and $`\beta =0.78`$ for Perseus. This is to be compared with the value of $`\beta =0.69`$ of the BSL’s velocity scaling for super–sonic turbulence.

The second basic assumption by She & Lévêque (1994) is:

$$F_{\mathrm{}}S_3^\gamma $$
(11)
which allows them to derive (5). The $`\gamma `$ test verifies directly the validity of (5) by first assuming the value of $`\beta `$ obtained from the $`\beta `$ test and then by plotting $`\eta (p)\chi (p,\beta )`$ versus $`p3\chi (p,\beta )`$, where $`\chi (p,\beta )=(1\beta ^p)/(1\beta ^3)`$. If the plot is a straight line the data pass the $`\gamma `$ test and the slope of the plot provides an estimate of the value of $`\gamma `$. The plot is shown in Figure 5, bottom panel. We are able to fit the plot to a straight line for orders $`p>11`$ for both Taurus and Perseus. Since $`\gamma `$ is a property of the very high order moments (see definition (11)), and since we have obtained a straight line for large moments, we can say that the Taurus and Perseus data pass the $`\gamma `$ test. We estimate $`\gamma =0.11`$ for Taurus and $`\gamma =0.06`$ for Perseus.

## 5 Conclusions

We have computed the structure functions of the integrated intensity images of two J=1-0 13CO maps of Taurus and Perseus. The structure functions scale as power laws within the range of scales 0.3–3 pc. The scaling exponents have been computed up to the 20th order and are statistically significant at least up to the 15th order in Taurus and 8th order in Perseus. They are found to follow very closely the velocity scaling of super–sonic turbulence proposed by Boldyrev (2002), within 5% and 10% accuracy for Taurus and Perseus respectively.

We have verified that the projected density field (or integrated intensity) of the Taurus and Perseus molecular cloud complexes can be described by a hierarchical model as the one proposed by She & Lévêque (1994) for the velocity structure functions of incompressible turbulence. We have done so by testing the validity of the two basic assumptions of the hierarchical model for our data. The validity of the assumptions of the hierarchical model means that the integrated intensity images we have analyzed (an approximate estimate of the projected density of the Taurus and Perseus molecular cloud complexes) are the result of a multiplicative process with Log–Poisson statistics (Dubrulle, 1994).

The complete derivation of the relation between the structure functions of velocity and projected density is the subject of future work. However, the close similarity of the structure functions of projected density in Taurus and Perseus with that of the velocity field of turbulence provides additional evidence that super–sonic turbulence is the major factor controlling the density field in the range of densities and scales sampled by the maps we have analyzed.

It is well established that super–sonic turbulence plays an important role in the dynamics of the cold ISM (Larson, 1981; Padoan et al., 1998b, 1999, 2001; Padoan & Nordlund, 1999) The statistical properties of this ISM turbulence need to be discovered and understood in order to elaborate a statistical theory of star formation (Padoan & Nordlund, 2002) We have shown in the present paper that this can be achieved with existing observational data, by studying the projected density field of molecular clouds. This type of study will be greatly improved by far infrared imaging of the dust thermal emission from turbulent ISM clouds, obtained by future satellite missions such as SIRTF.

The work of PP was performed while PP held a National Research Council Associateship Award at the Jet Propulsion Laboratory, California Institute of Technology.

Figure captions:

Figure 1: Structure functions of orders $`p=1`$ to to $`p=20`$ of the integrated intensity J=1-0 13CO maps of Taurus (Mizuno et al., 1995) and Perseus (Billawala et al., 1997). The solid lines are least square fits to the structure functions in the approximate range 0.3–3 pc, where the structure functions are well described as power laws.

Figure 2: Top panel: Structure function scaling exponents relative to the third order, $`\eta _3(p)=\eta (p)/\eta (3)`$, up to the 20th order, for Taurus (asterisk symbols), Perseus (diamond symbols), Boldyrev’s velocity scaling for super–sonic turbulence (solid line), She–Lévêque’s velocity scaling for incompressible turbulence (dashed line) and Kolmogorov’s velocity scaling (dotted line). The moments of the structure functions of the integrated intensity images of Taurus and Perseus follow very closely Boldyrev’s velocity scaling, within 5% accuracy up to the 20th order (and beyond) for Taurus and within 10% accuracy for Perseus. Bottom panel: Same as top panel, but relative to Boldyrev’s scaling, $`\eta _3(p)_{BSL}`$.

Figure 3: Left Panel: Probability density function (PDF) of the integrated intensity differences of the Taurus map. Differences are computed at $`\mathrm{\Delta }r=2.6`$ pc. Vertical dotted segments show the position of the peak of the integrand of (3) for orders $`p=3`$ to 15. The PDF defines well all moments up to the 15th order. Higher moments are converged due to the steep cutoff of the PDF. Right panel: Same as left panel, but for the Perseus map and for $`\mathrm{\Delta }r=2.8`$ pc. The PDF of Perseus has a very sharp cutoff at large values, which causes a fast convergence of high order moments. The PDF of Perseus is such that the peak of the integrand of (3) is not very pronounced and a large interval of values of the PDF give therefore a significant contribution to the integral.

Figure 4: Top Panel: The $`\beta `$ test (see § 4.1) that verifies the validity of the hierarchical structure model. Intensity of fluctuations from both Taurus and Perseus are well fit by a straight line in this log-log plot. The data therefore pass the $`\beta `$ test. The value of the parameter $`\beta `$ given by the slope of the least square fits is $`\beta =0.79`$ for Taurus and $`\beta =0.78`$ for Perseus. Bottom panel: The $`\gamma `$ test (see § 4.1). The data is well fit by a power law at high orders. The solid lines are least square fits for orders within the range $`p=1220`$. The data therefore pas the $`\gamma `$ test. The slope of the least square fits provides the values $`\gamma =0.11`$ for Taurus and $`gamma=0.06`$ for Perseus.


format: text

text:
# UWThPh-2002-21 One-loop corrections to the seesaw mechanism in the multi-Higgs-doublet Standard Model

## 1 Introduction

It is believed that the final confirmation of small but non-zero neutrino masses is around the corner, since neutrino oscillations  provide an excellent and natural solution for the atmospheric- and solar-neutrino problems. In the latter case, matter effects in the neutrino oscillations  seem to play a crucial role. For reviews, see for instance Ref. ; for the latest solar-neutrino results, see Ref.  and the references therein.

The results of the solar- and atmospheric-neutrino experiments have brought about an upsurge of model building for the lepton masses and mixings—for reviews of lepton-mass-matrix textures and models see Ref. . In this effort, the seesaw mechanism  figures prominently as a means for obtaining small Majorana masses for the neutrinos.

In this paper we discuss a very simple extension of the Standard Model (SM): we add to the SM right-handed neutrino singlets, and moreover admit an arbitrary number of Higgs doublets. The first addition allows to incorporate the seesaw mechanism in the SM. As shown in Ref. , this simple extension of the SM provides a framework in which interesting models can be constructed, which reproduce in a natural way at least part of the striking features of the neutrino masses and mixings. This extension of the SM has two large mass scales: the electroweak scale, of about 100 GeV, which is the order of magnitude of the $`Z`$ and $`W`$ boson masses $`m_Z`$ and $`m_W`$, respectively, and also of the scalar masses; and the seesaw scale $`m_R`$, much larger than $`m_Z`$, which is the scale of the singlet-neutrino masses. A small mass scale $`m_D`$, below the electroweak scale, gives the order of magnitude of the elements of the Dirac neutrino mass matrix $`M_D`$, in such a way that the light-neutrino masses end up being of order $`m_D^2/m_R`$.

It was noticed in Ref.  that, in this framework, the dominant radiative corrections to the seesaw mechanism are quadratic in the Yukawa couplings, just like the tree-level masses, and that the radiative corrections are smaller than the tree-level results solely because of the appearance of the factor $`\left(16\pi ^2\right)^1`$ in one-loop integrals. By using two Higgs doublets and only one right-handed neutrino singlet, this fact was later exploited in Ref.  for the radiative generation of the solar mass-squared difference; the ratio $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{\mathrm{atm}}^2`$ of solar over atmospheric mass-squared differences comes out of the correct order of magnitude for the large-mixing-angle MSW solution of the solar-neutrino problem (for the main features of that solution, with $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{\mathrm{atm}}^2`$ being about 0.02, see for instance Ref. ), with no other scale free to make adjustments.

We want to extend here the calculations of Ref. , where the focus was the radiative generation of neutrino masses vanishing at tree level. In the present work we perform a general calculation of the dominant one-loop corrections to the seesaw mechanism, with emphasis on the demonstration of the finiteness and gauge-independence (in $`R_\xi `$ gauges) of those radiative corrections. As a result of this effort we shall present a formula, easily applicable in our extension of the SM, whose purpose is twofold:

* To check the radiative stability of the seesaw mass matrix and of the relations derived therefrom, or else to compute the radiative corrections to those relations;
* To explore the possibilities of radiatively generating $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{\mathrm{atm}}^2`$.

Unfortunately, the second point appears quite difficult to us, when one wants to combine it with specific properties of the lepton mixing matrix $`U`$, like the small $`U_{e3}`$ and the practically maximal atmospheric mixing.

This paper is organized as follows. In Section 2 we review the seesaw mechanism and our formalism for the multi-Higgs-doublet SM. We discuss the calculation of the neutrino self-energies and write down their one-loop integrals in Section 3. The actual calculation of the dominant one-loop corrections to the seesaw mechanism and the demonstration of their finiteness and gauge-independence are carried out in Section 4. In Section 5 we give arguments in favour of the dominance of the corrections calculated in the previous section. In Section 6 we present the final results and a summary.

## 2 Framework

We consider the lepton sector of an extension of the Standard Model with $`n_L`$ left-handed doublets and $`n_L`$ right-handed charged singlets, i.e., $`n_L`$ families, plus $`n_R`$ right-handed neutrino fields and $`n_H`$ Higgs doublets. We denote the Higgs doublets by $`\varphi _k`$ ($`k=1,2,\mathrm{},n_H`$) and define $`\stackrel{~}{\varphi }_k=i\tau _2\varphi _k^{}`$. The vacuum expectation value (VEV) of the neutral component of $`\varphi _k`$ is $`v_k/\sqrt{2}`$. The Yukawa Lagrangian of the leptons is given by

$$_\mathrm{Y}=\underset{k=1}{\overset{n_H}{}}\left(\varphi _k^{}\overline{\mathrm{}}_R\mathrm{\Gamma }_k+\stackrel{~}{\varphi }_k^{}\overline{\nu }_R\mathrm{\Delta }_k\right)D_L+\mathrm{H}.\mathrm{c}.$$
(1)
We employ a vector and matrix notation where $`\mathrm{}_R`$, $`\nu _R`$, and $`D_L`$ are the vectors of the right-handed charged-lepton fields, of the right-handed neutrino singlets, and of the left-handed lepton doublets, respectively. The Yukawa coupling matrices $`\mathrm{\Gamma }_k`$ are $`n_L\times n_L`$, while the $`\mathrm{\Delta }_k`$ are $`n_R\times n_L`$. The charged-lepton mass matrix $`M_{\mathrm{}}`$ and the Dirac neutrino mass matrix $`M_D`$ are

$$M_{\mathrm{}}=\frac{1}{\sqrt{2}}\underset{k}{}v_k^{}\mathrm{\Gamma }_k\mathrm{and}M_D=\frac{1}{\sqrt{2}}\underset{k}{}v_k\mathrm{\Delta }_k,$$
(2)
respectively. Without loss of generality, we assume $`M_{\mathrm{}}`$ to be diagonal with real and positive diagonal elements: $`M_{\mathrm{}}=\mathrm{diag}(m_e,m_\mu ,m_\tau ,\mathrm{})`$. The mass terms for the neutrinos are

$$\overline{\nu }_RM_D\nu _L\frac{1}{2}\overline{\nu }_RCM_R\overline{\nu }_R^T+\mathrm{H}.\mathrm{c}.,$$
(3)
where $`C`$ is the charge-conjugation matrix and $`M_R`$ is non-singular and symmetric. It is well known  that Eq. (3) can be written in a compact form as a mass term with an $`(n_L+n_R)\times (n_L+n_R)`$ symmetric mass matrix

$$_{D+M}=\left(\begin{array}{cc}0& M_D^T\\ M_D& M_R\end{array}\right),$$
(4)
which may be diagonalized as

$$𝒰^T_{D+M}𝒰=\widehat{m}=\mathrm{diag}(m_1,m_2,\mathrm{},m_{n_L+n_R}),$$
(5)
where the $`m_i`$ are real and non-negative. In order to implement the seesaw mechanism  we assume that the elements of $`M_D`$ are of order $`m_D`$ and those of $`M_R`$ are of order $`m_R`$, with $`m_Dm_R`$. Then, the neutrino masses $`m_i`$ with $`i=1,2,\mathrm{},n_L`$ are of order $`m_D^2/m_R`$, while those with $`i=n_L+1,\mathrm{},n_L+n_R`$ are of order $`m_R`$. It is useful to decompose the $`(n_L+n_R)\times (n_L+n_R)`$ unitary matrix $`𝒰`$ as

$$𝒰=\left(\begin{array}{c}U_L\\ U_R^{}\end{array}\right),$$
(6)
where the submatrix $`U_L`$ is $`n_L\times (n_L+n_R)`$ and the submatrix $`U_R`$ is $`n_R\times (n_L+n_R)`$. With these submatrices, the left- and right-handed neutrinos are written as linear superpositions of the $`n_L+n_R`$ physical Majorana neutrino fields $`\chi _i`$:

$$\nu _L=U_L\gamma _L\chi \mathrm{and}\nu _R=U_R\gamma _R\chi ,$$
(7)
where $`\gamma _L=\left(1\gamma _5\right)/2`$ and $`\gamma _R=\left(1+\gamma _5\right)/2`$ are the projectors of chirality.

The unitarity relations for $`𝒰`$ are

$$U_LU_L^{}=\mathrm{𝟙}_{n_L},U_RU_R^{}=\mathrm{𝟙}_{n_R},U_LU_R^T=0_{n_L\times n_R},\mathrm{and}U_L^{}U_L+U_R^TU_R^{}=\mathrm{𝟙}_{n_L+n_R}.$$
(8)
From Eq. (5) we shall need the two relations

$$U_L^{}\widehat{m}U_L^{}=0\mathrm{and}U_R\widehat{m}U_L^{}=M_D.$$
(9)
It follows from Eqs. (8) and (9) that

$$U_R^{}M_D=U_R^{}U_R\widehat{m}U_L^{}=\left(\mathrm{𝟙}_{n_L+n_R}U_L^TU_L^{}\right)\widehat{m}U_L^{}=\widehat{m}U_L^{}.$$
(10)
The leptonic charged-current Lagrangian is

$$_{\mathrm{cc}}=\frac{g}{\sqrt{2}}W_\mu ^{}\overline{\mathrm{}}\gamma ^\mu \gamma _LU_L\chi +\mathrm{H}.\mathrm{c}.,$$
(11)
where $`g`$ is the $`SU(2)`$ gauge coupling constant. The interaction of the $`Z`$ boson with the neutrinos is given by

$$_{\mathrm{nc}}^{(\nu )}=\frac{g}{4c_w}Z_\mu \overline{\chi }\gamma ^\mu \left[\gamma _L\left(U_L^{}U_L\right)\gamma _R\left(U_L^TU_L^{}\right)\right]\chi ,$$
(12)
where $`c_w`$ is the cosine of the Weinberg angle.

A full account of our formalism for the scalar sector of the multi-Higgs-doublet SM is given in Ref. . We use the notation $`S_b^0`$ for the neutral-scalar mass eigenfields; the vectors $`b^{n_H}`$ characterize each of those fields—for their precise definitions see Refs. . The Yukawa couplings of the neutral scalars $`S_b^0`$ to the neutrinos are given by

$$_\mathrm{Y}^{(\nu )}\left(S^0\right)=\frac{1}{2\sqrt{2}}\underset{b}{}S_b^0\overline{\chi }\left[\left(U_R^{}\mathrm{\Delta }_bU_L+U_L^T\mathrm{\Delta }_b^TU_R^{}\right)\gamma _L+\left(U_L^{}\mathrm{\Delta }_b^{}U_R+U_R^T\mathrm{\Delta }_b^{}U_L^{}\right)\gamma _R\right]\chi ,$$
(13)
with $`\mathrm{\Delta }_b=_kb_k\mathrm{\Delta }_k`$. The neutral Goldstone boson $`G^0=S_{b_Z}^0`$ is given by the vector $`b_Z`$ with $`\left(b_Z\right)_k=iv_k/v`$ . Here, $`v=\left(\left|v_1\right|^2+\left|v_2\right|^2+\mathrm{}+\left|v_{n_H}\right|^2\right)^{1/2}=2m_W/g`$. We thus obtain

$$\mathrm{\Delta }_{b_Z}=\frac{ig}{\sqrt{2}m_W}M_D.$$
(14)
Similarly, the charged-scalar mass eigenfields $`S_a^\pm `$ are characterized by vectors $`a^{n_H}`$. The Yukawa couplings of the charged scalars are given by

$$_\mathrm{Y}\left(S^\pm \right)=\underset{a}{}S_a^{}\overline{\mathrm{}}\left[\gamma _R\left(\mathrm{\Delta }_a^{}U_R\right)\gamma _L\left(\mathrm{\Gamma }_aU_L\right)\right]\chi +\mathrm{H}.\mathrm{c}.,$$
(15)
where $`\mathrm{\Delta }_a=_ka_k\mathrm{\Delta }_k`$ and $`\mathrm{\Gamma }_a=_ka_k^{}\mathrm{\Gamma }_k`$. The charged Goldstone boson $`G^\pm =S_{a_W}^\pm `$ is given by the vector $`a_W`$ with $`\left(a_W\right)_k=v_k/v`$; therefore

$$\mathrm{\Delta }_{a_W}=\frac{g}{\sqrt{2}m_W}M_D\mathrm{and}\mathrm{\Gamma }_{a_W}=\frac{g}{\sqrt{2}m_W}M_{\mathrm{}}.$$
(16)
## 3 Neutrino self-energies

The seesaw mechanism  tells us that, at tree level, the mass matrix of the light neutrinos is given by

$$_\nu ^{\mathrm{tree}}=M_D^TM_R^1M_D,$$
(17)
where $`M_D`$ and $`M_R`$ appear in the full mass matrix $`_{D+M}`$ of Eq. (4). The one-loop corrections generate an $`n_L\times n_L`$ submatrix $`\delta M_L`$ in the upper left-hand corner of $`_{D+M}`$, where at tree level there is a zero submatrix. The submatrices $`M_D`$ and $`M_R`$ also receive one-loop corrections, denoted $`\delta M_D`$ and $`\delta M_R`$, respectively. We shall then have, for the effective light-neutrino mass matrix at one loop,

$$_\nu =_\nu ^{\mathrm{tree}}+\delta M_L\delta M_D^TM_R^1M_DM_D^TM_R^1\delta M_D+M_D^TM_R^1\delta M_RM_R^1M_D.$$
(18)
One-loop corrections to the neutrino masses originate in the one-loop neutrino self-energy $`\mathrm{\Sigma }(p)`$, where $`p`$ is the neutrino momentum. We use the decomposition

$$\mathrm{\Sigma }(p)=A_L\left(p^2\right)p/\gamma _L+A_R\left(p^2\right)p/\gamma _R+B_L\left(p^2\right)\gamma _L+B_R\left(p^2\right)\gamma _R.$$
(19)
The non-absorptive parts of $`A_{L,R}`$ are Hermitian, while those of $`B_{L,R}`$ are the Hermitian conjugates of each other. Furthermore, since the neutrino field vector $`\chi `$ consists of Majorana fields, the self-energy must fulfil the consistency condition

$$\mathrm{\Sigma }(p)=C\left[\mathrm{\Sigma }(p)\right]^TC^1,$$
(20)
and that consistency condition translates into

$$A_L=A_R^T,B_L=B_L^T,\mathrm{and}B_R=B_R^T.$$
(21)
After computing the neutrino self-energies one must renormalize the neutrino fields:

$$\chi _L=\left(1+\frac{1}{2}z_L\right)\chi _L^r,\mathrm{hence}\chi _R=\left(1+\frac{1}{2}z_L^{}\right)\chi _R^r,$$
(22)
since $`\chi _R=C\overline{\chi }_L^T`$. The neutrino mass matrices must also be renormalized. Mass counterterms are indicated by $`\delta ^c`$. The upper-left submatrix of $`_{D+M}`$, which we call $`M_L`$, has no counterterm since it vanishes at tree level: $`\delta ^cM_L=0`$. The operation $`\delta ^c`$ only acts on $`M_D`$ and on $`M_R`$, with, from the second Eq. (2),

$$\delta ^cM_D=\frac{1}{\sqrt{2}}\underset{k}{}\left[\left(\delta ^cv_k\right)\mathrm{\Delta }_k+v_k\left(\delta ^c\mathrm{\Delta }_k\right)\right].$$
(23)
Then,

$$\delta ^c\widehat{m}=𝒰^T\left(\begin{array}{cc}0& \left(\delta ^cM_D\right)^T\\ \delta ^cM_D& \delta ^cM_R\end{array}\right)𝒰.$$
(24)
The *renormalized* neutrino self-energy, in the basis where the tree-level neutrino mass matrix is diagonal, has the structure

$`i\mathrm{\Sigma }^r(p)`$ $`=`$ $`i\mathrm{\Sigma }(p)+i{\displaystyle \frac{1}{2}}\left(z_L+z_L^{}\right)p/\gamma _L+i{\displaystyle \frac{1}{2}}\left(z_L^{}+z_L^T\right)p/\gamma _R`$ (25)
$`i\left(\delta ^c\widehat{m}+{\displaystyle \frac{1}{2}}\widehat{m}z_L+{\displaystyle \frac{1}{2}}z_L^T\widehat{m}\right)\gamma _Li\left(\delta ^c\widehat{m}+{\displaystyle \frac{1}{2}}\widehat{m}z_L+{\displaystyle \frac{1}{2}}z_L^T\widehat{m}\right)^{}\gamma _R.`$
We indicate one-loop contributions to $`\mathrm{\Sigma }(p)`$ by $`\delta ^{1\text{-}\mathrm{loop}}`$. There are tadpole diagrams indicated by $`\delta ^{\mathrm{tadpole}}`$, but we do not include them in $`\mathrm{\Sigma }(p)`$.

Now we explain our strategy for calculating the dominant one-loop corrections to the seesaw mechanism:

1. We calculate $`\mathrm{\Sigma }(p)`$ in the basis where the tree-level neutrino mass matrix is diagonal.
2. We are especially interested in $`B_L\left(p^2\right)`$, because

$$\delta M_L=\delta ^{1\text{-}\mathrm{loop}}M_L=U_L^{}B_L\left(0\right)U_L^{}.$$
(26)
We have used the inverse of the relation (5) to transform back to the basis of $`_{D+M}`$.
3. Denoting by $`\mathrm{\Delta }`$ a typical coupling constant of the coupling matrices $`\mathrm{\Delta }_k`$, the light-neutrino masses are of order $`\mathrm{\Delta }^2`$—see Eq. (17). Therefore, we are allowed to evaluate $`B_L`$ at $`p^2=0`$ in Eq. (26). Note that $`\delta M_L`$ is identical with $`\delta ^{1\text{-}\mathrm{loop}}M_L`$, since $`\delta ^cM_L=0`$. For the same reason, $`\delta ^{1\text{-}\mathrm{loop}}M_L`$ must be finite. We shall see this explicitly in the next section, where we shall also demonstrate the gauge invariance of $`\delta ^{1\text{-}\mathrm{loop}}M_L`$.
4. On the other hand, the one-loop corrections to $`M_D`$ and to $`M_R`$ are given by

$$\delta ^{1\text{-}\mathrm{loop}}M_D=U_RB_L(0)U_L^{}\mathrm{and}\delta ^{1\text{-}\mathrm{loop}}M_R=U_RB_L(0)U_R^T,$$
(27)
respectively. We shall argue in Section 5 that one-loop corrections to $`M_R`$ are irrelevant, while those to $`M_D`$ lead to subdominant corrections in $`_\nu `$.

First we consider the contribution of the one-loop graph with $`W^\pm `$ to $`\mathrm{\Sigma }(p)`$:

$`\mathrm{\Sigma }_{ij}^{(W)}(p)`$ $`=`$ $`i{\displaystyle \frac{g^2}{2}}{\displaystyle \frac{d^dk}{\left(2\pi \right)^d}\left[S_{\mu \nu }^W\left(kp\right)\right]\underset{\mathrm{}}{}\frac{1}{k^2m_{\mathrm{}}^2}}`$ (28)
$`\times \left[\left(U_L^{}\right)_i\mathrm{}\left(U_L\right)_\mathrm{}j\gamma ^\mu k/\gamma ^\nu \gamma _L+\left(U_L^T\right)_i\mathrm{}\left(U_L^{}\right)_\mathrm{}j\gamma ^\mu k/\gamma ^\nu \gamma _R\right],`$
where $`iS_{\mu \nu }^W\left(kp\right)`$ is the propagator of a $`W^\pm `$ with momentum $`kp`$, and $`d`$ is the dimension of space–time. Clearly, $`\mathrm{\Sigma }_{ij}^{(W)}(p)`$ only contributes to $`A_L\left(p^2\right)`$ and $`A_R\left(p^2\right)`$.

Next we consider the contribution from the $`Z`$ boson. This yields, in a general $`R_\xi `$ gauge parameterized by the gauge parameter $`\xi _Z`$,

$`\left(B_L\right)_{ij}^{(Z)}\left(p^2\right)`$ $`=`$ $`i{\displaystyle \frac{g^2}{4c_w^2}}{\displaystyle \frac{d^dk}{\left(2\pi \right)^d}\underset{\iota }{}\frac{m_\iota }{k^2m_\iota ^2}\left(U_L^TU_L^{}\right)_{i\iota }\left(U_L^{}U_L\right)_{\iota j}}`$ (29)
$`\times \left\{{\displaystyle \frac{d}{\left(kp\right)^2m_Z^2}}+{\displaystyle \frac{1}{m_Z^2}}\left[{\displaystyle \frac{\left(kp\right)^2}{\left(kp\right)^2\xi _Zm_Z^2}}{\displaystyle \frac{\left(kp\right)^2}{\left(kp\right)^2m_Z^2}}\right]\right\}.`$
Now we consider the contribution of the one-loop graph with $`S_a^\pm `$. It is

$`\left(B_L\right)_{ij}^{(S^\pm )}\left(p^2\right)`$ $`=`$ $`i{\displaystyle \frac{d^dk}{\left(2\pi \right)^d}\underset{a}{}\frac{1}{\left(kp\right)^2m_a^2}\underset{\mathrm{}}{}\frac{m_{\mathrm{}}}{k^2m_{\mathrm{}}^2}}`$ (30)
$`\times \left[\left(U_R^{}\mathrm{\Delta }_a\right)_i\mathrm{}\left(\mathrm{\Gamma }_aU_L\right)_\mathrm{}j+\left(U_L^T\mathrm{\Gamma }_a^T\right)_i\mathrm{}\left(\mathrm{\Delta }_a^TU_R^{}\right)_\mathrm{}j\right].`$
When $`S_a^\pm `$ is the charged Goldstone boson $`G^\pm `$, one uses Eqs. (16) together with Eq. (10) to derive

$`\left(B_L\right)_{ij}^{(G^\pm )}\left(p^2\right)`$ $`=`$ $`i{\displaystyle \frac{g^2}{2m_W^2}}{\displaystyle \frac{d^dk}{\left(2\pi \right)^d}\frac{1}{\left(kp\right)^2\xi _Wm_W^2}\underset{\mathrm{}}{}\frac{m_{\mathrm{}}^2}{k^2m_{\mathrm{}}^2}}`$ (31)
$`\times \left[m_i\left(U_L^{}\right)_i\mathrm{}\left(U_L\right)_\mathrm{}j+\left(U_L^T\right)_i\mathrm{}\left(U_L^{}\right)_\mathrm{}jm_j\right].`$
In similar fashion, the neutral scalars $`S_b^0`$ contribute (neglecting tadpoles)

$`\left(B_L\right)_{ij}^{(S^0)}\left(p^2\right)`$ $`=`$ $`{\displaystyle \frac{i}{2}}{\displaystyle \frac{d^dk}{\left(2\pi \right)^d}\underset{b}{}\frac{1}{\left(kp\right)^2m_b^2}\underset{\iota }{}\frac{m_\iota }{k^2m_\iota ^2}}`$ (32)
$`\times \left(U_R^{}\mathrm{\Delta }_bU_L+U_L^T\mathrm{\Delta }_b^TU_R^{}\right)_{i\iota }\left(U_R^{}\mathrm{\Delta }_bU_L+U_L^T\mathrm{\Delta }_b^TU_R^{}\right)_{\iota j}.`$
If the scalar $`S_b^0`$ is the neutral Goldstone boson $`G^0`$, then, using Eqs. (14) and (10),

$`\left(B_L\right)_{ij}^{(G^0)}\left(p^2\right)`$ $`=`$ $`{\displaystyle \frac{ig^2}{4m_W^2}}{\displaystyle \frac{d^dk}{\left(2\pi \right)^d}\frac{1}{\left(kp\right)^2\xi _Zm_Z^2}\underset{\iota }{}\frac{m_\iota }{k^2m_\iota ^2}}`$ (33)
$`\times \left(\widehat{m}U_L^{}U_L+U_L^TU_L^{}\widehat{m}\right)_{i\iota }\left(\widehat{m}U_L^{}U_L+U_L^TU_L^{}\widehat{m}\right)_{\iota j}.`$
## 4 The computation of $`\delta M_L`$

In this section, we explicitly perform the calculation of $`\delta M_L`$, defined in Eq. (26).

It is clear from Eq. (30) and from the third Eq. (8) that the exchange of $`S^\pm `$ contributes neither to $`\delta ^{1\text{-}\mathrm{loop}}M_L`$ nor to $`\delta ^{1\text{-}\mathrm{loop}}M_R`$; it only contributes to $`\delta ^{1\text{-}\mathrm{loop}}M_D`$. We therefore have

$$\delta M_L=\delta M_L\left(Z\right)+\delta M_L\left(G^0\right)+\underset{bb_Z}{}\delta M_L\left(S_b^0\right).$$
(34)
The $`Z`$ boson part is given by

$$\delta M_L\left(Z\right)=i\frac{g^2}{4c_w^2}\frac{d^dk}{\left(2\pi \right)^d}\left[\frac{d}{k^2m_Z^2}+\frac{1}{m_Z^2}\left(\frac{k^2}{k^2\xi _Zm_Z^2}\frac{k^2}{k^2m_Z^2}\right)\right]U_L^{}\frac{\widehat{m}}{k^2\widehat{m}^2}U_L^{}.$$
(35)
The neutral Goldstone boson contributes

$$\delta M_L\left(G^0\right)=\frac{ig^2}{4m_W^2}\frac{d^dk}{\left(2\pi \right)^d}\frac{1}{k^2\xi _Zm_Z^2}U_L^{}\frac{\widehat{m}^3}{k^2\widehat{m}^2}U_L^{}.$$
(36)
Using $`m_W^2=c_w^2m_Z^2`$ and the relation

$$U_L^{}\frac{k^2\widehat{m}}{k^2\widehat{m}^2}U_L^{}=U_L^{}\frac{\widehat{m}^3}{k^2\widehat{m}^2}U_L^{},$$
(37)
which holds by virtue of $`U_L^{}\widehat{m}U_L^{}=0`$, it is clear that the $`\xi _Z`$-dependent terms in Eqs. (35) and (36) cancel and therefore *$`\delta M_L`$ is gauge-invariant*. One has

$$\delta M_L\left(Z\right)+\delta M_L\left(G^0\right)=i\frac{g^2}{4c_w^2}\frac{d^dk}{\left(2\pi \right)^d}\frac{1}{k^2m_Z^2}U_L^{}\left(d\frac{\widehat{m}^2}{m_Z^2}\right)\frac{\widehat{m}}{k^2\widehat{m}^2}U_L^{}.$$
(38)
The contribution to $`\delta M_L`$ from the exchange of all neutral scalars, except the neutral Goldstone boson, is

$$\delta M_L\left(S_b^0\right)=\frac{i}{2}\frac{d^dk}{\left(2\pi \right)^d}\frac{1}{k^2m_b^2}\mathrm{\Delta }_b^TU_R^{}\frac{\widehat{m}}{k^2\widehat{m}^2}U_R^{}\mathrm{\Delta }_b.$$
(39)
Now we perform the actual calculation of the integrals. Let the dimension of space–time be $`d=42ϵ`$, with $`ϵ0`$. We define the divergent quantity

$$k=ϵ^1+\gamma \mathrm{ln}\left(4\pi \right)1,$$
(40)
where $`\gamma `$ is Euler’s constant. Then, we obtain

$$\delta M_L\left(S_b^0\right)=\frac{1}{32\pi ^2}\mathrm{\Delta }_b^TU_R^{}\widehat{m}\left(k+\mathrm{ln}\widehat{m}^2+\frac{\mathrm{ln}r_b}{r_b1}\right)U_R^{}\mathrm{\Delta }_b,\mathrm{with}r_b=\frac{\widehat{m}^2}{m_b^2}.$$
(41)
Similarly,

$$\delta M_L(Z,1)=\frac{g^2}{16\pi ^2c_w^2}U_L^{}\widehat{m}\left(k+\frac{1}{2}+\mathrm{ln}\widehat{m}^2+\frac{\mathrm{ln}r_Z}{r_Z1}\right)U_L^{},\mathrm{with}r_Z=\frac{\widehat{m}^2}{m_Z^2},$$
(42)
is the contribution to $`\delta M_L\left(Z\right)+\delta M_L\left(G^0\right)`$ from the term with $`d`$ in Eq. (38); the term with $`\widehat{m}^2/m_Z^2`$ yields

$$\delta M_L(Z,2)=\frac{g^2}{64\pi ^2m_W^2}U_L^{}\widehat{m}^3\left(k+\mathrm{ln}\widehat{m}^2+\frac{\mathrm{ln}r_Z}{r_Z1}\right)U_L^{}.$$
(43)
In order to demonstrate the cancellation of infinities we need the orthogonality relation of the vectors $`b^{n_H}`$ :

$$\underset{bb_Z}{}b_jb_k+\left(b_Z\right)_j\left(b_Z\right)_k=0.$$
(44)
Remembering that $`\mathrm{\Delta }_b=_kb_k\mathrm{\Delta }_k`$, we find that all terms independent of the boson masses—in particular the infinities—drop out in the sum $`_{bb_Z}\delta M_L\left(S_b^0\right)+\delta M_L(Z,2)`$. We may thus write

$`\delta M_L^{}\left(S_b^0\right)`$ $`=`$ $`{\displaystyle \frac{1}{32\pi ^2}}\mathrm{\Delta }_b^TU_R^{}\widehat{m}{\displaystyle \frac{\mathrm{ln}r_b}{r_b1}}U_R^{}\mathrm{\Delta }_b,`$ (45)
$`\delta M_L^{}(Z,2)`$ $`=`$ $`{\displaystyle \frac{g^2}{64\pi ^2m_W^2}}U_L^{}\widehat{m}^3{\displaystyle \frac{\mathrm{ln}r_Z}{r_Z1}}U_L^{}.`$ (46)
Also, because of the first Eq. (9), the terms in $`\delta M_L(Z,1)`$ of Eq. (42) which are proportional to $`\widehat{m}`$ (but not the term with $`\widehat{m}\mathrm{ln}\widehat{m}^2`$), in particular the infinities, drop out. Thus,

$$\delta M_L^{}(Z,1)=\frac{g^2}{16\pi ^2c_w^2}U_L^{}\widehat{m}\frac{r_Z\mathrm{ln}r_Z}{r_Z1}U_L^{}.$$
(47)
It is obvious that $`\delta M_L^{}(Z,1)=4\delta M_L^{}(Z,2)`$ We thus we have the final result

$$\delta M_L=\underset{bb_Z}{}\frac{1}{32\pi ^2}\mathrm{\Delta }_b^TU_R^{}\widehat{m}\frac{\mathrm{ln}r_b}{r_b1}U_R^{}\mathrm{\Delta }_b+\frac{3g^2}{64\pi ^2m_W^2}U_L^{}\widehat{m}^3\frac{\mathrm{ln}r_Z}{r_Z1}U_L^{}.$$
(48)
Using Eq. (10), the last term of the previous relation reads

$$\frac{3g^2}{64\pi ^2m_W^2}M_D^TU_R^{}\widehat{m}\frac{\mathrm{ln}r_Z}{r_Z1}U_R^{}M_D.$$
(49)
This expression, together with $`m_W=gv/2`$, shows that the $`Z`$ contribution to $`\delta M_L`$ is of the same order of magnitude as the neutral-scalar contributions.

## 5 Dominance of $`\delta M_L`$ among the radiative corrections

We now want to give arguments for the dominance of $`\delta M_L`$ in the radiative corrections to $`_\nu ^{\mathrm{tree}}`$. The main point in our argumentation is the observation that the terms in $`\delta M_L`$ in Eq. (48) are smaller than the tree-level masses of the light neutrinos solely because of the factor $`\left(16\pi ^2\right)^1`$ from the one-loop integrals. In this context, we have to remember that, in our framework, it is natural to assume all scalar masses to be of the order of the electroweak scale.

The heavy-neutrino masses $`m_{n_L+1},\mathrm{},m_{n_L+n_R}`$ (or the elements of $`M_R`$) are free parameters of the theory; therefore, corrections to $`M_R`$ are irrelevant for us. We stress that this contrasts to the light-neutrino masses $`m_1,m_2,\mathrm{},m_{n_L}`$, which are calculable, i.e., they are functions of the parameters of the theory; this is reflected in the finiteness and gauge-invariance of $`\delta M_L`$.

Concentrating on $`\delta ^{1\text{-}\mathrm{loop}}M_D`$ as given by Eq. (27), we note that both $`S_a^\pm `$ and $`S_b^0`$ contribute to it, see Eqs. (30) and (32), respectively. Checking Eq. (29), we note the absence of a $`Z`$ contribution. In the integral of Eq. (30), for charged-scalar exchange, the large mass scale $`m_R`$ is absent; therefore, if $`Y`$ is a typical Yukawa coupling, then we must have $`\delta ^{1\text{-}\mathrm{loop}}M_D(S_a^\pm )Y^2m_{\mathrm{}}`$, where $`m_{\mathrm{}}`$ is a generic charged-lepton mass. With $`m_{\mathrm{}}m_D`$, we have $`Y^2m_D^2/m_R`$ as the order of magnitude of the charged-scalar contribution to $`_\nu `$, which is suppressed as compared to the $`\delta M_L`$ contribution. Turning now to $`S_b^0`$ exchange, with Eqs. (27) and (32) the result of the one-loop calculation is

$$\delta ^{1\text{-}\mathrm{loop}}M_D\left(S_b^0\right)=\frac{1}{32\pi ^2}\mathrm{\Delta }_bU_L\widehat{m}\left(k+\mathrm{ln}\widehat{m}^2+\frac{\mathrm{ln}r_b}{r_b1}\right)U_R^{}\mathrm{\Delta }_b.$$
(50)
Again, we arrive at $`\delta ^{1\text{-}\mathrm{loop}}M_D(S_b^0)Y^2m_D`$. Therefore, the contributions to $`_\nu `$ of Eq. (18) from $`\delta ^{1\text{-}\mathrm{loop}}M_D`$ are negligible when compared to the ones from $`\delta M_L`$.

Some remarks concerning the wave function renormalization matrix $`z_L`$ are in order. Usually, that matrix is defined in such a way that the renormalized fields remain in the mass basis . With Eqs. (17) and (26) we do not adopt this convention: the mass matrix of the light neutrinos, $`_\nu `$, is not diagonal at order $`Y^2`$. This is consistent with the fact that, at zeroth order in the Yukawa couplings, the light neutrinos are massless and, therefore, degenerate. This degeneracy is lifted only at order $`Y^2`$, and a large rotation, with angles of order one, will be needed, in general, for the diagonalization of $`_\nu `$. Since we assume that $`z_L`$ is “small” as compared to one, we do not include the large rotation in this matrix. In the same vein, concerning the terms $`\widehat{m}z_L+z_L^T\widehat{m}`$ in Eq. (25), we are allowed to set to zero the light-neutrino masses. Using the explicit form of $`z_L`$ for on-shell renormalization of the *heavy* neutrinos, as found, for instance, in Refs. , and taking into account that the small elements of $`𝒰`$ are of order $`m_D/m_R`$, we find that the wave function renormalization terms of Eq. (25) contribute at order $`Y^2m_D^2/m_R`$ or $`g^2m_D^2/m_R`$ to $`_\nu `$, and then we do not need to consider them further.

The tadpole graphs contribute only to $`\delta M_D`$. For their treatment and their relation to the $`\delta ^cv_k`$ terms in Eq. (23), see the extensive discussion in Ref. . They have no impact on our discussion.

## 6 Conclusions

We have argued in the preceding section that the dominant radiative corrections to $`_\nu ^{\mathrm{tree}}`$ are given by Eq. (26), which has contributions only from neutral-scalar and $`Z`$ exchange; equation (48) is the explicit form of Eq. (26) after carrying out the one-loop calculations. These radiative corrections are, like the tree-level mass term of Eq. (17), quadratic in the Yukawa couplings. In Eq. (48)—as compared to $`_\nu ^{\mathrm{tree}}`$—only the terms of first order in $`1/m_R`$ are relevant ($`m_R`$ is the scale of $`M_R`$). We may use the approximations

$$U_R(0,W),\mathrm{with}W^{}M_RW^{}\stackrel{~}{m}\mathrm{diag}(m_{n_L+1},\mathrm{},m_{n_L+n_R}),$$
(51)
where $`W`$ is a unitary $`n_R\times n_R`$ matrix whose elements are not suppressed by $`m_D/m_R`$ ($`m_D`$ is the scale of $`M_D`$). We thus have the final result

$$_\nu =M_D^TM_R^1M_D+\delta M_L,$$
(52)
with

$`\delta M_L`$ $`=`$ $`{\displaystyle \underset{bb_Z}{}}{\displaystyle \frac{m_b^2}{32\pi ^2}}\mathrm{\Delta }_b^TW^{}\left({\displaystyle \frac{1}{\stackrel{~}{m}}}\mathrm{ln}{\displaystyle \frac{\stackrel{~}{m}^2}{m_b^2}}\right)W^{}\mathrm{\Delta }_b`$ (53)
$`+{\displaystyle \frac{3g^2}{64\pi ^2c_w^2}}M_D^TW^{}\left({\displaystyle \frac{1}{\stackrel{~}{m}}}\mathrm{ln}{\displaystyle \frac{\stackrel{~}{m}^2}{m_Z^2}}\right)W^{}M_D.`$
The sum over the neutral scalars includes only the physical ones. The second line in Eq. (53) is the $`Z`$ contribution. Note that if accidentally $`_\nu ^{\mathrm{tree}}=M_D^TW^{}\left(1/\stackrel{~}{m}\right)W^{}M_D`$ vanishes, then, even for a single Higgs doublet, the one-loop $`Z`$ and Higgs contributions to $`\delta M_L`$ will in general be non-zero .

Having more than one Higgs doublet leads to flavour-changing neutral interactions. Neutral-scalar exchange at tree level then leads to processes like $`\mu ^\pm e^\pm e^+e^{}`$; however, such processes only restrict elements of the coupling matrices $`\mathrm{\Gamma }_k`$, which do not occur in $`_\nu `$ of Eq. (52). Moreover, with more than one Higgs doublet and at one-loop level, radiative decays like $`\mu ^\pm e^\pm \gamma `$ proceed via charged-scalar exchange. However, moderately large charged-scalar masses, still in the range of the electroweak scale, can render their decay rates compatible with experimental bounds .

Working in the physical basis of the charged leptons, Eqs. (52) and (53) determine the neutrino masses and mixings in the Standard Model with an arbitrary number of Higgs doublets. As we have shown, $`\delta M_L`$ is finite and gauge-invariant. Moreover, we have argued that all other one-loop contributions to $`_\nu `$ are of order $`Y^2m_D^2/m_R`$ (or $`g^2m_D^2/m_R`$), where $`Y`$ is a typical Yukawa coupling constant; consequently, $`\delta M_L`$ dominates the one-loop corrections.

The fact that the order of magnitude of the one-loop corrections in Eq. (53) is smaller than the tree-level mass term solely by a factor of order $`\left(16\pi ^2\right)^1\mathrm{ln}\left(m_R/m_0\right)`$, where $`m_0`$ is a typical neutral-scalar mass, can be exploited for generating the ratio $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{\mathrm{atm}}^2`$ radiatively, such that the order of magnitude of that ratio is precisely the one adequate for the large-mixing-angle MSW solution of the solar-neutrino problem . However, the question of whether this framework also allows to explain the specific features of lepton mixing, namely the small element $`U_{e3}`$ of the mixing matrix, and an atmospheric mixing angle close to $`45^{}`$, requires further investigation.


text:
# Global dynamics and stability limits for planetary systems around HD 12661, HD 38529, HD 37124 and HD 160691.

## 1 Introduction

The recent explosion in the number of newly detected extrasolar planets has brought the total number of such planets to about one hundred; it appears that the exact number changes too rapidly to be quoted. Some of the planets form multiple planetary systems around their parental stars. At the moment the total number of systems with two or more planetary-mass companions around main sequence stars has reached 11. Four planetary systems - HD 12661, HD 38529, HD 37124 and HD 160691 - have recently been elevated to the status of multiple, following the discovery of a second planetary companion in each of them. In most of these systems the strong dynamical interaction between planets makes planetary orbital parameters (see Table 1), found using standard two-body Keplerian fits, unreliable (Laughlin & Chambers 2001, 2002). There is also a great uncertainty in the determination of planetary masses. All those leave us a substantial available parameter space to be explored in order to exclude the initial conditions which lead to dynamically unstable configurations.

A classical method that allows one to distinguish between regular and chaotic dynamical states is the method of Lyapunov Characteristic Numbers (LCN). Let us note that chaotic in the Poincaré sense means that the dynamical behavior is not quasi-periodic, and does not necessarily mean that the system will disintegrate during any limited period of time. The estimation of LCN usually requires computations over long evolutionary time, sometimes much longer than the lifetime of the system studied.

In our previous papers (see for example, Goździewski et al. 2001a, 2001b) we showed that a new method developed by Cincotta & Simó (2000) and called MEGNO (the acronym of Mean Exponential Growth of Nearby Orbits), can be successfully applied to the studies of dynamical stability of extrasolar planetary systems. This method is based on the same ideas as LCN but converges about 100 times faster, and is more sensitive. For example its application to the Gliese 876 system (Goździewski et al. 2001b) clearly identified the exact location of the 2:1 mean motion resonance and its width. MEGNO helped to show (Goździewski & Maciejewski 2001) that while the orbital parameters of the HD 82943 system derived from the Keplerian fit lead to an unstable self-destructing system, there are some small changes in this fit which lead to stable configurations.

In this work we apply the technique to all four new planetary systems. We are especially interested in a comparative study of the global dynamics of these systems, as they are likely to represent different types of dynamical behavior. Taking into account the rather preliminary nature of all orbital fits, we hope that this paper will provide a useful guide for available stable orbital parameters (such as $`a`$, $`e`$ and $`\omega `$) for more sophisticated self-consistent fits (Laughlin & Chambers 2001, 2002; Marcy et al. 2002). It may also provide tighter constraints on the parameter space available to some planets due to unsufficient number of observations as in the case of the HD 37124c planet where any eccentricity $`e_c`$ between 0.3 and 0.8 fits the observational data within the velocity errors (Butler et al. 2002).

## 2 Results.

All our results obtained with MEGNO were confirmed by direct integration of the equations of motion using the Lie-series method (e.g. Hanslmeier & Dvorak, 1984) which has been already used for many numerical simulations of Solar System dynamics (Tsiganis et al. 2001) and dynamics of extrasolar planetary systems (Pilat-Lohinger & Dvorak 2002). This integration method is based on recurrence formulae and uses an automatic step-size control. The four systems under study were integrated over evolutionary times of up to $`10^7`$ yrs. In all our models we considered the Keplerian orbital elements from Table 1 in reference to the coordinate system centered at the central star. However we are aware that for multiple planet systems orbital elements derived from summed 2-body fits to a radial velocity data might be better interpreted as referring to Jacobi coordinates (see, for example, Rivera & Lissauer 2000). We will try this approach in our future studies. Our study revealed many interesting and rather unexpected dynamical features in all of these systems. It is remarkable that being discovered almost simultaneously they represent a real ‘dynamical zoo’ that we describe below. Let us note that because of the preliminary nature of the orbital parameters in all systems we discuss in this paper, conclusions about the exact dynamical status of a particular system may differ from results presented here as the determination of orbital parameters will continue to be improved. However, we expect that the method we use and the general conclusions about dynamical properties of the parameter space available for each system will have a longer term value.

### 2.1 HD 12661

In HD 12661 the two planets are rather close to each other. In fact, for this system the ratio of the periastron distance of the outer orbit to the apastron distance of the inner orbit is very close to the corresponding critical value for three-body dynamical stability derived by Eggleton & Kiseleva (1995) 111Note that in their Eq2 the second sign should be ”-” for less extreme mass ratios (triple stars). This places HD 12661 in a very interesting zone of critical parameters where small changes can lead to totally different dynamical regimes. Fig.1 shows the variations of the semi-major axis $`a_b`$ of the innermost planet in HD 12661 over 1 million years. In the lower panel, which represents the system with nominal initial orbital parameters (Table 1), $`a_b`$ displays a regular quasi-periodic behavior with no sign of unpredictable variations over time. The top panel shows totally different behavior of $`a_b`$ for the same planetary system with only a very small difference in initial value of $`a_c`$: $`a_c`$=2.58 AU instead of 2.56 AU. This $`\mathrm{\Delta }a_c=0.02`$ AU is well within the error of the $`a_c`$ determinaton (see, for example, Fischer et al. 2001). In this ‘new’ system the semi-major axis $`a_b`$ displays irregular and therefore unpredictable jumps from one ‘mean’ value to another. We call such an orbit chaotic. This phenomenon of unpredictibilty is well-known in dynamical systems since Poincaré, and it is a characteristic of chaos in such systems. These results are in a very good agreement with MEGNO tests: for the system with the nominal parameters of HD 12661 $`a_c`$=2.56 AU the MEGNO indicator (see Goździewski et al. 2001a,b) converges to the value 2 characteristic for regular systems, while for the system with $`a_c`$=2.58 AU $`<`$Y$`>`$ is $``$6, which indicates the presence of chaos in this system. Let us note we do not claim that such a system will disintegrate after any particular period of time. We just point out that its dynamical state over a long time is unpredictable. From the conventional point of view one would expect the system with larger $`a_c`$ to be more stable. However, this conventional wisdom fails for systems whose orbital parameters are located within regions of marginal stability in the phase space (see below).

Fig.2 shows the MEGNO stability map for initial $`a_b`$-$`a_c`$ of HD 12661 in a rather large neighbourhood of the nominal values of these orbital parameters. Filled circles mark stable regular systems where $`<`$Y$`>`$\[1.97.203\]. Small dots mark chaotic regions where $`<`$Y$`><`$ 1.95 or $`<`$Y$`>>`$ 2.05. Open circles with dots mark intermediate cases. One can see that the stability map itself looks rather chaotic, as it displays a number of chaotic initial orbital parameters (marked with open circles and with small dots) unsystematically located within fairly large regular zones (), marked with the filled circles. This particular structure of the phase space is typical for regions close to the border between stable and unstable motions (e.g. Dvorak et al. 1989). In such regions minor changes in each of the parameters may dramatically transfer the system from a regular into a chaotic dynamical state. This behavior is not observed in the motion of the planets in our Solar System (e.g. Murray & Holman 2001), but it is a well-known phenomenon for the dynamics of asteroids (Wisdom 1983). Kiseleva et al. (1994a, 1994b) found a similar phenomenon for simulated triple stars near the stability limit. The 9:2 mean motion resonance may be responsable for the rather wide strip of chaotic parameters.

### 2.2 HD 37124 and HD 38529

The orbits in both these systems are solidly stable with no chaotic features for a wide range of orbital parameters including high relative inclinations between the two planetary orbits. Let us note that all results presented in this paper were obtained for nearly coplanar systems unless it is indicated overwise. When applied to the planetary orbits of the HD 38529 system with different variations of its orbital parameters, the MEGNO indicator $`<`$Y$`>`$ almost always converges to 2, with a few exceptions for very large eccentricities $`e_c`$ of the outer planet (see Fig. 3c). Direct N-body integrations of the equations of motion performed for selected samples also showed a quasi-periodic behaviour for all orbital parameters considered. The same is true for HD 37124: this system becomes unstable when $`e_c`$ approaches 0.6 (see Fig. 3b). Butler et al. (2002) pointed out that with only 30 observations over more than 5 years the eccentricity of the outer planet $`e_c`$ in this system is poorly constrained and can be anything between 0.3 and 0.8. Our results are in very good agreement (Fig. 3c) with their critical for stability value $`e_c0.65`$ for co-planar configurations with $`\mathrm{sin}(i)`$=1.0. This is precisely the case when MEGNO can be used much more efficiently than standard numerical experiments to set tight dynamical constraints on the uncertain parameters of planets.

However, the orbital parameters of planets in both systems, especially the planetary eccentricities (which are the orbital parameter most sensitive to external perturbations), are not unaffected by the presence of the other planetary companion. For example, the eccentricity of the planet HD 37124b changes periodically between 0.1 and 0.43. Even in the HD 38529 system where the planets are not expected to have much influence on each other’s dynamics, the eccentricity of the very close to the star inner planet ($`P_b=14.3`$ days) fluctuates within the interval $`[0.292,0.368]`$ for a nearly coplanar model with nominal parameters. The massive outer planet with nearly substellar mass $`M_c^{min}=12.78M_J`$ is also subjected to small eccentricity fluctuations.

It can be interesting to consider a case where the two planetary planes in HD 38529 have high relative inclination. MEGNO tests showed that the orbits in this system remains stable for all relative inclinations, including orthogonal orbits with $`i_r=90^\mathrm{o}`$. Taking into account the sub-stellar minimal mass of planet C, completely independent formation scenarios for companions B and C are possible. Therefore a highly non-coplanar configuration would not be unlikely for this system. The innermost planet B is close enough to the star for tidal friction to become important when the ‘Kozai effect’ (Kozai 1962) induces large fluctuations on the orbital eccentricity ($`e_{max}=`$1 for $`i_r=90^\mathrm{o}`$). This may significantly affect the dynamics of this system (see, for example, Kiseleva et al. 1998; Eggleton et al. 1998; Eggleton & Kiseleva-Eggleton 2001a, 2001b; Wu & Murray, in preparation). We are going to investigate these effects in a future paper.

### 2.3 HD 160991

Our simulations show that for HD 160691 the initial orbital parameters in Table 1 lead to a highly unstable self-destructing system. However, the orbital parameters of both planetary orbits are rather speculative at the moment: even the existence of the second planetary companion is not yet fully confirmed (Jones et al. 2002). Therefore we looked for stable configurations among wide ranges of orbital parameters. We found a few stable configurations located not very far in the parameter space from the currently defined parameters of the system (see Fig. 3d). It looks as if a high eccentricity of the outermost planet $`e_c>0.7`$ is an important stabilizing factor in this system. Using MEGNO we were able to identify a few stability zones in a parameter space which includes together with the parameters listed above also parameters not determined from observations such as the relative inclination $`i_r`$ between the two planetary orbits. Most of these stable configurations are associated with the 2:1 mean motion resonance. We are going to discuss in detail the possible combinations of parameters resulting in stable orbits in HD 160991 in our forthcoming paper (Bois et al., in preparation).

## 3 Comparative stability analysis and conclusions

Fig. 3 presents a comparative visualisation of our MEGNO stability analysis for all four new planetary systems. For this visual presentation we choose to consider orbital stability as a function of both orbital eccentricities simply because the values of $`e_b`$ and $`e_c`$ can only be changed between 0 and 1, and therefore do not require any scaling (as for example would semi-major axes) for different systems. All other orbital parameters are the nominal ones from Table 1. It is easy to see the differences in dynamical status of the four systems and to identify the ranges of eccentricities which allow stable planetary orbits in each of them. We also produced simular maps for other pairs of orbital parameters. Such ranges of stable parameters can be very useful for the improvement of preliminary orbits in the new planetary systems discussed in this paper. In our future papers we plan to present detailed dynamical analyses of each system taking into account angular orbital parameters not constrained by observational data ($`i_r`$x, $`\mathrm{\Omega }`$), as well as $`\mathrm{sin}i`$ and the resulting different planetary masses.

This work was supported by the John Templeton Foundation (Agreement ID #398-COS272). RD thanks the Austrian FWF P14375-TPH for a support in travel expenses. LK-E thanks the University Bordeaux 1 for a research fellowship. LK-E and RD thank the Bordeaux Observatory for its hospitality. All authors thank the referee Greg Laughlin for very useful comments and suggestions.


text:
# Galaxy Clusters with Chandra11footnote 1Contribution to XIII Rencontres de Blois 2001, ed. L. M. Celnikier

## Galaxy Clusters with Chandra111Contribution to XIII Rencontres de Blois 2001, ed. L. M. Celnikier

#### 1) Smithsonian Astrophysical Observatory, Cambridge, MA, USA
2) MPI fur Astrophysik, Garching, Germany
3) Space Research Institute, Moscow, Russia

### 0.1 Cluster Formation at High Angular Resolution

For many years clusters were thought to be dynamically relaxed systems evolving slowly after an initial, short-lived episode of violent relaxation. However, while the dynamical timescale for the richest cluster cores is comfortably less than the Hubble time, other less dense clusters have dynamical timescales comparable to or longer than the age of the Universe . X-ray images, starting with Einstein and continuing with ROSAT and ASCA, and now with Chandra and XMM-Newton, provide a powerful technique to “map” the structure in the gravitational potential of massive systems containing hot gaseous atmospheres.

The X-ray observations supported the now prevalent idea that structure in the Universe has grown through gravitational amplification of small scale instabilities or hierarchical clustering. At one extreme, some clusters grow, in their final phase, through mergers of nearly equal mass components. Such mergers can be spectacular events involving kinetic energies as large as $`10^{64}`$ ergs, the most energetic events since the Big Bang. More common are smaller mergers and accretion of material from large scale filaments. The ROSAT image of A85 shows the relationship between large scale structure and cluster merging where small groups are detected infalling along a filament into the main cluster . Chandra and XMM with high angular resolution and high throughput provide new insights into cluster formation and evolution.

##### Cluster Cold Fronts

Prior to the launch of Chandra, sharp gas density discontinuities had been observed in the ROSAT images of A2142 and A3667 . Since both clusters exhibited characteristics of major mergers, these features were expected to be shock fronts. However, the first Chandra observations showed that these were not shocks, but a new kind of structure – cold fronts . Their study has provided new and detailed insights into the physics of the intracluster medium (ICM) .

Fig. 1a shows a striking example of a cold front as NGC1404 and its gaseous corona approaches the cluster center (to the northwest) . The image clearly shows the sharp edge of the surface brightness discontinuity, shaped by the ram pressure of the cluster gas. The temperature map (Fig. 1b) confirms that the infalling cloud is cold compared to the hotter Fornax ICM.

##### Physics of Cold Fronts

The best studied cold front is that in A3667 ($`z=0.055`$). Although expected to exhibit a shock front based on its ROSAT image , the sharp feature observed in the Chandra observation is the boundary of a dense cold cloud, a merger remnant . The surface brightness profile, converted to gas density, and the gas temperature distribution yield the gas pressure on both sides of the cold front. The difference between the two pressures is a measure of the ram pressure of the ICM on the moving cold front. Hence, the precise measurement of the gas parameters yields the cloud velocity. The factor of two difference in pressures between the free streaming region and the region immediately inside the cold front yields a Mach number of the cloud of $`1\pm 0.2`$ ($`1430\pm 290`$ km s-1.

In addition to the edge, a weak shock is detected. The distance between the cold front and the weak shock ($`350`$ kpc) and the observed gas density jump at the shock (a factor of 1.1-1.2) yield the shock’s propagation velocity, $`1600`$ km s-1, which is consistent with that derived independently from the pressure jump across the cold front .

The A3667 observation provides important information on the efficiency of transport processes in clusters. As the surface brightness profile shows (see ), the density “edge” is very sharp. Quantitatively, Vikhlinin et al. found that the width of the front was less than $`3.5^{\prime \prime }`$ (5 kpc). This sharp edge requires that transport processes across the edge be suppressed, presumably by magnetic fields. Without such suppression, the edge should be broader since the electron Coulomb mean free path is about 13 kpc, several times the width of the cold front . Furthermore, Vikhlinin et al. observed that the cold front appears sharp only over a sector of about $`\pm 30^{}`$ centered on the direction of motion, while at larger angles, the sharp boundary disappears . The disappearance can be explained by the onset of Kelvin-Helmholtz instabilities, as the ambient ICM gas flows past the moving cold front. To explain the observed extent of the sharp boundary, the instability must be partially suppressed, e.g., by a magnetic field parallel to the boundary with a strength of $`716\mu `$G with a corresponding pressure of only 10-20% of the thermal pressure .

##### ZW3146 – Multiple “Fronts”

ZW3146 is a moderately distant ($`z=0.2906`$; 5.74 kpc per arcsec) cluster with a remarkably high mass deposition rate that was estimated to exceed 1000 M yr-1 . The Chandra image further demonstrates the remarkable nature of this cluster – on scales from $`3^{\prime \prime }`$ to $`30^{\prime \prime }`$ ($`20`$ kpc to 170 kpc), three separate X-ray surface brightness edges are detected (see Fig. 2 and Forman et al. ). At the smallest radii, two edges are seen to the northwest and north of the center (see Fig.2a). The first, at a radius of $`3^{\prime \prime }`$ (17 kpc), spans an angle of nearly $`180^{}`$ with a surface brightness drop of almost a factor of 2. The second edge, at a radius of $`8^{\prime \prime }`$ (45 kpc) spans only $`90^{}`$ but has a surface brightness drop of almost a factor of 4. The third edge (see Fig. 2b) lies to the southeast, about $`35^{\prime \prime }`$ (200 kpc) from the cluster center, has a decrease of about a factor of 2, and, as with the first edge, extends over an angle of almost $`180^{}`$.

The variety of morphologies and scales exhibited by these sharp edges or cold fronts is quite remarkable. Possibly the edges may arise from moving cold gas clouds that are the remnants of merger activity as observed in A2142 and A3667 or as oscillations (or “sloshing”) of the cool gas at the center of the cluster potential as observed in A1795 . The extremely regular morphology of ZW3146 on large linear scales seems to exclude a recent merger and, hence, “sloshing” of the gas seems the more likely explanation for the observed edges. High resolution, large scale structure simulations show that dense halos, formed at very early epochs, survive cluster collapse . While most of the dark matter halos, having galaxy size masses, are associated with the sites of galaxy formation, the larger mass halos also may survive (without their gas) or may have fallen into the cluster only recently. Hence, we might expect to find a range of halo mass distributions moving within the cluster potential. We speculate that, as these halos move, the varying gravitational potential could accelerate the cool dense gas that has accumulated in the cluster core and could produce the “sloshing” needed to give rise to the multiple surface brightness edges observed in some clusters.

##### CL0657 – A Prototypical Cluster Shock Front

CL0657 ($`z=0.296`$) was discovered by Tucker et al. as part of a search for “failed” clusters, clusters that were X-ray bright but had few, if any, optical galaxies . From ASCA observations, this cluster was found to have a remarkably hot gas temperature of about 17 keV, making it the hottest cluster known .

The Chandra image of CL0657 shows the classic properties of a supersonic merger (see Markevitch et al. for a detailed discussion of this cluster ). We see a dense (cold) core moving to the west after having traversed, and disrupted, the core of the main cluster. Leading the cold, dense core is a density discontinuity which appears as a shock front (Mach cone). The spectral data show that the gas to the east (trailing the shock) has been heated by the passage of the shock. The detailed gas density parameters confirm that the “bullet” is moving to the west with a velocity of 3000–4000 km sec-1, approximately 2-3 times the sound speed of the ambient gas. CL0657 is the first clear example of a relatively strong shock arising from cluster mergers.

### 0.2 The Radio—X-ray Connection – or Bubbles, Bubbles Everywhere

Prior to the launch of Chandra, ROSAT observations of NGC1275 and M87 provided hints of complex interactions between radio emitting plasmas ejected from AGN within the nuclei of dominant, central cluster galaxies . With the launch of Chandra, the interaction between relativistic plasma, with and without detected radio emission, and the hot intracluster medium (ICM) has been observed in many systems and is now a major area of investigation .

##### M87 and the Evolution of Buoyant Plasma Bubbles

The 327 MHz high resolution, high dynamic range radio map (see Fig. 4) of M87 shows a well-defined torus-like eastern bubble and a less well-defined western bubble, both of which are connected to the central emission by a column, and two very faint almost circular emission regions northeast and southwest of the center . The correlation between X–ray and radio emitting features has been remarked by several authors .

Motivated by the similarity in appearance between M87 and hot bubbles rising in a gaseous atmosphere, Churazov et al. developed a simple model of the M87 bubbles which is generally applicable to the many bubble-like systems seen in the Chandra observations . Fig. 5 shows an initially buoyant bubble transform into a torus as it rises through a galaxy or cluster atmosphere.

By entraining cool gas as it rises, the bubble produces a characteristic “mushroom” appearance, similar to an atmospheric nuclear explosion. This qualitatively explains the correlation of the radio and X–ray emitting plasmas and naturally accounts for the thermal nature of the X-ray emission associated with the rising torus . Finally, in the last evolutionary phase, the bubble reaches a height at which the ambient gas density equals that of the bubble. The bubble then expands to form a thin layer (a “pancake”). The large low surface brightness features in the M87 radio map could be just such pancakes. In the simulations performed by Churazov et al. the buoyant bubbles behaved as expected and did produce the features observed in both X-rays and radio for M87. Although the exact form of the rising bubbles was sensitive to initial conditions, the toroidal structures were a common feature. Ambient gas was uplifted in the cluster atmosphere reducing the effects of cooling flows and producing the “stem” of the mushroom that is brighter than the surrounding regions . Note that the XMM-Newton observation shows that the bright X-ray columns are cool, exactly as expected in the buoyant bubble scenario . A sketch of the overall source structure of M87, based on the evolution of buoyant bubbles, is shown in Fig. 4 .

##### NGC507 - the central galaxy in a group

NGC507 is the central galaxy in a nearby ($`z=0.016`$) group that has been studied extensively in X-rays . The galaxy is the site of a weak B2 radio source (luminosity $`10^{37}`$ ergs s-1. The Chandra X-ray image, shown in Fig. 6a, covers only the central, high surface brightness emission of the group around NGC507. The 0.5-2.0 keV surface brightness distribution shows sharp edges to the southwest, southeast and north, reminiscent of those in the clusters A2142 and A3667. In addition to the edges, there are two X-ray peaks. The first, to the east, coincides with the nucleus of NGC507. A second peak, $`1^{}`$ to the west has no optical counterpart. However, comparing the X-ray and the radio map (Fig. 6a, c) shows that the western radio lobe lies precisely in the surface brightness trough between the nucleus and the peak to the west. Thus, it seems likely that the radio lobe, probably a buoyant bubble, has displaced X-ray emitting gas generating a trough in the X-ray surface brightness distribution.

The origin of the peculiar sharp surface brightness discontinuities around NGC507 is unclear. The bright emission is well fit by a thermal model with gas temperatures near 1 keV, consistent with the mean ASCA temperature of $`1.10\pm 0.05`$ keV . Detailed spectroscopic analyses were complicated by unexpected high background at high energies. A second Chandra observation is planned to allow detailed spectroscopy. The emission from the central region is resolved and hence the contribution from a central AGN is relatively small (see Forman et al.  for additional discussion of NGC507). Perhaps the X-ray surface brightness features arise either from motion of NGC507 and its dark halo within the larger group potential as suggested for the multiple edges in clusters .

### 0.3 Conclusions

We did not expect the rich variety of new structures seen in the Chandra high angular resolution observations of clusters and early type galaxies. Instead of confirming our prejudices, Chandra has brought us a wealth of new information on the interaction of radio sources with the hot gas in both galaxy and cluster atmospheres. We see “edges” in many systems with hot and cold gas in close proximity and have been able to extract important new parameters of the ICM. We have only barely begun to digest the import of the Chandra cluster and galaxy observations. We can only expect the unexpected as Chandra observations continue and as our understanding of how best to use this new observatory matures.

We acknowledge support from NASA contract NAS8 39073, NASA grants NAG5-3065 and NAG5-6749 and the Smithsonian Institution.

Note Added on ZW3146

One of the puzzles relating to the X-ray observation of ZW3146 was the apparent misalignment by a few arcseconds between the X-ray brightness peak and the optical center of the central, dominant cD galaxy. The solution to this apparent misalignment was resolved by the FIRST radio image which shows two radio sources with flux densities of $`2`$ mJy separated by $`14^{\prime \prime }`$ (see Fig. 7a). As Fig. 7b shows, the northern component of the pair of radio sources is precisely aligned with the optical cD galaxy at the cluster center. Fig. 7c shows that the radio source, and therefore the center of the cD galaxy, lies at the approximate center of irregular X-ray features. The bright X-ray structure northwest of the cD is elongated perpendicular to the direction to the galaxy center. A second X-ray bright feature lies southeast of the radio peak. We suggest that the radio emitting plasma formed a bubble that produced a trough in the X-ray surface brightness at the galaxy center. The bubble is surrounded by an irregular shell of X-ray emission. Thus, the brightest X-ray structures do not lie at the galaxy center, but instead form a partial shell surrounding the radio emitting plasma . Such structures are similar to the X-ray surface brightness enhancements seen around radio plasma bubbles in nearby galaxies e.g., M84 .