text:
# A Module for Radiation Hydrodynamic Calculations With ZEUS-2D Using Flux-Limited Diffusion

## 1 INTRODUCTION

In many astrophysical systems, radiation is the dominant energy transport mechanism. In some circumstances, including winds from Wolf-Rayet stars (Lucy & Abbott 1993) and asymptotic giant branch stars (Habing 1996), interiors of massive stars (Kippenhahn & Weigert 1990), supernova blast waves (Arnett et al. 1989), and compact object accretion flows (Klein et al. 1996) and disks (Shakura & Sunyaev 1973), radiation comprises a substantial fraction of the total energy density, momentum density, or pressure. The dynamics of such flows must be modeled using a radiation hydrodynamical (RHD) approach, in which energy and momentum conservation laws for the radiation field are solved along with those for the gas. The outcomes of these calculations depend on the exchange of internal energy and momentum between material and radiation. The length and time scales associated with the exchanges are determined by atomic processes, and may be orders of magnitude shorter or faster than dynamical scales, so that implicit time differencing is usually necessary. An accurate description of the angular dependence of the radiation field may be important when parts of the gas are optically thin, localized sources of radiation are present, or shadows are cast within the flow. It is possible to compute the angular dependence by solving the transfer equation along rays through the points of interest. However, the number of rays needed to adequately specify the angular dependence at every point is often large. Such a method has been implemented in two spatial dimensions in the ZEUS-2D magnetohydrodynamics (MHD) code (Stone & Norman 1992b) by Stone, Mihalas, & Norman (1992), but is complex and computationally intensive. Despite the promise of this full-transport method, it must be considered still under development. In this paper we describe a method for calculations using ZEUS-2D, in which the angular dependence of the radiation field is assumed to be given by the flux-limited diffusion (FLD) approximation. This method is simple, robust, and relatively cheap in computer time. We assume the gas is in local thermodynamic equilibrium (LTE) at a temperature which need not correspond to that of the radiation field. Frequency dependence of the opacities is neglected, though it may be included in a straightforward fashion. The RHD equations solved are discussed in §2, the FLD approximation in §3, and the numerical algorithm in §4. Results of test calculations are presented in §5, and the advantages and limitations of the method are summarized in §6.

## 2 EQUATIONS OF RADIATION HYDRODYNAMICS

Much as the gas dynamical equations are derived by taking velocity moments of the material particle kinetic equation, or Boltzmann equation, the conservation laws for the radiation field are generated by taking angular moments of the photon kinetic equation, or radiation transport equation. In a frame comoving with the radiating fluid, assuming LTE, and to order unity in $`v/c`$, the coupled equations of RHD are (Mihalas & Mihalas 1984)

$$\frac{D\rho }{Dt}+\rho 𝐯=0,$$
(1)
$$\rho \frac{D𝐯}{Dt}=p+\frac{1}{c}\chi _F𝐅,$$
(2)
$$\rho \frac{D}{Dt}\left(\frac{E}{\rho }\right)=𝐅𝐯:𝖯+4\pi \kappa _PBc\kappa _EE,$$
(3)
$$\rho \frac{D}{Dt}\left(\frac{e}{\rho }\right)=p𝐯4\pi \kappa _PB+c\kappa _EE,$$
(4)
and

$$\frac{\rho }{c^2}\frac{D}{Dt}\left(\frac{𝐅}{\rho }\right)=𝖯\frac{1}{c}\chi _F𝐅.$$
(5)
Here the convective derivative $`D/Dt/t+𝐯`$. The dependent quantities $`\rho `$, $`e`$, $`𝐯`$, and $`p`$ are the material mass density, energy density, velocity, and scalar isotropic pressure respectively, while $`E`$, $`𝐅`$, and $`𝖯`$ are the total frequency-integrated radiation energy density, momentum density or flux, and pressure tensor, respectively. The latter are the zeroth, first, and second angular moments of the radiation specific intensity. The specific intensity, $`I(𝐱,𝛀,\nu ,t)`$, in general varies with spatial position $`𝐱`$, viewing direction $`𝛀`$, frequency $`\nu `$, and time $`t`$. Complete description of $`I`$ in a three-dimensional space therefore requires seven scalar coordinates. The angular moments of $`I`$ appearing in the RHD equations above may be written

$$E(𝐱,t)=\frac{1}{c}_0^{\mathrm{}}𝑑\nu 𝑑\mathrm{\Omega }I(𝐱,𝛀,\nu ,t),$$
(6)
$$𝐅(𝐱,t)=_0^{\mathrm{}}𝑑\nu 𝑑\mathrm{\Omega }I(𝐱,𝛀,\nu ,t)𝐧,$$
(7)
and

$$𝖯(𝐱,t)=\frac{1}{c}_0^{\mathrm{}}𝑑\nu 𝑑\mathrm{\Omega }I(𝐱,𝛀,\nu ,t)\mathrm{𝐧𝐧}.$$
(8)
The assumption of LTE allows the source function specifying the rate of emission of radiation from the gas in equations (3) and (4) to be written as the Planck function $`B`$. Equations (2) through (5) have been integrated over frequency, leading to the flux mean total opacity $`\chi _F`$, and the Planck mean and energy mean absorption opacities, $`\kappa _P`$ and $`\kappa _E`$. In the present paper the opacities are assumed to be independent of frequency, so that $`\kappa _P=\kappa _E`$ and the subscripts may be omitted. The total opacity $`\chi `$ is the sum of components due to absorption, $`\kappa `$, and scattering, $`\sigma `$, all having dimensions of inverse length. In a pure hydrogen gas above $`10^4`$ K, simple forms for $`\kappa `$ and $`\sigma `$ are Kramers’ Law and the Thomson scattering coefficient. Equations (1) to (5) correctly describe the flow only in the comoving frame, where material properties are isotropic and the form of the material-radiation interaction terms is greatly simplified. In writing the equations, effects of self-gravity and magnetic fields have been ignored. The treatment of these effects in ZEUS-2D is independent of the algorithms for RHD described here.

The equations of RHD may be closed by the addition of constitutive relations for the Planck function and opacities, an equation of state specifying the gas pressure, and an assumption about the relationship between the angular moments of the radiation field. Since $`𝖯=\frac{1}{3}E`$ when the field is isotropic, one might choose to assume this relation holds everywhere. This is referred to as the Eddington approximation, and is analogous to assuming isotropy of the material particle distribution function in deriving the gas equation of state. The Eddington approximation implies that in steady-state, equation (5) becomes

$$𝐅=\frac{c}{3\chi }E.$$
(9)
This expression gives the correct flux in optically thick regions, where $`\chi `$ is large, and tends to infinity in optically thin regions where $`\chi 0`$. However the rate at which radiation transports energy is finite and, from equations (6) and (7), is limited to $`|𝐅|cE`$. This causality problem with equation (9) occurs because in optically thin regions photons travel freely, their mean free paths may exceed characteristic lengths in the system, and the radiation field may be anisotropic. In this paper we implement an extension of the Eddington approximation in which causality is preserved by assuming a particular form for the angular dependence of the radiation field in optically-thin regions. This extension is described in the next section.

## 3 THE FLUX-LIMITED DIFFUSION APPROXIMATION

FLD techniques were first used in astrophysics to solve the radiation transport equation, in calculations of accretion onto a neutron star (Alme & Wilson 1974). Subsequently, Levermore & Pomraning (1981; hereafter LP) generalized the method to approximately handle transport phenomena while preserving causality in regions where spatial variations occur over distances smaller than a mean free path. Pomraning (1983) included relativistic terms of order $`v/c`$ in a moving fluid. FLD methods have since been used in pure scattering media (Melia & Zylstra 1991), and for fully general-relativistic calculations (Anile & Romano 1992). The key ingredient in all cases is the assumption that the specific intensity is a slowly-varying function of space and time. In one space dimension, this assumption is valid in both the optically-thick diffusion limit and the optically-thin free-streaming limit. One hopes it holds approximately in the intermediate regime and in multi-dimensions, and this may be tested by comparison with results of full-transport calculations. Given slow spatial and time variation, analytic relations may be obtained between the angular moments of the radiation field. LP showed that the radiation flux may be written in the form of Fick’s law of diffusion, as

$$𝐅=DE,$$
(10)
with a diffusion coefficient $`D`$ given by

$$D=\frac{c\lambda }{\chi }.$$
(11)
The dimensionless function $`\lambda =\lambda (E)`$ is called the flux limiter. In this framework, the radiation pressure tensor may be expressed in terms of the radiation energy density via

$$𝖯=𝖿E,$$
(12)
where the components of the Eddington tensor $`𝖿`$ are given by

$$𝖿=\frac{1}{2}(1f)𝖨+\frac{1}{2}(3f1)\mathrm{𝐧𝐧},$$
(13)
$`𝐧=E/|E|`$ is the unit vector in the direction of the radiation energy density gradient, and the dimensionless scalar function $`f=f(E)`$ is called the Eddington factor. The flux limiter $`\lambda `$ and Eddington factor $`f`$ are related through implicit constraints between the moments $`𝐅`$ and $`𝖯`$, so that

$$f=\lambda +\lambda ^2R^2,$$
(14)
where $`R`$ is the dimensionless quantity $`R=|E|/(\chi E)`$.

Equations (10) through (14) close the equations of RHD, eliminating the need to solve the radiation momentum equation (5), and greatly simplifying the computation. Since the angular distribution of the radiation field is not explicitly computed, no formal solution of the transfer equation is needed. An important choice which must be made to achieve this simplification is the expression to be used for the flux limiter $`\lambda `$. Many expressions are possible which preserve causality, and are consistent with the assumption of smoothness in the radiation field. These are distinguished by different assumptions regarding the details of the angular dependence of the specific intensity. In one example, LP chose to constrain the normalized specific intensity $`\psi =I/(cE)`$ by

$$\frac{1}{c}\frac{\psi }{t}+𝛀\psi =0,$$
(15)
resulting in the flux limiter

$$\lambda (R)=\frac{2+R}{6+3R+R^2}.$$
(16)
Physically, the derivation of this result follows an analysis similar to the Chapman-Enskog approach for developing solutions to the Boltzmann equation in the kinetic theory of gases. In a second example of a choice of constraint on the anisotropy of the radiation field, Minerbo (1978) assumed a piecewise linear variation of the specific intensity with angle, and found the flux limiter

