text:
# Rayleigh-Bénard Convection with a Radial Ramp in Plate Separation

## I Introduction

Rayleigh-Bénard convection in a thin horizontal fluid layer heated from below is a canonical example of pattern formation in a continuous dissipative system far from equilibrium Cross and Hohenberg (1993). Under various conditions, wavenumber selection mechanisms have been identified Cross et al. (1986); Catton (1988); Cross and Hohenberg (1993); Getling (1998) that reduce the band of stable wavenumbers, sometimes to a single value. One such selection mechanism occurs when there is a one-dimensional spatial variation or ramping of the control parameter $`ϵ`$; where $`ϵ(RR_c)/R_c`$, $`R`$ is the Rayleigh number, and $`R_c`$ its critical value Eagles (1980); Kramer et al. (1982); Pomeau and Zaleski (1983); Kramer and Riecke (1985); Buell and Catton (1986); Bajaj et al. (1999). This can be accomplished by varying the plate separation $`d`$ such that $`ϵ`$ goes from $`ϵ=ϵ_o>0`$ in the bulk of the layer (i.e., the unramped region) to $`ϵ<0`$ as a lateral boundary is approached. It is expected that in the idealized case of an infinitely gradual one-dimensional ramp the wavenumber will equal $`k_c`$ at the position where the layer depth yields $`R_c`$. It has been shown under very general conditions, that as long as the layer becomes critical somewhere along the ramp this is sufficient to fix the wavenumber $`k_s`$ in the bulk and over the rest of the ramp Kramer et al. (1982). For slightly supercritical conditions it is expected that the selected wavenumber in the bulk can be expressed as

$$k_s=\stackrel{~}{k}_c+\alpha ϵ_o$$
(1)
where $`\stackrel{~}{k}_c=k_c=3.117`$ Bajaj et al. (1999) and $`\alpha `$ depends on the Prandtl number, $`\sigma `$, and the specifics of the ramp.

Recent Rayleigh-Bénard convection experiments Bajaj et al. (1999); Ahlers et al. (2001) in a cylindrical cell with a two-dimensional radial ramp in plate separation have generated intriguing results. Specifically, the plate separation as a function of radius used in experiment is, using the layer depth $`d`$ to nondimensionalize,