$$\lambda (R)=\{\begin{array}{cc}2/(3+\sqrt{9+12R^2})\hfill & \text{if }0R3/2\hfill \\ (1+R+\sqrt{1+2R})^1\hfill & \text{if }3/2<R<\mathrm{}\hfill \end{array}$$
(17)
In the optically thin limit $`R\mathrm{}`$, both the LP and Minerbo flux limiters give to first order in $`R^1`$

$$\underset{R\mathrm{}}{lim}\lambda (R)=\frac{1}{R},$$
(18)
so the magnitude of the flux approaches $`|𝐅|=c|E|/(\chi R)=cE`$, which obeys the causality constraint. In the optically thick or diffusion limit $`R0`$, both examples give to first order in $`R`$

$$\underset{R0}{lim}\lambda (R)=\frac{1}{3},$$
(19)
so the flux takes the value given by equation (9). For intermediate values of $`R`$, these two forms for the flux limiter differ substantially, as shown in figure 1. Thus a major uncertainty in FLD calculations is the appropriate choice for the form of the flux limiter. Some alternatives are described by Levermore (1984). The LP form has been used in studies of accretion onto black holes (Eggum, Coroniti, & Katz 1988) and protostellar collapse (Bodenheimer et al. 1990). Dynamics of radiation-dominated accretion disks computed using the LP and Minerbo forms were compared by Turner & Stone (2001). Kley (1989) derived a flux limiter appropriate for accretion disk boundary layers, using exact solutions of the transfer equation for spherically-symmetric stellar atmospheres. This limiter is quite similar to the Minerbo form over the whole range of $`R`$. A flux limiter appropriate for a scattering medium was derived and used in accretion disk corona models by Melia & Zylstra (1991). The LP and Minerbo flux limiters have been widely used, and bracket the range of values of most other flux limiters. Both are available in the ZEUS-2D implementation of FLD described here. Since the appropriate form of the flux limiter is not always clear beforehand, in some applications new limiters may have to be derived in order to achieve accurate results. Comparisons between FLD and full transport results remain important.

In summary, applying the FLD approximation involves computing the flux limiter from the value of $`R`$ in each zone, and then using equations (10) to (14) to obtain the radiation flux and radiation pressure from the radiation energy density.

## 4 NUMERICAL METHOD

As with the MHD equations in ZEUS-2D, the dynamical equations for the radiation are operator-split into source and transport terms. In the source step, the radiation and material energy densities are updated using finite difference approximations to

$$\frac{E}{t}=𝐅$$
(20)
and separately,

$$\frac{E}{t}=𝐯:𝖯+4\pi \kappa Bc\kappa E$$
(21)
and

$$\frac{e}{t}=p𝐯4\pi \kappa B+c\kappa E.$$
(22)
In the transport step, an integral formulation is used to generate a conservative differencing scheme for the advection terms,

$$\frac{d}{dt}_VE𝑑V=_{dV}E(𝐯𝐯_𝐠)\mathrm{𝐝𝐒},$$
(23)
where $`𝐯_𝐠`$ is an arbitrary coordinate velocity which allows grid motion. Equations (20) to (23) are evolved in the radiation module, as is the source term due to radiative acceleration in the gas momentum equation (2). The rest of the set of coupled equations is solved using the existing hydrodynamic algorithms. Centering of the radiation variables on the hydrodynamic grid is described in §4.1, update of the radiation flux term in §4.2 and the remaining source terms in §4.3, treatment of the compression terms in the case of strictly scattering opacity in §4.4, and the transport step in §4.5. Stability of the method and the choice of timestep are discussed in §4.6.

### 4.1 Centering of Variables

The dependent variables for the radiation are discretized onto the computational mesh as in ZEUS-2D (Stone & Norman 1992a, hereafter SN). Scalars, and tensors of even rank, are zone-centered, while tensors of odd rank are face-centered. The mesh is labeled by coordinate vectors $`x_1`$ and $`x_2`$, with the 3-direction taken to be orthogonal to the computational plane. The mesh whose grid lines mark zone edges is labeled by coordinates $`x1a_i`$ and $`x2a_j`$, while the mesh whose lines intersect at zone centers is labeled by coordinates $`x1b_i`$ and $`x2b_j`$, as shown in figure 2. On this mesh, the radiation variables are discretized according to

$`E(x_1,x_2)`$ $``$ $`E(x1b_i,x2b_j)=E_{i,j}`$
$`F_1(x_1,x_2)`$ $``$ $`F1(x1a_i,x2b_j)=F1_{i,j}`$
$`F_2(x_1,x_2)`$ $``$ $`F2(x1b_i,x2a_j)=F2_{i,j}`$
$`P_{11}(x_1,x_2)`$ $``$ $`P11(x1b_i,x2b_j)=P11_{i,j}`$
$`P_{22}(x_1,x_2)`$ $``$ $`P22(x1b_i,x2b_j)=P22_{i,j}`$
$`P_{33}(x_1,x_2)`$ $``$ $`P33(x1b_i,x2b_j)=P33_{i,j}`$
$`P_{12}(x_1,x_2)`$ $``$ $`P12(x1b_i,x2b_j)=P12_{i,j}`$
The components of the Eddington tensor $`𝖿`$ are centered identically to the components of $`𝖯`$. The opacity, absorption and scattering coefficients, and Planck function are all zone-centered. The radiation diffusion coefficient $`D`$ is computed on zone 1-faces ($`D1_{i,j}`$) and 2-faces ($`D2_{i,j}`$) along with the components of the flux.

### 4.2 Radiation Flux Divergence Term

Under the FLD approximation, the flux divergence term in the radiation energy equation (20) is a function of the gradient in radiation energy density. When implicitly differenced, it therefore depends on the time-advanced values of $`E`$ in adjacent grid zones, and may be updated by a matrix inversion. To simplify the form of the matrix, we operator-split this from the remaining source terms, and integrate

$$\frac{E}{t}=(DE)$$
(24)
with $`D`$ given by equation (11). The solution method is described here for a uniformly-spaced square Cartesian grid. Extension to the curvilinear coordinates used in ZEUS-2D is straightforward. The implicit difference equations to be solved are of the form

$`{\displaystyle \frac{(\mathrm{\Delta }x)^2}{\mathrm{\Delta }t}}(E_{i,j}^{n+1}E_{i,j}^n)`$ $`=`$ $`D1_{i+1,j}^n(E_{i+1,j}^{n+1}E_{i,j}^{n+1})D1_{i,j}^n(E_{i,j}^{n+1}E_{i1,j}^{n+1})+`$ (25)
$`D2_{i,j+1}^n(E_{i,j+1}^{n+1}E_{i,j}^{n+1})D2_{i,j}^n(E_{i,j}^{n+1}E_{i,j1}^{n+1}),`$
where $`\mathrm{\Delta }x`$ is the grid spacing, and the spatial indices span ranges $`i=1\mathrm{}N`$ and $`j=1\mathrm{}M`$. Radiation energy densities are to be determined at timestep $`n+1`$. The difference equations are kept linear by using diffusion coefficients from step $`n`$. This is a good approximation when $`|E|/(\chi E)`$ changes little each timestep, but may be less appropriate when the radiation field varies rapidly with time in regions of intermediate or low optical depth. In section 5.5 it is shown that a propagating radiation front can be more accurately evolved in this situation when the diffusion coefficients are recalculated every diffusion substep.

Considering $`E_{i,j}^{n+1}`$ as a single vector of length $`N\times M`$, the updated values may be determined by solving an equation with a matrix of $`(N\times M)^2`$ elements, using a sparse-matrix technique such as incomplete Cholesky – conjugate gradient (Press et al. 1992). These methods have proven useful for particular problems with fixed boundary conditions, but the form of the matrix and the sparse-matrix method required vary with the boundary conditions. Solutions of equation (25) may also be obtained by combining simple solvers with multi-grid acceleration techniques (Hackbusch 1985). However, here we choose an alternating-direction-implicit (ADI) method on a single grid. This may be less efficient than the best sparse-matrix and multi-grid algorithms, though the number of operations required is proportional to $`N\times M`$. In ADI, successive approximations to $`E^{n+1}`$ are computed by advancing towards a $`w`$-stationary state the equation

$$\frac{E_{i,j}^{n+1}}{w}=\frac{E_{i,j}^{n+1}E_{i,j}^n}{\mathrm{\Delta }t}+\left[(DE)\right]_{i,j}^{n+1}.$$
(26)
The new variable $`w`$ may be thought of as the pseudo-time. Each $`w`$-step is split into two parts. In the first, the update is $`w`$-implicit along the 1-direction, and explicit along the 2-direction. In the second, the differencing schemes for the two axes are exchanged. Labeling the approximate value of $`E_{i,j}^{n+1}`$ at pseudo-timestep $`m`$ by $`E_{i,j}^m`$, the difference equations solved are

$`E_{i,j}^{m+\frac{1}{2}}E_{i,j}^m`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }w}{2\mathrm{\Delta }t}}(E_{i,j}^nE_{i,j}^{m+\frac{1}{2}})+{\displaystyle \frac{\mathrm{\Delta }w}{2(\mathrm{\Delta }x)^2}}\times `$ (27)
$`\{D1_{i+1,j}^n(E_{i+1,j}^{m+\frac{1}{2}}E_{i,j}^{m+\frac{1}{2}})D1_{i,j}^n(E_{i,j}^{m+\frac{1}{2}}E_{i1,j}^{m+\frac{1}{2}})`$
$`+D2_{i,j+1}^n(E_{i,j+1}^mE_{i,j}^m)D2_{i,j}^n(E_{i,j}^mE_{i,j1}^m)\}`$
on sweeps which are implicit along the 1-direction, and

$`E_{i,j}^{m+1}E_{i,j}^{m+\frac{1}{2}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }w}{2\mathrm{\Delta }t}}(E_{i,j}^nE_{i,j}^{m+1})+{\displaystyle \frac{\mathrm{\Delta }w}{2(\mathrm{\Delta }x)^2}}\times `$ (28)
$`\{D1_{i+1,j}^n(E_{i+1,j}^{m+\frac{1}{2}}E_{i,j}^{m+\frac{1}{2}})D1_{i,j}^n(E_{i,j}^{m+\frac{1}{2}}E_{i1,j}^{m+\frac{1}{2}})`$
$`+D2_{i,j+1}^n(E_{i,j+1}^{m+1}E_{i,j}^{m+1})D2_{i,j}^n(E_{i,j}^{m+1}E_{i,j1}^{m+1})\}`$
on sweeps which are implicit along the 2-direction. Each sweep involves solving a tridiagonal matrix equation. For example, for boundary condition $`E=0`$, the matrix equation solved on the $`j`$th sweep implicit along the 1-direction is

$`\left[\begin{array}{ccccc}\beta hD1_{1,j}& hD1_{2,j}& & & \\ hD1_{2,j}& \beta & hD1_{3,j}& & \\ & hD1_{3,j}& \beta & hD1_{4,j}& \\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ & & hD1_{N1,j}& \beta & hD1_{N,j}\\ & & & hD1_{N,j}& \beta hD1_{N+1,j}\end{array}\right]\left[\begin{array}{c}E_{1,j}\\ \mathrm{}\\ E_{N,j}\end{array}\right]^{m+\frac{1}{2}}=𝐛,`$ (38)
where $`h=\frac{\mathrm{\Delta }w}{2(\mathrm{\Delta }x)^2}`$, $`\beta =1+\frac{\mathrm{\Delta }w}{2\mathrm{\Delta }t}+h(D1_{i+1,j}+D1_{i,j})`$, and the elements of vector $`𝐛`$ are

$$b_i=[1h(D2_{i,j+1}+D2_{i,j})]E_{i,j}^m+hD2_{i,j+1}E_{i,j+1}^m+hD2_{i,j}E_{i,j1}^m+\frac{\mathrm{\Delta }w}{2\mathrm{\Delta }t}E_{i,j}^n.$$
(39)
Advancing a pseudo-timestep involves solving $`M`$ such equations in $`N\times N`$ matrices for the 1-implicit sweeps, and $`N`$ similar $`M\times M`$ matrix equations for the 2-implicit sweeps. A standard method involving $`LU`$ decomposition and forward- and back-substitution is used. In the case of periodic boundaries, the lower left and upper right corner elements of the matrices are non-zero, and this condition is dealt with using the Sherman-Morrison formula (Press et al. 1992) to compute a correction to the solution of the tridiagonal part. The diffusion coefficients are held fixed during the entire real timestep. Numerical accuracy of the solution is improved by scaling the lengths, times, and diffusion coefficients in equation (38) by a common factor chosen so the diagonal matrix elements are near unity.

To bring the solution into pseudo-time steady-state rapidly on different size scales, the length of the pseudo-timestep $`\mathrm{\Delta }w`$ is increased exponentially with iteration number $`m=1\mathrm{}W`$, so that

$$\mathrm{\Delta }w=\mathrm{\Delta }w_0\left(\frac{\mathrm{\Delta }w_1}{\mathrm{\Delta }w_0}\right)^{\frac{m1}{W1}},$$
(40)
where the initial pseudo-timestep $`\mathrm{\Delta }w_0`$ is one-quarter of the square of the grid spacing, and the final pseudo-timestep $`\mathrm{\Delta }w_1`$ is one-quarter of the square of the total grid size (Black & Bodenheimer 1975). In cases where the number of zones is large, the rescaling of lengths and times applied to equation (38) may not be sufficient to keep the condition numbers of the matrices near unity for all pseudo-timesteps, and the reduced accuracy of solutions obtained using $`LU`$ decomposition may lead to slow convergence. Under these conditions it might be helpful to iteratively improve the initial solution, or to use an alternative decomposition such as $`QR`$ or Householder (Press et al. 1992). For the test problems discussed in section 5, reducing the condition numbers of the matrices by adjusting the rescaling factor did not speed convergence.

The error in the time-advanced solution $`E^{n+1}`$ obtained from the series of $`w`$-steps is checked by substituting in the original difference equation (25). The residual is normalized against the ratio of the time-centered $`E`$ to the hydrodynamical timestep. If the biggest error on the grid is too large, the method is applied again with larger $`W`$ and more closely-spaced pseudo-timesteps. If a sufficiently small error is not reached after many pseudo-timesteps, then the real timestep for the flux divergence term is halved, and the term is applied twice per hydrodynamical timestep. If a sufficiently small error is not obtained with the real timestep finely divided, the calculation is terminated. On the other hand, if the error obtained is very small, then fewer pseudo-timesteps are tried initially in the following diffusion substep.

### 4.3 Matter-Radiation Interaction and Compressive Heating Terms

The absorption, emission, and compressive heating source terms are approximated by implicit difference equations, to obtain numerical stability for timesteps which may be long compared with the matter-radiation interaction timescale. The solution is advanced in these terms from step $`n`$ to step $`n+1`$ by solving simultaneously

$$E_{i,j}^{n+1}E_{i,j}^n=\mathrm{\Delta }t\{(𝐯:𝖯)_{i,j}^{n+1}+4\pi \kappa B_{i,j}^{n+1}c\kappa E_{i,j}^{n+1}\}$$
(41)
and

$$e_{i,j}^{n+1}e_{i,j}^n=\mathrm{\Delta }t\left\{p_{i,j}^{n+1}(𝐯)_{i,j}4\pi \kappa B_{i,j}^{n+1}+c\kappa E_{i,j}^{n+1}\right\}.$$
(42)
These depend on the values of $`e`$ and $`E`$ in zone $`i,j`$ only, so the spatial coordinate subscripts are dropped in subsequent equations. The gas is assumed to be ideal, with ratio of specific heats $`\gamma `$, and equation of state $`p^{n+1}=(\gamma 1)e^{n+1}`$. The Planck function $`B^{n+1}=\frac{1}{\pi }\sigma _BT^4`$ is computed from the temperature $`T=(\gamma 1)\mu e^{n+1}/(\rho ^n)`$ using the density $`\rho ^n`$ at the preceding time. Here $`\mu `$ is the dimensionless mean particle mass, $``$ the gas constant, and $`\sigma _B`$ the Boltzmann constant. The rate $`(𝐯:𝖯)^{n+1}`$ at which work is done on the radiation by the flow is

$$\left\{𝐯_{11}f11+𝐯_{22}f22+𝐯_{33}(1f11f22)+(𝐯_{12}+𝐯_{21})f12\right\}E^{n+1}.$$
(43)
The Eddington tensor components $`f11`$, $`f22`$, and $`f12`$, the fluid velocity $`𝐯`$, and the absorption opacity $`\kappa `$ are evaluated at timestep $`n`$, and the components of the velocity gradient tensor are as in Stone, Mihalas, & Norman (1992) Appendix A. With these choices, the only unknowns in equations (41) and (42) are $`e^{n+1}`$ and $`E^{n+1}`$. Eliminating $`E^{n+1}`$ between the two equations yields a quartic polynomial with $`e^{n+1}`$ as one root. With $`x`$ in place of $`e^{n+1}`$, the quartic is

$$x^4+\frac{(1+a_4)(1+a_2+a_3)}{a_1(1+a_3)}x\frac{(1+a_2+a_3)e^n+a_2E^n}{a_1(1+a_3)}=0,$$
(44)
where

$$a_1=4\kappa \sigma _B\left\{\frac{\mu (\gamma 1)}{\rho ^n}\right\}^4\mathrm{\Delta }t,$$
(45)
$$a_2=c\kappa \mathrm{\Delta }t,$$
(46)
$$a_3=\frac{(𝐯:𝖯)^{n+1}}{E^{n+1}}\mathrm{\Delta }t,$$
(47)
and

$$a_4=(\gamma 1)(𝐯)\mathrm{\Delta }t.$$
(48)
The quartic typically has a single positive real root when the timestep satisfies the Courant condition of §4.6. Writing the coefficient of $`x`$ in equation (44) as $`c_1`$ and the constant coefficient as $`c_0`$, this single root lies between zero, and the smaller of $`|c_0/c_1|`$ and $`|c_0|^{1/4}`$. The updated gas energy density is obtained by finding the root using Newton-Raphson iteration on this interval, with bisection when the Newton-Raphson method fails. The updated radiation energy density is determined by substituting the updated gas energy density in equation (41).

### 4.4 The Case of Zero Absorption Opacity

When the absorption opacity $`\kappa `$ is zero and the material-radiation interaction terms vanish, equations (41) and (42) are no longer coupled. In this instance, numerical energy conservation may be improved by using time-centered values for the gas and radiation pressures. In this differencing scheme, the time-advanced gas energy density $`e`$ is

$$e^{n+1}=e^n\left[\frac{1\frac{1}{2}q\mathrm{\Delta }t}{1+\frac{1}{2}q\mathrm{\Delta }t}\right],$$
(49)
with $`q=(\gamma 1)(𝐯)^n`$. For the radiation energy density update, $`q`$ is replaced by the quantity multiplying $`E^{n+1}`$ in equation (43).

### 4.5 Transport Step

Advection of the radiation energy proceeds exactly as described by SN for the hydrodynamic variables. Having expressed the advection terms in integral form in equation (23), we apply the conservative differencing scheme used for the hydrodynamic variables. In this scheme, the flux of advected radiation energy is computed from the mass flux in order to reduce the relative numerical diffusion of radiation with respect to the gas. The flux across every zone interface is computed using an interpolation method which may be selected as either donor cell, van Leer, or piecewise parabolic. These fluxes are then used to update the radiation energy density in a directionally split fashion.

Thus, the flux of radiation energy density along the 1-direction is constructed using

$$_{i,j}^1=(E/d)_{1,i,j}^{}\dot{M}_{i,j}^1g2a_i^{n+\frac{1}{2}}g31a_i^{n+\frac{1}{2}}dvl2a_j^n,$$
(50)
where the mass fluxes in the 1-direction, $`\dot{M}_{i,j}^1`$, are defined by equation  in SN. These fluxes are used to update the radiation energy under 1-transport via

$$E_{i,j}^{n+1}dvl1a_i^{n+1}dvl2a_j^nE_{i,j}^ndvl1a_i^ndvl2a_j^n=\mathrm{\Delta }t(_{i,j+1}^1_{i,j}^1).$$
(51)
Advection fluxes of radiation energy density in the 2-direction are computed by

$$_{i,j}^2=(E/d)_{2,i,j}^{}\dot{M}_{i,j}^2g31b_i^ndx1b_i^ng32a_j^{n+\frac{1}{2}},$$
(52)
where the mass fluxes in the 2-direction, $`\dot{M}_{i,j}^2`$, are defined by equation  in SN. These fluxes are then used to update the radiation energy due to advection in the 2-direction via

$$E_{i,j}^{n+1}dvl1a_i^ndvl2a_j^{n+1}E_{i,j}^ndvl1a_i^ndvl2a_j^n=\mathrm{\Delta }t(_{i,j+1}^2_{i,j}^2).$$
(53)
To reduce the systematic effects of the directional splitting discussed in section 5.3, the order in which the two advection substeps are applied is reversed each timestep.

### 4.6 Stability and the Timestep Criterion

Since the advection terms in the RHD equations are treated time-explicitly, the calculations may be numerically unstable unless the timestep is less than the time for waves or advection to carry energy across a grid zone (Richtmyer & Morton 1957). The wave family with the largest group speed may be either adiabatic sound waves or radiation acoustic waves, depending on the ratio of gas and radiation energy densities. The sound speed for purposes of computing the timestep is therefore chosen to be

$$c_s=\left[\mathrm{max}(\gamma ,\frac{4}{3})\frac{P_{tot}}{\rho }\right]^{1/2},$$
(54)
where $`P_{tot}`$ is the sum of the gas pressure and the largest component of the radiation pressure tensor. The remainder of the calculation of the timestep proceeds as in the hydrodynamic portion of ZEUS-2D (SN). Though resulting steps may be longer than the time for radiation to diffuse across a zone $`\frac{(\mathrm{\Delta }x)^2}{D}`$, the implicit differencing of the radiation diffusion term (§4.2) ensures stability of the method.

## 5 TEST CALCULATIONS

The problems used to test the radiation module fall into three categories. Calculations outlined in the first three sections below test in isolation the heating and cooling terms (§5.1), the radiation flux divergence term (§5.2), and the transport terms (§5.3). Linear RHD wave tests involving all the radiation terms are described in §5.4. Calculations which exercise radiation terms in the non-linear regime are propagation of optically-thin radiation fronts (§5.5), structure of radiating gas pressure dominated shocks of moderate optical depth (§5.6), and formation of a steady optically-thick shock in which radiation pressure dominates downstream (§5.7).

### 5.1 Heating and Cooling Terms

In a static uniform absorbing fluid initially out of thermal balance, an analytic solution for the time evolution of the gas energy may be obtained for the case where radiation energy dominates the total. The radiation energy emitted or absorbed in reaching equilibrium is a small fraction of the initial value, so radiation energy density $`E`$ is assumed constant. The time evolution of the gas energy density $`e`$ is obtained by solving the ordinary differential equation

$$\frac{de}{dt}=c\kappa E4\pi \kappa B(e).$$
(55)
Solutions are plotted together with the numerical results in figure 3, for a density of $`10^7`$ g cm<sup>-3</sup>, opacity $`4\times 10^8`$ cm<sup>-1</sup>, mean dimensionless mass per particle 0.6, and index in the equation of state $`\gamma =5/3`$. Two cases are shown. In the first, the initial thermal energy density is less than the equilibrium value. The agreement of the computed result with the analytic solution in this case is excellent. In the second case, $`e`$ is initially larger than the equilibrium value, lags the expected decline for the first few steps, and reaches equilibrium in the correct time. These results are consistent with the use of a fully-implicit method and a timestep of a few times $`10^{11}`$ seconds, much longer than the cooling time. When the timestep is held to less than the cooling time, the agreement between the calculation and the analytic solution is improved, as shown by the dotted line in figure 3. The sum of thermal and radiation energies is conserved in these three cases to about the floating-point precision.

### 5.2 Flux Divergence Term

An analytic solution of the diffusion equation (24) on the unit square, with periodic boundaries and unit diffusion coefficient, is (Lamb 1995)

$$E(x,y,t)=2+e^{8\pi ^2t}\mathrm{sin}2\pi x\mathrm{sin}2\pi y.$$
(56)
This function at $`t=0`$ is discretized onto a grid with 100 zones along each side, and the evolution followed using the radiation diffusion term alone. The timestep chosen is a factor of a hundred longer than the time $`\frac{(\mathrm{\Delta }x)^2}{D}`$ for radiation to diffuse across a zone. The largest zone-by-zone difference between the numerical and analytic solutions during the evolution towards equilibrium is 5.3%, and the largest root-mean-square difference is 2.5%. When the timestep is shortened by a factor ten, the corresponding residuals are 0.7% and 0.3%. The implicit numerical solution relaxes towards equilibrium more slowly than the analytic solution in each case.

### 5.3 Transport Terms

When the directionally-split advection substeps described in section 4.5 are applied in the same order each timestep, the splitting leads to a systematic delay in transport along the first axis, and an advance in transport along the second axis. As a test of this effect, the source terms were switched off, and a circular region of increased radiation energy density initially at the center of the domain was carried twice diagonally across a grid of $`32^2`$ zones with periodic boundaries. About 2000 timesteps were taken, and the advection substeps were applied always in the order 1, 2. At the end of this calculation, the centroid of the region lay off the grid diagonal by 0.02 zones along each axis. When the calculation was repeated with the order of the substeps reversed each timestep, the centroid remained within $`3\times 10^7`$ zones of the diagonal. We have chosen to reverse the order of advection each timestep in all further calculations discussed in this article, in order to minimize the effects of the directional splitting of the transport terms.

### 5.4 Linear RHD Waves

Waves in a radiating fluid may be damped by escape of photons. Damping is often rapid for oscillations with frequency near the rate at which material cools by radiating, and also for oscillations with wavelength near the distance radiation diffuses in a wave period. Mihalas & Mihalas (1983) obtained a linear dispersion relation for radiation and acoustic plane waves, driven by a boundary into a uniform medium with opacity due to absorption. They made the Eddington approximation, and assumed a uniform background in thermal equilibrium. The coefficients of the dispersion relation may be specified by two dimensionless parameters. The Boltzmann number $`\text{Bo}4\gamma c_ae/(cE)`$ is a ratio of typical rates of energy transport due to material advection and radiation in the background state, and $`rE/(4\gamma e)`$ fixes the ratio of the energy densities. In these definitions, $`\gamma `$ is the index in the equation of state and $`c_a`$ is the adiabatic sound speed. For the tests discussed below, $`\text{Bo}=10^3`$ and $`r=0.1`$. Three frequencies are considered, selected for optical depth per wavelength $`\tau _\lambda =10^3`$, unity, and 1000. The first of these frequencies is near the radiative cooling rate, and the last mode has wavelength near the radiative diffusion distance. In the test calculations, absorption opacity is set to 0.4 cm<sup>2</sup> g<sup>-1</sup> times the density. The domain is initially stationary and uniform, and the chosen Bo and $`r`$ correspond to density $`3.216\times 10^9`$ g cm<sup>-3</sup>, gas energy density 26020 erg cm<sup>-3</sup>, and radiation energy density 17340 erg cm<sup>-3</sup>. A region one wavelength long at the left of the domain is updated each timestep with the time and space variation found from the linear analysis. A wave propagates to the right from this driving region. Late-time results for waves of the three optical depths are shown in figures 4 to 6. There are ten zones per wavelength in each case. Timesteps were chosen according to the Courant condition of §4.6 for the calculations in the left-hand panel of each figure. The resulting steps, though having Courant numbers of 0.5, are longer than the time $`(\mathrm{\Delta }x)^2/D`$ for radiation to diffuse across a zone. The fully-implicit method used remains stable for these long timesteps, but loses accuracy, and the waves are damped more strongly than expected. Similar damping occurs when the Eddington approximation is made, indicating the effect is not due to the explicit differencing of the radiation diffusion coefficient in equation (25). In the right-hand panels of figures 4 to 6 are results from calculations in which the timesteps were reduced so the method adequately reproduced the linear analytic results. The wavelengths obtained are one to two percent shorter than the predicted wavelengths. Similar agreement with the linear predictions was observed for $`\text{Bo}=10^6`$ and $`r=1000`$.

### 5.5 Propagating Radiation Fronts in Optically Thin Media

In optically-thin regions, changes in radiation energy density propagate at the speed of light. To examine how well the method follows such changes, a domain 1 cm long is placed in thermal equilibrium with a radiation energy density of $`10^{22}`$ erg cm<sup>-3</sup>. The mass density is set to 0.025 g cm<sup>-3</sup> and the opacity to 0.4 cm<sup>2</sup> g<sup>-1</sup>, so optical depth across the domain is 0.01. One hundred grid zones are used. At time $`t=0`$, the radiation energy density in zones to the left of $`x=0.1`$ is raised to unity, and the increase is allowed to propagate to the right. The left-hand boundary condition is inflow, the right-hand outflow. Results using the LP flux limiter and scattering opacity are shown in figure 7. The front moves at approximately the speed of light. It is spread over a large number of zones, indicating the method is quite diffusive in this situation. For problems in which details of such fronts are important, it may be better to solve the time-dependent transfer equation using steps set by the radiation flow time. With timesteps longer than $`\mathrm{\Delta }x/(2c)`$, the front is wider still, however the method remains stable. The width is less when the diffusion coefficients are recalculated every diffusion substep, rather than once per hydrodynamic timestep. Where the leading edge of the front reaches the end of the grid, a region of uniform lower radiation energy density is seen. The effect is greatest when the timestep is long and the front is broad, as in the curve shown by dots on the bottom panel of figure 7. Radiation energy accumulates in this region because the boundary condition $`E=0`$ forces the radiative flux to zero at the boundary. Other conditions on radiation energy density which may be appropriate across outflow boundaries are fixed values of $`E`$, and $`𝐅=0`$.

The off-diagonal components of the Eddington tensor (equation (13)) may be non-zero in radiation fronts inclined with respect to the coordinate axes. As a test of the method in such a situation, the domain used in the one-dimensional radiation front calculation is expanded to a 1 cm square of $`100\times 100`$ zones. The high radiation energy density is placed initially in zones within 0.1 cm of grid center. The boundaries are made periodic in both directions, and the timestep is fixed at $`\mathrm{\Delta }x/(2c)`$. The distribution of radiation after $`10^{11}`$ seconds is shown in figure 8. The lowest contours on the figure are non-circular. After the leading edge of the front reaches the midpoints of the grid boundaries, the gradient in radiation energy density there is directed along the boundaries, rather than away from grid center. Since in the FLD approximation the flux is assumed parallel to $`E`$, an incorrect result is obtained in these regions. A similar error will result in FLD calculations whenever inclined radiation fronts meet. After the front has passed the corners of the grid, at $`5\times 10^{11}`$ sec, total energy differs from that at the start of the calculation by about one part in $`10^{14}`$, indicating that energy is conserved to high precision in this instance.

### 5.6 Subcritical and Supercritical Shocks

The structure of a shock in radiating fluid is sometimes altered by photons propagating upstream to preheat approaching material. If material reaching the shock front remains cooler than the shocked gas, the shock is termed subcritical, and the main effect is a slightly raised post-shock temperature, which declines downstream as radiation is emitted upstream. In stronger, supercritical shocks, upstream material is heated almost to the temperature of the post-shock gas, and material immediately upstream is in thermal equilibrium with the radiation emitted through the front. Gas reaching the front then overshoots the final temperature significantly, and cools to the post-shock temperature by radiating from a layer with optical depth less than unity. In this section we compare structures calculated for subcritical and supercritical shocks against analytic estimates by Zel’dovich & Raizer (1967), and numerical results obtained by Sincell, Gehmeyr, & Mihalas (1999) using a radiation hydrodynamics code with adaptive refinement of a one-dimensional grid. The problem is set up as by Sincell et al. (1999). The gas is given absorption opacity 0.4 cm<sup>2</sup> g<sup>-1</sup>, mean molecular weight 0.5, and uniform density $`7.78\times 10^{10}`$ g cm<sup>-3</sup>. Gas temperature $`10+75x/(7\times 10^{10}`$ cm$`)`$ K increases linearly with distance $`x`$ from the origin, and the radiation energy density is chosen for equilibrium. Grid spacing is set to $`1.4\times 10^8`$ cm, so initially the optical depth per zone is 0.044. A piston moves into the gas from the origin at uniform speed, and shocked material accumulates in front of the piston.

Results for a subcritical shock, with piston speed 6 km s<sup>-1</sup>, are shown in figure 9. The time is $`10^4`$ seconds. Temperature declines from the post-shock maximum more rapidly than expected, and the flux is too large in several zones downstream from the front. These effects may be due to an incorrect angular distribution of specific intensity. They are similar when the numerical resolution is doubled, suggesting they are not due to the finite optical depth of the front imposed by the artificial hydrodynamic viscosity. Sincell et al. (1999) found, by directly solving the transfer equation, that the Eddington factor near the shock front is slightly less than $`\frac{1}{3}`$ because the emitting layer has greatest line-of-sight optical depth when seen at grazing angles. In the FLD approximation, angular variation is assumed to be either smooth (LP, eq. ) or stepwise (Minerbo), and this feature is missed. Upstream, flux is higher with the LP than with the Minerbo flux limiter. As may be seen from figure 1, the flux is reduced below the Eddington value at larger optical depths with the Minerbo limiter.

The structure calculated for a supercritical shock with piston speed 16 km s<sup>-1</sup> is shown in figure 10. The time is again $`10^4`$ seconds. The temperature spike at $`\tau =0`$ is about one-fifth as tall as found by Sincell et al. (1999), and about five times as wide, due to the limited spatial resolution. The peak temperature in the spike is 314 K above the downstream temperature $`T_1`$, whereas Zel’dovich & Raizer estimate that the spike amplitude is $`(\frac{4}{\gamma +1}1)T_1=1619`$ K. With spatial resolution doubled, the spike is narrower and its height is larger at 506 K. The feature has approximately the expected degree of asymmetry, extending to a larger optical depth downstream than upstream. The flux in the equilibrium preheated region is slightly overestimated, but away from the front varies with optical depth in a manner close to that predicted by Zel’dovich & Raizer (1967). The flux upstream from the equilibrium region depends on the choice of limiter, with the Minerbo limiter producing smaller flux and lower Eddington factor, as in the subcritical case.

In both subcritical and supercritical shock calculations, the outcome is independent of the timestep for steps smaller than about ten zone radiation diffusion times $`\frac{(\mathrm{\Delta }x)^2}{D}`$. We conclude from these results that on a uniformly spaced grid, the FLD approximation is a poor choice for determining detailed structure of shocks of intermediate optical depth. It may be sufficiently accurate for multi-dimensional calculations in which such shocks must be followed approximately as part of a more complex flow.

### 5.7 Radiation-Dominated Shock

If material upstream and downstream from a radiating shock have sufficient optical depth, radiation is trapped near the interface, diffusing a limited distance upstream. In this section we examine whether the code adequately matches the jump conditions and thickness of such a shock. Parameters are chosen so the energy density is mostly in the gas upstream and in the radiation downstream. Opacity is assumed to be due to absorption as in §5.6, and material and radiation are assumed to be in thermal equilibrium far from the shock. Initially, the left of the grid is set to density 0.01 g cm<sup>-3</sup>, temperature $`10^4`$ K, and speed $`10^9`$ cm s<sup>-1</sup>, yielding a Mach number of 658. The right is set to density 0.0685847 g cm<sup>-3</sup>, temperature $`4.239\times 10^7`$ K, and speed $`1.458\times 10^8`$ cm s<sup>-1</sup>, as computed from the jump conditions. The boundary conditions are inflow on the left, and outflow on the right. After a brief transient, a steady shock is established. The situation after 50 flow crossing times is shown in figure 11. The approximate expected thickness of the shock is the distance $`l`$ for which the time to diffuse upstream, $`l^2/D`$, matches the time to sweep downstream, $`l/u_1`$ (Mihalas & Mihalas 1984). Here $`D`$ is the radiation diffusion coefficient and $`u_1`$ the upstream gas speed. The shock thickness obtained in the calculation is consistent with the resulting estimate of 7500 cm. The values downstream in the calculation differ by a maximum of 1.0% from those predicted using the jump conditions. When zone spacing is decreased a factor four, the discrepancies are reduced to 0.2%, indicating that mass, momentum, and energy are adequately conserved in material passing through the shock.

## 6 SUMMARY

We have described a set of subroutines which may be used with the ZEUS-2D magnetohydrodynamics code for radiation MHD calculations in the flux-limited diffusion approximation. Results were obtained typically in three to ten times more floating-point operations per timestep than required for similar calculations without radiation, with the greater number needed in the optically-thin tests. The primary limitation of flux-limited diffusion is that the flux is not evolved as an independent variable, but is assumed anti-parallel to the gradient in radiation energy density. This is incorrect in flows where shadows are cast or inclined radiation fronts cross. The incorporation in the flux limiter of a specified angular dependence for the radiation field has greatest effect in regions of intermediate to low optical depth, and may contribute to the broadening of radiation fronts. Flux limiters derived for angular dependences occurring in particular applications may be useful in obtaining more accurate solutions.

The limitations of FLD revealed in the tests mean that care must be taken in applying the module to problems with highly anisotropic radiation fields in optically thin regions. In such cases, full transport methods may be more appropriate. However, in optically thick flows, FLD is more efficient. The module has proven useful in studies of the internal dynamics of optically thick, radiation-dominated accretion disks (Agol et al. 2001; Turner & Stone 2001). Although the method remains stable for long timesteps, it shares with other implicit schemes the property that fluctuations with periods less than the timestep may be rapidly and unphysically damped, as seen in the linear wave and sub- and supercritical shock tests. When the diffusion timestep constraint of §4.6 is violated, it is advisable to check that results are independent of the timestep.

The radiation terms included in the equations solved are those of order unity in $`v/c`$ in at least one of the free-streaming, static diffusion, and dynamic diffusion regimes (Mihalas & Mihalas 1984; Stone, Mihalas, & Norman 1992). Adding higher-order terms would allow study of special relativistic effects such as radiation drag and photon viscosity. In the implementation described here, the radiation source terms are operator-split, with the flux divergence term evolved separately from the remainder. Linear waves and shocks to which both groups of terms contribute were evolved adequately, but situations may be encountered where the operator splitting leads to inaccuracies. We have assumed the gas is in local thermodynamic equilibrium, and relaxing this assumption may in the future allow a great variety of additional problems to be addressed.

Financial support of this work from the United States Department of Energy is gratefully acknowledged. We thank Eric Agol and Julian Krolik for their contributions to the tests of the module, and the anonymous referee for helpful comments on the presentation and the importance of the matrix condition numbers.


format: text

text:
# Large Mixing Angle MSW Solution in an 𝑆𝑈(3)_𝐿×𝑈(1)_𝑁 Gauge Model with Two-loop Radiative Mechanism

## Abstract

An $`SU(3)_L\times U(1)_N`$ gauge model, where $`E^i`$ ($`i`$=1,2,3) as heavy leptons are placed in the third member of lepton triplets, is shown to provide the large mixing angle (LMA) MSW solution to the solar neutrino problem. By introducing approximate conservation of $`L_eL_\mu L_\tau `$ ($`L^{}`$) into the model, we have estimated $`\mathrm{\Delta }m_{}^2`$ $``$ $`\sqrt{2}ϵ\mathrm{\Delta }m_{atm}^2m_{E^1}`$/$`(m_{E^2}+m_{E^3})`$, where $`ϵ`$ represents $`L^{}`$-breaking effects and $`m_{E^i}`$ is a mass of $`E^i`$. Since $`\mathrm{\Delta }m_{atm}^23\times 10^3`$ eV2, the LMA solution is obtained for $`ϵ0.010.1`$ in the case of $`m_{E^1}`$ $``$ $`m_{E^{2,3}}`$. This mechanism utilizes two-loop radiative mechanism, which contains a triplet version of the singly charged Zee boson, $`h^+`$, given by $`\xi `$ = ($`\xi ^{++}`$, $`\overline{\xi }^+`$, $`h^+`$)T. Almost bimaximal mixing is dynamically realized by the approximate equality between the $`\xi `$-couplings of the first family to the second and third families.

The recent analysis of solar neutrino oscillation data done by the Super-Kamiokande collaboration has so far shown no evidence for a distortion of the energy spectrum, nor for Earth regeneration effects, nor for seasonal variations of the neutrino flux . These observations suggest that solutions with large mixing angles are favored while solutions with small mixing angles as well as the vacuum oscillation (VO) solution are disfavored at the 95% confidence level. However, see also other analyses in Ref. for saving the disfavored solutions. Solar neutrino oscillations  characterized by large mixing angles are controlled by $`\mathrm{\Delta }m_{}^210^510^4`$ eV2 for the large mixing angle (LMA) MSW solution, $`\mathrm{\Delta }m_{}^210^7`$ eV2 for the LOW solution with low mass and low probability, $`\mathrm{\Delta }m_{}^210^9`$ eV2 for the quasi-VO (QVO) solution similar to VO but with small matter corrections taken into account  and $`\mathrm{\Delta }m_{}^210^{10}`$ eV2 for the VO solution. Such neutrino oscillations with large mixing angles are also observed for atmospheric neutrinos . The preliminary analysis on long-baseline neutrino oscillations done by K2K  has shown that $`\mathrm{\Delta }m_{atm}^23\times 10^3`$ eV2 well reproduces their data. Therefore, there are two mass scales, $`\mathrm{\Delta }m_{atm}^2=𝒪`$($`10^3`$) eV2 and $`\mathrm{\Delta }m_{}^2𝒪`$($`10^4`$) eV2, which indicate the hierarchy of $`\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2`$. Neutrino masses  are implied to be about several $`\times `$ 10-2 eV and such tiny neutrino masses can be generated by the seesaw mechanism or the radiative mechanism .

Since both oscillations are now characterized by large mixing angles, the observed pattern of the neutrino oscillations is compatible with bimaximal mixing scheme proposed earlier . One of the underlying physics behind bimaximal mixing is the presence of a new symmetry based on a lepton number of $`L_eL_\mu L_\tau `$ ($``$ $`L^{}`$) . Especially, in radiative mechanism of creating tiny neutrino masses , the use of $`U(1)_L^{}`$ based on $`L^{}`$ has recently been advocated . The $`U(1)_L^{}`$ symmetry applied to three active Majorana neutrino mass terms guarantees the hierarchy of $`\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2`$ with the maximal mixing for solar neutrinos. However, $`\mathrm{\Delta }m_{}^2`$ = 0 is also derived and is experimentally unacceptable. Finite but tiny amount of $`\mathrm{\Delta }m_{}^2`$ can be induced if tiny $`U(1)_L^{}`$-breaking interactions are present. There is a scenario that naturally contains tiny $`U(1)_L^{}`$-breaking interactions in the radiative mechanism, where one-loop radiative mechanism is responsible for $`\mathrm{\Delta }m_{atm}^2`$ while two-loop radiative mechanism accompanied by $`U(1)_L`$-breaking effects creates $`\mathrm{\Delta }m_{}^2`$ . Since two-loop radiative effects are much more suppressed than one-loop radiative ones, it is naturally expected to obtain $`\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2`$.

In the Zee’s radiative mechanism , a singly charged $`SU(2)_L`$-singlet Higgs scalar, $`h^+`$, which couples to charged lepton - neutrino pairs, is a key ingredient. It has been claimed that, in $`SU(3)_L\times U(1)_N`$ gauge models , $`h^+`$ together with the standard Higgs, $`\varphi `$, can be unified into a Higgs triplet either with ($`\varphi ^+`$, $`\varphi ^0`$, $`h^+`$) or with ($`\varphi ^0`$, $`\varphi ^{}`$, $`h^+`$) and that $`\mathrm{\Delta }m_{atm}^2`$ from one-loop radiative mechanism  as well as $`\mathrm{\Delta }m_{}^2`$ from two-loop radiative mechanism  are successfully generated. However, we have ended up with the VO and QVO solutions but not with the most favorable LMA solution . The reason lies in the fact that two-loop effects are too suppressed to obtain $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{atm}^2=𝒪`$(10-2) for the LMA solution. Roughly speaking, we estimated $`m^{(1)}`$ $`\stackrel{>}{}`$ $`fm_\tau ^2/(16\pi ^2v)`$ as a one-loop mass scale and $`m^{(2)}`$ $`\stackrel{<}{}`$ $`ϵf^2v/(16\pi ^2)^2`$ as a two-loop mass scale, where $`f`$ represents an $`h^+`$-coupling to a neutrino-charged lepton (or neutrino-heavy lepton) pair, $`v`$ $``$ $`v_{weak}`$ = $`(2\sqrt{2}G_F)^{1/2}`$ (= 174 GeV) and $`ϵ`$, which is presumably much smaller than 1, measures $`L^{}`$-breaking effects. This estimate gives $`\mathrm{\Delta }m_{atm}^2`$ $``$ $`m^{(1)2}`$ and $`\mathrm{\Delta }m_{}^2`$ $``$ $`m^{(1)}m^{(2)}`$ because of $`m^{(1)}`$ $``$ $`m^{(2)}`$, leading to $`\mathrm{\Delta }m_{}^2`$/$`\mathrm{\Delta }m_{atm}^2`$ $``$ $`m^{(2)}`$/$`m^{(1)}`$ $`\stackrel{<}{}`$ $`ϵf/16\pi ^2\times 10^4`$ ($`<`$ $`f/16\pi ^2\times 10^4`$), which must include the magnitude of $`𝒪`$(10-2) for the LMA solution. Therefore, the LMA solution at least requires $`f`$ $`\stackrel{>}{}`$ $`10^4`$, which is, however, not the case because $`\mathrm{\Delta }m_{atm}^2`$ $``$ $`3\times 10^3`$ eV2 fixes $`f`$ to be about 10-7, thereby, leading to $`\mathrm{\Delta }m_{}^2`$ $`\stackrel{<}{}`$10-8 eV2, which covers the VO and QVO solutions as mentioned.

In the present article, we seek other possibilities to be consistent with the LMA solution within the $`SU(3)_L\times U(1)_N`$ frameworks. A model with heavy leptons with the same charges as those of charged leptons  turns out to be relevant to provide the LMA solution. Since lepton triplets consist of ($`\nu _L^i`$, $`\mathrm{}_L^i`$, $`E_L^i`$) with $`i`$ = 1, 2, 3, this charge assignment requires the scalar of $`h^+`$ together with a new $`SU(2)_L`$-doublet Higgs ($`\xi ^{++}`$, $`\overline{\xi }^+`$) to be included in a Higgs triplet of ($`\xi ^{++}`$, $`\overline{\xi }^+`$, $`h^+`$). It is demonstrated that all neutrino masses are generated by two-loop radiative mechanism but not by one-loop radiative mechanism . The suppression of $`\mathrm{\Delta }m_{}^2`$ is essentially ascribed to tiny breaking of $`U(1)_L^{}`$ without the extra loop-suppression. The model is found to yield $`\mathrm{\Delta }m_{}^2`$ $``$ $`\sqrt{2}ϵm_{E^1}`$/$`(m_{E^2}+m_{E^3})`$ $`\mathrm{\Delta }m_{atm}^2`$, where $`m_{E^i}`$ is a mass of the $`i`$-th heavy lepton, which reproduces the LMA solution to solar neutrino problem for $`ϵ(0.010.1)`$ with $`m_{E^1}m_{E^{2,3}}`$.