$$h(r)=\{\begin{array}{cc}1,\hfill & r<r_0\hfill \\ 1\delta _r\left[1\mathrm{cos}\left(\frac{rr_0}{r_1r_0}\pi \right)\right],\hfill & rr_0\hfill \end{array}$$
(2)
where $`\delta _r`$ is a constant and the radius values $`r_0`$ and $`r_1`$ are the locations where the ramp begins and ends, respectively. The ramp always extends to the sidewall. A schematic of the cosine ramp is shown in Fig. 1; note that $`r_0`$ and $`r_1`$ are geometric constants but $`r_c`$, the location where the plate separation yields $`R_c`$, is a function of the ramp shape and $`ϵ_0`$ such that, for a ramp given by Eq. (2), $`r_c/ϵ_0>0`$.

Experiments using the cosine ramp defined by Eq. (2) have yielded unexpected results for the wavenumber Ahlers et al. (2001); Bajaj et al. (1999). Mean pattern wavenumber measurements (using the Fourier methods discussed in Morris et al. (1993)) yielded $`\stackrel{~}{k}_c=2.97<k_c`$. Additionally, measurements of the local wavenumber defined at each position in space (method discussed in Egolf et al. (1998)) displayed interesting variation as $`ϵ_0`$ is increased. For the time independent patterns the bulk region, $`rr_0`$, contains approximately straight parallel rolls. Near threshold, $`ϵ_00.048`$, a centered egg-shaped domain of convection rolls with small wavenumber, $`k<k_c`$, extends through the convection cell with the long-axis parallel to the roll axes. A dramatic roll expansion from $`k3.6`$ at the edge where the ramp begins to $`k2.6`$ in the center of the domain is observed. As $`ϵ_0`$ increases, the wavenumber field evolves into a domain characterized by large wavenumbers extending through the layer with the long-axis perpendicular to the roll axes (see Fig. 3 of Bajaj et al. (1999)). Similar experiments without a ramp do not exhibit this wavenumber behavior Hu et al. (1993, 1994, 1995).

In the experiments with a radial ramp in plate separation, for $`ϵ_00.03`$, time dependent states were found through the repeated formation of defects via an Eckhaus mechanism consistent with the local wavenumbers crossing the Eckhaus stability boundary for an ideal infinite layer of two-dimensional rolls. Also in the ramped experiments, for $`ϵ_00.18`$, defects were formed via a skewed varicose mechanism consistent with the local wavenumber exceeding the skewed varicose stability boundary for an ideal infinite layer of two-dimensional rolls similar to what has been observed in experiments in unramped cylindrical domains with rigid sidewalls Croquette and Williams (1989).

It has been suggested that these features may be the result of the interaction of the convective roll pattern and weak large scale flows Bajaj et al. (1999). The visualization and quantification of these large scale flows is not possible in the current generation of experiments. However, we are able to make these measurements by performing full numerical simulations with a new spectral element code (discussed further in Section III). We utilize the complete knowledge of the flow field together with analytical results valid near threshold to explore this further.

## II Large Scale Flows

In this work the terminology large scale flows is used to describe flows that extend over distances larger than that of the convection roll scale. We would like to distinguish between two different large scale flows: large scale counterflow and mean flow.

### II.1 Large Scale Counterflow

In the presence of a spatial ramp in plate separation a large scale counterflow is present for all values of the bulk control parameter $`ϵ_0`$, including $`ϵ_0<0`$. Warm fluid ascends the ramp eventually reaching the sidewall and is forced to flow back toward the center of the domain over the cold top wall causing it to descend resulting in a large zone of circulation in the vertical plane over the ramp. The magnitude of the large scale counterflow depends upon the specifics of the ramp and is roughly independent of $`ϵ_0`$ and $`\sigma `$ Walton (1982). Figure 2 illustrates this with a vertical slice from a three-dimensional numerical simulation where the entire convection layer is subcritical. As shown for this subcritical case the fluid motion of the large scale counterflow generates axisymmetric convection near the base of the ramp which extends a couple of roll widths toward the center of the domain. For the gradual ramp used in experiment the large scale counterflow is small in magnitude and has not been measured.

### II.2 Mean Flow

In discussing the mean flow it will be convenient to first present the governing heat and fluid equations. The velocity $`\stackrel{}{u}`$, temperature $`T`$, and pressure $`p`$, evolve according to the Boussinesq equations,

$`\sigma ^1\left(_t+\stackrel{}{u}\stackrel{}{}\stackrel{}{}\right)\stackrel{}{u}`$ $`=`$ $`\stackrel{}{}p+RT\widehat{z}+^2\stackrel{}{u},`$ (3)
$`\left(_t+\stackrel{}{u}\stackrel{}{}\stackrel{}{}\right)T`$ $`=`$ $`^2T,`$ (4)
$`\stackrel{}{}\stackrel{}{}\stackrel{}{u}`$ $`=`$ $`0,`$ (5)
where $`_t`$ indicates time differentiation, and $`\widehat{z}`$ is a unit vector in the vertical direction opposite to gravity. The equations are nondimensionalized in the standard manner using the layer depth $`d`$, the vertical diffusion time for heat $`\tau _vd^2/\kappa `$ where $`\kappa `$ is the thermal diffusivity, and $`\mathrm{\Delta }T`$ the temperature difference between the top and bottom surfaces, as the length, time, and temperature scales, respectively. The lower and upper surfaces are no-slip and are held at constant temperature. The sidewalls are also no-slip and, unless otherwise noted, are insulating.

The mean flow field, $`\stackrel{}{U}(x,y)`$, is the horizontal velocity integrated over the depth and originates from the Reynolds stress induced by pattern distortions. As illustrated by the fluid equations, Eqs. (3) and (5), it is evident that the pressure is not an independent dynamic variable. The pressure is determined implicitly to enforce incompressibility,

$$^2p=\sigma ^1\stackrel{}{}\stackrel{}{}\left[\left(\stackrel{}{u}\stackrel{}{}\stackrel{}{}\right)\stackrel{}{u}\right]+R_zT.$$
(6)
Focussing on the nonlinear Reynolds stress term and rewriting the pressure as $`p=p_o(x,y)+\overline{p}(x,y,z)`$ yields,

$$p_o(x,y)\sigma ^1𝑑x^{}𝑑y^{}\mathrm{ln}\left(1/\left|rr^{}\right|\right)\stackrel{}{}^{}\stackrel{}{}\left[\left(\stackrel{}{u}\stackrel{}{}\stackrel{}{}\right)\stackrel{}{u}\right]_z,$$
(7)
where $`_z`$ represents an average in the z direction. In Eq. (7) the $`\mathrm{ln}(1/|rr^{}|)`$ is not exact, in order to be more precise the finite system Green’s function would be required; however the long range behavior persists. This gives a contribution to the pressure that depends on distant parts of the convection pattern. The Poiseuille-like flow driven by this pressure field subtracts from the Reynolds stress induced flow leading to a divergence free horizontal flow that can be described in terms of a vertical vorticity.

Near threshold an explicit expression for the mean flow is Cross and Newell (1984)

$$\stackrel{}{U}(x,y)=\gamma \stackrel{}{k}\stackrel{}{}_{}\stackrel{}{}\left(\stackrel{}{k}A^2\right)\stackrel{}{}_{}p_o(x,y)$$
(8)
where $`\gamma `$ is a coupling constant given by $`\gamma =0.42\sigma ^1(\sigma +0.34)(\sigma +0.51)^1`$, $`A^2`$ is the convection amplitude normalized so that the convective heat flow per unit area relative to the conducted heat flow at $`R_c`$ is $`|A|^2R/R_c`$, $`p_o`$ is a slowly varying pressure, see Eq. (7), and $`\stackrel{}{}_{}`$ is the horizontal gradient operator \[see Siggia and Zippelius (1981); Newell et al. (1990) for the complete analysis and more details\]. The mean flow is important not because of its strength; under most conditions the magnitude of the mean flow is substantially smaller than the magnitude of the roll flow making it extremely difficult to quantify experimentally. The mean flow is important because it is a nonlocal effect acting over large distances (many roll widths) and changes important general predictions of the phase equation Cross and Newell (1984). The mean flow is driven by roll curvature, roll compression and gradients in the convection amplitude. The resulting mean flow advects the pattern giving an additional slow time dependence. It is important to note, that unlike the long range counterflow, the magnitude of the mean flow vanishes when the convection layer becomes critical, $`|\stackrel{}{U}|ϵ_0`$ for $`ϵ_01`$.

## III Numerical Simulation

We have performed full numerical simulations of the governing fluid and heat equations, Eqs. (3)-(5), in a cylindrical geometry with a radial ramp in plate separation using a parallel spectral element algorithm (described in detail elsewhere Fischer (1997), see Paul et al. (2001, submitted March 2002) for related applications).

For discussion purposes it will be convenient to define Cartesian $`(x,y)`$ and polar $`(r,\theta )`$ coordinates centered on a mid-depth horizontal slice of a cylindrical convection layer containing a field of straight parallel x-rolls with wavevector $`\stackrel{}{k}=k_o\widehat{x}`$ as shown in Fig. 3. The x-axis is perpendicular to the roll axes, the y-axis is parallel to the roll axes and $`\theta `$ measures the angle from the positive x-axis.

We have investigated the results found in experiment by performing simulations on a ramped cylindrical convection layer for a variety of scenarios and initial conditions. Simulations were performed over the range $`ϵ_00.2`$ and for simulation times of $`t_f>\tau _h`$ where $`\tau _hr_0^2`$ is the time required for heat to diffuse horizontally across the bulk region of the layer which has been suggested as the earliest time scale for the flow field to reach equilibrium Cross and Newell (1984).

The mean flow present in the simulation flow fields, $`\stackrel{}{U}_s(x,y)`$, is investigated by calculating the depth averaged horizontal velocity,

$$\stackrel{}{U}_s(x,y)=_0^1\stackrel{}{u}_{}𝑑z$$
(9)
where $`\stackrel{}{u}_{}`$ is the horizontal velocity field. Furthermore it will be convenient to work with the vorticity potential, $`\zeta `$, defined as

$$_{}^2\zeta =\widehat{z}\stackrel{}{}\left(_{}\times \stackrel{}{U}_s\right)=\omega _z$$
(10)
where $`\omega _z`$ is the vertical vorticity and $`_{}^2`$ is the horizontal Laplacian.

## IV Analytical Development

Near threshold, assuming straight parallel rolls, it is possible to approximately determine $`\omega _z`$ analytically. It will be convenient to start from $`\omega _z`$ given by the vertical component of the curl of Eq. (8),

$$\omega _z=\widehat{z}\stackrel{}{}\left(\stackrel{}{}_{}\times \stackrel{}{U}\right)=\gamma \widehat{z}\stackrel{}{}\stackrel{}{}_{}\times \left[\stackrel{}{k}\stackrel{}{}_{}\stackrel{}{}\left(\stackrel{}{k}|A|^2\right)\right].$$
(11)
Consider a cylindrical convection layer with a radial ramp in plate separation containing a field of x-rolls given by $`\stackrel{}{k}=k_o\widehat{x}`$. The amplitude can be represented for large $`ϵ_0`$, using an adiabatic approximation, as $`|A|^2=ϵ(r)/g_o`$ for $`ϵ>0`$ and $`|A|^2=0`$ for $`ϵ(r)<0`$, making the amplitude a function of radius only $`|A|^2=f(r)`$. This approximation is good except for the kink at $`r_c`$ where $`ϵ=0`$. Inserting $`|A|^2=f(r)`$ into Eq. (11) yields, after some manipulation, the following expression for the vertical vorticity,

$$\omega _z=\frac{\gamma k_{o}^{}{}_{}{}^{2}}{2}\left[\frac{d^2|A|^2}{dr^2}\frac{1}{r}\frac{d|A|^2}{dr}\right]\mathrm{sin}2\theta .$$
(12)
To correct for nonadiabaticity and to smooth $`|A(r)|^2`$ near $`r_c`$, the one-dimensional time independent amplitude equation is solved,

$$0=ϵ(r)A+\xi _{o}^{}{}_{}{}^{2}\mathrm{cos}^2\theta \frac{^2A}{r^2}g_o|A|^2A,$$
(13)
where angular derivatives have been assumed small, $`\xi _{o}^{}{}_{}{}^{2}=0.148`$, $`g_o=0.69950.0047\sigma ^1+0.0083\sigma ^2`$ and $`ϵ(r)`$ is determined by

$$ϵ(r)=\{\begin{array}{cc}ϵ_0,\hfill & r<r_0\hfill \\ ϵ_0(h^3h_c^3)/(1h_c^3),\hfill & rr_0\hfill \end{array}$$
(14)
where $`h_c=h(r_c)`$. In Eq. (11) $`\omega _z`$ is dominated by radial derivatives so Eq. (12) is still a good approximation for the case $`|A|^2=f(r,\theta )`$, $`d/dr/r`$. Due to the angular dependence of the nonadiabaticity $`|A|^2`$ is now $`\theta `$ dependent which will induce higher angular harmonics in $`\omega _z`$. We will neglect these higher harmonics, assume a $`\mathrm{sin}2\theta `$ dependence and approximately evaluate the magnitude using Eqs. (12) and (13) at $`\theta =\pi /4`$. Equation (13) can be rewritten in a more convenient form as,

$$0=\frac{ϵ(\overline{r})}{ϵ_0}\overline{A}+\frac{^2\overline{A}}{\overline{r}^2}\overline{A}^3$$
(15)
where $`r=\left(ϵ_{0}^{}{}_{}{}^{1/2}\xi _o\mathrm{cos}\theta \right)\overline{r}`$ and $`A=(ϵ_0/g_o)^{1/2}\overline{A}`$ Cross and Newell (1984) (since the amplitude goes to zero it can be shown that $`A`$ is now a real quantity Cross et al. (1983)). Equation (15) is solved numerically using the boundary conditions $`_{\overline{r}}\overline{A}=0`$ at $`\overline{r}=0`$, and $`\overline{A}=0`$ at $`\overline{r}=\overline{r}_1`$.

These analytical results are now used to investigate the vertical vorticity generation in a large radially ramped cylindrical convection layer. We start by looking at the configuration used in the recent experiments with input parameters: $`r_0=42.29`$, $`r_1=101.33`$, $`\delta _r=0.036`$ and $`\sigma =0.87`$.

When $`ϵ_{0}^{}{}_{}{}^{1/2}\xi _{o}^{}{}_{}{}^{1}(r_cr_0)1`$ the amplitude $`A^2(r)`$ is unable to adiabatically follow the ramp, this nonadiabaticity results in a considerable deviation from $`ϵ(r)/g_o`$ as shown in Fig. 4. However, when $`ϵ_{0}^{}{}_{}{}^{1/2}\xi _{o}^{}{}_{}{}^{1}(r_cr_0)1`$ the amplitude $`A^2(r)`$ follows $`ϵ(r)/g_o`$ adiabatically almost over the entire ramp except for the small kink at $`r_c`$ as shown in Fig. 5. The structure of $`\omega _z`$ depends upon this adiabaticity and is shown for various values of $`ϵ_0`$ in Fig. 6 where the $`\mathrm{sin}2\theta `$ dependence has been removed by choosing $`\theta =\pi /4`$.

If the $`\theta `$ dependence is included it is evident from Fig. 6 that the vertical vorticity has a quadrupole angular structure for small $`ϵ_0`$, i.e. four lobes of alternating positive and negative vorticity with one lobe per quadrant, and makes a transition to an octupole angular dependence for larger $`ϵ_0`$, octupole in the sense of an inner, $`rr_c`$, and outer, $`rr_c`$, quadrupole. In addition, since $`r_c/ϵ_0>0`$ there is a radial shift of the vorticity curves as $`ϵ_0`$ is increased.

In all cases the amplitude $`A^2(r)`$ decreases monotonically with $`r`$ and as a result $`r^1d|A|^2/dr0`$ thus generating only positive vorticity. However, the term $`d^2|A|^2/dr^2`$ can be of either sign and is responsible for the quadrupole and octupole angular structure in $`\omega _z`$.

As $`ϵ_0`$ approaches zero the nonadiabaticity of $`A^2(r)`$ increases until $`A^2(r)`$ exhibits a quadratic fall-off with $`r`$ for $`rr_0`$ resulting in $`\omega _z(r)0`$ for all $`r`$.

The mean flow generated by these vorticity distributions is determined by solving Eq. (11) with the boundary condition $`\zeta (r_1)=0`$. The vorticity potential is related to the mean flow in polar coordinates by $`(U_r,U_\theta )=(r^1_\theta \zeta ,_r\zeta )`$. The vorticity potential is expanded radially in second order Bessel functions while maintaining the $`\mathrm{sin}2\theta `$ angular dependence. Of particular interest is the mean flow perpendicular to the convection rolls, $`U_r(\theta =0)`$ or equivalently $`U_x(y=0)`$, which is shown in Fig. 7.

As expected, regions of negative and positive vorticity yield corresponding negative and positive values of the mean flow. As $`ϵ_0`$ vanishes $`U0`$ for all $`r`$ providing a mechanism for roll expansion in the bulk. For larger $`ϵ_0`$ the mean flow becomes larger in magnitude and increasingly negative for $`rr_0`$ providing a mechanism for roll compression.

To make the connection between mean flow and wavenumber quantitative it is noted that the wavenumber variation resulting from a mean flow across a field on x-rolls can be determined from the one-dimensional phase equation,

$$U_x\varphi =D_{}_{xx}\varphi $$
(16)
where the wavenumber is the gradient of the phase, $`k=_x\varphi `$, $`D_{}=\xi _o^2\tau _o^1`$, and $`\tau _o^1=19.65\sigma (\sigma +0.5117)^1`$ Cross and Hohenberg (1993). Assuming that for $`ϵ_01`$ the wavenumber is approximately $`kk_c`$ everywhere and that the rolls are exposed to a constant mean flow the wavenumber change over the bulk can be expressed as

$$\mathrm{\Delta }k=k(r_0)k(r=0)=Uk_cD_{}^1r_0.$$
(17)
For example, for curve (c) in Fig. 7 the maximum value of the mean flow is $`U=4.92\times 10^5`$ which yields a small roll expansion of $`\mathrm{\Delta }k=0.0035`$. If the mean flow were solely responsible for the dramatic roll expansion seen in experiment of $`\mathrm{\Delta }k1.0`$ ($`ϵ_0=0.012`$ see Fig. 3(a) of Bajaj et al. (1999)) a mean flow of $`U0.014`$ would be required, which is not found in the analytic results.

## V Discussion

The large scale flows discussed in Section II cannot be measured in current experiments placing us in a unique position to use the complete flow field information from our full numerical simulations of Eqs. (3)-(5) together with the analytical results of Section IV to investigate how the mean flow and the large scale counterflow induce wavenumber distortions and the variation of this distortion with Rayleigh number.

It is computationally expensive to perform full three-dimensional numerical simulations for the very large system used in experiment. We have, however, performed a variety of simulations for radially ramped cylindrical convection layers. The full three-dimensional simulations are of smaller spatial extent with the precise ramp defined by Eq. (2) and the specific input parameters: $`r_0=11.31`$, $`r_1=20`$, $`\delta _r=0.036`$, and $`\sigma =0.87`$. Two-dimensional simulations of a vertical slice of a three-dimensional domain (see Fig. 1) were also performed for both the large experimental configuration and the smaller computational domain just described. Three-dimensional simulations were also conducted without a large scale counterflow by a specific choice of ramp parameters that will be discussed below.

Initially we consider prescribed x-roll initial conditions given by $`\stackrel{}{k}=k_c\widehat{x}`$. Other initial conditions such as random thermal perturbations or initial x-rolls of varying wavenumbers were also investigated and found not to affect the final pattern wavenumber or any of the conclusions drawn. Simulations were performed for $`ϵ_0=0.025`$, 0.054, 0.113, and 0.171. Figure 8 compares the wavenumbers found in these simulations with recent experiments and will be discussed in detail below.

The final patterns in the simulations maintain the x-roll configuration imposed by the initial conditions. Figure 9 displays the final pattern observed for three-dimensional simulations with $`ϵ_0=0.025`$ in panel (a) and $`ϵ_0=0.171`$ in panel (b). Figure 9a illustrates that near threshold the convection rolls exhibit very little curvature indicating that the assumption of straight parallel x-rolls in Section IV is valid. There is more roll curvature apparent in Fig. 9b as would be expected for larger $`ϵ_0`$. Figure 9 also illustrates the decreasing size of the subcritical region as the supercriticality of the bulk increases. All simulations settled to a time independent state.

It is illustrative to compare the analytical results of Section IV with the results of simulation. Figure 10a displays $`A^2(r)`$ for the case $`ϵ_0=0.025`$, as determined by Eq. (15). A significant nonadiabaticity is present for this case as shown by the deviation of $`A^2(r)`$ from $`ϵ(r)/g_o`$. For the ramped domain used in simulation the distance $`r_c(ϵ_0)r_0`$ is smaller than in the larger domain with a more shallow ramp used in experiment. This results in the presence of more nonadiabaticity in the simulations when compared to experimental results at the same control parameter. This is beneficial because this allows the exploration of highly nonadiabatic situations without having to perform the task of simulating near the convective threshold, which becomes computationally difficult because of the diverging time scales.

A comparison between theory and simulation of the vertical vorticity and the resulting mean flow is shown in Fig. 10b and c. The theoretical predictions are based on the amplitude variation caused when straight parallel convection rolls encounter a radial ramp in plate separation as discussed in Section IV. For both the vertical vorticity and the mean flow the comparison is made in the absence of any adjustable parameters. For the vertical vorticity calculated in simulation an angular average, weighted by $`\mathrm{sin}2\theta `$, is used for the comparison. The agreement between theory and simulation is quite good. This illustrates quantitatively that the major source of vertical vorticity and mean flow is indeed the variation in the convective amplitude caused by the radial ramp in plate separation. Over the bulk of the domain the mean flow is negative and very small in magnitude with a maximum value of $`U_x(y=0)=7.0\times 10^4`$ and by Eq. (17) the wavenumber variation would be extremely small in agreement with the near constant bulk wavenumbers found in simulation.

A similar comparison between theory and simulation is made in Fig. 11 for $`ϵ_0=0.171`$. As shown in Fig. 11a $`A^2(r)`$ is much more able to follow the ramp, $`ϵ(r)/g_o`$, and exhibits very little nonadiabaticity except for the kink near $`r_c`$. This results in a much stronger negative vertical vorticity in the bulk which in turn yields a larger negative mean flow as shown if Fig. 11b and c. The agreement between theory and simulation for the vertical vorticity is still quite good. The discrepancy in the mean flow comparison may be due to the fact that as $`ϵ_0`$ increases other mean flow sources such as roll curvature, see Fig. 9b, become important.

Figure 12 illustrates the octupole structure in the vorticity potential in panel (a) and the roll compression occurring in the bulk by plotting contours of the local wavenumber in panel (b) for $`ϵ_0=0.171`$. As illustrated in panel (a) the mean flow has significant structure over the ramped region as well as extending into the subcritical region of the layer, $`r>r_c`$. It has also been suggested that the mean flow extends into a subcritical region in related experiments implementing “finned” boundaries Pocheau and Daviaud (1997).

The vorticity potential displays an octupole structure containing a pair of counter-rotating vortices in each quadrant. The inner quadrupole is localized around $`r_0`$ where gradients in the amplitude of convection occur as the ramp in plate separation begins. The direction of rotation of the inner quadrupole causes a focussing of the mean flow into the bulk region of the domain and is responsible for the larger wavenumbers found as $`ϵ_0`$ is increased as shown by the ($``$) curve in Fig. 8.

To make the connection between mean flow and wavenumber quantitative Eq. (16) is applied to the simulation results in the form,

$$U=D_{}k_c^1_xk.$$
(18)
Figure 13a illustrates the wavenumber variation, $`k(r)`$ for $`rr_0`$, found in simulation by simply measuring the distance between roll boundaries and makes evident the roll compression, $`k(r=0)>k(r_0)`$. Figure 13b compares the mean flow calculated from simulation with the predicted value of the mean flow required to produce the wavenumber variation shown in Fig. 13a using Eq. (18). The agreement is good and the discrepancy near $`r_0`$, which is contained within one roll wavelength from where the ramp begins, is expected because the influence of the ramp was not included in Eq. (16). This illustrates quantitatively that the mean flow compresses the rolls in the bulk of the domain.

As mentioned earlier, the mean flow vanishes as $`ϵ_0`$ approaches critical whereas the large scale counterflow is present for all $`ϵ_0`$ and therefore could play a role near threshold in the determination of the final convection pattern. In order to gain further insight into this possibility a radial ramp was constructed that did not drive a large scale counterflow. This was accomplished by setting the temperature of the ramped surface, $`T_b(r)`$, to the value of the linear conduction profile at that height, $`T_b(r)=h(r)`$. This ramp, therefore, does not bend the isotherms which is the source of the large scale counterflow. The wavenumber variation for these simulations, see curve labelled with ($``$) in Fig. 8, does not differ strongly from the simulations with a ramp producing large scale counterflow, see the curve labelled with ($``$). The similarity in wavenumber results is strongest for small $`ϵ_0`$ suggesting that the large scale counterflow is not responsible for the shift of the critical wavenumber to smaller values as seen in experiment.

To study the large scale counterflow further two-dimensional simulations were also performed, corresponding to a vertical slice of the domain considered thus far, and in addition to a more spatially extended domain as used in experiment. In two dimensions the mean flow is absent however the large scale counterflow persists. As shown by the ($``$) and ($`\mathrm{}`$) curves in Fig. 8 the wavenumbers measured in the two-dimensional simulations are not compressed to the same extent as $`ϵ_0`$ increases as in the three-dimensional simulations with both mean flow and large scale counterflow present. As expected, the wavenumbers found in the two-dimensional simulations are also independent of aspect ratio. Additionally, for small $`ϵ_0`$ the wavenumber found in simulation does not deviate markedly from its critical value suggesting that the large scale counterflow is not responsible for the wavenumber shift observed in experiment near threshold regardless of the spatial extent of the domain.

We also investigated the possibility that the vertical large scale counterflow could bifurcate into a horizontal flow similar to the mean flow in the presence of a slight spatial asymmetry. This was accomplished by giving the ramped domain used in the full three-dimensional simulations an eccentricity of $`e0.8`$ for a variety of ramps $`0.036\delta _r0.25`$ and simulating over a range of subcritical and supercritical conditions. For all of the scenarios tested the large scale counterflow remained vertical and did not undergo any significant changes.

Lastly, the possibility of wavenumber pinning was studied by varying the aspect ratio in increments of less than half of a roll width for both the two- and three-dimensional domains. In all of the scenarios tested the final pattern wavenumbers were not appreciably affected by these small changes in aspect ratio.

## VI Conclusion

We have analytically and numerically investigated pattern formation in a cylindrical convection layer with a radial ramp in plate separation. In particular, we have studied quantitatively the effects of two large scale flows; large scale counterflow and mean flow. These large scale flows are important theoretically yet are extremely difficult to measure experimentally.

Our results suggest that the mean flow plays an important role in the observed pattern wavenumber and is generated in a novel way, by the spatial variation of $`|A|^2`$ driven by the variation of $`ϵ`$ rather than the more usual variations in roll curvature and wavenumber. The mean flow sources are quantified analytically and agreement is found with numerical results. The geometric structure and magnitude of the mean flow is used to explain quantitatively the wavenumber variation found in the simulations.

The large scale counterflow is investigated numerically and our results indicate a small roll compression effect away from threshold. In particular, although the large scale counterflow is present at and near threshold it does not appear responsible for the dramatic wavenumber shift to values less than critical as seen in experiment.

Although it is too expensive computationally to simulate the very large systems used in the experiments we can use our quantitative understanding of the ramp-generated mean flow, validated by the simulations at smaller aspect ratio, to extrapolate our results to these larger systems. Furthermore our two-dimensional simulations in sizes equal to the experimental ones allow us to estimate the effect of the large scale counterflow on the wavenumber distribution. Despite these exhaustive efforts, we are unable to reproduce the large shift to smaller wavenumbers observed near threshold in the experiments, and the physical origin of these results remains a mystery.

## ACKNOWLEDGEMENTS

We are grateful to H. S. Greenside and G. Ahlers for helpful discussions. This research was supported by the U.S. Department of Energy, Grant DE-FT02-98ER14892, and the Mathematical, Information, and Computational Sciences Division subprogram of the Office of Advanced Scientific Computing Research, U.S. Department of Energy, under Contract W-31-109-Eng-38. We also acknowledge the Caltech Center for Advanced Computing Research and the North Carolina Supercomputing Center.


format: text

text:
# The Shape of Nereid’s Phase Function Changed From 1998 to 1999 and 2000

## 1 Introduction

Nereid has an unusual orbit around Neptune ($`e=0.75`$, $`P=360`$ days, and $`i=28\mathrm{°}`$) which suggests that it might be a captured Kuiper Belt Object (Schaefer & Schaefer, 1995, 2000). Nereid also has the most unusual known photometric history of all objects in the Solar System (Schaefer & Schaefer, 1988, 2000; Schaefer & Tourtellotte, 2001). From 1987 to 1991, several groups found that Nereid had fast photometric variations with amplitudes up to around one magnitude (Schaefer & Schaefer, 1988, 2000; Williams, Jones, & Taylor, 1991; Bus, & Larson, 1989). From 1991 to 1997, Nereid’s amplitude was smaller at $`0.4`$ mag (Schaefer & Schaefer, 2000), with highly significant evidence that it showed fast variations (Buratti, Goguen, & Mosher, 1997; Brown, & Webster, 1998). In 1998, we observed Nereid on 52 nights and found that Nereid displayed no significant variability other than a very large opposition surge (Schaefer & Tourtellotte, 2001).

This unique behavior of large-and-small variations that change from year to year could be caused by chaotic rotation (Dobrovolskis, 1995), much like Hyperion (Klavetter, 1989). This possibility is appealing because Nereid has the necessary ingredients of a $`>1\%`$ out-of-round shape (due to its small size) and a highly eccentric orbit around a planet with a large J2 component in its gravity. Unfortunately, Dobrovolskis has shown that Nereid’s rotation will be chaotic only if its rotational period is longer than roughly two weeks, and this requirement is not consistent with the observed brightness changes on time scales of one day. Thus, the cause of Nereid’s variations is not known (Schaefer & Tourtellotte, 2001).

The 1998 light curve (Schaefer & Tourtellotte, 2001, see Figure 1) displayed no significant brightness changes other than those associated with the normal changing of the solar phase angle, $`\alpha `$, and the opposition surge. The phase function (see Figure 2) shows an opposition surge of 0.52 magnitudes over a range of 2$`\mathrm{°}`$ in phase and this is among the largest known in our Solar System. The shape of the phase function is definitely not linear, while a broken line fits well.

Given Nereid’s long history of changes in its variability, it is prudent to continue photometric monitoring. The Yale 1-m telescope on Cerro Tololo is operated in a queue mode by a resident operator, and hence is perfect for the long-term synoptic study of Nereid. So again we have used this telescope to monitor Nereid, this time during its 1999 and 2000 oppositions. We find that Nereid’s phase function has significantly changed its shape from 1998 to 1999 and 2000.

## 2 Photometry

All photometry was taken with the Yale 1-m telescope (operated by the YALO consortium) with the ANDICAM CCD camera which has a $`2048\times 2048`$ array of $`0.29\mathrm{}`$ pixels. All exposures were through a standard V-band filter for 900 seconds of exposure. A complete description of our analysis procedures are presented in Schaefer & Tourtellotte (2001).

In 1999, we obtained a total of 54 usable images on 54 separate nights in a 130 day interval between 20 June 1999 (JD 2451350) and 28 October 1999 (JD 2451480). Nereid passed through opposition on 27 July 1999 (JD 2451386) and was at maximum phase angle on 26 April 1999 (JD 2451294) and 25 October 1999 (JD 2451476).

In 2000, we obtained a total of 14 usable images on 14 separate nights in a 102 day interval between 22 July 2000 (JD 2451748) and 1 November 2000 (JD 2451850). Nereid passed through opposition on 28 July 2000 (JD 2451753) and was at maximum phase angle on 28 April 2000 (JD 2451662) and 26 October 2000 (JD 2451843).

All our V-band photometry is presented in Table 1. The first column gives the heliocentric Julian Date of the start of the exposure. The second column gives solar phase angle, $`\alpha `$, in degrees as given by the JPL Horizons program (http://ssd.jpl.nasa.gov/cgi-bin/eph). The third column gives our observed V magnitude and one-sigma uncertainty. The fourth column gives the V magnitude corrected to the 1998 opposition distance of 29.123 AU as $`V_{OPP}=V5\mathrm{log}(\mathrm{\Delta }/29.123)`$ where $`\mathrm{\Delta }`$ is the Earth-Nereid distance in units of AU, to allow direct comparisons with the 1998 phase function. The error bars for $`V_{OPP}`$ are the same as for $`V`$ in the previous column. (The total range of variation of the Sun-Nereid distance over the observed time intervals in 1998-2000 was only 0.08 AU, thus resulting in a negligible correction of under 0.003 mag from the mean.) The light curves are displayed in Figure 1.

We notice that the one measurement with $`\alpha <0.11\mathrm{°}`$ in 1999 (on JD 2451386.829) is bright by $``$ 0.2 mag when compared to observations taken a few days before and afterwards. We have carefully examined this image for bad columns, measured Nereid’s image shape for evidence of cosmic rays, checked other deeper images for background stars, and checked Nereid’s position on the CCD chip with other images for hot pixels, all with the result that we have no grounds for impeaching this ’bright’ magnitude. (Actually, similar checks were performed on all our Nereid images.) The coincidence of this bright Nereid with phase angle zero suggests that this might be evidence of a narrow spike in Nereid’s phase function (much like that for four of the Uranian moons (Buratti, Gibson, & Mosher, 1992). Yet such a spike does not appear in the well-sampled 1998 light curve, so the ’spike’ would have to be time variable if it exists.

## 3 Phase Function

In 1999, the observed light curve comes to a peak around 26 July 1999 and a minimum in October 1999 with a symmetric light curve around the peak. In 2000, the observed light curve comes to a peak around 26 July 2000 and a minimum in October 2000. The coincidence of the observed maximum and minimum dates with the opposition and quadrature dates as well as the symmetric light curves are the hallmarks of brightness variations caused by the opposition surge. Thus, it appears that (as in 1998), most of the variations arises from simple phase effects. The phase functions for 1999 and 2000 are presented in Figure 2.

We first tried to fit the 1999 and 2000 data to linear and broken-line functional forms. For the 1999 data, the broken-line form is strongly preferred over the linear form. The best fit is $`V_{OPP}=18.97+1.09\alpha `$ for $`\alpha <0.29\mathrm{°}`$ and $`V_{OPP}=19.25+0.13\alpha `$ for $`\alpha >0.29\mathrm{°}`$ with a chi-square of 40.3 for 50 degrees of freedom. The 1999 data is strongly inconsistent with the 1998 broken-line fit (Schaefer & Tourtellotte, 2001). The 1999 and 2000 data have nearly identical fits.

Fits to lines and broken-lines have the virtue of simplicity, but they have no physical motivation. Hapke (1993) developed a physical model that predicts the opposition surge component of the phase function to vary as $`V_{OPP}=V_02.5\mathrm{log}\{[1+B_0/(1+0.5\alpha /h)]/[1+B_0]\}`$ in the small angle limit. Here, $`V_0`$ is the V-band magnitude at zero phase angle, $`B_0`$ is a measure of the amplitude of the opposition surge, and $`h`$ is the angular width of the opposition surge. The use of this physical model has the disadvantage that it might not be a complete description of the physics (see discussion in the next section). Also, the physical model does not describe the observed variations over small ranges of phase angles, in particular at low phase angles.

The best fit for the 1999 data has $`V_0=19.01`$ mag, $`B_0=0.72`$, and $`h=0.19\mathrm{°}`$ for a chi-square of 48.8 with 51 degrees of freedom. The best fit for the 2000 data has $`V_0=19.07`$ mag, $`B_0=0.78`$, and $`h=0.34\mathrm{°}`$ for a chi-square of 9.6 with 11 degrees of freedom, although the 2000 errors are relatively large and correlated. The fit parameters are similar for 1999 and 2000. The best way to show this is by evaluating the chi-squares for each year with the other year’s fit. The 1999 data have a chi-square of 52.1 for the 2000 best fit; while the 2000 data have a chi-square of 16.4 for the 1999 best fit. The reduced chi-square values are near unity and F-Tests show that the cross-fits are not significantly different from each other. A joint fit of the 1999+2000 data reveal a best fit with $`V_0=19.04`$ mag, $`B_0=0.73`$, and $`h=0.25\mathrm{°}`$ for a chi-square of $`49.4+12.0=61.4`$ with 65 degrees of freedom. Again, the differences in chi-square are not significant. Hence, it appears that the phase functions in 1999 and 2000 were similar if not identical. So we will adopt the best fit to the Hapke model for the joint 1999+2000 data as a good description of the phase function for both years.

Is the 1998 phase function (Schaefer & Tourtellotte, 2001) the same as what we are here reporting for 1999 and 2000? The 1998 data were fitted to the Hapke model with $`V_0=19.07`$ mag, $`B_0=1.32`$, and $`h=0.52\mathrm{°}`$ for a chi-square of 70.0 with 54 degrees of freedom. Both the 1998 and the 1999+2000 model phase functions are drawn in all three panels of Figure 2. In examining this figure, the answer is obvious that the observed phase function has changed from 1998 to 1999+2000. In particular, the observed 1999+2000 phase function appears to be more peaked at low phases and substantially brighter at large phases when compared to the observed 1998 phase function.

The 1998 phase function and light curve is flat at low phase angles, while the 1999+2000 phase function and light curve is strongly peaked to opposition. To be quantitative, in bins of 0.05$`\mathrm{°}`$ phase from 0.00$`\mathrm{°}`$ to 0.30$`\mathrm{°}`$, the 1998 observed values are $`19.143\pm 0.026`$, $`19.137\pm 0.021`$, $`19.105\pm 0.020`$, $`19.14\pm 0.05`$, $`19.118\pm 0.027`$, and $`19.199\pm 0.028`$ mag. Similarly the 1999+2000 binned phase function is $`19.064\pm 0.022`$ (or $`19.087\pm 0.024`$ if the point nearest opposition is arbitrarily ignored), $`19.12\pm 0.03`$, $`19.155\pm 0.019`$, $`19.180\pm 0.018`$, $`19.173\pm 0.026`$, and $`19.253\pm 0.032`$ mag. All 1998 points with $`0.00\mathrm{°}<\alpha <0.26\mathrm{°}`$ are consistent with a perfectly flat phase curve at $`19.126\pm 0.011`$ mag. The 1999+2000 magnitudes form a simple linear decline that is not consistent with a flat phase function. Formal fits for the slopes yield $`0.12\pm 0.12`$ mag deg-1 for 1998 and $`0.73\pm 0.14`$ mag deg-1 for 1999+2000 ($`0.63\pm 0.14`$ mag deg-1 if the near opposition datum is excluded despite the lack of grounds).

At large phases with $`\alpha >0.85\mathrm{°}`$, all 1998 points lie below the 1999+2000 fit and all 1999+2000 points lie above the 1998 fit. (Both statements have a single exception which in both cases is the one point with the largest error bar.) Taking points with $`\alpha >1.5\mathrm{°}`$, the weighted averages for 1998 and 1999+2000 are $`19.551\pm 0.013`$ and $`19.485\pm 0.010`$ mag respectively, for a difference of $`0.066\pm 0.016`$ mag which is significant at the 4.0 sigma level. Hence, the observed phase function at large phases changed significantly after 1998.

Overall, the change in the shape of the phase function can be tested by comparing the fits to the Hapke model. The 1998 phase function when applied to the 1999+2000 data yields a chi-square of $`119.3+19.9=139.2`$ for 68 degrees of freedom (versus a chi-square of 61.4 for the 1999+2000 phase function), and this strongly rejects the 1998 function for the 1999+2000 data. Similarly, the 1999+2000 phase function when applied to the 1998 data yields a chi-square of 120.7 for 57 degrees of freedom (versus a chi-square of 70.0 for the 1998 phase function), and this strongly rejects the 1999+2000 phase function for the 1998 data. Thus, Nereid’s phase function changed from 1998 to 1999+2000.

This is an unprecedented result, so it is prudent to consider whether there is any realistic possibility of observational error. For example, what if we adopted magnitudes for the 1998 comparison stars used at high phase angles that are too faint by 0.07 mag due to some calibration error? (This would still not resolve the difference in the slope at low phases, so a separate error would have to be postulated.) A 0.07 mag error is far larger than any reasonable statistical error, and our repeated independent derivation of all the calibrations argues against numerical errors. The same telescope, CCD chip, filters, analysis procedures, and many of the same Landolt standard stars were used throughout each of the three years. The presence of undetected light cirrus clouds that might dim the comparison star images in comparison to the standard star images is ruled out since the $`\alpha >1.3\mathrm{°}`$ comparison star data in 1998 was independently calibrated against Landolt standard stars on three separate nights and these all agree to within $``$ 0.015 mag. Similarly, the comparison stars for the 1999 data with $`\alpha >1.3\mathrm{°}`$ were independently calibrated on three nights with the YALO telescope and one night with the McDonald 2.1-m telescope. Nereid crossed its own path on nights around 20 June 1998 and 20 October 1999, for which our independent calibrations of the comparison stars are in agreement to within $`0.03\pm 0.03`$ mag. In all, we have strong reasons to be confident in our photometry and thus we conclude that Nereid’s phase function did indeed change from 1998 to 1999+2000.

Nereid has a long history of photometric variations, so could the phase function changes actually be small amplitude long-term changes? Conceptually, it is difficult to know how the two cases can be distinguished, especially when we do not know the cause of the variability. Certainly we could postulate changes in Nereid’s overall brightness (but not its surface texture) that happen to make it faint in September and October 1998 (or bright in September and October of both 1999 and 2000) as well as fade Nereid by $``$ 0.2 mag in the week centered on the 1998 opposition (or brightened Nereid in the week centered on both the 1999 and 2000 oppositions). But such a postulate would require two rather improbable coincidences. The first coincidence is that the fading in 1998 is centered on the opposition date, since this time is only special from an Earthly viewpoint. (At this time, Nereid is far from Neptune and its 360 day orbital period means that it will be at roughly the same orbital phase for the oppositions from 1998-2000.) The second coincidence is that Nereid’s brightness always returns to the same magnitude the five times when $`0.4\mathrm{°}<\alpha <0.8\mathrm{°}`$ (early July 1998, middle August 1998, early July 1999, middle August 1999, and middle August 2000). The point to these coincidences is that Nereid’s brightness appears to be correlated with the position of the Earth. This correlation only makes sense if the observed variations are caused by a changing phase function which is explicitly a function of the Earth’s position. So again, we are concluding that Nereid’s phase function changed from 1998 to 1999+2000.

## 4 Implications

The best fit for the 1998 data has $`B_0`$ significantly larger than unity, whereas the model of Hapke requires $`B_01`$. This is actually a fairly broad problem, as many objects throughout the Solar System have this same difficulty (Helfenstein, Veverka, & Hillier, 1997). The likely solution is that the true situation is more complex than given in the basic model, perhaps involving a combination of coherent-backscattering plus strongly backscattering single-particle phase functions.

The narrow width of the opposition surge for Nereid has precedents from four Uranian moons (Miranda, Ariel, Titania, and Oberon) with high albedos (Buratti, Gibson, & Mosher, 1992) as well as for Europa (Doningue et al., 1991; Helfenstein et al., 1998). In general, simple shadow-hiding cannot account for such narrow widths without presuming a very low filling factor for objects in the regolith, with $`1\%`$ being a typical value (cf. Hapke 1993, equations 8.84 and 8.85). Such a low value seems implausible, although Hapke (1983) has proposed that scatterers might be concentrated towards the surface of a relatively deep and transparent medium. Nevertheless, the narrow component in the opposition surge is generally thought to arise from coherent-backscattering (Hapke, 1993). In this case, $`h\lambda /(4\pi L)`$ where $`\lambda `$ is the wavelength of the light and $`L`$ is the photon mean free path in the medium (Hapke 1993, equation 8.88). For visible light ($`\lambda =0.55\mathrm{\mu m}`$) and $`h=0.0043`$ rad ($`0.25\mathrm{°}`$), $`L=10.2\mathrm{\mu m}`$.

We interpret the changes displayed in Figure 2 as due to a variable phase function, but other explanations are possible. It is always possible to postulate subtle observational errors that can explain systematic offsets of $``$ 0.1 mag. Nevertheless, in our case, the stability of our equipment and procedures along with our many independent calibrations and cross checks gives us strong confidence that our results are reliable. It is always possible to postulate some arbitrary fading and brightening caused by whatever mechanism has previously made Nereid a variable, such that the light curves in Figure 1 are produced without changing the surface texture. Nevertheless, such postulated changes appear to be correlated with the Earth’s position due to the low-phase differences being centered on opposition and due to the return to the same magnitude on five occasions when the phase is $`0.4\mathrm{°}<\alpha <0.8\mathrm{°}`$. This connection with the Earth argues that the light curve changes arise from phase function changes. These arguments give us the grounds to identify the light curve changes as arising from changes in the phase function of Nereid.

A central mystery of Nereid remains to explain its large-and-small variations, only now there is the associated mystery of how it can change its phase function on a time scale of one year. Presumably, as the phase function depends on the texture of the surface, then Nereid’s texture must have changed. Could this have been caused by the deposition of particles or the formation of ice crystals on the surface? If so, then what could have caused this surface change? How is this resurfacing mechanism related to the brightness changes seen by so many groups? Or could chaotic rotation of Nereid have moved a different hemisphere into view, with each hemisphere having different textures? Unfortunately, we do not have answers to any of these questions.

We know of no precedents for significant changes in the shape of the phase function of Solar System objects. But this statement is not as strong as might be expected, since few objects in our Solar System have been examined with well-sampled phase functions over more than one opposition with measures that can be directly compared. Likely, Nereid’s variable phase function is related somehow to its previously reported variations, but the causes of both types of changes are unknown. To solve this mystery, a series of well-sampled light curves is certainly needed, ideally in several different bands and accompanied by spectroscopy.


text:
# Oases in the Desert: Three New Proposals

## 1 Introduction

The title of the talk I gave at Beyond the Desert 2002 was “New Unexpected Gauge Extensions of the Standard Model”. At the request of the editor of these Proceedings, I am including some very recent work on the prediction of the Majorana neutrino mass matrix, hence the present title.

I start with a brief review of the symmetries of the Standard Model and previous (trivial) gauge extensions. I then discuss two newly discovered nontrivial U(1) gauge extensions , based on the particle content of the Standard Model and its supersymmetric generalization respectively. Finally I switch gear and present a new natural understanding of the Majorana neutrino mass matrix, which automatically gives the observed pattern of solar  and atmospheric  neutrino oscillations and a prediction of neutrinoless double beta decay at the level of 0.4 eV .

## 2 Anomaly Cancellation in the Standard Model

The gauge group of the Standard Model is $`SU(3)_C\times SU(2)_L\times U(1)_Y`$, under which each family of quarks and leptons transforms as follows:

$`(u,d)_L(3,2,1/6),u_R(3,1,2/3),d_R(3,1,1/3);`$ (1)
$`(\nu ,e)_L(1,2,1/2),e_R(1,1,1).`$ (2)
The axial-vector triangle anomaly  is absent because

$`[SU(3)]^2Y:\left({\displaystyle \frac{1}{2}}\right)\left[2\left({\displaystyle \frac{1}{6}}\right){\displaystyle \frac{2}{3}}\left({\displaystyle \frac{1}{3}}\right)\right]=0,`$ (3)
$`[SU(2)]^2Y:\left({\displaystyle \frac{1}{2}}\right)\left[3\left({\displaystyle \frac{1}{6}}\right)+\left({\displaystyle \frac{1}{2}}\right)\right]=0,`$ (4)
$`Y^3:6\left({\displaystyle \frac{1}{6}}\right)^33\left({\displaystyle \frac{2}{3}}\right)^33\left({\displaystyle \frac{1}{3}}\right)^3+2\left({\displaystyle \frac{1}{2}}\right)^3(1)^3=0.`$ (5)
Note the nontrivial cancellation between quarks and leptons in the last two equations. The original 1967 model of leptons by Weinberg  was thus anomalous.

## 3 Automatic Symmetries of the Standard Model

Given the minimal particle content of the Standard Model, there are four automatic global symmetries: baryon number $`B`$ and the 3 lepton numbers $`L_e`$, $`L_\mu `$, $`L_\tau `$. Each is anomalous but the combination $`a_BB+a_eL_e+a_\mu L_\mu +a_\tau L_\tau `$ is not for the following cases.

(1) If there is no $`N_R(1,1,0)`$, then the only solution is $`a_B=0`$ and $`(a_e,a_\mu ,a_\tau )=(1,1,0)`$, $`(1,0,1)`$, or $`(0,1,1)`$. This allows $`L_iL_j`$ to be gauged .

(2) If only one family has an additional $`N_R`$, then the only solution is $`a_B=1`$ and $`(a_e,a_\mu ,a_\tau )=(3,0,0)`$, $`(0,3,0)`$, or $`(0,0,3)`$. This allows $`B3L_i`$ to be gauged .

(3) If two $`N_R`$’s are added, then $`a_B=1`$ and $`a_e=0`$, $`a_\mu +a_\tau =3`$, or $`a_\mu =0`$, $`a_e+a_\tau =3`$, or $`a_\tau =0`$, $`a_e+a_\mu =3`$ is a solution. This allows for example $`B(3/2)(L_\mu +L_\tau )`$ to be gauged .

(4) If there are three $`N_R`$’s, then $`a_B=1`$ and $`a_e+a_\mu +a_\tau =3`$ is a solution. The well-known example of $`a_e=a_\mu =a_\tau =1`$ allows $`BL`$ to be gauged . Another solution is $`a_B=0`$ and $`a_e+a_\mu +a_\tau =0`$. This means for example that $`2L_eL_\mu L_\tau `$ may be gauged.

## 4 Neutrino Mass and a New U(1) Gauge Symmetry

The addition of $`N_R`$ allows the doublet neutrinos $`\nu _i`$ to acquire small Majorana masses via the famous seesaw mechanism. However, if $`N_R`$ is replaced by a heavy fermion triplet

$$(\mathrm{\Sigma }^+,\mathrm{\Sigma }^0,\mathrm{\Sigma }^{})_R(1,3,0),$$
(6)
neutrinos get seesaw masses just as effectively .

Since the addition of one $`N_R`$ per family leads to the well-known U(1) gauge symmetry $`BL`$, the same question may be raised as to the addition of one $`\mathrm{\Sigma }_R`$ triplet per family. It has recently been shown  that indeed such an $`U(1)_X`$ exists, under which

$`(u,d)_Ln_1,u_Rn_2={\displaystyle \frac{1}{4}}(7n_13n_4),d_Rn_3={\displaystyle \frac{1}{4}}(n_1+3n_4),`$ (7)
$`(\nu ,e)_Ln_4,e_Rn_5={\displaystyle \frac{1}{4}}(9n_1+5n_4),\mathrm{\Sigma }_Rn_6={\displaystyle \frac{1}{4}}(3n_1+n_4).`$ (8)
Note that this does not correspond to any combination of known quantum-number assignments, such as $`Q`$, $`Y`$, $`B`$ or $`L`$.

To show that $`U(1)_X`$ has no anomalies, consider first

$`[SU(3)]^2X:2n_1n_2n_3=0,\mathrm{and}Y^2X:`$ (9)
$`6\left({\displaystyle \frac{1}{6}}\right)^2n_13\left({\displaystyle \frac{2}{3}}\right)^2n_23\left({\displaystyle \frac{1}{3}}\right)^2n_3+2\left({\displaystyle \frac{1}{2}}\right)^2n_4(1)^2n_5=0,`$ (10)
$`\mathrm{and}[SU(2)]^2X:\left({\displaystyle \frac{1}{2}}\right)(3n_1+n_4)(2)n_6=0.`$ (11)
These imply

$$n_3=2n_1n_2,n_5=\frac{1}{2}n_1n_2+\frac{1}{2}n_4,n_6=\frac{1}{4}(3n_1+n_4).$$
(12)
Consider next $`YX^2:`$

$$6\left(\frac{1}{6}\right)n_1^23\left(\frac{2}{3}\right)n_2^23\left(\frac{1}{3}\right)n_3^2+2\left(\frac{1}{2}\right)n_4^2(1)n_5^2=0.$$
(13)
Using Eq. (12), this implies

$$\frac{1}{4}(3n_1+n_4)(7n_14n_23n_4)=0.$$
(14)
If $`3n_1+n_4=0`$, then $`n_6=0`$ as well, and $`U(1)_X`$ is proportional to $`U(1)_Y`$, i.e. no new gauge symmetry has been discovered. On the other hand, if $`7n_14n_23n_4=0`$ is chosen instead, the solution given by Eqs. (7) and (8) is obtained. Nevertheless, there is still one more condition to be checked, i.e. the sum of the cubes of all $`U(1)_X`$ charges. Remarkably,

$$X^3:6n_1^33n_2^33n_3^3+2n_4^3n_5^33n_6^3=0$$
(15)
automatically. Furthermore, the sum of all $`U(1)_X`$ charges themselves,

$$X:6n_13n_23n_3+2n_4n_53n_6=0$$
(16)
as well. Hence the mixed gravitational-gauge anomaly  is also absent automatically. These are highly nontrivial results.

The Higgs sector of this gauge extension requires two doublets, one with $`U(1)_X`$ charge $`(9n_1n_4)/4`$ to give mass to the charged leptons, and the other with charge $`3(n_1n_4)/4`$ for the other fermions.

In general, for the fermion multiplet $`(1,2p+1,0;n_6)_R`$, there are 3 conditions to be satisfied, i.e. the analogs of Eqs. (11), (15) and (16):

$`{\displaystyle \frac{1}{2}}(3n_1+n_4)={\displaystyle \frac{1}{3}}p(p+1)(2p+1)n_6,`$ (17)
$`6n_1^33n_2^33n_3^3+2n_4^3n_5^3=(2p+1)n_6^3,`$ (18)
$`6n_13n_23n_3+2n_4n_5=(2p+1)n_6.`$ (19)
For $`p0`$ and $`3n_1+n_40`$, these imply

$$\frac{4n_6}{3n_1+n_4}=\frac{6}{p(p+1)(2p+1)}=\frac{3}{2p+1}=\left(\frac{3}{2p+1}\right)^{1/3}.$$
(20)
This determines $`p=1`$ (or $`p=2`$ which is the same as $`p=1`$ with $`n_6n_6`$ and $`\mathrm{\Sigma }_R\mathrm{\Sigma }_L`$.) In other words, the solution I have found is nontrivial and unique.

## 5 NuTeV Discrepancy

Since the gauge boson $`X`$ couples to quarks and leptons according to Eqs. (7) and (8), it may have measurable effects in precision electroweak data, such as the deep inelastic scattering of $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ on nucleons in the NuTeV experiment . A discrepancy has been reported in the value of the effective $`\mathrm{sin}^2\theta _W`$, i.e. $`0.2277\pm 0.0013\pm 0.0009`$ versus the expected $`0.2227\pm 0.00037`$ of the Standard Model. To see how this may be explained in the $`U(1)_X`$ model, consider the effective Hamiltonian of neutrino interactions:

$$_{int}=\frac{G_F}{\sqrt{2}}\overline{\nu }\gamma ^\mu (1\gamma _5)\nu [ϵ_L^q\overline{q}\gamma _\mu (1\gamma _5)q+ϵ_R^q\overline{q}\gamma _\mu (1+\gamma _5)q].$$
(21)
Assume no $`ZX`$ mixing (so that the precision $`Z`$-pole measurements are not affected), then

$`ϵ_L^u={\displaystyle \frac{1}{2}}{\displaystyle \frac{2}{3}}\mathrm{sin}^2\theta _W+n_1\zeta ,`$ (22)
$`ϵ_L^d={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{3}}\mathrm{sin}^2\theta _W+n_1\zeta ,`$ (23)
$`ϵ_R^u={\displaystyle \frac{2}{3}}\mathrm{sin}^2\theta _W+n_2\zeta ,`$ (24)
$`ϵ_R^d={\displaystyle \frac{1}{3}}\mathrm{sin}^2\theta _W+n_3\zeta ,`$ (25)
where

$$\zeta =2n_4\left(\frac{g_X^2}{M_X^2}\right)\left(\frac{M_Z^2}{g_Z^2}\right).$$
(26)
The NuTeV resluts versus the Standard Model predictions are:

$`(ϵ_L^u)^2+(ϵ_L^d)^2=0.3005\pm 0.0014\mathrm{versus}0.3042,`$ (27)
$`(ϵ_R^u)^2+(ϵ_R^d)^2=0.0310\pm 0.0011\mathrm{versus}0.0301.`$ (28)
Consider $`n_1=1`$, $`n_4=4/3`$, then $`n_2=3/4`$ and $`n_3=5/4`$. The discrepancies in the above are thus

$`\mathrm{\Delta }_L={\displaystyle \frac{2}{3}}\mathrm{sin}^2\theta _W\zeta +2\zeta ^2,`$ (29)
$`\mathrm{\Delta }_R={\displaystyle \frac{1}{6}}\mathrm{sin}^2\theta _W\zeta +{\displaystyle \frac{17}{8}}\zeta ^2.`$ (30)
Let $`\zeta =\mathrm{sin}^2\theta _W/6`$, then a very good fit is obtained, i.e.

$`\mathrm{\Delta }_L=0.0028\mathrm{versus}0.0037\pm 0.0014,`$ (31)
$`\mathrm{\Delta }_R=+0.0016\mathrm{versus}+0.0009\pm 0.0011.`$ (32)
This choice also implies that $`G_X/G_F=\mathrm{sin}^2\theta _W/16=0.014`$. Note that this solution assumes that $`X`$ couples to $`\mu `$ and the $`u`$ and $`d`$ quarks. If the same couplings were used for the electron, then the constraints from atomic parity violation would be grossly violated.

## 6 New U(1) Gauge Extension of the Supersymmetric Standard Model

In extending the Minimal Standard Model to include supersymmetry, one old problem and two new problems have to be faced. The old problem is neutrino mass. This may be solved again by adding heavy singlet neutral superfields $`N^c`$ (analog of $`N_R`$). The two new problems are rapid proton decay and the value of the $`\mu `$ term. It has been shown recently  that a nontrivial U(1) gauge extension exists which cures all three problems.

Consider the following superfields with their $`SU(3)_C\times SU(2)_L\times U(1)_Y\times U(1)_X`$ assignments:

$`Q=(u,d)(3,2,1/6;n_1),u^c(3^{},1,2/3;n_2),`$ (33)
$`d^c(3^{},1,1/3;n_3),L=(\nu ,e)(1,2,1/2;n_4),`$ (34)
$`e^c(1,1,1;n_5),N^c(1,1,0;n_6),`$ (35)
$`\varphi _1(1,2,1/2;n_1n_3),\varphi _2(1,2,1/2;n_1n_2).`$ (36)
Without $`U(1)_X`$, the trilinear terms $`LLe^c`$, $`LQd^c`$, and $`u^cd^cd^c`$ are allowed in the superpotential, thus causing rapid proton decay. The usual solution is to impose $`R`$ parity, i.e. $`R(1)^{3B+L+2j}`$, to forbid these terms. Even so, dimension-5 terms such as $`\stackrel{~}{q}\stackrel{~}{q}ql`$ are still allowed and may well be too big. Secondly, the term $`\mu \varphi _1\varphi _2`$ is allowed by supersymmetry, so there is no understanding of why $`\mu `$ should be of the order 1 TeV and not some very large mass such as the Planck scale or the string scale.

With $`U(1)_X`$, the $`\mu `$ problem is solved by replacing it by a singlet neutral superfield

$$\chi (1,1,0;2n_1+n_2+n_3),$$
(37)
where $`2n_1+n_2+n_30`$. The subsequent spontaneous breaking of $`U(1)_X`$ at the TeV scale generates the desirable value of the effective $`\mu `$ parameter. The solution of the proton decay problem now depends on finding an anomaly-free set of $`n_i`$’s which also forbids the undesirable terms mentioned above. With the superfield content as it is, there is no solution. However, if 2 sets of

$$U(3,1,2/3;n_7),U^c(3^{},1,2/3;n_8),$$
(38)
and 1 set of

$$D(3,1,1/3;n_7),D^c(3^{},1,1/3;n_8)$$
(39)
are added, a solution again appears. This will involve the remarkable exact factorization of the sum of 11 cubic terms as shown below.

There are 8 $`n_i`$’s, but 3 conditions are imposed:

$$n_1+n_3=n_4+n_5,n_1+n_2=n_4+n_6,n_7+n_8=2n_1n_2n_3.$$
(40)
These allow $`\varphi _1`$ to give mass to the $`d`$ quarks and the charged leptons, $`\varphi _2`$ to give mass to the $`u`$ quarks and the neutrinos (i.e. $`\nu N^c`$), and $`\chi `$ to give mass to the $`U`$ and $`D`$ quarks. Consider now the anomaly-free conditions linear in $`X`$:

$`X[SU(3)]^2:2n_1+n_2+n_3+n_7+n_8=0,`$ (41)
$`X[SU(2)]^2:n_2+n_3=7n_1+3n_4,`$ (42)
$`X[Y]^2:n_2+n_3=7n_1+3n_4.`$ (43)
The first condition is automatically satisfied, while the other two conditions are identical and eliminate just one additional $`n_i`$. Let the remaining 4 independent $`n_i`$’s be $`n_1`$, $`n_2`$, $`n_4`$, and $`n_7`$, then $`X^2[Y]`$ implies

$`3n_1^22n_2^2+n_3^2n_4^2+n_5^2+3n_7^23n_8^2(n_1+n_3)^2+(n_1+n_2)^2`$
$`=6(3n_1+n_4)(2n_14n_23n_7)=0.`$ (44)
The condition $`3n_1+n_4=0`$ contradicts $`2n_1+n_2+n_30`$ and must be discarded. Thus $`2n_14n_23n_7=0`$ is required and only 3 independent $`n_i`$’s are left, which are chosen finally to be $`n_1`$, $`n_4`$, and $`n_6`$. The most nontrivial condition, i.e. $`X^3`$, is then

$`18n_1^3+9n_2^3+9n_3^3+6n_4^3+3n_5^3+3n_6^3+9n_7^3+9n_8^3`$
$`2(n_1+n_3)^32(n_1+n_2)^3+(2n_1+n_2+n_3)^3=0.`$ (45)
Amazingly, this exactly factorizes to

$$36(3n_1+n_4)(9n_1+n_42n_6)(6n_1n_4n_6)=0.$$
(46)
Whereas the first factor cannot be zero, either of the other two can be chosen to be zero, and two solutions have been found. They are summarized in Table 1.

In Table 2 the various conditions of forbidding $`B`$ and $`L`$ violation, etc. in the superpotential are listed. The condition $`3n_1+n_40`$ which forbids the $`\mu `$ term also forbids the higher-dimensional terms $`QQQL`$ and $`u^cu^cd^ce^c`$, which are allowed by $`R`$ parity. Thus proton decay is much more effectively suppressed in this case. Note also that solutions (A) and (B) are the same for $`n_1=n_4`$.

## 7 Naturally Small Dirac Neutrino Masses

The sum of $`X`$ charges is

$`3(6n_1+3n_2+3n_3+2n_4+n_5+n_6)+3(3n_7+3n_8)+2(n_1n_3)`$
$`+2(n_1n_2)+(2n_1+n_2+n_3)=6(3n_1+n_4)0.`$ (47)
To get rid of this mixed gravitational-gauge anomaly, add the following singlet superfields in units of $`3n_1+n_4`$: one with charge 3, three ($`S^c`$) with charge $`2`$, and three ($`N`$) with charge $`1`$. Then since

$$1(3)+3(2)+3(1)=6,1(3)^3+3(2)^3+3(1)^3=0,$$
(48)
the mixed gravitational-gauge anomaly is canceled without affecting Eq. (45).

If now $`n_6=3n_1+n_4`$ is assumed, then there are 3 copies each of $`N^c(n_6)`$, $`N(n_6)`$, and $`S^c(2n_6)`$. Thus $`NN^c`$ forms a large invariant mass $`M`$, $`NS^c\chi `$ implies a mass $`m_2`$ proportional to $`\chi `$, and $`\nu N^c\varphi _2^0`$ implies a mass $`m_1`$ proportional to $`\varphi _2^0`$. The $`12\times 12`$ mass matrix spanning $`(\nu ,S^c,N,N^c)`$ is then of the form

$$=\left(\begin{array}{cccc}0& 0& 0& m_1\\ 0& 0& m_2& 0\\ 0& m_2& 0& M\\ m_1& 0& M& 0\end{array}\right).$$
(49)
With $`m_110^2`$ GeV, $`m_210^3`$ GeV, and $`M10^{16}`$ GeV, this implies that neutrinos have naturally small Dirac masses of order

$$m_\nu =\frac{m_1m_2}{M}10^2\mathrm{eV}.$$
(50)
## 8 Nearly Degenerate Majorana Neutrino Masses

I now switch gear and consider the derivation of the observed pattern of solar  and atmospheric  neutrino oscillations in the context of nearly degenerate Majorana neutrino masses. This is based  on the non-Abelian discrete symmetry $`A_4`$ (i.e. the symmetry group of a regular tetrahedron, the simplest of the 5 perfect geometric solids). It will be shown  that radiative corrections automatically generate the desired neutrino mass matrix and if these come from softly broken supersymmetry, then the effective mass measured in neutrinoless double beta decay  should not be much smaller than about 0.4 eV.

Suppose that at some high energy scale, the charged lepton mass matrix and the Majorana neutrino mass matrix are such that after diagonalizing the former, i.e.

$$_l=\left(\begin{array}{ccc}m_e& 0& 0\\ 0& m_\mu & 0\\ 0& 0& m_\tau \end{array}\right),$$
(51)
the latter is of the form

$$_\nu =\left(\begin{array}{ccc}m_0& 0& 0\\ 0& 0& m_0\\ 0& m_0& 0\end{array}\right).$$
(52)
From the high scale to the electroweak scale, one-loop radiative corrections will change $`_\nu `$ as follows:

$$(_\nu )_{ij}(_\nu )_{ij}+R_{ik}(_\nu )_{kj}+(_\nu )_{ik}R_{kj},$$
(53)
where the radiative correction matrix is assumed to be of the most general form, i.e.

$$R=\left(\begin{array}{ccc}r_{ee}& r_{e\mu }& r_{e\tau }\\ r_{e\mu }^{}& r_{\mu \mu }& r_{\mu \tau }\\ r_{e\tau }^{}& r_{\mu \tau }^{}& r_{\tau \tau }\end{array}\right).$$
(54)
Thus the observed neutrino mass matrix is given by

$$_\nu =m_0\left(\begin{array}{ccc}1+2r_{ee}& r_{e\tau }+r_{e\mu }^{}& r_{e\mu }+r_{e\tau }^{}\\ r_{e\mu }^{}+r_{e\tau }& 2r_{\mu \tau }& 1+r_{\mu \mu }+r_{\tau \tau }\\ r_{e\tau }^{}+r_{e\mu }& 1+r_{\mu \mu }+r_{\tau \tau }& 2r_{\mu \tau }^{}\end{array}\right).$$
(55)
Consider first the case where all the parameters are real. Redefine them as follows:

$`\delta _0r_{\mu \mu }+r_{\tau \tau }2r_{\mu \tau },`$ (56)
$`\delta 2r_{\mu \tau },`$ (57)
$`\delta ^{}r_{ee}{\displaystyle \frac{1}{2}}r_{\mu \mu }{\displaystyle \frac{1}{2}}r_{\tau \tau }r_{\mu \tau },`$ (58)
$`\delta ^{\prime \prime }r_{e\mu }+r_{e\tau }.`$ (59)
Then the neutrino mass matrix becomes

$$_\nu =m_0\left(\begin{array}{ccc}1+\delta _0+2\delta +2\delta ^{}& \delta ^{\prime \prime }& \delta ^{\prime \prime }\\ \delta ^{\prime \prime }& \delta & 1+\delta _0+\delta \\ \delta ^{\prime \prime }& 1+\delta _0+\delta & \delta \end{array}\right).$$
(60)
This matrix is $`exactly`$ diagonalized with the eigenvalues

$`m_1=m_0(1+2\delta +\delta ^{}\sqrt{\delta ^2+2\delta ^{\prime \prime 2}}),`$ (61)
$`m_2=m_0(1+2\delta +\delta ^{}+\sqrt{\delta ^2+2\delta ^{\prime \prime 2}}),`$ (62)
$`m_3=m_0,`$ (63)
where $`\delta _0`$ has been set equal to zero by a trivial rescaling of the other parameters, and the neutrino mixing matrix is given by

$$\left(\begin{array}{c}\nu _e\\ \nu _\mu \\ \nu _\tau \end{array}\right)=\left(\begin{array}{ccc}\mathrm{cos}\theta & \mathrm{sin}\theta & 0\\ \mathrm{sin}\theta /\sqrt{2}& \mathrm{cos}\theta /\sqrt{2}& 1/\sqrt{2}\\ \mathrm{sin}\theta /\sqrt{2}& \mathrm{cos}\theta /\sqrt{2}& 1/\sqrt{2}\end{array}\right)\left(\begin{array}{c}\nu _1\\ \nu _2\\ \nu _3\end{array}\right),$$
(64)
where

$$\mathrm{tan}\theta =\frac{\sqrt{2}\delta ^{\prime \prime }}{\sqrt{\delta ^2+2\delta ^{\prime \prime 2}}\delta ^{}},(\delta ^{}<0).$$
(65)
This is exactly the right description of present data on solar and atmospheric neutrino oscillations, with $`\mathrm{sin}^22\theta _{atm}=1`$ and $`\mathrm{tan}^2\theta _{sol}=0.38`$ if $`|\delta ^{\prime \prime }/\delta ^{}|=\sqrt{2}`$ for example. It also predicts that the common mass of the three neutrinos, i.e. $`m_0`$, is what is measured in neutrinoless double beta decay.

## 9 Discrete $`A_4`$ Symmetry

The successful derivation of Eq. (64) depends on having Eqs. (51) and (52). To be sensible theoretically, they should be maintained by a symmetry. At first sight, it appears impossible that there can be a symmetry which allows them to coexist. Here is where the non-Abelian discrete symmetry $`A_4`$ comes into play. The key is that $`A_4`$ has three inequivalent one-dimensional representations 1, 1, 1<sup>′′</sup>, and one three-dimensional reprsentation 3, with the decomposition

$$\underset{¯}{3}\times \underset{¯}{3}=\underset{¯}{1}+\underset{¯}{1}^{}+\underset{¯}{1}^{\prime \prime }+\underset{¯}{3}+\underset{¯}{3}.$$
(66)
This allows the following natural assignments of quarks and leptons:

$`(u_i,d_i)_L,(\nu _i,e_i)_L\underset{¯}{3},`$ (67)
$`u_{1R},d_{1R},e_{1R}\underset{¯}{1},`$ (68)
$`u_{2R},d_{2R},e_{2R}\underset{¯}{1}^{},`$ (69)
$`u_{3R},d_{3R},e_{3R}\underset{¯}{1}^{\prime \prime }.`$ (70)
Heavy fermion singlets are then added:

$$U_{iL(R)},D_{iL(R)},E_{iL(R)},N_{iR}\underset{¯}{3},$$
(71)
together with the usual Higgs doublet and new heavy singlets:

$$(\varphi ^+,\varphi ^0)\underset{¯}{1},\chi _i^0\underset{¯}{3}.$$
(72)
With this structure, charged leptons acquire an effective Yukawa coupling matrix $`\overline{e}_{iL}e_{jR}\varphi ^0`$ which has 3 arbitrary eigenvalues (because of the 3 independent couplings to the 3 inequivalent one-dimensional representations) and for the case of equal vacuum expectation values of $`\chi _i`$, i.e.

$$\chi _1=\chi _2=\chi _3=u,$$
(73)
the unitary transformation $`U_L`$ which diagonalizes $`_l`$ is given by

$$U_L=\frac{1}{\sqrt{3}}\left(\begin{array}{ccc}1& 1& 1\\ 1& \omega & \omega ^2\\ 1& \omega ^2& \omega \end{array}\right),$$
(74)
where $`\omega =e^{2\pi i/3}`$. This implies that the effective neutrino mass operator, i.e. $`\nu _i\nu _j\varphi ^0\varphi ^0`$, is proportional to

$$U_L^TU_L=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right),$$
(75)
exactly as desired .

## 10 Softly Broken Supersymmetry

To derive Eq. (75), the validity of Eq. (73) has to be proved. This is naturally accomplished in the context of supersymmetry. Let $`\widehat{\chi }_i`$ be superfields, then its superpotential is given by

$$\widehat{W}=\frac{1}{2}M_\chi (\widehat{\chi }_1\widehat{\chi }_1+\widehat{\chi }_2\widehat{\chi }_2+\widehat{\chi }_3\widehat{\chi }_3)+h_\chi \widehat{\chi }_1\widehat{\chi }_2\widehat{\chi }_3.$$
(76)
Note that the $`h_\chi `$ term is invariant under $`A_4`$, a property not found in $`SU(2)`$ or $`SU(3)`$. The resulting scalar potential is

$$V=|M_\chi +h_\chi \chi _2\chi _3|^2+|M_\chi \chi _2+h_\chi \chi _3\chi _1|^2+|M_\chi \chi _3+h_\chi \chi _1\chi _2|^2.$$
(77)
Thus a supersymmetric vacuum $`(V=0)`$ exists for

$$\chi _1=\chi _2=\chi _3=u=M_\chi /h_\chi ,$$
(78)
proving Eq. (73), with the important additional result that the spontaneous breaking of $`A_4`$ at the high scale $`u`$ does not break the supersymmetry.

To generate the proper radiative corrections which will result in a realistic Majorana neutrino mass matrix, $`A_4`$ is assumed broken also by the soft supersymmetry breaking terms. In particular, the mass-squared matrix of the left sleptons will be assumed to be arbitrary. This allows $`r_{\mu \tau }`$ to be nonzero through $`\stackrel{~}{\mu }_L\stackrel{~}{\tau }_L`$ mixing, from which the parameter $`\delta `$ may be evaluated. For illustration, using the approximation that $`\stackrel{~}{m}_1^2>>\mu ^2>>M_{1,2}^2=\stackrel{~}{m}_2^2`$, where $`\mu `$ is the Higgsino mass and $`M_{1,2}`$ are gaugino masses, I find

$$\delta =\frac{\mathrm{sin}\theta \mathrm{cos}\theta }{16\pi ^2}\left[(3g_2^2g_1^2)\mathrm{ln}\frac{\stackrel{~}{m}_1^2}{\mu ^2}\frac{1}{4}(3g_2^2+g_1^2)\left(\mathrm{ln}\frac{\stackrel{~}{m}_1^2}{\stackrel{~}{m}_2^2}\frac{1}{2}\right)\right].$$
(79)
Using $`\mathrm{\Delta }m_{32}^2=2.5\times 10^3`$ eV<sup>2</sup> from the atmospheric neutrino data, this implies that

$$\left[\mathrm{ln}\frac{\stackrel{~}{m}_1^2}{\mu ^2}0.3\left(\mathrm{ln}\frac{\stackrel{~}{m}_1^2}{\stackrel{~}{m}_2^2}\frac{1}{2}\right)\right]\mathrm{sin}\theta \mathrm{cos}\theta 0.535\left(\frac{0.4\mathrm{eV}}{m_0}\right)^2.$$
(80)
To the extent that the factor on the left cannot be much greater than unity, this means that $`m_0`$ cannot be much smaller than about 0.4 eV .

## 11 Nonzero $`U_{e3}`$ and CP Violation

The matrix of Eq. (55) has complex phases. It can easily been shown that only one phase remains after all possible redefinitions. As a result, the most general Majorana neutrino mass matrix derivable from $`A_4`$ is actually of the form

$$_\nu =m_0\left(\begin{array}{ccc}1+2\delta +2\delta ^{}& \delta ^{\prime \prime }& \delta ^{\prime \prime }\\ \delta ^{\prime \prime }& \delta & 1+\delta \\ \delta ^{\prime \prime }& 1+\delta & \delta \end{array}\right),$$
(81)
where only $`\delta ^{\prime \prime }`$ is complex. In this case, $`U_{e3}`$ becomes nonzero, i.e.

$$U_{e3}\frac{iIm\delta ^{\prime \prime }}{\sqrt{2}\delta }.$$
(82)
Note the important result that $`U_{e3}`$ is purely imaginary. Thus CP violation in neutrino oscillations is predicted to be maximal, which is desirable for future long-baseline neutrino experiments.

In the presence of $`Im\delta ^{\prime \prime }`$, the previous expressions are still approximately valid with the replacement of $`\delta ^{}`$ by $`\delta ^{}+(Im\delta ^{\prime \prime })^2/2\delta `$ and of $`\delta ^{\prime \prime }`$ by $`Re\delta ^{\prime \prime }`$. There is also the relationship

$$\left[\frac{\mathrm{\Delta }m_{12}^2}{\mathrm{\Delta }m_{32}^2}\right]^2\left[\frac{\delta ^{}}{\delta }+|U_{e3}|^2\right]^2+\left[\frac{Re\delta ^{\prime \prime }}{\delta }\right]^2.$$
(83)
Using $`\mathrm{\Delta }m_{12}^25\times 10^5`$ eV<sup>2</sup> from solar neutrino data and $`|U_{e3}|<0.16`$ from reactor neutrino data , I find

$`Im\delta ^{\prime \prime }<8.8\times 10^4(0.4\mathrm{eV}/m_0)^2,`$ (84)
$`Re\delta ^{\prime \prime }<7.8\times 10^5(0.4\mathrm{eV}/m_0)^2.`$ (85)
## 12 Conclusions

There are undoubtedly oases in the desert beyond the Standard Model. Two first examples of nontrivial U(1) gauge extensions of the Standard Model and its supersymmetric generalization have been discovered. In the first case , a heavy lepton triplet $`(\mathrm{\Sigma }^+,\mathrm{\Sigma }^0,\mathrm{\Sigma }^{})`$ per family is added to the Standard Model as the anchor for a naturally small seesaw Majorana neutrino mass (instead of using the canonical singlet $`N_R`$). This allows for the remarkable existence of a new U(1) gauge symmetry, with possible phenomenological consequences already at the electroweak scale, such as the NuTeV discrepancy. In the second case , with the addition of new superfields to the Supersymmetric Standard Model, a new U(1) gauge symmetry is again discovered which forbids proton decay, explains the size of the $`\mu `$ term, and allows for naturally small Dirac neutrino masses. The existence of this new gauge symmetry involves the highly nontrivial exact factorization of the sum of 11 cubic terms.

A third possible oasis is that of an underlying $`A_4`$ symmetry at some high energy scale, which allows the observed Majorana neutrino mass matrix to be derived from radiative corrections. It has been shown  that this automatically leads to $`\mathrm{sin}^22\theta _{atm}=1`$ and a large (but not maximal) solar mixing angle. Using neutrino oscillation data, and assuming radiative corrections from soft supersymmetry breaking, the effective mass measured in neutrinoless double beta decay is predicted to be not much less than 0.4 eV.

## Acknowledgements

I thank Hans Klapdor, Juha Peltoniemi, and the other organizers of Beyond 2002 for their great hospitality at Oulu. This work was supported in part by the U. S. Department of Energy under Grant No. DE-FG03-94ER40837.

## References


text:
# Wave packets of a harmonic oscillator with various degrees of rigidity

## I  Introduction

In recent years, the time evolution of localized wave packets of various quantum-mechanical systems has been widely discussed in the literature (see Ref. Bluhm et al. (1996) for a review). As for wave packets of the harmonic oscillator, the subject is as old as quantum mechanics itself. Though some aspects of the Gaussian wave packets (displaced or squeezed ground state) are still being discussed in the current literature de Castro and da Cruz (1999); Waldenstrøm and Naqvi (1999), the problem was essentially solved in the very early years Schrödinger (1926); Kennard (1927); Schiff (1968). In the 1950’s, several authors studied more general displaced number states and found that they keep their shape unchanged while their center oscillates like a classical particle Senitzky (1954); Plebanski (1956, 1954); Epstein (1959). Displaced and squeezed number states of special forms were also considered in some of these papers Plebanski (1956, 1954); it turns out that their width is also oscillating, so the shapes of such wave packets change with time apparently. Displaced and squeezed number states of the most general form were studied in a recent paper Nieto (1997). This may represent the most general case where the time dependence of the wave packets could be worked out explicitly.

The questions we are concerned with are: First, if the time dependence of the wave packet cannot be obtained explicitly, what conclusions could be made in regard to its time evolution? Second, if the shape of the wave packet changes with time, how can we describe the level of changes for different cases? For the first question, it is known that the center of the wave packet moves like a classical particle while its width pulsates. However, one can find wave packets whose width keeps unchanged but whose shape still changes with time. To describe the change of the shape in the general case we introduce general-order moments of $`x`$ about its center, defined as the mean value of $`(x\overline{x}_t)^K`$ and denoted by $`Q_K(t)`$, where $`\overline{x}_t`$ is the mean value of $`x`$ (which represents the center of the wave packet) and $`K`$ a natural number. Differential equations for these moments are established for general $`K`$, and are solved for the most general initial conditions up to $`K=4`$. They are essentially oscillating, but higher-order ones involve higher frequencies and thus their time dependence is more complicated. Wave packets could be found whose $`Q_K(t)`$ are time independent up to $`K=4`$. Obviously the shape of such wave packets changes with time less than those that only keep their width \[represented by $`\sqrt{Q_2(t)}`$\] unchanged. At this stage the answer to the second question is becoming clear. We introduce a natural number called the “degree of rigidity” for a wave packet: If $`Q_K(t)`$ is time independent for $`K=2,3,\mathrm{},2N`$ \[note that $`Q_0(t)1`$ and $`Q_1(t)0`$\] but not for $`K=2N+2`$ (for $`K=2N+1`$ it may or may not be time independent), we say that the wave packet has degree of rigidity $`N`$. It is obvious that the larger $`N`$ is, the less the shape changes with time. The shape of a wave packet with a time-dependent width changes with time apparently and thus has no rigidity. One with a constant width but with all higher-order moments time dependent has the “ground” degree of rigidity 1. On the other hand, if $`Q_K(t)`$ is time independent for all $`K`$, then the degree of rigidity for such a wave packet is infinity. In other words, its shape is perfectly rigid. Typical examples of such wave packets are displaced number states. Two classes of wave packets with an arbitrarily given degree of rigidity are presented.

## II  Some general results

Consider localized wave packets of the harmonic oscillator in one dimension, whose time evolution is governed by the Hamiltonian

$$H=\frac{p^2}{2\mu }+\frac{1}{2}\mu \omega ^2x^2.$$
(1)
The normalized initial wave function $`\psi (x,0)`$ at $`t=0`$ will be denoted by $`\psi _0`$ for convenience. Localization means that the mean value of $`x^kp^l`$ in $`\psi _0`$ is finite for any nonnegative $`k`$ and $`l`$. The wave function $`\psi (x,t)`$ at time $`t`$ will be denoted by $`\psi _t`$. The mean values of $`x`$ and $`p`$ at time $`t`$ are

$$\overline{x}_t=(\psi _t,x\psi _t),\overline{p}_t=(\psi _t,p\psi _t).$$
(2)
Using the Schrödinger equation it is easy to show that

$$\dot{\overline{x}}_t=\overline{p}_t/\mu ,\dot{\overline{p}}_t=\mu \omega ^2\overline{x}_t.$$
(3)
The solution of these equations is

$$\overline{x}_t=\overline{x}_0\mathrm{cos}\omega t+\frac{\overline{p}_0}{\mu \omega }\mathrm{sin}\omega t,\overline{p}_t=\overline{p}_0\mathrm{cos}\omega t\mu \omega \overline{x}_0\mathrm{sin}\omega t.$$
(4)
The solution means that the center of an arbitrary wave packet moves like a classical particle. This is a well known result.

To study the change with time of the shape of the wave packet we consider the mean value of $`(x\overline{x}_t)^K`$ in $`\psi _t`$ where $`K`$ is a natural number. These mean values are denoted by $`Q_K(t)`$ in the following. If the $`Q_K(t)`$’s are time independent for all $`K`$, then the wave packet obviously keeps its shape unchanged while it is moving, and this is the case for displaced number states Senitzky (1954); Plebanski (1956, 1954); Epstein (1959). In order to examine their time dependence for a general wave packet, we define the real quantities

$$R_{kl}(t)=\frac{1}{2}\mathbf{(}\psi _t,\{(x\overline{x}_t)^k,(p\overline{p}_t)^l\}\psi _t\mathbf{)},S_{kl}(t)=\frac{1}{2\mathrm{i}}\mathbf{(}\psi _t,[(x\overline{x}_t)^k,(p\overline{p}_t)^l]\psi _t\mathbf{)},$$
(5)
where $`k`$ and $`l`$ are nonnegative integers, and $`\{F,G\}=FG+GF`$, $`[F,G]=FGGF`$ for any operator $`F`$ and $`G`$. Obviously $`R_{K0}(t)=Q_K(t)`$. In the following we also denote $`R_{0K}(t)`$ by $`P_K(t)`$. It is straightforward to show that they satisfy the following equations

$$\dot{R}_{kl}(t)=\frac{k}{\mu }R_{k1,l+1}(t)l\mu \omega ^2R_{k+1,l1}(t)+\frac{\mathrm{}}{2\mu }k(k1)S_{k2,l}(t)\frac{\mathrm{}\mu \omega ^2}{2}l(l1)S_{k,l2}(t),$$
(6a)
$$\dot{S}_{kl}(t)=\frac{k}{\mu }S_{k1,l+1}(t)l\mu \omega ^2S_{k+1,l1}(t)\frac{\mathrm{}}{2\mu }k(k1)R_{k2,l}(t)+\frac{\mathrm{}\mu \omega ^2}{2}l(l1)R_{k,l2}(t).$$
(6b)
Note that by definition we have $`S_{k0}(t)=S_{0k}(t)=0`$ for all $`k`$. We are mainly interested in $`R_{kl}(t)`$, especially $`Q_K(t)`$. It is easy to realize that the equations for the subset $`\{R_{kl}(t)|k+l=K\}`$ close among themselves provided that the lower-order subset $`\{S_{kl}(t)|k+l=K2\}`$ has been obtained. The equations for the subset $`\{S_{kl}(t)|k+l=K\}`$ also close among themselves, provided that the lower-order subset $`\{R_{kl}(t)|k+l=K2\}`$ is known. Now for $`k=l=0`$ we have $`R_{00}(t)=1`$ and $`S_{00}(t)=0`$ by definition. Similarly for $`k+l=1`$ we have $`R_{10}(t)=R_{01}(t)=S_{10}(t)=S_{01}(t)=0`$. These enable us to solve the cases with $`k+l2`$.

The case $`k+l=2`$ is of essential interest, because it involves the width of the wave packet. By definition we have $`S_{20}(t)=S_{02}(t)=0`$, and $`S_{11}(t)=\mathrm{}/2`$. These will be useful when the case $`k+l=4`$ is considered. For the moment we are interested in the remaining three. It is not difficult to find that

$`Q_2(t)={\displaystyle \frac{\mu ^2\omega ^2Q_2(0)+P_2(0)}{2\mu ^2\omega ^2}}+{\displaystyle \frac{\mu ^2\omega ^2Q_2(0)P_2(0)}{2\mu ^2\omega ^2}}\mathrm{cos}2\omega t+{\displaystyle \frac{R_{11}(0)}{\mu \omega }}\mathrm{sin}2\omega t,`$ (7a)
$`P_2(t)={\displaystyle \frac{\mu ^2\omega ^2Q_2(0)+P_2(0)}{2}}{\displaystyle \frac{\mu ^2\omega ^2Q_2(0)P_2(0)}{2}}\mathrm{cos}2\omega t\mu \omega R_{11}(0)\mathrm{sin}2\omega t,`$ (7b)
$`R_{11}(t)=R_{11}(0)\mathrm{cos}2\omega t{\displaystyle \frac{\mu ^2\omega ^2Q_2(0)P_2(0)}{2\mu \omega }}\mathrm{sin}2\omega t.`$ (7c)
The width of the wave packet at time $`t`$ is characterized by the quantity $`\mathrm{\Delta }_tx=\sqrt{Q_2(t)}`$, and that in the momentum space by $`\mathrm{\Delta }_tp=\sqrt{P_2(t)}`$. Thus the width of the wave packet oscillates with frequency $`2\omega `$, just like that for a squeezed ground state. It can be shown that both $`Q_2(t)`$ and $`P_2(t)`$ are positive as they should be. The uncertainty product $`\mathrm{\Delta }_tx\mathrm{\Delta }_tp`$ is also oscillating. These conclusions (for a general wave packet) have been discussed in some different way previously Plebanski (1956, 1954), so we will not go into more details. A relation that seems not to be emphasized in the literature is

$$\mu ^2\omega ^2(\mathrm{\Delta }_tx)^2+(\mathrm{\Delta }_tp)^2=\mu ^2\omega ^2(\mathrm{\Delta }_0x)^2+(\mathrm{\Delta }_0p)^2.$$
(8)
According to this result, $`\mathrm{\Delta }_tx`$ reaches its minimum when $`\mathrm{\Delta }_tp`$ reaches its maximum and vice versa.

When $`k+l=3`$, all the $`S_{kl}(t)`$ can be found to be zero. The equations for the $`R_{kl}(t)`$ can be solved without much difficulty. They are all linear combinations of $`\{\mathrm{sin}\omega t,\mathrm{cos}\omega t,\mathrm{sin}3\omega t,\mathrm{cos}3\omega t\}`$. Since the results are lengthy and not important for further discussions we will not write them down. We just point out that they all vanish if the corresponding initial values $`R_{kl}(0)`$ are all zero.

When $`k+l=4`$, it can be found that $`S_{40}(t)=S_{04}(t)=0`$, $`S_{31}(t)=3\mathrm{}Q_2(t)/2`$, $`S_{13}(t)=3\mathrm{}P_2(t)/2`$, and $`S_{22}(t)=2\mathrm{}R_{11}(t)`$. However, the results for the $`R_{kl}(t)`$ are rather lengthy, they are linear combinations of $`\{1,\mathrm{sin}2\omega t,\mathrm{cos}2\omega t,\mathrm{sin}4\omega t,\mathrm{cos}4\omega t\}`$. We only write down one of them here:

$`Q_4(t)`$ $`=`$ $`{\displaystyle \frac{3\mu ^4\omega ^4Q_4(0)+3P_4(0)+6\mu ^2\omega ^2R_{22}(0)+3\mathrm{}^2\mu ^2\omega ^2}{8\mu ^4\omega ^4}}`$ (9)
$`+`$ $`{\displaystyle \frac{\mu ^4\omega ^4Q_4(0)P_4(0)}{2\mu ^4\omega ^4}}\mathrm{cos}2\omega t+{\displaystyle \frac{R_{13}(0)+\mu ^2\omega ^2R_{31}(0)}{\mu ^3\omega ^3}}\mathrm{sin}2\omega t`$
$`+`$ $`{\displaystyle \frac{\mu ^4\omega ^4Q_4(0)+P_4(0)6\mu ^2\omega ^2R_{22}(0)3\mathrm{}^2\mu ^2\omega ^2}{8\mu ^4\omega ^4}}\mathrm{cos}4\omega t`$
$``$ $`{\displaystyle \frac{R_{13}(0)\mu ^2\omega ^2R_{31}(0)}{2\mu ^3\omega ^3}}\mathrm{sin}4\omega t.`$
We see that higher-order moments about the center have more complicated time dependence, but they are essentially oscillating and hence are periodic functions. This is also true for still higher ones. Indeed, $`\overline{x}_t`$, $`\overline{p}_t`$ and $`\psi _t`$ are all periodic, and so are $`R_{kl}(t)`$ and $`S_{kl}(t)`$.

For larger values of $`k+l`$, the solutions are more difficult to find. We will not proceed further in this respect.

## III  Simplified results for special cases

The results obtained in Sec. II are rather complicated. In the following sections we will discuss the results for special cases with the initial condition

$$\psi _0=\psi (x,0)=\phi (xx_0)\mathrm{e}^{\mathrm{i}p_0x/\mathrm{}},$$
(10a)
where $`x_0`$ and $`p_0`$ are real constants, and $`\phi (x)`$ has definite parity, namely
$$\phi (x)=\pm \phi (x).$$
(10b)
In this initial state it is not difficult to show that

$$\overline{x}_0=x_0,\overline{p}_0=p_0,$$
(11)
and

$$R_{kl}(0)=0,S_{kl}(0)=0,k+l=1,3,5,\mathrm{}.$$
(12)
The latter could be derived from the useful relations

$`\mathbf{(}\psi _0,(xx_0)^k(pp_0)^l\psi _0\mathbf{)}=\mathbf{(}\phi (x),x^kp^l\phi (x)\mathbf{)},`$
$`\mathbf{(}\psi _0,(pp_0)^l(xx_0)^k\psi _0\mathbf{)}=\mathbf{(}\phi (x),p^lx^k\phi (x)\mathbf{)},`$
$`k,l=0,1,2,\mathrm{},`$ (13)
which are easy to show.

The shape of the initial wave packet, $`|\psi _0|^2`$, is obviously symmetric about the center $`x_0`$. This symmetry will be kept at later times. To prove this it is sufficient to show that $`Q_{2K1}(t)=0`$ for all natural numbers $`K`$. Indeed, we will show that

$$R_{kl}(t)=0,S_{kl}(t)=0,k+l=1,3,5,\mathrm{}.$$
(14)
For $`k+l=1`$ it is true as given in Sec. II. Now suppose that it is true for $`k+l=2K1`$, and consider the case with $`k+l=2K+1`$. It is easy to realize that the equations for the subset $`\{R_{kl}(t)|k+l=2K+1\}`$ close among themselves and are all homogeneous because all of the $`S_{kl}(t)`$ in the subset $`\{S_{kl}(t)|k+l=2K1\}`$ vanish (according to the assumption). Since the initial conditions are all homogeneous too, the solutions are obviously $`R_{kl}(t)=0`$ ($`k+l=2K+1`$). Similarly, $`S_{kl}(t)=0`$ for $`k+l=2K+1`$ as well. The case with $`k+l=3`$ has been explicitly calculated, and the result is consistent with this general conclusion, as pointed out in Sec. II.

If $`\phi (x)`$ is real, we have another useful consequence:

$$R_{kl}(0)=0,k,l=1,3,5,\mathrm{}.$$
(15)
As a result the solutions $`R_{kl}(t)`$ with $`k+l=2K`$ are also simplified. For example, the sine terms in Eqs. (7a), (7b), and (9) all vanish. However, that $`\phi (x)`$ is real is just a sufficient condition, not a necessary one for Eq. (15).

Now if $`\phi (x)`$ in the initial state satisfies the condition

$$\mathbf{(}\phi (x),(xp+px)\phi (x)\mathbf{)}=0,\mathbf{(}\phi (x),\mu ^2\omega ^2x^2\phi (x)\mathbf{)}=\mathbf{(}\phi (x),p^2\phi (x)\mathbf{)},$$
(16)
we have \[cf. Eq. (III)\]

$$R_{11}(0)=0,\mu ^2\omega ^2Q_2(0)=P_2(0),$$

and from Eq. (7) we obtain

$$Q_2(t)=Q_2(0),P_2(t)=P_2(0),R_{11}(t)=0.$$

That is, the width of the wave packet keeps unchanged. The second condition in Eq. (16) means that the kinetic energy and the potential energy have the same mean value in $`\phi (x)`$. It is easy to realize that all the number states $`\phi _n(x)`$ satisfy the conditions (10b) and (16). However, there exist many other functions that satisfy these conditions. For example,

$$\phi ^{\text{even}}(x)=\underset{i=1}{\overset{\mathrm{}}{}}a_i\phi _{2n_i}(x),0n_1<n_2<\mathrm{},\underset{i=1}{\overset{\mathrm{}}{}}|a_i|^2=1$$
(17a)
and
$$\phi ^{\text{odd}}(x)=\underset{i=1}{\overset{\mathrm{}}{}}b_i\phi _{2n_i+1}(x),0n_1<n_2<\mathrm{},\underset{i=1}{\overset{\mathrm{}}{}}|b_i|^2=1$$
(17b)
all have definite parity. And as long as $`n_{i+1}n_i2`$ for all $`i`$, they also satisfy the conditions in Eq. (16).

However, the shapes of such wave packets change with time, though their width keeps unchanged. The reason is that $`Q_4(t)`$, and in general the higher-order moments, are still oscillating. The conditions for $`Q_4(t)`$, $`P_4(t)`$ etc. to be time independent are

$$R_{13}(0)=R_{31}(0)=0,\mu ^4\omega ^4Q_4(0)=P_4(0),2\mu ^2\omega ^2Q_4(0)6R_{22}(0)=3\mathrm{}^2.$$
(18)
If $`\phi (x)=\phi _n(x)`$, a number state, then it can be shown that these conditions are satisfied, as expected. However, there exist many other functions that satisfy these conditions. In fact, the two classes of functions given in Eq. (17) do if $`n_{i+1}n_i3`$ for all $`i`$. In this case Eq. (16) is of course satisfied too. Since now $`Q_4(t)`$ is also time independent, the shapes of the wave packets changes with time less than the case where only the width is kept unchanged. In other words, their shapes are more rigid.

In general, the higher-order moments, say $`Q_6(t)`$, of the above wave packets will still be time dependent. It will become more and more difficult to discuss the problem in the above manner. Since we are not going to find the most general results for these higher-order moments, we will proceed in a different way in the next section.

## IV  Wave packets with various degrees of rigidity

In the last section we have found wave packets whose width keeps unchanged with time, and ones whose $`Q_K(t)`$ up to $`K=4`$ all keep unchanged with time. Obviously the shape of the latter changes with time less than the former. Now we introduce a natural number called the “degree of rigidity” for a wave packet: If $`Q_K(t)`$ is time independent for $`K=2,3,\mathrm{},2N`$ \[note that $`Q_0(t)1`$ and $`Q_1(t)0`$\] but not for $`K=2N+2`$ (the situation for $`K=2N+1`$ is not important in this definition), we say that the wave packet has degree of rigidity $`N`$. Thus the two cases mentioned above have degrees of rigidity 1 and 2, respectively. It is obvious that the larger $`N`$ is, the less the shape changes with time. If the width of a wave packet changes with time, then its shape changes with time apparently and thus has no rigidity. On the other hand, if $`Q_K(t)`$ is time independent for all $`K`$, then the degree of rigidity for such a wave packet is infinity. In other words, its shape is perfectly rigid. Typical examples of such wave packets are displaced number states.

In order to find wave packets with a given degree of rigidity, we consider the quantity

$$W_{kl}(t)=\mathbf{(}\psi _t,(x\overline{x}_t)^k(p\overline{p}_t)^l\psi _t\mathbf{)}=R_{kl}(t)+\mathrm{i}S_{kl}(t).$$
(19)
Because $`\psi _t=\mathrm{e}^{\mathrm{i}Ht/\mathrm{}}\psi _0`$, we have

$$W_{kl}(t)=\mathbf{(}\psi _0,(x_t\overline{x}_t)^k(p_t\overline{p}_t)^l\psi _0\mathbf{)},$$
(20)
where

$$x_t=\mathrm{e}^{\mathrm{i}Ht/\mathrm{}}x\mathrm{e}^{\mathrm{i}Ht/\mathrm{}}=x\mathrm{cos}\omega t+\frac{p}{\mu \omega }\mathrm{sin}\omega t,$$
(21a)
$$p_t=\mathrm{e}^{\mathrm{i}Ht/\mathrm{}}p\mathrm{e}^{\mathrm{i}Ht/\mathrm{}}=p\mathrm{cos}\omega t\mu \omega x\mathrm{sin}\omega t.$$
(21b)
Now we confine ourselves to initial states of the form (10). On account of Eqs. (4) and (11), we have

$$W_{kl}(t)=\mathbf{(}\psi _0,[(xx_0)\mathrm{cos}\omega t+(\mu \omega )^1(pp_0)\mathrm{sin}\omega t]^k[(pp_0)\mathrm{cos}\omega t\mu \omega (xx_0)\mathrm{sin}\omega t]^l\psi _0\mathbf{)}.$$
(22)
Note that for any function $`f(x)`$ we have $`(pp_0)[\mathrm{e}^{\mathrm{i}p_0x/\mathrm{}}f(x)]=\mathrm{e}^{\mathrm{i}p_0x/\mathrm{}}pf(x)`$, and that the inner product is an integral over $`x`$, the above equation can be simplified as

$$W_{kl}(t)=\mathbf{(}\phi (x),[x\mathrm{cos}\omega t+(\mu \omega )^1p\mathrm{sin}\omega t]^k(p\mathrm{cos}\omega t\mu \omega x\mathrm{sin}\omega t)^l\phi (x)\mathbf{)}.$$
(23)
Using Eqs. (21) again this becomes

$$W_{kl}(t)=\mathbf{(}\phi _t(x),x^kp^l\phi _t(x)\mathbf{)},$$
(24)
where

$$\phi _t(x)=\mathrm{e}^{\mathrm{i}Ht/\mathrm{}}\phi (x).$$
(25)
If $`\phi (x)=\phi _n(x)`$, a number state, then $`\phi _t(x)=\mathrm{e}^{\mathrm{i}E_nt/\mathrm{}}\phi _n(x)`$, and

$$W_{kl}(t)=\mathbf{(}\phi _n(x),x^kp^l\phi _n(x)\mathbf{)}=W_{kl}(0),k,l=0,1,2,\mathrm{}.$$
(26)
Especially, $`Q_K(t)=W_{K0}(t)=Q_K(0)`$ for all $`K`$, so that the shape of a displaced number state is perfectly rigid. This is a well-known conclusion. Here we obtain the conclusion in a different and also simple way.

Consider the function $`\phi ^{\text{even}}(x)`$ given in Eq. (17a), we have

$$\phi _t^{\text{even}}(x)=\underset{i=1}{\overset{\mathrm{}}{}}a_i\mathrm{exp}(\mathrm{i}E_{2n_i}t/\mathrm{})\phi _{2n_i}(x).$$
(27)
Assuming that $`n_{i+1}n_iN+1`$. Because $`x^kp^l\phi _{2n_i}(x)`$ is a linear combination of $`\{\phi _{2n_i+k+l}(x),\phi _{2n_i+k+l2}(x),\mathrm{},\phi _{2n_i(k+l)}(x)\}`$, we have for $`k+l2N`$ that

$$W_{kl}(t)=\mathbf{(}\phi _t^{\text{even}}(x),x^kp^l\phi _t^{\text{even}}(x)\mathbf{)}=\underset{i=1}{\overset{\mathrm{}}{}}|a_i|^2\mathbf{(}\phi _{2n_i}(x),x^kp^l\phi _{2n_i}(x)\mathbf{)}=W_{kl}(0).$$
(28)
Especially

$$Q_K(t)=Q_K(0),K=2,3,\mathrm{},2N.$$
(29)
The same result can be obtained for the function $`\phi ^{\text{odd}}(x)`$ given in Eq. (17b). \[Remember that $`Q_K(t)=0`$ for all odd $`K`$ as proved in Sec. III.\] Thus, if $`n_{i+1}n_iN+1`$, the two classes of functions given in Eq. (17) lead to wave packets \[whose initial wave functions are given by Eq. (10a)\] that have degree of rigidity not less than $`N`$. If some of the differences $`n_{i+1}n_i`$ equals $`N+1`$, they will in general have degree of rigidity $`N`$. Otherwise the degree of rigidity will be greater than $`N`$. For large $`N`$, the shapes of these wave packets are almost unchanged with time.

## V  Summary

In this paper we studied the time evolution of a general wave packet in the harmonic oscillator potential by examining the time dependence of the various moments $`Q_K(t)`$ of $`x`$ about the center of the wave packet. The differential equations for these objects are derived for general $`K`$, and are solved for the most general initial conditions up to $`K=4`$. Wave packets with constant width, and ones whose $`Q_K(t)`$ up to $`K=4`$ are all constants, are discussed. These include displaced number states as special cases. In general the shapes of these wave packets still change with time because the higher order moments are oscillating. A natural number $`N`$, called the degree of rigidity, is introduced to describe qualitatively how much the shape is changed with time. The larger $`N`$ is, the less the shape is changed with time. Displaced number states are perfectly rigid and have $`N=\mathrm{}`$. Two classes of wave packets with an arbitrarily given degree of rigidity are given explicitly.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of the People’s Republic of China, and by the Foundation of the Advanced Research Center of Sun Yat-Sen University.