In $`SU(3)_L\times U(1)_N`$ models, the minimal set of Higgs scalars consists of three triplets to be denoted by $`\eta `$, $`\rho `$ and $`\chi `$ in Eq.(5), which are minimal because, without these scalars, masses of triplet quarks would not be generated. In addition, our model uses a triplet version of the Zee scalar , $`\xi `$, and a doubly charged scalar of the Zee-Babu type , $`k^{++}`$, that initiates two-loop radiative mechanism. The particle content is specified by the $`U(1)_N`$ quantum number, $`N/2`$, related to the hypercharge, $`Y`$, as $`Y`$=$`\lambda ^8/\sqrt{3}+N`$, leading to $`Q_{em}=(\lambda ^3+Y)/2`$, where $`\lambda ^a`$ is the $`SU(3)`$ generator with Tr$`(\lambda ^a\lambda ^b)=2\delta ^{ab}`$ ($`a,b`$=1$``$8). Summarized as follows is the whole particle content in our model, where the quantum numbers are specified by ($`SU(3)_L`$, $`U(1)_N`$):

$`\psi _L^{i=1,2,3}=(\nu ^i,\mathrm{}^i,E^i)_L^T:(\text{3},2/3),\mathrm{}_R^{1,2,3}:(\text{1},1),E_R^{1,2,3}:(\text{1},1),`$ (1)
as leptons, where we have denoted heavy leptons $`E^i`$ by $`E^i`$ ($`i`$=1, 2, 3),

$`Q_L^1=(u^1,d^1,d^1)_L^T:(\text{3},0),Q_L^{i=2,3}=(d^i,u^i,u^i)_L^T:(\text{3}^{},1/3),`$ (2)
$`u_R^{1,2,3}:(\text{1},2/3),d_R^{1,2,3}:(\text{1},1/3),u_R^{2,3}:(\text{1},2/3),d_R^1:(\text{1},1/3),`$ (3)
as quarks, and

$`\eta =(\eta ^0,\eta ^{},\overline{\eta }^{})^T:(\text{3},2/3),\rho =(\rho ^+,\rho ^0,\overline{\rho }^0)^T:(\text{3},1/3),`$ (4)
$`\chi =(\chi ^+,\overline{\chi }^0,\chi ^0)^T:(\text{3},1/3),\xi =(\xi ^{++},\overline{\xi }^+,\xi ^+)^T:(\text{3},4/3),`$ (5)
as triplet Higgs scalars and

$$k^{++}:(\text{1},2),$$
(6)
as a singlet Higgs scalar. The Higgs scalars have the following vacuum expectation values (VEV’s):

$$0|\eta |0=(v_\eta ,0,0)^T,0|\rho |0=(0,v_\rho ,0)^T,0|\chi |0=(0,0,v_\chi )^T,$$
(7)
and quarks and leptons will acquire masses via these VEV’s, where the orthogonal choice of these VEV’s is ensured by appropriate Higgs interactions to be shown in Eq.(12). Quarks and leptons are collaborating to cancel gauge anomalies , where the obvious cancellation of the pure $`SU(3)_L`$ anomaly is done by two triplets and two anti-triplets.

The lepton number, $`L`$, associated with $`\psi _L^i`$ is assigned to other particles such that their possible interactions, which are $`SU(3)_L\times U(1)_N`$-singlets, maximally conserve the lepton number. For our newly introduced Higgs scalars, the possible interactions are $`\psi _L^i\psi _L^j\xi `$, $`\mathrm{}_R^i\mathrm{}_R^jk^{++}`$ and $`E_R^iE_R^jk^{++}`$, which can be used to determine $`L`$ = $``$2. Other $`L`$-assignment is shown in Table I. As a result, all Yukawa interactions can be described by the $`L`$-conserving ones.

In order to explain the observed behavior in neutrino oscillations, we also introduce another lepton number $`L_eL_\mu L_\tau `$ (=$`L^{}`$) and require its (approximate) conservation in our interactions. It is the tiny $`L^{}`$-breaking interactions that lead to the hierarchy of $`\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2`$. In order to import $`L^{}`$-breaking interactions, we employ the duplicate of $`k^{++}`$:

$$k^{++}:(\text{1},2),$$
(8)
but with the different $`L^{}`$ of $`L^{}`$ = $``$2 so as to develop possible Yukawa interactions, which will respect all relevant symmetries. The $`L^{}`$-assignment to other participating particles can been read off from Table I.

In addition to these global symmetries, another discrete symmetry based on $`Z_2`$ is used to suppress dangerous flavor-changing-neutral-currents (FCNC) interactions. Since our model contains quarks with the same charge, whose mass terms can be generated by $`\rho `$ and $`\chi `$ between $`Q_L^1`$ and down-type quarks and by $`\rho ^{}`$ and $`\chi ^{}`$ between $`Q_L^{2,3}`$ and up-type quarks, FCNC is induced at the phenomenologically unacceptable level . To avoid these interactions, Yukawa interactions must be constrained such that a triplet quark flavor gains a mass from only one Higgs field . The lepton sector also contains the similar FCNC problem because $`\mathrm{}^i`$ and $`E^i`$ ($`i`$=1, 2, 3) has the same charge. It turns out that such a constraint on FCNC is satisfied by introducing the following $`Z_2`$ symmetry into the model:

$`\psi _L^{1,2,3}\psi _L^{1,2,3},\mathrm{}_R^{1,2,3}\mathrm{}_R^{1,2,3},E_R^{1,2,3}E_R^{1,2,3},`$ (9)
$`Q_L^{1,2,3}Q_L^{1,2,3},u_R^{1,2,3}u_R^{1,2,3},d_R^{1,2,3}d_R^{1,2,3},u_R^{2,3}u_R^{2,3},d_R^1d_R^1,`$ (10)
$`\eta \eta ,\rho \rho ,\chi \chi ,\xi \xi ,k^{++}k^{++},k^{++}k^{++}.`$ (11)
Our Higgs interactions are, now, given by self-Hermitian terms composed of $`\varphi _\alpha \varphi _\beta ^{}`$ ($`\varphi `$ = $`\rho `$, $`\eta `$, $`\chi `$, $`\xi `$, $`k^{++}`$, $`k^{++}`$), which include the potential term of $`V_{\eta \rho \chi }`$

$`V_{\eta \rho \chi }=\lambda _1|\eta \times \rho |^2+\lambda _2|\rho \times \chi |^2+\lambda _3|\chi \times \eta |^2`$ (12)
with the definition of $`a\times b`$ $``$ $`ϵ^{\alpha \beta \gamma }a_\beta b_\gamma `$ and by non self-Hermitian Higgs potentials, $`V_0`$ and $`V_b`$:

$`V_0`$ $`=`$ $`\lambda _4(\eta ^{}\chi )(\xi ^{}\chi )+\lambda _5(\eta ^{}\rho )(\xi ^{}\rho )+\mu _0\xi ^{}\eta k^{++}+(h.c.),`$ (13)
$`V_b`$ $`=`$ $`\mu _b\xi ^{}\eta k^{++}+(h.c.),`$ (14)
where $`\lambda _i`$ ($`i=15`$) stands for a coupling constant and $`\mu _0`$ and $`\mu _b`$, respectively, denote $`L^{}`$-conserving and $`L^{}`$-breaking mass scale. The Higgs potential of $`V_0`$ conserves $`L^{}`$ and is responsible for the emergence of the bimaximal mixing while $`V_b`$ breaks the $`L^{}`$-conservation and is responsible for the emergence of the hierarchy of $`\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2`$. The Yukawa interactions are given by the following lagrangian:

$`_Y`$ $`=`$ $`{\displaystyle \underset{i=2,3}{}}f_{[1i]}ϵ^{\alpha \beta \gamma }\overline{\left(\psi _{\alpha L}^1\right)^c}\psi _{\beta L}^i\xi _\gamma +{\displaystyle \underset{i=1,2,3}{}}\overline{\psi _L^i}\left(f_{\mathrm{}}^i\rho \mathrm{}_R^i+f_E^i\chi E_R^i\right)`$ (17)
$`+{\displaystyle \underset{i=2,3}{}}\left[f_{\{i1\}}^{\mathrm{}}\overline{(\mathrm{}_R^i)^c}\mathrm{}_R^1+f_{\{i1\}}^E\overline{(E_R^i)^c}E_R^1\right]k^{++}+{\displaystyle \frac{1}{2}}\left[f_{\{11\}}^{\mathrm{}}\overline{(\mathrm{}_R^1)^c}\mathrm{}_R^1+f_{\{11\}}^E\overline{(E_R^1)^c}E_R^1\right]k^{++}`$
$`+\overline{Q_L^1}(\eta U_R^1+\rho D_R^1+\chi D_R^1)+{\displaystyle \underset{i=2,3}{}}\overline{Q_L^i}(\eta ^{}D_R^i+\rho ^{}U_R^i+\chi ^{}U_R^i)+(h.c.),`$
where $`f`$’s are Yukawa couplings with $`f_{[ij]}=f_{[ji]}`$ and $`f_{\{ij\}}^{\mathrm{},E}=f_{\{ji\}}^{\mathrm{},E}`$, the lepton mass terms are assumed to be diagonal for simplicity and right-handed quarks are denoted by $`U_R^i=_{j=1}^3f_{uj}^iu_R^j`$, $`D_R^i=_{j=1}^3f_{dj}^id_R^j`$, $`U_R^i=f_{u^{}2,3}^iu_R^{2,3}`$ and $`D_R^1=f_{d^{}1}^id_R^1`$.

We here note four comments on physical aspects of our interactions:

1. the suppression of FCNC is ensured by the absence of $`\chi D_R^1`$ and $`\rho D_R^1`$ for $`Q_L^1`$, of $`\chi ^{}U_R^i`$ and $`\rho ^{}U_R^i`$ for $`Q_L^{2,3}`$ and of $`\chi \mathrm{}_R^i`$ and $`\rho E_R^i`$ for $`\psi _L^i`$;
2. the orthogonal choice of VEV’s of $`\eta `$, $`\rho `$ and $`\chi `$ as in Eq.(7) is supported by $`V_{\eta \rho \chi }`$ if all $`\lambda `$’s are negative. It is because $`V_{\eta \rho \chi }`$ gets lowered if $`\eta `$, $`\rho `$ and $`\chi `$ develop VEV’s. So, one can choose VEV’s such that $`0|\eta _1|0`$ $``$ 0, $`0|\rho _2|0`$ $``$ 0 and $`0||\chi _3|0`$ $``$ 0;
3. the breaking of $`L`$-conservation is supplied by $`(\eta ^{}\chi )(\xi ^{}\chi )`$ and $`(\eta ^{}\rho )(\xi ^{}\rho )`$ in $`V_0`$; and
4. the breaking of $`L^{}`$-conservation is supplied by $`\xi ^{}\chi k^{++}`$ in $`V_b`$.

The additional effects of applying the $`Z_2`$ symmetry to the model include the absence of interactions given by

1. $`(\mathrm{}_R^1E_R^i+\mathrm{}_R^iE_R^1)k^{++}`$ ($`i`$=2,3) and $`\mathrm{}_R^1E_R^1k^{++}`$ that would give extra contributions to neutrino masses,
2. $`\eta \rho \chi `$ that would enhance the orthogonal choice of VEV’s, and
3. $`(\eta ^{}\rho )(\xi ^{}\chi )`$ that would induce one-loop radiative interactions.

Therefore, the suppression of FCNC simultaneously allows us to formulate our $`SU(3)_L\times U(1)_N`$ model with two-loop radiative mechanism but without one-loop radiative mechanism. The absence of one-loop radiative mechanism is essential for our discussions to obtain the LMA solution to the solar neutrino problem.

Now, let us explain how our $`SU(3)_L\times U(1)_N`$ model yields the LMA solution to the solar neutrino problem. Tiny neutrino masses are generated by two-loop radiative mechanism, where the interplay between $`\xi `$ and $`k^{++}`$ provides mixing masses for $`\nu _{eL}`$ and $`\nu _{\mu L,\tau L}`$ while the one between $`\xi `$ and $`k^{++}`$ provides tiny mixing masses among $`\nu _{\mu L,\tau L}`$. Indeed, these mixings are generated by interactions corresponding to $`L^{}`$-conserving diagrams as in Fig.1(a) and to $`L^{}`$-violating diagrams as in Fig.1(b), which contain Higgs couplings described by two types of non self-Hermitian Higgs potentials in Eq.(14). These interactions can be used to derive $`\mathrm{\Delta }m_{}^2`$ $``$ $`\sqrt{2}ϵm_{E^1}`$/$`(m_{E^2}+m_{E^3})`$ $`\mathrm{\Delta }m_{atm}^2`$, where $`ϵ`$ represents $`L^{}`$-breaking effects and $`m_{E^i}`$ is the mass of $`E^i`$. The LMA solution will correspond to $`ϵ0.010.1`$ for $`m_{E^1}`$ $``$ $`m_{E^{2,3}}`$.

From two-loop calculations, we obtain the following neutrino mass matrix parameterized by $`m_{1i}`$ calculated from $`L^{}`$-conserving two-loop diagrams and by $`m_{ij}^{}`$ calculated from $`L^{}`$-violating diagrams ($`i,j`$=2,3):

$$M_\nu =\left(\begin{array}{ccc}0& m_{12}& m_{13}\\ m_{12}& m_{22}^{}& m_{23}^{}\\ m_{13}& m_{23}^{}& m_{33}^{}\end{array}\right),$$
(18)
where

$`m_{1i}=2f_{[1i]}c_0\mu _0,m_{ij}^{}=f_{[1i]}f_{[1j]}c_b\mu _b`$ (19)
with

$`c_0`$ $`=`$ $`{\displaystyle \underset{i=2,3}{}}f_{[1i]}[\lambda _4f_{\{i1\}}^{\mathrm{}}m_\mathrm{}^im_ev_\chi ^2I(m_\mathrm{}^i^2,m_e^2,m_{\xi ^+}^2,m_{\overline{\eta }^{}}^2,m_{\xi ^+}^2,m_k^2)`$ (21)
$`+\lambda _5f_{\{i1\}}^Em_{E^i}m_{E^1}v_\rho ^2I(m_{E^i}^2,m_{E^1}^2,m_{\overline{\xi }^+}^2,m_\eta ^{}^2,m_{\overline{\xi }^+}^2,m_k^2)],`$
$`c_b`$ $`=`$ $`\lambda _4f_{\{11\}}^{\mathrm{}}m_e^2v_\chi ^2I(m_e^2,m_e^2,m_{\xi ^+}^2,m_{\overline{\eta }^{}}^2,m_{\xi ^+}^2,m_k^{}^2)`$ (23)
$`+\lambda _5f_{\{11\}}^Em_{E^1}^2v_\rho ^2I(m_{E^1}^2,m_{E^1}^2,m_{\overline{\xi }^+}^2,m_\eta ^{}^2,m_{\overline{\xi }^+}^2,m_k^{}^2).`$
Here, the mass of the $`i`$-th charged lepton and the mass of the $`i`$-th heavy lepton are, respectively, denoted by $`m_\mathrm{}i`$ ($`f_\mathrm{}i^iv_\rho `$) and $`m_{E^i}`$ ($`f_{E^i}^iv_\chi `$) and masses of Higgs scalars are denoted by the subscripts in terms of their fields . From the Appendix of Ref., we obtain the 2-loop integral, $`I`$, estimated as follows:

$`I(m_c^2,m_d^2.m_1^2,m_2^2,m_3^2,m_4^2)`$ $`=`$ $`{\displaystyle \frac{J(m_c^2,m_d^2,m_1^2,m_2^2,m_4^2)J(m_c^2,m_d^2,m_1^2,m_3^2,m_4^2)}{m_2^2m_3^2}},`$ (24)
$`J(m_c^2,m_d^2,m_1^2,m^2,m_4^2)`$ $`=`$ $`{\displaystyle \frac{1}{m_4^2}}G(m_c^2,m_1^2,m_4^2)G(m_d^2,m^2,m_4^2),`$ (25)
$`G(x,y,z)`$ $`=`$ $`{\displaystyle \frac{1}{16\pi ^2}}{\displaystyle \frac{x\mathrm{ln}(x/z)y\mathrm{ln}(y/z)}{xy}}.`$ (26)
The integral of Eq.(24) can be approximated by the chief contributions from $`m_{k,k^{}}^2`$ in Eq.(26) for $`m_{k,k^{}}^2`$ (other masses squared). Therefore, it can be taken to be universal for $`c_0`$ and $`c_b`$ and is to be denoted by $`I_{univ}`$. Hereafter, the dominance of $`m_{k,k^{}}^2`$ in masses of the model is assumed for the following analyses. We also assume that $`\lambda _4`$ $``$ $`\lambda _5`$.

The form of $`M_\nu `$ with the anticipation of $`m_{1i}m_{ij}`$ gives the following expressions of $`\mathrm{\Delta }m_{atm}^2`$ and $`\mathrm{\Delta }m_{}^2`$:

$`\mathrm{\Delta }m_{atm}^2=m_{12}^2+m_{13}^2(m_\nu ^2),`$ (27)
$`\mathrm{\Delta }m_{}^2=2\left(m_{22}^{}\mathrm{cos}^2\vartheta +2m_{23}^{}\mathrm{sin}\vartheta \mathrm{cos}\vartheta +m_{33}^{}\mathrm{sin}^2\vartheta \right)m_\nu ,`$ (28)
where $`\vartheta `$ is the mixing angle for solar neutrinos and is defined by $`\mathrm{cos}\vartheta =m_{12}/m_\nu `$ and $`\mathrm{sin}\vartheta =m_{13}/m_\nu `$. The maximal atmospheric neutrino mixing can be achieved by setting $`m_{12}`$ = $`m_{13}`$, thus requiring

$$f_{[12]}=f_{[13]},$$
(29)
as can been seen from Eq.(19). In the case that the heavy-lepton-contributions dominate the charged-lepton-contributions, which will be the case, our estimates of Eq.(23) yield

$`c_0f\lambda _5m_{E^1}v_\rho ^2\left(f_{\{12\}}^Em_{E^2}+f_{\{13\}}^Em_{E^3}\right)I_{univ},c_bf_{\{11\}}^E\lambda _5m_{E^1}^2v_\rho ^2I_{univ}`$ (30)
with $`f`$ = $`f_{[12]}`$ = $`f_{[13]}`$, which in turn give a relation of

$$\mathrm{\Delta }m_{}^2ϵ\frac{\sqrt{2}c_b}{c_0}\mathrm{\Delta }m_{atm}^2ϵ\frac{\sqrt{2}f_{\{11\}}^Em_{E^1}}{f_{\{12\}}^Em_{E^2}+f_{\{13\}}^Em_{E^3}}\mathrm{\Delta }m_{atm}^2,$$
(31)
where $`ϵ`$ ($``$ $`\mu _b/\mu _0`$) describes the $`U(1)_L^{}`$-breaking effects. The hierarchy of $`\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2`$ is simply ascribed to the breaking of $`U(1)_L^{}`$ due to $`ϵ`$ unless the couplings of $`f_{\{1i\}}^E`$ ($`i`$=2,3) and $`f_{\{11\}}^E`$ as well as the masses of $`E^1`$ and $`E^{2,3}`$ develop large hierarchy. In the case that $`f_{\{1i\}}^E`$ $``$ $`f_{\{11\}}^E`$, we reach a simpler relation expressed in terms of heavy lepton masses:

$$\mathrm{\Delta }m_{}^2ϵ\frac{\sqrt{2}m_{E^1}}{m_{E^2}+m_{E^3}}\mathrm{\Delta }m_{atm}^2.$$
(32)
Since $`\mathrm{\Delta }m_{atm}^2`$ $``$ 3$`\times `$10-3 eV2, we finally obtain $`\mathrm{\Delta }m_{}^2`$ $``$ $`ϵ\mathrm{\Delta }m_{atm}^2`$ for $`m_{E^1}`$ $``$ $`m_{E^{2,3}}`$, which, thus, covers the LMA solution for $`ϵ0.010.1`$, the LOW solution for $`ϵ`$ $``$ $`10^4`$, the QVO solution for $`ϵ`$ $``$ $`10^6`$ and the VO solution for $`ϵ`$ $``$ $`10^7`$.

More accurate numerical analysis is done by Gaussian integral method to evaluate the integral of Eq.(24). For the various parameters in this model, where the scale of the model is set by $`v_\chi `$ that breaks $`SU(3)_L\times U(1)_N`$ down to $`SU(2)_L\times U(1)_Y`$, we make the following assumptions on relevant free parameters to compute $`\mathrm{\Delta }m_{atm}^2`$ and $`\mathrm{\Delta }m_{}^2`$ in Eq.(28):

\- $`v_{\eta ,\rho }=v_{weak}/\sqrt{2}`$ to set $`(v_\eta ^2+v_\rho ^2)^{1/2}`$ = $`v_{weak}`$ for the weak boson masses, where $`v_{weak}`$ = $`(2\sqrt{2}G_F)^{1/2}`$ = 174 GeV,

\- $`m_{\eta ,\rho }=v_{weak}`$ and $`m_\xi =3v_{weak}`$ to suppress the $`\xi `$-contributions in low-energy phenomenology,

\- $`v_\chi v_{weak}`$ to suppress interactions for heavy leptons, exotic quarks and exotic gauge bosons, from which $`v_\chi `$ $`=`$ $`10v_{weak}`$, $`m_\chi =v_\chi `$ and $`m_{E^1}m_{E^2}m_{E^3}ev_\chi `$ are taken, where $`e`$ stands for the electromagnetic coupling,

\- $`m_{k,k^{}}=v_\chi `$ to suppress the $`k^{++}`$\- and $`k^{++}`$-contributions in low-energy phenomenology,

\- $`\lambda _5=2\times 10^6`$ (with $`\lambda _4\lambda _5`$) to reproduce $`\mathrm{\Delta }m_{atm}^2`$ and $`\mu _bv_\chi /10`$ to cover the LMA solution, where $`\lambda _51`$ and $`\mu _bv_\chi `$ can be considered to be natural because these couplings explicitly break the $`L`$\- or $`L^{}`$-conservation  while $`f_{[1i]}`$ $`=`$ $`e`$, $`f_{\{1i\}}^E`$ $`=`$ $`f_{\{11\}}^E`$ $`=`$ 1 and $`\mu _0=v_\chi `$ because the absence of these couplings does not enhance symmetries of the model.

The magnitude of $`\lambda _5`$ $``$ $`10^6`$ can been seen from the approximation of $`\mathrm{\Delta }m_{atm}^2`$ $`3.4\times 10^4\lambda _5^2(v_\chi /1\mathrm{GeV})^2`$ eV2 using $`I_{univ}`$ $``$ $`1/[(16\pi ^2)^2v_\chi ^4`$\], which is set equal to $`3\times 10^3`$ eV2. Our numerical analysis finally reproduces $`\mathrm{\Delta }m_{atm}^2`$ = $`3.03\times 10^3`$ eV2 for $`\lambda _5=2\times 10^6`$ and leads to $`\mathrm{\Delta }m_{}^2`$ = ($`0.212.14)\times 10^4`$ eV2 as the LMA solution for $`ϵ0.010.1`$.

It is useful to estimate $`U_{e3}`$, the matrix element that connects $`\nu _e`$ to $`\nu _3`$, since it can be used to classify various models for neutrino oscillations. This element of $`U_{e3}`$ is subject to the experimental constraint of $`|U_{e3}|^2\stackrel{<}{}0.0150.05`$ . In the present mass matrix of Eq.(18). this element is calculated to be

$$U_{e3}=\frac{1}{m_\nu }\left[\left(\mathrm{cos}^2\vartheta \mathrm{sin}^2\vartheta \right)m_{23}^{}+\mathrm{cos}\vartheta \mathrm{sin}\vartheta \left(m_{33}^{}m_{22}^{}\right)\right],$$
(33)
which turns out to vanish since $`m_{ij}^{}f_{[1i]}f_{[1j]}`$ ($`i`$=2,3) as in Eq.(19) and $`\mathrm{sin}\vartheta /\mathrm{cos}\vartheta =f_{[13]}/f_{[12]}`$ specific to the present two-loop kinematics. One can further check that, up to the third order of $`m_{ij}^{}`$, $`U_{e3}`$ is proportional to Eq.(33) as shown in the Appendix so that $`U_{e3}`$ vanishes up to the third order. Therefore, we expect that $`|U_{e3}|\stackrel{<}{}(m_{ij}^{}/m_\nu )^4`$, which is about $`(\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{atm}^2)^4\stackrel{<}{}10^4`$. It should be noted that this property of the vanishing $`U_{e3}`$ up to the third order does not rely upon Eq.(29) for the baimaximal mixing.

Two comments are in order. The first one concerns the choice of $`L^{}`$ for $`k^{++}`$. If $`L^{}`$ = 2 instead of $`L^{}`$ = $`2`$, $`k^{++}`$ couples to the 2nd and 3rd families. As a result, we reach the mass matrix with the (1,1)-entry and without other entries. Then, $`m_{11}^{}`$ is roughly proportional to $`(m_{E^2}+m_{E^3})^2`$, which results in $`\mathrm{\Delta }m_{}^2`$ $``$ $`ϵ\mathrm{\Delta }m_{atm}^2(m_{E^2}+m_{E^3})/\sqrt{2}m_{E^1}`$ estimated to be the similar magnitude to the one by Eq.(32) for $`m_{E^1}`$ $``$ $`m_{E^{2,3}}`$. Therefore, the same conclusion comes out. The second one is about phenomenological constraints on masses and couplings of $`\xi ^+`$, $`k^{++}`$ and $`k^{++}`$ because these contributions to leptonic processes of $`\mu `$\- and $`\tau `$-decays, $`e^{}e^{}e^{}e^{}`$ and $`\nu _ee^{}\nu _ee^{}`$ would disturb the well-established low-energy phenomenology. They are given by,

1. for $`\mu ^{}e^{}e^{}e^+`$ and $`\mu ^{}e^{}\gamma `$  induced by the collaboration of $`k^{++}`$ and $`k^{++}`$, $`\kappa f_{\{12\}}^{\mathrm{}}f_{\{11\}}^{\mathrm{}}/\overline{m}_k^2`$ $`\stackrel{<}{}`$ $`1.2\times 10^{10}`$ GeV-2 from $`B(\mu ^{}e^{}e^{}e^+)`$ $`<`$ $`10^{12}`$  and $`\stackrel{<}{}`$ $`2.4\times 10^8`$ GeV-2 from $`B(\mu ^{}e^{}\gamma )`$ $`<`$ $`1.2\times 10^{11}`$ , where $`\overline{m}_k`$ $``$ $`m_k`$ $``$ $`m_k^{}`$ and $`\kappa `$ estimated to be about $`\mu _b\mu _0/[16\pi ^2\overline{m}_k^2]`$ ($``$ 1) reads the suppression due to the approximate $`L^{}`$-conservation accompanied by the loop factor of $`(16\pi ^2)^1`$, which arises from the $`h^+`$-loop for the mixing of $`k^{++}`$ and $`k^{++}`$,
2. for $`\tau ^{}\mu ^{}e^{}e^+`$ and $`\tau ^{}\mu ^{}\gamma `$ mediated by $`k^{++}`$ and by $`\xi ^+`$, $`\left|f_{\{13\}}^{\mathrm{}}f_{\{12\}}^{\mathrm{}}/\overline{m}_k^2\right|`$ $`\stackrel{<}{}`$ $`2.1\times 10^7`$ GeV-2 from $`B(\tau ^{}\mu ^{}e^{}e^+)`$ $`<`$ $`1.7\times 10^6`$  and $`\stackrel{<}{}`$ $`4.2\times 10^6`$ GeV-2 from $`B(\tau ^{}\mu ^{}\gamma )`$ $`<`$ $`1.1\times 10^6`$  and $`\left|f_{[13]}f_{[12]}/m_{\xi ^+}^2\right|`$ $`\stackrel{<}{}`$ $`4.2\times 10^6`$ GeV-2 also from $`B(\tau ^{}\mu ^{}\gamma )`$,
3. for $`e^{}e^{}e^{}e^{}`$ mediated by $`k^{++}`$ , $`\left|f_{\{11\}}^{\mathrm{}}/m_k^{}\right|^2`$ $`\stackrel{<}{}`$ $`4.8\times 10^5`$ GeV-2 and,
4. for $`\nu _\mu e^{}\nu _\mu e^{}`$ mediated by $`\xi ^+`$ , $`\left|f_{[12]}/m_{\xi ^+}\right|^2`$ $`\stackrel{<}{}`$ $`1.7\times 10^6`$ GeV-2.

It should be noted that the leading contribution of $`\xi ^+`$ to $`\mu ^{}e^{}\gamma `$, which gives the most stringent constraint on $`\xi ^+`$, is forbidden by the $`U(1)_L^{}`$-invariant coupling structure. Our choice of $`f_{[1i]}`$ $``$ $`e`$ with $`m_{\xi ^+}`$ $``$ $`3v_{weak}`$ is consistent with the constraints. As far as these constraints are satisfied, the charged lepton-exchanges in Eq.(23) turn out to yield well suppressed contributions compared with those from the heavy-lepton exchanges. §§§ The charged-lepton contributions are suppressed as far as $`\lambda _4`$ $`\stackrel{<}{}`$ $`\lambda _5m_{E^1}m_{E^3}v_\rho ^2/(m_em_\tau v_\chi ^2)`$, roughly, leading to $`\lambda _4`$ $`\stackrel{<}{}`$ 0.8.

In summary, we have demonstrated that the solar neutrino mixing based on the currently most favorable LMA solution is compatible with the one discussed in the specific class of $`SU(3)_L\times U(1)_N`$ gauge models that involves $`\psi _L=(\nu ^i,\mathrm{}^i,E^i)_L^T`$ as the lepton triplets and $`\xi =(\xi ^{++},\overline{\xi }^+,\xi ^+)^T`$ as the Higgs triplet. Neutrinos acquire Majorana masses through the two-loop radiative mechanism combined with the approximate conservation of $`L^{}`$ (=$`L_eL_\mu L_\tau `$). Our two-loop radiative mechanism utilizes $`\xi `$ that couples to $`\psi _L^1\psi _L^{2,3}`$ and $`k^{++}`$ with $`L^{}`$ = 0 that couples to $`e_R\mathrm{}_R^i`$ and to $`E_R^1E_R^i`$ ($`i`$=2,3) as well as $`k^{++}`$ with $`L^{}`$ = $`2`$ that couples to $`e_Re_R`$ and to $`E_R^1E_R^1`$. As a result, the neutrino mass matrix exhibits the bimaximal structure if the approximate equality between $`f_{[12]}`$ and $`f_{[13]}`$ is realized. The relation of $`\mathrm{\Delta }m_{}^2`$ $``$ $`\sqrt{2}ϵ\mathrm{\Delta }m_{atm}^2m_{E^1}/(m_{E^2}+m_{E^3})`$ successfully covers the LMA solution.

The work of M.Y. is supported by the Grant-in-Aid for Scientific Research No 12047223 from the Ministry of Education, Science, Sports and Culture, Japan.

Appendix

The matrix element $`U_{e3}`$ can be defined by

$$U_{e3}=\frac{1}{\sqrt{2}}\left(c_{31}^{(1)}+c_{31}^{(2)}+c_{31}^{(3)}c_{32}^{(1)}c_{32}^{(2)}c_{32}^{(3)}\right),$$
(34)
where $`c_{3i}^{(1,2,3)}`$ ($`i`$=1,2) are the coefficients in the perturbative expansion of $`|\nu _3`$ as the eigenstate of $`\nu _3`$ for $`M_\nu `$ of Eq.(18) in terms of $`m_{ij}^{}`$ ($`i,j`$=2,3). In fact, $`U_{e3}`$ is expressed in terms of unperturbed state of $`|i^{(0)}`$ ($`i`$=1,2,3):

$`|\nu _3=\left(\begin{array}{ccc}U_{e3}& & \\ & & \\ & & \end{array}\right)=|3^{(0)}+{\displaystyle \underset{i=1,2,3}{}}\left(c_{3i}^{(1)}+c_{3i}^{(2)}+c_{3i}^{(3)}\right)|i^{(0)},`$ (38)
up to the third order. Since $`|3^{(0)}`$ has no relevant entry for $`U_{e3}`$, $`|3^{(0)}`$ and $`c_{33}^{(1,2,3)}`$ can be neglected. The unperturbed states are defined by $`|1^{(0)}=(1,c,s)^T/\sqrt{2}`$, $`|2^{(0)}=(1,c,s)^T/\sqrt{2}`$ and $`|3^{(0)}=(0,s,c)^T`$ for $`c=\mathrm{cos}\vartheta `$ and $`s=\mathrm{sin}\vartheta `$, whose masses are given by ($`m_1^{(0)}`$, $`m_2^{(0)}`$, $`m_3^{(0)}`$) = ($`m_\nu `$, $`m_\nu `$, 0), respectively.

The coefficients can be determined by the well-known formula:

$`c_{ij}^{(1)}=j^{(0)}|M_\nu ^{}|i^{(0)}/\left(m_j^{(0)}m_i^{(0)}\right),c_{ij}^{(2)}=\left({\displaystyle \underset{k}{}}c_{ik}^{(1)}j^{(0)}|M_\nu ^{}|k^{(0)}m_i^{(1)}c_{ij}^{(1)}\right)/\left(m_j^{(0)}m_i^{(0)}\right),`$ (39)
$`c_{ij}^{(3)}=\left({\displaystyle \underset{k}{}}c_{ik}^{(2)}j^{(0)}|M_\nu ^{}|k^{(0)}m_i^{(1)}c_{ij}^{(2)}m_i^{(2)}c_{ij}^{(1)}\right)/\left(m_j^{(0)}m_i^{(0)}\right),`$ (40)
$`c_{ii}^{(1)}=0,c_{ii}^{(2)}={\displaystyle \underset{j}{}}|c_{ij}^{(1)}|^2/2,c_{ii}^{(3)}={\displaystyle \underset{j}{}}\left(c_{ij}^{(2)}c_{ij}^{(1)}+c_{ij}^{(1)}c_{ij}^{(2)}\right)/2,`$ (41)
where

$`m_i^{(1)}=i^{(0)}|M_\nu ^{}|i^{(0)},m_i^{(2)}={\displaystyle \underset{j}{}}c_{ij}^{(1)}i^{(0)}|M_\nu ^{}|j^{(0)},`$ (42)
and $`M_\nu ^{}`$ is the $`3\times 3`$ matrix only with $`m_{ij}^{}`$ in the same entries as in $`M_\nu `$. These equations lead to

$`U_{e3}`$ $`=`$ $`\sqrt{2}c_{31}^{(1)}\sqrt{2}\left[\left(1^{(0)}|M_\nu ^{}|2^{(0)}+m_1^{(1)}m_3^{(1)}\right)c_{31}^{(2)}+1^{(0)}|M_\nu ^{}|3^{(0)}c_{33}^{(2)}\right]/m_\nu `$ (43)
with

$`c_{31}^{(1)}`$ $`=`$ $`\left[\left(c^2s^2\right)m_{23}^{}+cs\left(m_{33}^{}m_{22}^{}\right)\right]/\sqrt{2}m_\nu ,c_{33}^{(2)}=c_{31}^{(1)2},`$ (44)
$`c_{31}^{(2)}`$ $`=`$ $`\left(s^2m_{22}^{}2csm_{23}^{}+c^2m_{33}^{}\right)c_{31}^{(1)}/m_\nu ,`$ (45)
where we have used various relations such as $`c_{31}^{(1)}=c_{32}^{(1)}`$, $`c_{31}^{(2)}=c_{32}^{(2)}`$, $`m_1^{(1)}=m_2^{(1)}`$ and so on. Especially, the relation of $`c_{31}^{(2)}=c_{32}^{(2)}`$ makes the second order $`U_{e3}`$ to vanish. We, then, find that $`U_{e3}`$ is proportional to $`c_{31}^{(1)}`$, which vanishes because of the proportionality of $`m_{ij}^{}f_{[1i]}f_{[1j]}`$ ($`i,j`$=2,3) indicated by Eq.(19) with $`f_{[12]}c`$ and $`f_{[13]}s`$.

Table Captions

TABLE I: $`L`$ and $`L^{}`$ quantum numbers.

Figure Captions


text:
# Hole carrier in MgB2 characterized by Hall Measurements

\[

## Abstract

The longitudinal resistivity $`\rho _{xx}`$ and Hall coefficient $`R_H`$ were measured for MgB2 sintered under high pressure. We found that $`R_H`$ is positive like cuprate high-$`T_c`$ superconductors, and decreases as temperature increases for 40 K $`<T<`$ 300 K. The cotangent of Hall angle was found to follow $`a+b`$$`T`$2 behavior from $`T_c`$ to 300 K. At $`T`$ = 100 K, $`R_H=4.1\times 10^{11}`$ m3/C from which hole carrier density was determined to be $`1.5\times 10^{23}`$/cm3. This carrier density is 2 - 3 orders of magnitude larger than those of Nb3Sn and optimally doped YBa2Cu3Oy superconductors.

\]

Recently, MgB2 was found to be metallic superconductor with transition temperature ($`T_c`$) of about 40 K., and has provided great scientific interest. Several thermodynamic parameters have been estimated, such as a upper critical field $`H_{c2}=`$ 13 - 18 T, a Ginzburg-Landau parameter $`\kappa `$ 26, and the critical supercurrent density $`J_c(0)`$ 105 A/cm2. In order to probe the nature of gap, tunneling spectroscopy measurements have been reported , and they observed superconducting energy gap ($`\mathrm{\Delta }`$) of 5 - 7 meV in the framework of the BCS model. The conventional BCS electron-phonon interaction was proposed as the origin of the superconductivity based on a band calculation. The possible origin of the enhanced $`T_c`$ is suggested to originate from a strong electron-phonon interaction and a enhanced phonon frequency due to the light boron mass in MgB2. Most of the charge carrier density at the Fermi level comes from the boron band. Indeed, the boron isotope effect has been reported with an exponent of $`\alpha _B0.26`$.  No experimental study on the electronic structure has been reported yet.

Another interesting feature concerning the normal-state Hall effect in high-$`T_c`$ cuprates is the universal temperature dependence of the cotangent of the Hall angle (cot$`\theta _H`$). Anderson and Chien et al. have proposed that the charge transport is governed by two different scattering times with different temperature dependences. In this model, the cot$`\theta _H`$ should be proportional to $`T^2`$ since the Hall angle is proportional to the inverse of the Hall scattering time $`\tau _H`$ ($`T^2`$), and has been observed for most high-$`T_c`$ superconductors. In the mixed-state, the flux-flow Hall effect is also quite interesting. A puzzling sign anomaly has been observed in some conventional superconductors as well as in most of high-$`T_c`$ superconductors. Even double or triple sign changes have been observed in some high-$`T_c`$ superconductors. Furthermore, a universal scaling behavior between the Hall resistivity and the longitudinal resistivity has attracted much experimental  and theoretical interest. However, these Hall effects in the mixed state are not well understood.

To understand the superconductivity in MgB2, it is essential to know the type of charge carrier and it’s density, but these have not been reported yet . Theoretically, Hirsch  proposed that the 40 K superconductivity of MgB2 originates mainly from the hole carriers with boron planes acting like the CuO2 planes in cuprate high-temperature superconductors. He proposed that pairing of hole carriers leads to hole undressing, which is driven by Coulomb interactions. To the best of our knowledge, the Hall coefficient ($`R_H`$) for MgB2 has not been reported. For other metal diborides (MB2, where M = Y, Nb, Ti, V, Cr, Zr, Mo, Ta, or W), the sign of $`R_H`$ was observed to be negative, and no observation of superconductivity has been reported. Therefore, the measurement of the Hall effect for MgB2 has received much scientific attention.

To obtain reliable results from transport measurements by using polycrystalline samples, one must make the sample strong and dense. In this case, samples sintered under high pressure suitable. In our previous report,  we showed that the mechanical properties, as well as the superconducting properties, were vastly enhanced for samples sintered at 950 C under high pressure (3 GPa range).

In this paper, we report the first measurement of the $`R_H`$ of MgB2, which was carried out using carefully prepared sample with a thin bar shape. We found that the sign of $`R_H`$ is positive like those of cuprate high-$`T_c`$ superconductors. This is contrary to most other metal diborides with the same structure as MgB2. Also the $`R_H`$ decreased as the temperature increased, and the cotangent of the Hall angle follows $`a+bT`$1.8 for most of the measured temperature region from 40 to 300 K. At T = 100 K, $`R_H=4.1\times 10^{11}`$ m3/C, and the calculated hole carrier density is $`1.5\times 10^{23}`$/cm3.

The polycrystalline samples (4.5 $`mm`$ in diameter and 3.3 $`mm`$ in height) used in this study were sintered at 950 $`{}_{}{}^{}C`$ under 3 GPa. The fabrication method was reported in detail by Jung et al..  The sample purity was more than 99% as determined by X-ray diffraction analysis. No grain boundries were observed using the scanning electron microscopy. In order to obtain a higher Hall voltage signal, we cut the sample into a bar shape with a length of 4 $`mm`$ and a width of 2.4 $`mm`$, and then mechanically polished it until it was very thin ( 50 - 100 $`\mu m`$). The standard photolithography technique was adopted to align the electrical pads shown in the upper inset of Fig. 1. To obtain good ohmic contacts ($`<`$1 $`\mathrm{\Omega }`$), we coated Au film on contact pads after cleaning the sample surface with an Ar ion beam. This process was done $`insitu`$ in a high vacuum chamber. The voltage noise, which is detrimental to precise measurements, was successfully reduced to a lower level by preparing a very thin and optically clean specimen polished from a strong, dense samples. After installing a low noise preamplifier (N11, EM Electronics ) prior to the nanovoltmeter (HP 34420A), we achieved a voltage resolution of below 1 nV under a vias current of 50 -100 mA. Fine temperature control was crucial since the Hall signal was very small. The longitudinal and Hall voltages were measured simultaneously by using the standard dc 6-probe method. The magnetic field was applied perpendicular to the sample surface by using a superconducting magnet system (Maglab2000 Oxford Ltd.) and the applied current density was $`42`$ A/cm2. The Hall voltage was extracted from the antisymmetric parts of the transverse voltages measured under opposite directions to remove the longitudinal component due to the misalignment of the Hall voltage pads. The Hall voltage was found to be linear in both the current and the magnetic field.

Fig. 1 shows the temperature dependence of the longitudinal resistivity, $`\rho _{xx}`$. The low-field magnetization in the zero-field-cooling state for the original bulk sample is shown in the lower inset of Fig. 1. The diamagnetism is 100% to almost $`T_c`$; thus we normalized to the value at lower temperature. The superconducting transition temperature is 38.4 K with a narrow transition width of $``$0.6 K, as judged from the 10 to 90% superconducting transition. The resistivity value of $`\rho `$70 $`\mu \mathrm{\Omega }`$cm at 300 K is comparable to that of single crystalline intermetallic superconductors. As reported in our previous work, the normal-state $`\rho _{xx}`$ follows roughly a $`T^2`$ behavior rather than a $`T^3`$ behavior, for the entire temperature region below room temperature. No magnetoresistance was observed from $`T_c`$ to 300 K , which is consistent with the previous results by Jung et al.  and Takano et al., but different from the data by Finnemore et al..

The temperature dependence of the Hall coefficient is shown in Fig. 2. The two curves in the inset represent the Hall voltage measured at 100 K for opposite magnetic fields up to 5 T. The clearly symmetric and linear shape demonstrates that the signal to noise ratio for our measurement is high. The Hall coefficient was positive for all temperatures above $`T_c`$. At 100 K, $`R_H=4.1\times 10^{11}`$ m3/C, and the hole carrier density was calculated to be $`1.5\times 10^{23}`$ /cm3. The absolute value of the hole carrier density is two orders of magnitude larger than that of Nb3Sn superconductors and nearly three orders of magnitude larger than that of optimally doped YBa2Cu3Oy.

Hirsch offers an explanation based on a universal mechanism by assuming that superconductivity in MgB2 is similar to that in cuprate superconductors and is driven by pairing of heavily dressed hole carriers in a band that is almost full, whereby they gain enough kinetic energy to overcome the Coulomb energy. Based on this assumption, he claimed that the type of the charge carrier is positive.

In Fig. 3, we show the temperature dependence of cot$`\theta _H`$ at 5 T. A good linear fit to $`a+bT^{1.8}`$, rather than $`a+bT^2`$, is observed for the temperature range from $`T_c`$ to 300 K, In high-$`T_c`$ cuprates, the charge transport is governed by two different scattering times with different temperature dependences According to this two-scattering-time model, the longitudinal conductivity ($`\sigma _{xx})`$ is governed by the transport scattering time $`\tau _{tr}`$, which is proportional to $`1/T`$, whereas the Hall conductivity ($`\sigma _{xy}\tau _H\tau _{tr}`$) follows $`1/T^3`$ since the Hall relaxation rate is proportional to $`1/T^2`$. As a result, the cot$`\theta _H`$ (= $`\sigma _{xy}/\sigma _{xx}`$) should follow an $`T^{2\text{ }}`$law. Our data is in fair agreement with a cot$`\theta _HT^{2\text{ }}`$law, which is consistent with the observations in most high-Tc superconductors. However, this observation cannot be explained by the above model because the data show $`\rho _{xx}T^2`$, as shown in Fig. 1.

In summary, we report the temperature dependence of $`R_H`$ for the recently discovered binary superconductor MgB2 which has a remarkably high transition temperature. We find that $`R_H`$ is positive like those for cuprate high-$`T_c`$ superconductors and that the cot$`\theta _H`$ follows $`a+b`$$`T`$2 from $`T_c`$ to 300 K. At T = 100 K, $`R_H=4.1\times 10^{11}`$ m3/C, and the hole carrier density was $`1.5\times 10^{23}`$/cm2. We discussed the implication of the hole superconductivity based on a recent model.

###### Acknowledgements.

We appreciate valuable discussion with J. L. Tallon, This work is supported by the Ministry of Science and Technology of Korea through the Creative Research Initiative Program.