text:
# References

BIHEP-TH-2003-38

Noncommutative Theory in Light of Neutrino Oscillation

Shao-Xia Chen111ruxanna@mail.ihep.ac.cn and Zhao-Yu Yang222yangzhy@mail.ihep.ac.cn

Institute of High Energy Physics, Chinese Academy of Sciences

P.O.Box 918(4), 100039 Beijing, China

## Abstract

Solar neutrino problem and atmospheric neutrino anomaly which are both long-standing issues studied intensively by physicists in the past several decades, are reckoned to be able to be solved simultaneously in the framework of the assumption of the neutrino oscillation. For the presence of the Lorentz invariance in the Standard Model, the massless neutrino can’t have flavor mixing and oscillation. However, we exploit the $`q`$-deformed noncommutative theory to derive a general modified dispersion relation, which implies some violation of the Lorentz invariance. Then it is found that the application of the $`q`$-deformed dispersion relation to the neutrino oscillation can provide a sound explanation for the current data from the reactor and long baseline experiments.

PACS numbers: 14.60.Pq, 11.10.Nx, 11.30.Cp.

The current knowledge concerning high energy neutrinos mainly comes from solar, atmospheric and long baseline neutrino experiments, which all present some inconsistence with the electro-weak Standard Model. In the Super-Kamiokande (S-K) experiment , the results of directly detecting $`\nu _ee`$ scattering via the observation of the Cherenkov light in the large water Cherenkov detector, and the detection in the SNO - using a heavy water Cherenkov detector have significantly confirmed that the number of solar neutrinos detected in these experiments is less than that predicted by the theories of the standard solar model. Some experiments detecting atmospheric neutrino, including several water Cerenkov experiments(Kamiokande , S-K and IMB ) together with the iron calorimeter Soudan 2 experiment , observing atmospheric neutrinos over energies varying from sub-GeV to tens of GeV have shown obvious discrepancy between the expected and measured ratio of the numbers of muon and electron events $`R_{\mu /e}`$. As one of the new generation of long baseline neutrino oscillation searches, KamLAND  has published its first results recently, which give strong evidence for the disappearance of neutrinos travelling from a power reactor to a far detector. Meanwhile, the results of long baseline K2K experiment  imply a reduction of $`\nu _\mu `$ flux together with a distortion of the energy spectrum. Thanks to all the experiments above, it is now convincing that the long-standing solar neutrino $`\nu _e`$ deficit and the atmospheric neutrino $`\nu _\mu `$ anomaly are both from neutrino oscillations between the flavor eigenstates.

Neutrino oscillations may arise, in general, when some terms are added to the neutrino sector of the Standard Model Lagrangian in the way generating the discordance of the eigenstates of the total Hamiltonian with the neutrino flavor eigenstates. In addition to the general scenario of neutrino oscillations - from the mismatch between the diagonalization of the charged lepton mass matrix and that of the neutrino mass matrix in an arbitrary flavor basis, there have appeared considerable papers - discussing the neutrino oscillation in the context of the premise that the oscillations may be the effects of the disagreement of the flavor eigenstates with the velocity eigenstates which are defined as the energy eigenstates of the massless neutrino and much has already been achieved. In the perspective of the violation of the Lorentz invariance several kinds of outcomes have appeared in the literature: $`\lambda E^0`$ , $`\lambda E^1`$ , $`\lambda E^2`$  and $`\lambda E^3`$ , here $`\lambda `$ is the neutrino oscillation length which will be defined in the subsequent discussion and $`E`$ is the mean value of observed neutrino energy.

In this Letter, we will start from a modified dispersion relation induced by the $`q`$-deformed noncommutative theory which also indicates the violation of the Lorentz invariance, to investigate the oscillations between the different flavor eigenstates of the massless neutrinos.

The $`q`$-deformed noncommutative theory, as one of the significant application of the quantum group - to the modern physics, plays a useful role in the study of the high energy physical phenomena. In the $`q`$-deformed noncommutative theory physical quantity $`[x]`$ is defined in terms of a parameter $`q`$ (taken to be real for simplicity):

$$[x]=\frac{q^xq^x}{qq^1}.$$

The parameter $`q`$ measures the degree of the deviation of the considered system from the usual commutative case. Now that neutrino oscillation is a physical phenomenon related to the GeV scale energy, it is natural to make use of the $`q`$-deformed noncommutative theory to investigate it. The dispersion relation for a fermion in the $`q`$-deformed case is of the form :

$$E=\sqrt{m^2+p^2}+\frac{\sqrt{m^2+p^2}(4m^2+4p^2\omega ^2)}{24\omega ^2}(q1)^2.$$
(1)
For the purpose of subsequent discussion, we derive the approximated form of Eq. (1) for massless neutrinos:

$$E=p\left(1+\frac{(q1)^2p^2}{6\omega ^2}\right)=p(1+\alpha p^2),$$
(2)
here the notation $`\alpha `$=$`{\displaystyle \frac{(q1)^2}{6\omega ^2}}`$ has been introduced and characterizes the Lorentz invariance violation in $`q`$-deformed noncommutative theory.

As argued in , the Lorentz invariance violation would indicate that neutrinos may differ in their maximal attainable velocities, which disagree with the light velocity $`c`$ in vacua. In this senario, neutrino oscillations can occur if the neutrino flavor eigenstates decided by the weak interactions are coherent superposition of the neutrino velocity eigenstates. In the case of two neutrino flavors, the flavor eigenstates can be expressed in terms of the velocity eigenstates:

$$\nu _\mu =\nu _1\mathrm{cos}\theta +\nu _2\mathrm{sin}\theta ,\nu _e=\nu _2\mathrm{cos}\theta \nu _1\mathrm{sin}\theta .$$
(3)
here $`\theta `$ is the mixing angle between different flavor eigenstates. The probability that neutrino oscillates from $`\alpha `$ flavor to $`\beta `$ flavor is:

$$P_{\alpha \beta }=\mathrm{sin}^22\theta \mathrm{sin}^2\left(\frac{\pi L}{\lambda }\right),$$
(4)
where $`L`$ is the cosmological distance travelled by neutrinos between the emission and the detection, and the oscillation length $`\lambda `$ is defined as $`\lambda ={\displaystyle \frac{2\pi }{E_1E_2}}`$. The mixing angle $`\theta `$ only determines the amplitude of the oscillation but not affect the oscillation length. It should be noted that, although here we consider merely the simplified case of two neutrino flavors, however, in the case of more than two flavors, Eq. (4) will take a more complicated form but not contain more intrinsical physical contents, so we will concentrate only on the two flavors case.

During their evolution, neutrinos propagate as a linear superposition of their velocity eigenstates whose energy eigenvalues are $`E_a`$ ($`a=1,2`$). Exploiting Eq. (2), one obtains that the energy difference is

$$E_1E_2=\delta \alpha E^3,$$
(5)
where $`\delta \alpha `$=$`\alpha _1\alpha _2`$. Combining the definition of the oscillation length $`\lambda `$, (5) and (4), one can acquire

$$\lambda =\frac{2\pi }{\delta \alpha E^3},P_{\alpha \beta }=\mathrm{sin}^22\theta \mathrm{sin}^2\left(\frac{LE^3\delta \alpha }{2}\right).$$
(6)
It is obvious that the oscillation length $`\lambda `$ is of the form $`\lambda E^3`$, which is identical with that appeared in the literature . For a ulterior discussion, we can plot several curves in Fig.1 describing the dependence of the neutrino oscillation length $`\lambda `$ on $`\delta \alpha `$ for different energies. From the curves we find that only a tiny deviation of $`q`$ from 1 can unravel the existent data, which means that in the phenomenon of neutrino oscillation, at least on the current energy scale of detection, the degree of noncommutativity is very small. Of course, we must note that here we only take into account the neutrino oscillation from the disagreement of the neutrino flavor eigenstates with velocity eigenstates and omit the conventional oscillation mechanism in which neutrino oscillation comes from the difference between neutrino mass eigenstates and flavor eigenstates. However, it is evident that in the presence of neutrino mass the requisite deviation of $`q`$ from 1 will be smaller than that for massless neutrinos. In general, if neutrinos have mass, they may have simultaneous velocity and mass oscillations. Then exploiting the argument in , one can further discuss the situation in the presence of two oscillations.

Experimental Data

| EXP. | STATUS | $`E`$(GeV) | L(km) |
| --- | --- | --- | --- |
| CHORUS | closed 1997 | 26 | 0.85 |
| NOMAD | closed 1999 | 24 | 0.94 |
| SK | operating | 1.3 | $`1010^4`$ |
| K2K | operating | 1.3 | 250 |
| SNO | operatng | 0.008 | $`10^8`$ |
| MINOS | starting 2003 | 15 | 730 |
| CNGS | starting 2005 | 17 | 732 |

Table 1. Shown for each experiment are its operation status, mean value of observed neutrino energy and typical neutrino flight distance $`L`$. The table is quoted from the paper .

To summarize, we exhibit in this Letter the relation between the neutrino oscillation length $`\lambda `$ and the Lorentz invariance violation parameter $`\delta \alpha `$, which can be verified by the existent and future neutrino oscillation experiments. Especially, when the data from long baseline experiments and neutrino factories come to be available, Lorentz invariance violation in neutrino oscillation can be probed to new and significant levels. Once the oscillation length are measured, $`\delta \alpha `$ will be derived as additional contribution apart from the observations of ultrahigh energy cosmic rays  to determine the $`q`$-deformed noncommutativity parameter.

Acknowledgements

The work was supported in part by the Natural Science Foundation of China.


format: text

text:
# The CMB Anisotropy Power Spectrum from the Background Emission Anisotropy Scanning Telescope Experiment

## 1 Introduction

Understanding the mechanisms of structure formation in the early universe ($`10<z<1000`$) is one of the most important and active areas in Cosmology today and measurements of the Cosmic Microwave Background (CMB) anisotropy play a pivotal role in this field. In the framework of the standard cosmological model, the CMB radiation is interpreted as the blackbody radiation associated with a hot dense phase of the Universe, when matter and radiation were in thermal equilibrium (e.g. Peebles, 1993). On large angular scales the CMB radiation traces the primordial power spectrum set by physical processes during the first instants after the Big Bang. On smaller angular scales, CMB anisotropies are influenced by factors that control the expansion rate of the Universe and formation of large-scale structure, such as the cosmological constant, the matter density and the existence and nature of dark matter (e.g. Kolb and Turner, 1994). By measuring the angular power spectrum of CMB fluctuations, one can discriminate among various competing theories that predict the primordial mass distribution (e.g., inflation, cosmic strings and textures, primordial isocurvature baryonic perturbations) and understand the gravitational collapse that ultimately brought about the formation of galaxies. Since the fluctuation amplitudes at angular scales of a few degrees and smaller are also sensitive to the free electron distribution, CMB measurements can also be used to determine the ionization history of the universe.

After the release of the WMAP full-sky data (Bennett et al., 2003), sub-orbital CMB anisotropy experiments are still of high scientific interest as they can improve angular resolution and sensitivity over limited sky regions. The Background Emission Anisotropy Scanning Telescope (BEAST) is the only project currently on-going which is probing a frequency range overlapping with that of WMAP, with improved angular resolution (up to 0.38 degrees at $`40`$ GHz) and potentially better sensitivity over approximately 5% of the sky. The experiment is installed in a conventionally accessible, high altitude site and it has so far accomplished three observing campaigns, on which this paper is based. In this paper we discuss the constraints BEAST places on the power spectrum of CMB anisotropies and its consistency with data taken from a subset of previous experiments (MAXIMA1 (Hanany et al., 2000), TOCO (Miller et al., 1999), BOOMERANG02 (Ruhl et al., 2002), DASI01 (Halverson et al., 2001), VSA1 (Grainge et al., 2002), ACBAR1 (Kuo et al., 2002), CBI (Padin et al., 2001), WMAP (Bennett et al., 2003)).

We present a brief overview of the experiment in §2 and an overview of the estimator in §3. §4 details our implementation of the estimator for the BEAST data and §5 presents the power spectrum and the parameter estimation. We summarize the results in §6.

## 2 The BEAST Experiment

BEAST is a 2.2 meter off axis telescope, currently configured with an 8 element mixed Q (38-45 GHz) and Ka (26-36 GHz) focal plane, and a modulating flat mirror. BEAST was designed as a high altitude balloon system and had two flights: May 20-21, 2000 and October 16, 2000. Subsequent to the second flight BEAST was reconfigured to take advantage of the UC White Mountain Research Station, Barcroft Station at an altitude of 3.8 km in the Eastern Sierra of California. The instrument was fully installed and operational at Barcroft in July, 2001, and took data nearly continuously until December 2001 (except for weather and several equipment failures due to power surges and lightning). Two more weeks of data were obtained in February 2002. A second data taking campaign proceeded in August and September of 2002. The data used for determining the power spectrum presented in this paper are taken from all three of these campaigns.

The data presented in this paper were gathered using the BEAST telescope in a fixed elevation mode. The telescope is kept at a fixed elevation near 90 degrees and the rotation of the Earth provides the map scanning. This strategy results in a sky coverage which forms an annulus centered on the NCP. The annulus is 9 degrees wide and is located between 33 and 42 degrees in declination.

Other aspects of the BEAST experiment are described in the following papers: The instrument is described in Childers et al. (2003) and a more detailed discussion of the optics can be found in Figueiredo et al. (2003). The map-making procedure is described in Meinhold et al. (2003) and constraints on galactic foregrounds in Mejía et al. (2003).

## 3 The MASTER Method

We extract the CMB power spectrum from the BEAST data using the MASTER method, a binned pseudo-$`C_{\mathrm{}}`$ estimator (Wandelt, Hivon & Górski, 2001; Hivon et al., 2002). We chose this estimator for its ease of implementation and the flexibility it offers, which allows testing the analysis with a number of cuts and filtering schemes designed to remove galactic, terrestrial and instrumental foregrounds.

The MASTER method is a de-biasing scheme calibrated against Monte Carlo simulations. Pseudo-$`C_{\mathrm{}}`$ are calculated on the noisy maps over the observed region on the sky with no corrections made for the effect of this cut in terms of the couplings introduced between spherical harmonic modes. The expectation values of these Pseudo-$`C_{\mathrm{}}`$ are modeled in terms of an ansatz which involves, as parameters, an instrumental transfer function $`F_{\mathrm{}}`$ and a noise bias term $`N_{\mathrm{}}`$. These terms are estimated from Monte Carlo simulations of CMB signal and of experimental noise.

The signal and noise are simulated by taking separate random realizations of pure CMB signal and realistic simulations of experimental noise and subjecting them separately to exactly the same data processing (such as beam smoothing, scanning, cuts in the time-ordered data, filtering, template removal and map-making) as the real data.

The power spectra of the resulting signal and noise maps are averaged over the Monte Carlo runs to produce expectation values of the signal-only and noise-only power spectra. These are used to compute the transfer function and noise bias terms in the pseudo-$`C_{\mathrm{}}`$ estimator.

To the extent to which the MASTER ansatz models the expectation values of the pseudo-$`C_{\mathrm{}}`$ and to which our Monte Carlo procedure mimics the acquisition of the real data, we are guaranteed an unbiased power spectrum result.

The experimental data is now passed through the data processing pipeline and the pseudo-$`C_{\mathrm{}}`$ are calculated. Since the experiment covers only a fraction of the sky, a coupling is introduced when performing the spherical harmonic transforms to calculate the power spectra. By calculating the mode-mode coupling kernel for the observed unmasked region on the sky, it is possible to correct for this effect.

Lastly, a binning scheme is chosen in $`\mathrm{}`$ for the final power spectrum and a number of Monte Carlo simulations containing both signal and noise are performed. The covariance matrix of the estimates is calculated by computing the pseudo-$`C_{\mathrm{}}`$ estimator on these simulations. The diagonal elements of the binned covariance matrix are the variances of the binned power spectrum.

## 4 Implementation of MASTER for BEAST

In order to produce an accurate CMB power spectrum from the BEAST data, a detailed knowledge of the experimental beam shape and pointing is required. A residual $`\chi ^2`$ fit of a smoothed delta function to maps of Cygnus A and a best fit to the flux from Cygnus A lead us to characterize the beam as circularly symmetric, with an effective FWHM of $`23^{}\pm 1^{}`$. We use the pointing information reconstructed from a pointing model, which is included in the raw data files, to project our simulations onto the sky in the same manner that the real data is scanned.

A total of 682 individual hours of experimental data are used for the analysis. The data are naturally divided into 55 minute sections by our hourly calibration cycles. These hourly sections are a useful size for several reasons. In addition to the natural delineation by calibrations, 55 minutes is a very manageable size for manipulation in the *IDL* software package on a desktop computer. Also, sky rotation over one hour at our observing angle provides redundant scanning over a nearly symmetric sky patch. The most important effect of this choice of time slices is on our ’template removal’ described below. We tested the sensitivity of our results to varying the timescale of our template removal from the fiducial hour down to a minimum (set by sky rotation) of 600 seconds, and observed no significant changes.

The data has been inspected and spurious signal events, e.g. due to aircraft, have been removed. The data includes both the signal measured by the experiment and the experimental pointing at that instant. This information is used to construct a sky map of the observed signal. For all the maps created in the data analysis we use the HEALPix 111http://www.eso.org/science/healpix/ (Górski, Hivon & Wandelt, 1999) pixelization scheme with an *nside* parameter of 512. This results in a map containing 3,145,728 pixels. Given the size of the experimental beam and the high sampling frequency which is possible with a ground based instrument (450Hz for BEAST), the effects of pixel smoothing are negligible and are ignored here. For the experimental data we create a HEALPix map and calculate the CMB power spectrum using the HEALPix *anafast* package. Further details of the map making process can be found in Meinhold et al. (2003). Fig. 1 shows an overview of the steps in the BEAST simulation and analysis pipeline.

The WMAP (Bennett et al., 2003) best-fit theoretical power spectrum is used to create random realizations of the pure CMB sky. We tested the BEAST pipeline with the power spectra from two fiducial cosmological models and found the final power spectrum to be unchanged by this choice. The first model was a set of reasonable current estimates for cosmological parameters prior to the WMAP data release and the second was the best-fit power spectrum published by the WMAP team.

We scan these signal maps using our experimental pointing strategy read from the time-ordered data (TOD) files. We expect the time-averaged atmospheric contributions to the data to vary with elevation. To remove this foreground we fit a function of elevation angle to the TOD for each hour and subtract it from the TOD samples. Subsequently a 10Hz high pass filter is used. The simulation has now been subjected to exactly the same scanning and filtering as the real BEAST data and we project this simulated data back onto a sky map.

A foreground mask is applied to remove the Galaxy and point source contamination from known sources. We remove from the analysis all pixels with latitude $`b17.5^{}`$. We tested the analysis pipeline with a range of galactic latitude cuts and found that below $`b=17.5^{}`$ there was significant galactic foreground contamination. In addition to this, a separate analysis of the Galactic foregrounds for the BEAST experiment (Mejía et al., 2003) showed that $`b17.5^{}`$ gives an optimal compromise between maximizing the sky fraction observed by the experiment and minimizing the amount of foreground contamination. In this work it was also found that residual Galactic foregrounds outside the mask are small and they are ignored here.

Finally a power spectrum is generated from each signal map and these power spectra are averaged to produce an average signal-only power spectrum.

To construct noise-only maps we subtract our signal estimate for the map from each sample in the experimental TOD and assume that each hourly segment of experimental data is now noise-dominated. We further assume the noise to be piecewise stationary over one hour sections of data and that each one hour noise chunk is independent. We estimate noise power spectra using a windowed FFT on each hourly segment (Press et al., 1986). We are then able to generate synthetic noise simulations which have the same power spectrum as the actual noise from the experiment. We filter the simulated noise TOD in the same manner as for the data and signal simulations and project the noise onto a sky map, then calculate the average noise power spectrum. Comparisons of the data map and the maps created in the simulation pipeline are shown in Fig.2.

Since we have all of the pointing information, we can also create the experimental window function on the sky. This is a simple geometrical construction which is 1 for any HEALPix pixel which the experiment observes and 0 elsewhere. We use this window function to calculate the mode-mode coupling kernel, $`M_{\mathrm{}\mathrm{}^{}}`$, which depends only on the geometry of the observed region of sky. We use the ansatz for the expected pseudo-$`C_{\mathrm{}}`$ which was proposed in Hivon et al. (2002). From the signal-only simulations we can calibrate the transfer function

$$F_{\mathrm{}}=M_{\mathrm{}\mathrm{}^{}}^1C_{\mathrm{}}^{}{}_{s}{}^{}C_{\mathrm{}}^1(B_{\mathrm{}}^2)^1,$$

where $`C_{\mathrm{}}^{}{}_{s}{}^{}`$ are the signal-only pseudo $`C_{\mathrm{}}`$ and $`C_{\mathrm{}}`$ are the best-fit theory $`C_{\mathrm{}}`$ from the WMAP experiment. $`B_{\mathrm{}}`$ is the experimental beam, a Gaussian with FWHM of 23 arcmin in this case. Since the coupling kernel is ill conditioned we use an iterative approach for computing $`M_{\mathrm{}\mathrm{}^{}}^1C_{\mathrm{}}^{}{}_{s}{}^{}`$. The transfer function for the BEAST experiment is shown in Fig. 3.

Now our $`C_{\mathrm{}}`$ estimate is given by

$$\widehat{C_{\mathrm{}}}=\frac{M_{\mathrm{}\mathrm{}^{}}^1\stackrel{~}{C_{\mathrm{}}}\stackrel{~}{N_{\mathrm{}}}}{F_{\mathrm{}}B_{\mathrm{}}^2}$$

where $`\stackrel{~}{N_{\mathrm{}}}`$ are the pseudo-$`C_{\mathrm{}}`$ from the noise Monte Carlo simulations and $`\stackrel{~}{C_{\mathrm{}}}`$ are the pseudo-$`c_{\mathrm{}}`$ from the data.

In practice we use the binned version of the above equation as given in Hivon et al. (2002). The binned mode-mode coupling kernel is shown in Fig. 4.

By averaging the power spectrum over bins in $`\mathrm{}`$ we effectively reduce correlations between the $`C_{\mathrm{}}`$ bins which were introduced by the sky cut and we also reduce the errors on the resulting power spectrum estimator. We have tried different binning schemes and choose a bin width of $`\mathrm{\Delta }\mathrm{}=55`$.

Finally, we create sky simulations by adding the signal and the noise maps, produced as described above. The covariance matrix $`C_{bb^{}}`$ of the binned power spectrum is calculated from these simulations and the diagonal elements give us the error bars on the binned power spectrum estimator. The power spectrum obtained from this process is discussed in the next section.

The code for the BEAST analysis pipeline was written and executed on an IBM SP RS/6000 (*Seaborg*) at the National Energy Research Scientific Computing Center. The code was parallelized using MPI and ran on 640 processors. In order to obtain a stable PS estimate and to estimate our error bars to $``$20% accuracy we required 40 Monte Carlo runs. The operation count for our analysis pipeline scales approximately as $`N_{tod}log(N_{tod})`$ with a large prefactor, where $`N_{tod}`$ is the number of samples in the TOD.

In order to minimize the computational time, we modified the Healpix routines *synfast* (which makes a sky map from a power spectrum) and *anafast* (which calculates the power spectrum from a sky map) so that they only performed analysis on the portion of the sky where BEAST scans. Since the data set read in for the BEAST simulations is $``$ 80GB and the output maps for 40 MC runs are $``$1.7TB, we also implemented compression algorithms for storing the output maps on disk.

## 5 Power Spectrum and Parameter Estimation

The CMB power spectrum extracted from the BEAST data is shown in Fig. 5. The values of the power spectrum are shown in Table 1 . The 1-$`\sigma `$ error bars shown in the figure should be interpreted with some caution. 40 Monte Carlo simulations allow us to calculate these error bars to within 20%, which is sufficient for our purposes here, but more simulations would lead to more accurate error bars. In addition, we use the Monte Carlo simulations to calculate the transfer function ($`F_{\mathrm{}}`$), which is then used to produce the $`C_{\mathrm{}}`$ estimates and we use these same simulations to calculate the error bar on these estimates. Therefore, our estimate of the error bars on the power spectrum is not unbiased and we underestimate the size of these error bars. In calculating our 8 binned $`C_{\mathrm{}}`$ estimates, we effectively compute a binned transfer function $`T_b`$ and a binned noise estimate $`N_b`$ for each bin. We use 40 Monte Carlo simulations of noise to estimate $`N_b`$ and 40 signal simulations to estimate $`T_b`$. Based on the number of degrees of freedom used to produce these 8 binned $`N_b`$ and $`T_b`$, we estimate the bias in the error bar to be approximately 15%, of the same order as our Monte Carlo uncertainty in the errors. However, since this latter effect is a systematic bias, the comparison of the BEAST power spectrum estimates and the resulting parameter estimates to WMAP should be taken as ”worst-case” consistency checks.

A $`\chi ^2`$ comparison of the BEAST data and the WMAP data was performed. For this comparison the WMAP data was assumed to have zero error. We find a $`\chi ^2`$ parameter of 15.02. With 9 degrees of freedom this means a larger value of $`\chi ^2`$ would occur approximately 10% of the time, so the BEAST power spectrum is marginally consistent with the WMAP result. We show the BEAST power spectrum overplotted with the power spectra from several recent experiments in Fig. 6.

After the mean power spectrum was determined, its likelihood was sampled 40 times, producing 40 sample binned power spectra. The likelihood around the power spectrum is not, in general, Gaussian distributed, but through a change of variables - to the log-offset-normal variables of Bond, Jaffe and Knox (BJK parameterization (Knox et al., 1998)) - the distribution can be mapped into one that is much more nearly Gaussian. However, it was found that 40 samples of the power spectrum distribution was too few for a reliable determination of the BJK parameters, and thus it was decided that the power spectrum likelihood would be approximated as Gaussian-distributed. We then calculate the Likelihood $`L`$ of a theoretical power spectrum, $`D_i^{th}`$, as follows:

$`\chi ^2=_{ij}(D_i^{th}D_i^{ob})M_{ij}(D_j^{th}D_j^{ob})`$

$`L=\mathrm{exp}(\chi ^2/2)`$

$`D_i^{ob}C_i^{ob}l(l+1)/2\pi `$

where $`C_i^{ob}`$ is the observed band-power of the i-th bin, and $`M_{ij}`$ is the covariance matrix.

We determined the best-fit (maximally likely) points in parameter space for:

1) BEAST data + WMAP + MAXIMA1,MAT98,BOOMERANG02,DASI01,VSA1,ACBAR1,CBI + Hubble Key Project + Big Bang Nucleosynthesis relation between $`\mathrm{\Omega }_Bh^2`$ and $`He^4`$ mass-fraction ($`Y_p`$) (Huey, Cyburt & Wandelt, 2003) over the parameter space: $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`h`$, $`n_s`$, $`\mathrm{\Omega }_Bh^2`$, $`Y_p`$, $`\tau `$, $`n_t`$, $`r`$

2) BEAST data alone.

For 1) BEAST data + other recent cosmological data, we found the parameter values and errors via a Markov chain approach. Starting from a 30000 point Markov chain previously run with the experiments: WMAP + MAXIMA1,MAT98,BOOMERANG02,DASI01,VSA1,ACBAR1,CBI + Hubble Key Project + BBN $`\mathrm{\Omega }_Bh^2`$-$`He^4`$ relation, the Markov chain was thinned by discarding 99 out of every 100 points. Each point was then weighted by the Beast likelihood. From this weighted point distribution the parameter means and variance matrix were determined. The parameter estimates were taken to be the means, and the parameter errors were taken as square roots of the diagonal elements of the variance matrix.

For 2) BEAST data alone, cosmic parameter space was searched for the maximally likely point by first trying several candidate points, and then applying the Numerical Recipes Amoeba algorithm (Press et al., 1986) to minimize the trial $`\chi ^2`$. The Amoeba algorithm has no inherent minimum scale (similar to adaptive mesh refinement, the resolution increases as necessary, with the precision limited only by the machine floating point arithmetic), and makes no assumptions about the shape of likelihood function.

Once the BEAST-alone best-fit cosmic parameters have been found, we determined the errors in these values. Ideally a method that, again, does not depend on the parameter likelihood function having a particular shape (ie: Gaussian, for example), such as a Markov Chain algorithm, would be used. In this case however, a less computationally costly method can be employed. We determined the errors in the best-fit parameter values by fitting the likelihood function around that point to a multivariate Gaussian. The resulting estimate of the errors is crude, but sufficient to give an overall measure of the dispersion. To the extent to which the likelihoods are approximately Gaussian in the narrowly constrained case (1) we expect these errors to be more accurate. The results of the joint parameter estimation for BEAST plus other experiments are shown in Table 2. We found the BEAST alone parameters to be consistent with these values, although much less well constrained. For example we found $`\mathrm{\Omega }_k1\mathrm{\Omega }_{tot}=0.074\pm 0.070`$ for BEAST alone compared with $`0.014\pm 0.012`$ for the joint estimate.

In order to examine possible future directions for the BEAST experiment we analyzed the effect of increased quantities of data on the power spectrum error bars. A two- and four-fold reduction in the simulated noise were considered, equating to four and sixteen times more data respectively, assuming no improvement in radiometer sensitivity. We found that over the first peak in the power spectrum there was not a significant improvement in the error bars with more data (see Fig 7). This is expected, since in this region we are sample variance limited by the relatively small patch of sky we observe. However, at larger $`\mathrm{}`$ we do see a significant improvement in the power spectrum error bars, up to the point where the experimental beam cuts off around an $`\mathrm{}`$ of 600, when the error bars become large regardless of the amount of data.

## 6 Conclusions

We have presented the angular power spectrum of the cosmic microwave background as measured by the BEAST experiment. We have demonstrated that it is possible to extract cosmological signal from an easily accessible, ground based CMB experiment which is dominated by correlated noise and that the resulting power spectrum and parameter estimation is consistent with previous results.

The MASTER method was successfully implemented and although this method is approximate, it proved to be flexible and robust and produced a power spectrum with less than 1000 CPU hours of computational time. We believe the BEAST CMB dataset to be one of the largest TOD’s analyzed to date and this proved feasible within the MASTER framework. This suggests that the analysis of future, larger CMB datasets (e.g. Planck) should be computationally feasible.

We also analyzed how additional observing time would improve the power spectrum errors and found that significant improvements could be made between $`250l600`$ with additional time. We note that the atmospheric conditions at White Mountain allow for a better than 50% ’good observing’ fraction over the year and that the 26 days of data presented here were limited by funding constraints.

This work was partially supported by the University of Illinois at Champaign-Urbana. This work has been partially supported by the National Computational Science Alliance under grant number AST020003N. This work was funded by NASA grants NAG5-4078, NAG5-9073 , and NAG5-4185 , and by NSF grants 9813920 and 0118297. In addition we were supported by the White Mountain Research Station, the California Space Institute (CalSpace), and the UCSB Office of Research. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. The research described in this paper was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. J.M is supported by FAPESP grants 01/13235-9 and 02/08471-1. T.V. and C.A.W were partially supported by FAPESP grant 00/06770-2. T.V was partially supported by CNPq grants 466184/00-0 and 302266/88-7-FA. CAW was partially supported by CNPq grant 300409/97-4-FA and FAPESP grant 96/06501-4. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is provided by the NASA Office of Space Science. N.F and A.P. were partially supported by CNPq grant number 470531/2001-0. B.D.W. acknowledges the 2003/4 NCSA Faculty Fellowship. Some of the results in this paper have been derived using HEALPix (Górski, Hivon & Wandelt, 1999). We would like to thank Julian Borrill at NERSC for valuable discussions on the computational aspects of this project. We also thank members of the Planck community for stimulating discussions.


text:
# Weather forecasts, Weather derivatives, Black-Scholes, Feynmann-Kac and Fokker-Planck

## 1 Introduction

Weather forecasting is all about predicting what the weather is going to do next, or over the next 10 days. Weather derivatives are a way of insuring oneself against adverse weather conditions. The Black-Scholes equation is a financial model for the pricing of certain kinds of financial contracts. Feynmann-Kac theory is an esoteric result from the study of partial differential equations and stochastic calculus. The Fokker-Planck equation is a result from statistical physics concerning the probability density of the location of a particle. We describe, in an informal way, how these five topics are related.

## 2 Temperature forecasts as Brownian motion

Today is Monday, and we start by considering a weather forecast that attempts to predict the average temperature for next Saturday i.e. a 5 day forecast. To be specific, we will define average temperature, as is convention, to be the midpoint of the minimum and maximum temperatures during the 24 hour period of Saturday. We will consider a single-valued forecast, rather than a forecast of the whole distribution of possible temperatures. How should such a forecast be interpreted in statistical terms? The most obvious interpretation is that the forecast represents the *expectation*, or *mean*, of the distribution of possible temperatures.

Tomorrow we will have a new forecast for Saturday’s temperature, which will then be a 4 day forecast. The new forecast will most likely have a different value from today’s forecast, and will probably be more accurate, at least on average over many trials, since it is predicting less far into the future. Come Wednesday we will have a new forecast again, and so on, until the final forecast is delivered on Saturday, and the final measurements become available at the end of the day Saturday and we can look back and judge how well our forecasts have performed.

We have considered a situation in which we get new forecasts once a day. One could also imagine getting new forecasts more frequently, and could move towards a limit in which new forecast information is arriving continuously. Some of our mathematics will use this limit to simplify the notation, but our intuition will be based mostly on the idea of daily updates.

How do the forecasts we are considering change from day to day? It seems natural that they could go either up or down. Furthermore, if we were sure, on Tuesday, that the Wednesday forecast was going to be higher, then we should really incorporate that information into our Tuesday forecast. Having incorporated all the information we have into the Tuesday forecast, it would seem likely that the Wednesday forecast could be either higher or lower, but would be the same on average. This suggests we could model the changes in the forecast as a random change with mean zero. In statistics such a stochastic process that doesn’t go up or down on average, but jumps randomly, is known as a *martingale*. We offer three explanations for why we think that it is reasonable to model weather forecasts as martingales:

* The intuitive explanation given above: if they were *not* martingales then they would have a predictable component, which contains information that could have been included in the current forecast. We assume that the forecasts are produced by a rational and efficient process, and hence that this information *is* included in the forecast and so it *is* a martingale. We have previously called this assumption the ’efficient forecast hypothesis’, as a parody on the ’efficient market hypothesis’, which is a similar assumption used to justify the modelling of jumps in share prices as being random.
* A more general intuitive explanation that any expectation, of anything, can only change randomly, otherwise it cannot be an expectation.
* A mathematical proof known variously as the ’tower law’ (Baxter and Rennie, 1996) or the ’law of iterated expectations’(Bjork (1998), page 33) that shows that *all* expectations are martingales.

Temperature is very often reasonably close to normally distributed, and so are temperature forecasts. If temperature forecasts are normally distributed, then *changes* in temperature forecasts are also normally distributed. This makes our model for temperature forecasts more specific: they are *gaussian* martingales. In fact, there is only one gaussian martingale and that is Brownian motion, sometimes known as a Wiener process or a continuous random walk. So now we have arrived at the conclusion that weather forecasts can be modelled using Brownian motion. We will write this as the stochastic differential equation (SDE):

$$d\mu =\sigma dB$$
(1)
where $`\mu `$ is the weather forecast, $`B`$ is Brownian motion and $`\sigma `$ is a volatility.

If we knew $`\sigma `$ we could integrate this equation forward in time as an ensemble to give possible future values for forecasts between now and Saturday, and also the possible future values for temperature on Saturday.

So the next question is: what is $`\sigma `$, and how does it vary in time?

## 3 Weather forecast volatility

We start by considering the volatility between Monday and Tuesday. Our *first* model for this 5-day forecast to 4-day forecast volatility is that it is constant at all times of year. However, a cursory investigation of forecasts shows that they show a larger variance in winter, corresponding to the larger variance of temperatures at that time. As a result of this, the *changes* in forecasts from day to day show a larger variance, and so the volatility is larger. This is an entirely predictable effect. Our *second* model for the volatility is therefore that the volatility is deterministic, with a seasonal cycle corresponding to the seasonal cycle of variance of observed temperatures. We have used this model previously in Jewson (2003f).

We will now question our second model of volatility. Is the volatility really completely predictable, as we have assumed above, or does it change from day to day? (to be clear, at this point we mean: does the 5 day forecast to 4 day forecast volatility change as the starting and target day of the forecast move in time…we address the question of whether the 5 day forecast to 4 day forecast volatility differs from the 4 day forecast to 3 day forecast volatility later).

Are not some weather situations more predictable than others, and doesn’t that lead to periods of lower volatility in the forecast? It has been shown that weather forecasts can predict, with some skill, that the variance of the conditional distribution of temperature varies with time (for example, see Palmer and Tibaldi (1988)). Another way of saying this is that the variance of forecast errors varies in time. We have extended this to show that it is possible to predict the volatility with skill in advance (Jewson and Ziehmann, 2003). So from this we see that the volatility shows days to day variations in size which overlay the seasonal changes, and also that part of these variations are predictable.

So where does this leave us? We have a volatility that varies seasonally, in a way that is predictable infinitely far in advance. It also varies on shorter timescales in a way that could be considered completely unpredictable from a long distance in advance, but becomes partly predictable over short time scales using weather forecast models. In fact, however, our own research seems to show that the predictable short time-scale variations in forecast uncertainty are rather small, and do not add much, if anything, to the skill of temperature forecasts (see Jewson et al. (2003) and Jewson (2003e)). As a result, short term changes in the volatility of the forecasts are also small. Because they are small, and for the sake of being able to make some progress with simple mathematics, we are going to ignore them. This means we can write the volatility simply as a fixed function of time.

We now consider variations in the size of the volatility between Monday and Saturday i.e. the variation with lead time (now we are comparing the 5 day forecast to 4 day forecast volatility with the 4 day forecast to 3 day forecast volatility). Are the sizes of the changes in forecasts between Monday and Tuesday any different from the sizes of changes between Friday and Saturday? In other words, does the volatility change as we approach the target day? It is interesting to pose this question to a few meteorologists: we have done so, and found no consistency at all between the answers that we received. Our (albeit brief) data analysis, however, seems to show, rather interestingly, that these volatilities are roughly constant with lead time (Jewson, 2002a). Maybe there is some underlying reason for this, we don’t know. Either way, it certainly simplifies the modelling. Combined with the model of the changes in the forecasts as being Brownian motion, this gives us a linear variation in forecast mean square error, which has been noted many times (and is also more or less unexplained).

To summarise this section: we now have a model that says that a temperature forecast for a single target day follows a Brownian motion with deterministic volatility as we approach the target day. The volatility varies with forecast day (because it depends on the time of year) but does not vary with lead time.

Our arguments in this section have been somewhat ad-hoc: as an alternative one could consider trying to develop models for the volatility based on a detailed analysis of the statistics of the variability of daily temperatures and daily forecasts, and their dependencies. We have not tried this. It may be possible, but we suspect it will be very hard indeed, given how hard it is to model just temperatures alone (see our attempts to solve this much simpler problem in Caballero et al. (2002) and Jewson and Caballero (2003)).

## 4 Forecasts for the monthly mean

We will now consider a slightly different weather forecast. It is Monday 30th November. Tomorrow is the first day of December. We are interested in a forecast for the average temperature in December (defined, for the sake of being a little precise, as the arithmetic mean of the daily averages). Today’s estimate of that average temperature is based on a single valued forecast for the expected temperature that extends out for the next 10 days, followed by estimates of the expected temperature for the remainder of the month based on historical data. Tomorrow (December 1st) we will update this forecast: we will have one more day of forecasts relevant for the month of December, and will need to use one fewer day of historical data. The day after tomorrow we will get an observed value for December 1st (it won’t be the final, quality controlled value, as so is still subject to change, but we will ignore that). We will then use 10 days of forecast from the 2nd to the 11th, and historical data from the 12th to the 31st. The day after that we will get one more day of historical data, our forecast will move forward one more step, and so on. Then at some point the forecast starts to drop off the end, and a few days after that we get the final observed data point for December 31st.

We can now think about how such a forecast changes in time. In fact, the changes in the monthly forecast are made up of sums of changes in daily forecasts, which we have already considered. Since a sum of Brownian motions is a Brownian motion, we conclude that we can model the *monthly* forecast as a Brownian motion too.

The shape of the volatility for this monthly forecast Brownian motion, is, however, a little more interesting than before. When the forecast first starts to impinge upon December the volatility is going to be low. The forecast does not initially have much effect on the monthly mean, and so the day to day changes in the monthly forecast are small. As more and more forecast takes effect the volatility grows. During the bulk of the contract all of the forecast is relevant, and the volatility has a relatively constant value. Finally, at the end of the contract, the volatility starts to reduce as less and less of the forecast is used. What model could we use to represent this ramping up, constant level, and then ramping down, of the volatility? The actual shape of the ramps depends on the sizes of the volatilities for daily forecasts, and of the correlations between changes of the forecasts at different leads. To keep things simple we will make two assumptions: firstly, that the sizes of the volatilities for single day forecasts are constant…the same assumption we made in section 3. This is partly justified by analysis of data. Secondly, we assume that these changes are uncorrelated. This is not really true…the changes in forecasts at different leads do show some, albeit weak, correlations (see Jewson (2002a)). But it is not too far wrong, at least, and allows us to make some progress. In particular, it means that the ramping up and down is given by a straight line in volatility squared. We will call this the ’trapezium’ model. We have described this model in detail in Jewson (2002b) and Jewson (2003f), and used it in Jewson (2003c) and Jewson (2003d).

## 5 Statistical properties of weather forecasts

We have concluded that both daily and monthly forecasts can be modelled using Brownian motion. In the daily case the volatility is constant with lead, but varies seasonally. In the monthly case the volatility squared follows a trapezium shape. We can now derive a number of further results from these models. We start with a consideration of *transition probabilities*.

We first ask the question: if our forecast says $`\mu (t)`$ today (at time $`t`$) then what is the distribution of the different things it might say at some point in the future (at time $`T`$)? This can be solved by integrating equation 1 forward in time. This gives:

$$\mu N(\mu (t),_t^T\sigma ^2𝑑t)$$
(2)
i.e. the distribution of future values of our forecast is a normal distribution with mean $`\mu (t)`$ and variance given by the integral of the daily volatilities squared. Alternatively, the probabilities themselves as a function of the value of the forecast $`x`$ and time $`t`$ are given by solutions of the Fokker-Planck equation (also known as the Kolmogorov forward equation):

$$\frac{p}{t}=\frac{^2}{x^2}(\sigma ^2p)$$
(3)
along with a boundary condition based on what the forecast currently says. This equation is a diffusion type of equation. The probabilities diffuse outwards as time moves forwards.

Another related question is: if our forecast is going to say $`x`$ in the future, what are the possible things it could be saying today? These probabilities are given by the *backward* Fokker-Planck equation (also known as the Kolmogorov backward equation):

$$\frac{p}{t}=\sigma ^2\frac{^2p}{x^2}$$
(4)
along with a boundary condition that specifies the final condition. This equation is diffusion in reverse i.e. starting at the final forecasts, the probabilities diffuse outwards as time moves backwards.

For more details on these equations, see, for example, Gardiner (1985) or Bjork (1998).

## 6 The Feynmann-Kac formula

We now take a diversion from discussing weather forecasts, and present the Feynmann-Kac formula. In subsequent sections we will then see that this formula is related to the fair price of weather options.

Consider the following partial differential equation (PDE) for $`F(x,t)`$:

$$\frac{F}{t}+\mu \frac{F}{x}+\frac{1}{2}\sigma ^2\frac{^2F}{x^2}=0$$
(5)
with the boundary condition:

$$F(x,T)=\mathrm{\Phi }(x)$$
(6)
This PDE is very ”special”: because of the particular nature of the terms in the equation we can solve it by taking expectations of a function of Brownian motion. The Brownian Motion we must use is fixed by the coefficients in the PDE:

$$dX=\mu dt+\sigma dB$$
(7)
with the initial condition:

$$X=x$$
(8)
and the expectation we must calculate is

$$F(x,t)=E[\mathrm{\Phi }(X_T)]$$
(9)
What we see from this is that one interpretation of what the solution of the PDE actually does is to somehow take expectations of the final condition, at different points in time, and according to some stochastic process. The proof of this relation is easy, and is given in, for example, Bjork (1998).

## 7 Weather derivative fair prices are martingales

We will now link our discussion of weather forecasts to the pricing of weather derivatives. Weather derivatives are financial contracts that allow entities to insure themselves against adverse weather. They are based on a weather index, such as monthly mean temperature in December, and have a payoff function which converts the index into a financial amount that is to be paid from the seller of the contract to the buyer of the contract. The payoff function can in principle have any form, but a certain number of piecewise linear forms known as swaps, puts and calls are most common (see, for example, Jewson (2003a)).

The fair price of a weather derivative is defined, by convention, as the expectation of the distribution of possible payoffs from the contract. The fair price can form the basis of negotiations of what the actual price should be. Typically the seller of the contract (who sells such contracts as a business) would want to charge more than the fair price, otherwise they could not expect to make money on average. The buyer may be willing to pay above the fair price because they are grateful for the chance to insure themselves against potentially dire weather situations.

We can calculate the fair price using whatever relevant historical data and forecasts are available to estimate the final distribution of outcomes of the contract. This is the usual approach (details are given in Jewson et al. (2002) and Brix et al. (2002)). Alternatively, we can calculate the fair price by integrating equation 1 foreward to the end of the contract to give us the final distribution of settlement indices.

As we move through a contract, our estimate of the fair price will change in time as the available historical data and forecasts change. The fair price is an expectation, by definition, and we have argued in section 2 that expectations are martingales, so the fair price should change as a martingale.

## 8 Deriving the fair price PDE

It is useful to think of the conditional estimate of the fair price at some point during a weather contract as being dependent on the conditional estimate of the expected index, rather than, say, on the temperature so far during the contract. The former is much more convenient in many situations because the expected index is a Brownian motion while the temperature is very highly autocorrelated and hard to model. Thus we have:

$$V=V(\mu ,\sigma _x)$$
(10)
where $`\mu `$ and $`\sigma _x`$ are the conditional mean and standard deviation of the settlement index.

Differentiating this with respect to time gives:

$`dV`$ $`=`$ $`{\displaystyle \frac{V}{\mu }}d\mu +{\displaystyle \frac{V}{\sigma }}d\sigma +{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2V}{\mu ^2}}d\mu ^2+\mathrm{}`$
$`=`$ $`\delta d\mu +\zeta d\sigma +{\displaystyle \frac{1}{2}}\gamma d\mu ^2+\mathrm{}`$
$`=`$ $`\delta \sigma dB+\theta dt+{\displaystyle \frac{1}{2}}\gamma \sigma ^2dB^2+\mathrm{}`$
$`=`$ $`\delta \sigma dB+\theta dt+{\displaystyle \frac{1}{2}}\gamma \sigma ^2dt+\mathrm{}`$
$`=`$ $`\delta \sigma dB+dt(\theta +{\displaystyle \frac{1}{2}}\sigma ^2\gamma )+\mathrm{}`$
$``$ $`\delta \sigma dB+dt(\theta +{\displaystyle \frac{1}{2}}\sigma ^2\gamma )`$
This last expression is the Ito derivative, and we have used the standard definitions $`\delta =\frac{V}{\mu }`$, $`\zeta =\frac{V}{\sigma }`$, $`\theta =\frac{V}{t}`$ and $`\gamma =\frac{^2V}{\mu ^2}`$ (see Jewson (2003b)). We see that changes in $`V`$ are driven by stochastic jumps (the $`dB`$ term) and by a deterministic drift (the $`dt`$ term).

But we have already seen in section 7 that $`V`$ is a martingale, and hence that there can be no drift term. Thus we must have that the coefficient of $`dt`$ in this equation is zero, giving:

$$\theta +\frac{1}{2}\sigma ^2\gamma =0$$
(12)
or, re-expanding in terms of the full notation:

$$\frac{V}{t}+\frac{1}{2}\sigma ^2\frac{^2V}{\mu ^2}=0$$
(13)
We conclude that the fair price of a weather option satisfies a PDE, which is a backwards diffusion-type of equation. The diffusion coefficient comes from the volatility of weather forecasts.

### 8.1 Relation to Feynmann-Kac

There is a close relation to what we have just derived, and the Feynmann-Kac theory discussed in section 6. Equation 13 is a particular example of equation 5. Applying the Feynmann-Kac theorem, we see that we can solve this equation by integrating a stochastic process and taking an expectation. The stochastic process is given by equation 1. We have come full circle.

## 9 The stochastic process for the fair price and the VaR

Given the PDE for the fair price (equation 13) equation 8 now simplifies to:

$$dV=\delta \sigma dB$$
(14)
In other words, changes in the fair price of a weather derivative over short time horizons are normally distributed around the current fair price, and have a volatility given simply in terms of the $`\delta `$ of the contract, and the $`\sigma `$ of the underlying expected index.

Holders of financial contracts often like to know how much they could lose, and how quickly. Minimising, or limiting, the amount one could lose is known as risk management. One of the most common measures of risk used in risk management is the market value at risk, or market VaR. Market VaR is defined to be the 5% level of the distribution of possible market prices of a contract at a specified future time. So far we haven’t considered market prices, just fair value, so we will define *actuarial VaR* as the 5% level of the distribution of possible *fair* values of a contract at a specified future time.

Equation 14 actually gives us the actuarial VaR over short time horizons <sup>2</sup><sup>2</sup>2this is the useful result we mention in the abstract. The distribution of changes in the fair value is given by:

$$dVN(0,\delta \sigma )$$
(15)
and the VaR is a quantile from this distribution.

This equation doesn’t apply over finite time horizons because both $`\delta `$ and $`\sigma `$ change with time. We could integrate equation 14 forward in time to get around this and derive the distribution of $`dV`$ over finite time horizons. However, very small time steps are needed because $`\delta `$ can change very rapidly in certain situations. It is generally easier to calculate actuarial VaR over longer horizons by integrating the underlying process (equation 1) to give a distribution of future values of the expected index (which can be done analytically) and *then* calculating the distribution of values of $`V`$. These issues are discussed in more details in Jewson (2003c).

## 10 Arbitrage pricing, Black-Scholes and Black

Equity options are contracts that have a payoff dependent on the future level of some share price. By trading the shares themselves very frequently one can more or less replicate the payoff structure of an equity option. The cost of doing this replication tell us what the value of the option should be. This is the Black-Scholes theory of option pricing (Black and Scholes, 1973). If the replication is done using forward contracts on the equity rather than the equity itself then the Black-Scholes theory must be modified as described by Black (1976) (although if we set interest rates to zero, as we will, the two theories are the same).

We can apply the Black-Scholes argument to weather options, as described in Jewson and Zervos (2003). <sup>3</sup><sup>3</sup>3thanks to Anna-Maria Velioti for pointing out that equation 28 in this paper has a sign error in the last term We imagine hedging the weather option using weather forwards. The final equation is slightly different from the Black equation, because the underlying process we use to model the expectation of the weather index is slightly different from the processes that are used to model share prices. In the case where interest rates are zero we find that the Black equation for weather is the same as the PDE for the fair price, already given in equation 13.

To summarise: in an actuarial pricing world equation 13 gives us the fair price for weather options. In a Black-Scholes world (where frequent trading is possible) equation 13 gives us *both* the fair price *and* the market price.

## 11 Notes on the relation between Black-Scholes and Weather Black-Scholes

In the standard Black-Scholes world one way of expressing the price of an option is as the discounted expected payoff with expectations calculated *under the risk neutral measure*. A curious feature of the weather derivative version of the Black-Scholes model is that the change of measure is not necessary, since the underlying process (given by equation 1) does not have any drift. The objective measure is the risk neutral measure.

## 12 Summary

We have argued that temperature forecasts for a fixed target day can be modelled as Brownian motion with deterministic volatility. We have also argued that temperature forecasts for a fixed target month can be modelled as a Brownian motion, but with a more complex volatility structure. Standard results for Brownian motion thus apply such as the forward and backward Fokker-Planck equations, and equations for transition probabilities.

We have also argued that the fair price of a weather derivative must be a martingale. As a result, the fair price satisfies both a PDE and an SDE, and the SDE can be used to estimate actuarial VaR over short time horizons.

Finally, under additional assumptions about the weather market the PDE for the fair price is the weather derivative equivalent of the Black-Scholes and Black equations.

## 13 Legal statement

The author was employed by RMS at the time that this article was written.

However, neither the research behind this article nor the writing of this article were in the course of his employment, (where ’in the course of his employment’ is within the meaning of the Copyright, Designs and Patents Act 1988, Section 11), nor were they in the course of his normal duties, or in the course of duties falling outside his normal duties but specifically assigned to him (where ’in the course of his normal duties’ and ’in the course of duties falling outside his normal duties’ are within the meanings of the Patents Act 1977, Section 39). Furthermore the article does not contain any proprietary information or trade secrets of RMS. As a result, the author is the owner of all the intellectual property rights (including, but not limited to, copyright, moral rights, design rights and rights to inventions) associated with and arising from this article. The author reserves all these rights. No-one may reproduce, store or transmit, in any form or by any means, any part of this article without the author’s prior written permission. The moral rights of the author have been asserted.


text:
# Pressure dependence of magnetism in URu2Si2

(29 August 2003; printed )

## Abstract

Neutron-scattering and specific-heat measurements of the heavy-fermion superconductor URu2Si2 under hydrostatic pressure and with Rh-doping \[U(Ru0.98Rh0.02)2Si2\] show the existence of two magnetic phase transitions. At the second-order phase transition $`T_m`$17.5 K, a tiny ordered moment is established, while at $`T_M<T_m`$, a first-order phase transition (under pressure or doping) gives rise to a large moment. The results can be understood in terms of a hidden OP $`\psi `$ coupled to the ordered moment $`m`$, where $`m`$ and $`\psi `$ have the same symmetry.

The phase transition at $`T_m`$=17.5 K in URu2Si2 is one of the most puzzling features in heavy-fermion systems. Most bulk properties display anomalies at $`T_m`$, and in particular there is a huge jump of $`\mathrm{\Delta }C/T_m`$=0.3 JK-2mol-1 in the specific heat palstra . These anomalies cannot be explained by the tiny ordered moment of $``$0.03 $`\mu _B`$, which also develops at $`T_m`$. The antiferromagnetic (AFM) structure, observed both by neutron broholm ; walker ; fak ; mason ; santini ; FRL  and x-ray isaacs  scattering, is characterized by a propagation vector k=(001) and dipolar moments along the $`c`$ axis of the body-centered tetragonal structure (I4/mmm).

In a recent tour de force experiment combining neutron elastic scattering with high magnetic fields FRL , it was found that the magnetic moment in URu2Si2 did not disappear at 15 T as suggested by the extrapolation of earlier low-field data mason ; santini , but remains finite at 17 T with an inflection point near 7 T . The particular field dependence of the ordered moment, $`m(H)`$, agrees with the predictions of a model containing two order parameters (OPs) with the same symmetry shah , based on the theoretical work of Ref. agterberg .

In this Letter, we present new experimental data where hydrostatic or chemical pressure is used to tune the interactions in the system. The main result is the existence of two distinct phases in the $`pT`$ phase diagram of URu2Si2 for pressures $`p>p_M`$, where $`p_M4.9`$ kbar (see Fig. 1), or for small Rh doping. The existence of two phases have also been inferred from thermal expansion measurements sato , but, as will be discussed later, their phase diagram is fundamentally different from ours. We call the two magnetic phases for small-moment (SMAF) and large-moment (LMAF) antiferromagnetic phases. The SMAF phase, observed below the second-order phase transition at $`T_m17.5`$ K ($`T_m(p)`$ in Fig. 1 is determined from resistivity measurements under pressure, Ref. schmidt ), is the same as that observed at $`p`$=0 and is characterized by a small moment $`m`$ and large anomalies in the macroscopic properties. The low-T LMAF phase, observed below the transition $`T_M`$, which seems to be a first-order transition, is characterized by a large moment $`m`$ and small anomalies in the macroscopic properties. The SMAF and LMAF phases have the same magnetic structures.

Many phenomenological theories for URu2Si2 introduce two physical quantities or “order parameters”: an AFM parameter $`m`$ that describes ordering of localized dipolar 5$`f`$ moments and a hidden (unknown) order $`\psi `$ responsible for the $`T_m`$=17.5 K transition shah ; agterberg . It has been suggested that $`\psi `$ could originate from triple-spin correlators of 5$`f`$ electrons agterberg  or an exotic spin-density wave. The Landau free energy for the two OP’s (up to fourth-order terms) can be written as

$$F=\alpha _\psi \psi ^2+\alpha _mm^2+2\gamma \psi m+\beta _\psi \psi ^4+\beta _mm^4+2\beta _i\psi ^2m^2$$
(1)
The coupling term $`\psi m`$ is allowed only if $`m`$ and $`\psi `$ transform according to the same irreducible representation, otherwise $`\gamma 0`$. Coupling terms of type $`\psi ^3m`$ and $`\psi m^3`$ can be removed from Eq. (1) by linear transformation. Shah et alshah  studied the above free energy at zero pressure and finite magnetic fields $`H`$ and predicted an inflection point for $`m(H)`$ only in the case $`\gamma 0`$. This was recently observed experimentally FRL . Motivated by this agreement, we investigate here the functional Eq. (1) as a function of pressure.

For $`\gamma 0`$, the two OP’s appear simultaneously below the transition line from the paramagnetic state $`T_m(p)`$, given by $`\alpha _\psi \alpha _m=\gamma ^2`$. At $`p`$=0, the AFM OP behaves as $`m(\gamma /\alpha _m)\psi `$ (for $`|\alpha _\psi |,|\gamma |\alpha _m`$). A small coefficient $`\gamma /\alpha _m`$ implies weak ordered moments, while $`\psi `$ gives rise to a large anomaly in the specific heat at $`T_m`$, as in the SMAF phase. The presence and nature of extra transitions below $`T_m(p)`$ depend on quartic terms in Eq. (1). Since the transition at $`T_M`$ appears to be first order, we consider the case $`\beta _i>\sqrt{\beta _\psi \beta _m}`$.

The $`p`$$`T`$ phase diagram for $`\gamma `$=0 is sketched in the left panel of Fig. 2. The line of first-order phase transitions $`T_M(p)`$ between states with $`m`$=0 and $`\psi `$=0 emerges from the kink point on $`T_m(p)`$. The equation for $`T_M(p)`$ is $`\alpha _\psi ^2\beta _m=\alpha _m^2\beta _\psi `$. Such a phase diagram was recently considered by Chandra et alchandra  in their interpretation of small AFM moments as a phase separation effect.

The recently observed inflection point in $`m(H)`$ FRL , however, points to an intrinsic origin of the tiny ordered moment and a finite $`\gamma `$. For $`\gamma 0`$, the first-order transition line $`T_M(p)`$ is stable and its position in the $`p`$$`T`$ plane is given by the same equation as for $`\gamma `$=0. However, $`T_M(p)`$ splits from the line of second order transitions $`T_m(p)`$ and terminates at a critical point $`p_c`$ below $`T_m`$ (right panel of Fig. 2). The location of the critical point is given by $`\alpha _\psi ^c=2|\gamma |(\beta _\psi ^3\beta _m)^{1/4}/(\beta _i\sqrt{\beta _\psi \beta _m})`$. The two states above and below $`T_M(p)`$ are phases with large $`\psi _L`$ and small $`m_S`$ (SMAF) and with small $`\psi _S`$ and large $`m_L`$ (LMAF), respectively. A relative jump of the ordered AFM moments across $`T_M(p)`$ is given by $`(m_Lm_S)/(m_L+m_S)=\sqrt{(\alpha _\psi \alpha _\psi ^c)/(\alpha _\psi +\alpha _\psi ^c)}`$. The size of the jump varies continuously along $`T_M(p)`$ and vanishes at $`p=p_c`$. The corresponding phase diagram, shown schematically in the right panel of Fig. 2, resembles the vapor-liquid phase diagram.

The first evidence for the existence of two phase transitions were obtained from neutron scattering and specific-heat measurements on a Rh doped sample: U(Ru0.98Rh0.02)2Si2these . The $`T`$ dependence of the ordered magnetic moment (see inset of Fig. 3) clearly shows two transitions at $`T_M`$=8.3 K and $`T_m`$=13.2 K, with saturation moments (extrapolated to $`T0`$) of 0.25 and 0.04 $`\mu _B`$, respectively. The transition at $`T_M`$ has a first order shape (although no hysteresis was observed) while that at $`T_m`$ is second order. These measurements were performed on the DN3 neutron diffractometer at the Siloé reactor using a wavelength of 1.54 Å and on the thermal triple-axis spectrometer 2T at the Orphée reactor with a final wave vector $`k_f`$= 2.662 Å-1 on a 1.2 g annealed Czochralski-grown single crystal of U(Ru0.98Rh0.02)2Si2.

Specific heat measurements (see Fig. 3) show that while the large jump in the specific heat is associated with a small moment, the small jump is associated with a large moment estimated at 0.16 $`\mu _B`$ calcp . This shows that while there is no relation between the small moment $`m_S`$ and $`\mathrm{\Delta }C/T_m`$, the large moment $`m_L`$ is consistent with the specific heat anomaly at $`T_M`$. The entropy jump of the LMAF transition deduced from these measurements corresponds to $`\mathrm{\Delta }S(T_M)`$ = 23 mJK-1mol-1. The effect of Rh doping is more subtle than a simple (negative) pressure, since Rh is not isoelectronic with Ru. The added electron modifies the band structure, and while the volume increases slightly, the low-$`T`$ lattice parameter $`a`$ decreases and $`c`$ increases, as shown by X-ray measurements on the ID31 beam line at the ESRF id20 . Using the Clausius-Clapeyron relation with $`\mathrm{\Delta }V`$ deduced from the $`pure`$ URu2Si2sato , a large slope $`dT_M/dp`$ is found, characteristic of a first-order phase transition.

Experiments under hydrostatic pressure on the pure compound show strong similarities with the Rh-doped sample. Neutron-scattering measurements were performed at the Institut Laue-Langevin (ILL) on the D15 diffractometer using a wavelength of 0.85 Å and on the IN22 triple-axis spectrometer using a final wave vector of 2.662 Å-1. A 4.8-mm diameter crystal was cut from the same batch used in the high-field measurements of Ref. FRL  with the $`a`$ axis along the 4-mm long cylinder axis. Standard ILL CuBe and steel clamp pressure cells of inner diameter of 5 or 6 mm were used for most measurement with fluorinert as pressure transmitter. Complementary measurements were made with an ILL He-gas pressure cell. On D15, the pressure was determined from the lattice spacing of a NaCl crystal mounted inside the pressure cells.

The observed temperature dependence of the magnetic intensity ($`m^2`$) of URu2Si2 is shown in Fig. 4 for different pressures. The increase of the temperature slope, $`dm^2(T)/dT`$, with pressure near $`T_M`$ indicates that the phase transition is first order. This behavior is opposite to that expected for a pressure gradient if the transition is second order, and is thus an intrinsic effect. Fitting the expression $`m^2=m_0^2[1(T/T_M)^\alpha ]`$ fak  leads to $`\alpha `$2.5 for pressures below 5 kbar, while a much larger value, $`\alpha `$8.5, is obtained at higher pressures, suggestive of a first order transition. The phase diagram presented in Fig. 1 is constructed from the data shown in Fig. 4, where $`T_M`$ is taken as the midpoint of the transition, which is the correct estimate for a first-order transition in the presence of pressure gradients. The presence of pressure gradients, usually estimated to $`\mathrm{\Delta }p/p5`$% for the used clamp cells, is clear from the strong intensity observed at 4.5 kbar using a clamp cell (Fig. 4), since measurements in a He-gas pressure cell at $`p=4.45\pm 0.15`$ kbar show no (large) ordered moment on D15, where the detection limit is 0.06 $`\mu _B`$. We also note that the magnetic intensity saturates at a value of 0.33 $`\mu _B`$ for pressures $`p6.4`$ kbar. The critical pressure for the onset of the large moment is estimated to $`p_M=4.9\pm 0.2`$ kbar. As for the doped Rh sample (see inset of Fig. 3), no hystereris was observed in $`m^2(T)`$. However, the vertical slope of $`dT/dp`$ in Fig. 1 is a strong evidence of a first-order transition. The absence of magnetic scattering at the (0,0,2$`l`$+1) reflections shows that the ordered moments are along the $`c`$-axis at all pressures, i.e. the SMAF and LMAF phases have the same AFM structure. We finally note that while the AFM Bragg peaks are resolution limited in the LMAF phase amiprl ; id20 , they have a finite correlation length of a few hundred Å in the SMAF phase broholm ; fak ; amiprl ; id20 .

The longitudinally polarized magnetic excitations were studied under pressure on IN22 at the two wave vectors Q1=(1,0,0) and Q2=(1.4,0,0), which correspond to the minima in the $`p`$=0 dispersion curve of URu2Si2 broholm . Spectra for energy transfers up to 7 meV were performed with $`k_f`$=1.55 Å-1 and a 5-cm long Be filter in the scattered beam to reduce higher-order contamination; the energy resolution was 0.3 meV. Complementary spectra for energy transfers up to 16 meV were performed with $`k_f`$=2.662 Å-1 and a 4-cm long pyrolytic graphite filter in the scattered beam; the energy resolution was 0.9 meV.

The measurements show pronounced differences in the excitation spectra between the LMAF and SMAF phase. In the SMAF phase, $`T_M<T<T_m`$, the two excitations at a finite pressures of 5.0 kbar are very similar to those at zero pressure, with the same gap energies of $`\mathrm{\Delta }_1`$=1.6 and $`\mathrm{\Delta }_2`$=4.5 meV at Q1 and Q2, respectively \[see Figs. 5(a-b)\]. In the LMAF phase, however, the excitation at Q1 is absent \[see Fig. 5(c)\]. There is no clear feature at this wave vector in the measured spectra for energy transfers between 1 and 16 meV, in agreement with Ref. metoki . This suggests that the Q1 excitation is closely related to $`\psi `$, which is small in the LMAF phase. The excitation at Q2, on the other hand, is readily seen. Its energy depends on pressure: we found $`\mathrm{\Delta }_2`$=7.7, 8.6, and 9.2 meV at pressures of 5, 9, and 10.5 kbar, respectively.

Our neutron scattering measurements on the 2% Rh-doped sample reveal that the excitations (not shown) are very similar to the pure compound under pressure. In the SMAF phase, the excitations at both Q1 and Q2 are present. In the LMAF phase, the excitation at Q1 is absent, as in URu2Si2 under pressure, while the excitation at Q2 is still present.

We will now discuss our results in the context of other work. Recent thermal-expansion measurements also observe a first-order transition under pressure sato . The authors construct a phase diagram where the first-order transition line appears to join the second-order line around a pressure of 11 kbar. This would imply that the coefficient $`\gamma `$ in Eq. (1) is zero, in clear contradiction with high-field measurements FRL . Our data, cf. Fig. 1, clearly show that the first and second order lines are nearly parallel at these pressures, and the most likely scenario is that the first-order line has an endpoint (cf. right panel of Fig. 2), and hence $`\gamma 0`$.

NMR matsuda  and $`\mu `$SR amiuSR  measurements have also observed a large moment under pressure, i.e. in the LMAF phase. The coexistence of an AFM signal with a “paramagnetic” signal in the NMR measurements under pressure and a magnetic field of a few Tesla on an aligned but crushed powder sample by Matsuda et al. matsuda  was interpreted in terms of a spatially inhomogenous AFM ordering. There are several arguments against such a scenario. First, it would be a curious coincident if the AFM and the hidden order would have the same transition temperature in a phase-separated scenario, since $`m`$ and $`\psi `$ are not coupled as they occur in different parts of the sample. Secondly, one would not expect the Q1 excitation to disappear at a pressure where the two phases still coexist, just a smooth intensity change as the population of the two phases changes with pressure. Thirdly, a phase-separated scenario is inconsistent with a number of other experiments, such as the recent de Haas-van Alphen measurements under pressure nakashima  and the NMR measurements by Bernal et al. bernal . Since the temperature dependence observed in Ref. matsuda  is suggestive of either a large pressure gradient or strong internal strains, it would be important to repeat 29Si-NMR measurements under pressure on a well-characterized single crystal. We recall here that the magnetism of URu2Si2 is very sensitive to the sample quality fak . Indeed, the neutron scattering measurement under pressure on an annealed but otherwise uncharacterized crystal (closely related to that of Ref. matsuda ) presented in Ref. amiprl  are partly inconsistent with our neutron scattering experiments, the thermal-expansion measurements of Ref. sato , and their own $`\mu `$SR work amiuSR .

In conclusion, both the $`p`$-$`T`$ phase diagram (this work) and the magnetic-field dependence of the ordered moment FRL  in URu2Si2 are in excellent agreement with a scenario where a hidden OP is linearly coupled to the ordered moment and hence breaks time-reversal symmetry.

We acknowledge discussions with M.B. Walker and I. Fomin and experimental support from A.A. Menovsky, A.D. Huxley, E. Ressouche, J.L. Laborier, F. Fauth, L. Paolasini, A. Bombardi, F. Yakhou and N. Pyka. We thank G. Motoyama and N.K. Sato for sending us complementary thermal expansion data.


text:
# Order parameters with higher dimensionful composite fields

## Abstract

We discuss the possibility of the spontaneous symmetry breaking characterized by order parameters with higher dimensionful composite fields. By analyzing general Ginzburg-Landau potential for a complex scalar field $`\varphi =\varphi _1+\mathrm{i}\varphi _2`$ with $`O(2)`$ symmetry, we demonstrate that a phase characterized by $`\varphi _1^2\varphi _2^20`$ with $`\varphi _1=\varphi _2=0`$ is realized in a certain parameter region. To clarify the driving force to favor this phase, we study the $`O(2)`$ $`\varphi ^6`$ theory in three dimensions.

preprint: MIT-CTP 3457

The spontaneous symmetry breaking (SSB) should be enumerated as one of the key concepts in modern physics, in particular in condensed matter and elementary particle physics nam61 . Even if some global symmetry is manifest in the Hamiltonian, the ground state and excited spectra do not have to reflect the symmetry manifestly. The order parameter, which is a measure of SSB, is defined by the ground state expectation value of an operator being variant under the symmetry transformation. The choice of the operator is constrained by the symmetry breaking pattern. The purpose of this Letter is to draw attention to the issue of SSB through such order parameters with higher canonical dimensions. As we will argue, unbroken discrete groups generally require operators with higher canonical dimensions.

Let us first give a summary to construct generic order parameters. Suppose $`X`$ and $`Q_a`$ be a symmetry group and the generators of the symmetry transformation, respectively. We take an operator $`𝒪`$ which transforms non-trivially under $`Q_a`$; $`[\mathrm{i}Q_a,𝒪]0`$. SSB takes place if we have such an $`𝒪`$ that $`0|[\mathrm{i}Q_b,𝒪]|00`$ for some $`b`$ 111Precisely speaking, $`0|[\mathrm{i}Q_b,𝒪]|0`$ is ill-defined as it is when SSB occurs. Nevertheless, $`dx0|[\mathrm{i}j_b,𝒪]|0`$ with $`Q_b=dxj_b`$ can be finite and make sense.. If $`Q_b`$ is a generator of $`X`$ but is outside of its subgroup $`Y`$, the symmetry is spontaneously broken from $`X`$ to $`Y`$ with an order parameter $`0|[\mathrm{i}Q_b,𝒪]|0`$. $`𝒪`$ can be chosen to be any operators as long as it is invariant under $`Y`$.

In particular, if $`Y`$ contains a discrete subgroup of $`X`$, SSB generally requires higher dimensionful composite operators for $`𝒪`$ because the elementary operators are not invariant in $`Y`$. We pursue such exotic realization of SSB not only for academic interest but for phenomenological importance. In Quantum Chromodynamics (QCD) we sometimes encounter possibilities of higher dimensionful operators. A well-known example is the four-quark operator introduced for a non-standard chiral symmetry breaking pattern in $`N_\text{f}`$-flavors ste98  where $`X=SU_\text{L}(N_\text{f})\times SU_\text{R}(N_\text{f})`$ with $`Y=SU_\text{V}(N_\text{f})\times Z_\text{A}(N_\text{f})`$. Although this is proven to be forbidden in QCD at zero and finite temperature, it may be realized at finite baryon chemical potential. Interestingly enough, the remaining $`Z_\text{A}(N_\text{f})`$ symmetry could change the hadron spectrum significantly gir01 . Exotic hadron spectrum in a similar natureson00  is indeed known to be realized in the color-flavor-locking (CFL) phase in color superconductivity where $`U_\text{B}(1)`$ is broken down to $`Z(2)`$ subgroup alf99 . In particular, the existence of light kaonic excitations associated with this SSB opens a possibility for the kaon condensation in the CFL phase bed02 . Therefore, in connection with these phenomenological concerns, it is an important issue to probe SSB with an unbroken discrete subgroup and with higher dimensionful order parameter from the generic point of view.

In order to look into the problem in a concrete setting, let us take an $`O(2)`$ symmetric model with one complex scalar field $`\varphi =\varphi _1+\mathrm{i}\varphi _2`$. If we take $`𝒪=\mathrm{Im}\varphi =\varphi _2`$ which is not symmetric in any of the $`O(2)`$ transformation, we have $`\delta \varphi _2\mathrm{Re}\varphi =\varphi _1`$. It follows that the order parameter for the symmetry breaking pattern $`X=O(2)\text{none}`$ is $`\varphi _10|\varphi _1|0`$. This is, however, not a unique symmetry breaking pattern possible. Consider $`X=O(2)`$ and $`Y=Z_2`$ with $`Z_2`$ being a discrete subgroup with the operation; $`\varphi _{1,2}\varphi _{1,2}`$. In this case $`𝒪=\mathrm{Im}\varphi ^2=2\varphi _1\varphi _2`$, which is $`Y`$ invariant, leads to an order parameter,

$$\mathrm{Re}\varphi ^2=\varphi _1^2\varphi _2^2.$$
(1)
Note that the difference between $`\varphi _1^2`$ and $`\varphi _2^2`$ is an appropriate order parameter instead of the individuals because the common constant cancels out in the former. One can construct $`k`$-th composite order parameters for $`O(2)Z_k`$ in the same way;

$$\mathrm{Re}\varphi ^k.$$
(2)
If $`\mathrm{Re}\varphi `$ has a non-vanishing value, the $`O(2)`$ symmetry is totally broken so that the higher dimensionful order parameters, $`\mathrm{Re}\varphi ^k`$, are non-zero for all $`k1`$. The SSB with $`O(2)Z_k`$ can be characterized by the condition; $`\mathrm{Re}\varphi =\mathrm{}=\mathrm{Re}\varphi ^{k1}=0`$ and yet $`\mathrm{Re}\varphi ^k0`$.

Application of the above argument to $`O(N)`$ symmetric models with $`N>2`$ is straightforward. For the SSB pattern, $`O(N)O(2)Z_k`$, the higher dimensionful order parameter can be defined in the same way as above with respect to the $`O(2)`$ subgroup.

Now we shall go on to investigate whether there is a situation where the above SSB pattern, $`O(2)Z_k`$, is indeed realized. We will focus on a simplest case, $`k=2`$, namely the phase characterized by $`\varphi _1^2\varphi _2^2`$. For later purpose, we introduce complex variables $`W`$ and $`H`$, and a real variable $`F`$ as follows,

$$\begin{array}{cc}& W=\varphi ^2=\varphi _1^2\varphi _2^2+2\mathrm{i}\varphi _1\varphi _2,\hfill \\ & F=\varphi ^{}\varphi =\varphi _1^2+\varphi _2^2,\hfill \\ & H=\varphi =\varphi _1+\mathrm{i}\varphi _2.\hfill \end{array}$$
(3)
To clarify the phase structure in the $`W`$-$`F`$-$`H`$ space, let us consider a Ginzburg-Landau type effective potential with $`O(2)`$ symmetry up to terms of $`\varphi ^6`$ order which are at least necessary to induce a phase with $`W0`$;

$`V_{\text{eff}}[W,F,H]`$
$`=(d_1+d_2F+d_3F^2)|H|^2+(d_4+d_5F)|H|^4+d_6|H|^6`$
$`+(c_1+c_2F)(W^{}H^2+H^2W)`$
$`+a_1F+a_2F^2+a_3F^3+(b_1+b_2F)|W|^2.`$ (4)
The coefficients $`a_i`$, $`b_i`$, $`c_i`$, and $`d_i`$ are considered to be arbitrary coupling constants at the present stage. Taking a specific model with $`O(2)`$ symmetry, one may determine their magnitudes and relations. For more than two dimensions, one needs to subtract the short distant singularities in $`F`$ and $`W`$ to make them finite. This will be discussed in more detail later and we assume a simple cutoff in the short distant part here.

The standard phase corresponding to the SSB pattern, $`O(2)\text{none}`$, is characterized by $`H0`$ (and thus $`W0`$), while the non-trivial phase corresponds to $`W0`$ with $`H=0`$. Note that $`F0`$ in both cases. For convenience, we call the latter phase as “the WFH phase” after the notation in Eq. (3). To find the condition for the WFH phase, we choose $`d_i`$ appropriately so as to realize $`H=0`$. Then the GL potential is reduced to only the terms containing $`a_i`$ and $`b_i`$ in Eq. (4). Among these coupling constants, we can eliminate two of them by rescaling $`F`$ (and $`W`$) and coupling constants. Then one finds

$$V_{\text{eff}}[W,F]=a_1^{}F\pm F^2+F^3+(b_1^{}+b_2^{}F)W^2,$$
(5)
where $`W`$ is chosen to be real and positive without loss of generality and as a consequence we have a constraint, $`FW0`$. Note that the sign of the $`F^3`$ term should be positive to guarantee the stability of the potential, while it is not necessary for the coefficient of the $`F^2`$ term.

An instability toward $`W0`$ is induced by the negative $`W^2`$ term in Eq. (5). Therefore, once the instability to the WFH phase occurs, $`W=F`$ is realized because of the constraint $`FW0`$. To see the competition between the symmetric phase ($`W=0`$) and the WFH phase ($`W=F0`$), we compare the potential energies at two possible global minima as

$`V_{\text{eff}}[W=0,F]`$ $`=aF\pm F^2+F^3,`$ (6)
$`V_{\text{eff}}[W=F,F]`$ $`=aF+bF^2+cF^3,`$ (7)
where $`a=a_1^{}`$, $`b=\pm 1+b_1^{}`$, and $`c=1+b_2^{}`$.

The gap equation for $`F`$ is obtained from $`V_{\text{eff}}/F=0`$. Eq. (6) has a global minimum at $`F_0>0`$, i.e., $`V_{\text{eff}}[F_0,0]<V_{\text{eff}}[0,0]=0`$ is satisfied for $`a<0`$ (for $`+F^2`$) and for $`a<1/4`$ (for $`F^2`$). On the other hand, Eq. (7) has a global minimum at $`F_1>0`$ either for ($`a>0`$, $`b<0`$, $`b^2>4ac`$) or for $`a<0`$. Then the condition $`V_{\text{eff}}[F_1,F_1]<V_{\text{eff}}[0,F_0]`$ yields an inequality;

$$\begin{array}{cc}& \left\{(13a)^{3/2}+\frac{(b^23ac)^{3/2}}{c^2}\right\}\hfill \\ & +\left(\pm 1\frac{b^3}{c^2}\right)\frac{9}{2}\left(\pm a\frac{ab}{c}\right)>0.\hfill \end{array}$$
(8)
Shown in Fig. 1 is the above condition for the negative $`F^2`$ term in Eqs. (5) and (6). The region below the surface corresponds to the parameters where the WFH phase can exist. When the $`F^2`$ term is positive, it is sufficient to consider the potential only up to the $`F^2`$ term in order to guarantee the stability of the potential, and the WFH phase is realized for $`a<0`$ and $`0<b<1`$. The message here is that there is always a wide parameter region in the Ginzburg-Landau type effective potential so that the WFH phase is realized.

To study microscopic mechanism to induce the instability toward $`W0`$ in field theoretical models, let us consider an $`O(2)`$ $`\varphi ^6`$ model in three spatial dimensions with the Hamiltonian density,

$$=\frac{1}{2}|\varphi |^2+\frac{g_2}{2}|\varphi |^2+\frac{g_4}{4}|\varphi |^4+\frac{g_6}{6}|\varphi |^6.$$
(9)
We assume $`g_6>0`$ for stability of the system. If the symmetry group is $`O(4)`$ instead of $`O(2)`$, the model is considered as the 3d effective theory to describe the tricritical behavior of the chiral phase transition in QCD at finite temperature and baryon density step98 .

We treat the model in the Cornwall-Jackiw-Tomboulis (CJT) formalism cor74  which is best suited for studying the system with composite order parameters. The effective action is given in general by

$$\mathrm{\Gamma }[G,H]=I_0[H]+\frac{1}{2}\mathrm{Tr}\mathrm{ln}G^1+\frac{1}{2}\mathrm{Tr}G_0^1[H]G+\mathrm{\Gamma }_2[G,H],$$
(10)
where $`I_0[H]`$ and $`G_0[H]`$ are the tree-level potential and propagator, respectively. $`G`$ is the full propagator satisfying $`\delta \mathrm{\Gamma }/\delta G=0`$. $`\mathrm{\Gamma }_2`$ represents the contributions from 2 particle irreducible (2PI) diagrams with $`G`$.

Taking into account the terms up to leading order 2PI diagrams, which is equivalent to the Hartree-Fock (HF) approximation in many-body theories, $`G`$ can be written by the dynamical mass $`m_a`$ of the field $`\varphi _a`$. Then the expectation value of the local composite operator $`G_{ab}(x,x)=\varphi _a\varphi _b`$ reads

$`S_{ab}`$ $`G_{ab}(x,x)=\delta _{ab}{\displaystyle ^\mathrm{\Lambda }}{\displaystyle \frac{\mathrm{d}^3p}{(2\pi )^3}}{\displaystyle \frac{1}{p^2+m_a^2}}`$
$`=\delta _{ab}{\displaystyle \frac{1}{2\pi ^2}}\left\{\mathrm{\Lambda }\right|m_a|\mathrm{arctan}\left({\displaystyle \frac{\mathrm{\Lambda }}{|m_a|}}\right)\},`$ (11)
where $`\mathrm{\Lambda }`$ is the three dimensional cutoff of the loop integral. The term proportional to $`\mathrm{\Lambda }`$ in Eq. (11) corresponds to the short distance part independent of $`m_a`$. Thus, the WFH phase, in which $`\varphi _1^2\varphi _2^2=S_{11}S_{22}0`$ is realized, can be characterized by $`m_1m_2`$ within the present approximation.

From now on we take $`H=0`$ by hand as before to focus on the phase with $`W0`$. The dynamical masses, $`m_a`$, are determined by the gap equation derived from the variation of the CJT potential or from the self-consistent HF equation;

$`{\displaystyle \frac{m_a^2}{2}}={\displaystyle \frac{g_2}{2}}+{\displaystyle \frac{g_4}{4}}\left({\displaystyle \underset{b}{}}S_b+2S_a\right)`$
$`+{\displaystyle \frac{g_6}{6}}\left\{\left({\displaystyle \underset{b}{}}S_b\right)^2+4S_a{\displaystyle \underset{b}{}}S_b+2{\displaystyle \underset{b}{}}S_b^2+8S_a^2\right\},`$ (12)
where $`S_aS_{aa}`$ (see Eq. (11)). Each term in Eq. (12) has the graphical representation as shown in Fig. 2.

The genuine Hartree terms in the right hand side of Eq. (12) or the fully disconnected graphs in Fig. 2 are index-blind and give equivalent contributions to $`m_1`$ and $`m_2`$. Hence, if we take only Hartree terms, the standard solution, $`m_1=m_2`$, is obtained. On the other hand, the Fock terms in the right hand side of Eq. (12) are connected to the external lines and depend on the external index $`a`$. This leads to a possibility for $`m_1m_2`$.

Of course, the above argument does not guarantee the realization of an asymmetric solution, because $`m_1=m_2`$ is also a solution even if we have Fock terms. To study which solution is more stable, we have calculated the CJT effective potential written in terms of $`f=|m_1|+|m_2|`$ and $`w=|m_1||m_2|`$. For small $`m_a`$, they are related to $`F=S_{11}+S_{22}`$ and $`W=S_{11}S_{22}`$ as

$$\begin{array}{cc}& F=\frac{1}{2\pi ^2}\left(2\mathrm{\Lambda }\frac{\pi }{2}f+O(|m_a|^2)\right),\hfill \\ & W=\frac{1}{2\pi ^2}\left(\frac{\pi }{2}w+O(|m_a|^2)\right).\hfill \end{array}$$
(13)
For convenience, we use $`w`$ and $`f`$ as basic variables in the following.

Assuming that $`m_a`$ are small enough compared to $`\mathrm{\Lambda }`$, we expand the CJT potential in terms of $`f`$ and $`w`$ up to the same order with Eq. (5);

$`\overline{V}_{\text{eff}}[w,f]`$ $`={\displaystyle \frac{\overline{g}_2}{2}}f+{\displaystyle \frac{\overline{g}_4}{2}}f^2+(1\overline{g}_6)f^3`$
$`+{\displaystyle \frac{1}{2}}\left\{{\displaystyle \frac{\overline{g}_4}{2}}+3(2\overline{g}_6)f\right\}w^2.`$ (14)
Note that $`w^2`$ terms come from the Fock contributions and the bare coupling constants and the potential are redefined as

$$\begin{array}{cc}\hfill \overline{g}_6& =\frac{3}{2\pi ^2}g_6,\overline{g}_4=\frac{24}{\pi }g_4+\frac{36}{\pi ^3}g_6,\hfill \\ \hfill \overline{g}_2& =24g_2\frac{48}{\pi ^2}g_4\frac{144}{\pi ^4}g_6,\hfill \\ \hfill \overline{V}_{\text{eff}}& =96\pi V_{\text{eff}}\hfill \end{array}$$
(15)
with the cutoff $`\mathrm{\Lambda }`$ being set to 1. The stability of the system at large $`f`$ leads to a condition, $`\overline{g}_6<1`$ (or equivalently $`g_6<2\pi ^2/3`$ bar84 ). Although the definition of the order parameter is slightly different from the general analysis in Eq. (5), obtained effective potential has the same structure. The correspondence becomes even transparent if we rescale $`\overline{V}_{\text{eff}}`$ and $`f`$ (and $`w`$) to eliminate two of the coefficients in Eq. (14) to obtain

$$V_{\text{eff}}[w,f]=a_1f\pm f^2+f^3+\left(\pm \frac{1}{2}+b_2f\right)w^2,$$
(16)
where $`b_2>3`$, which comes from $`\overline{g}_6<1`$. In the following, we choose the negative signs in the $`f^2`$ and $`w^2`$ terms in Eq. (16) because the $`\frac{1}{2}w^2`$ term, which originates from the Fock contribution, induces the instability toward the WFH phase.

The condition that $`V_{\text{eff}}[w=m_1,f=m_1]`$ (the WFH phase) has at least a local minimum with respect to $`m_1`$ leads to $`a_1(b_2+1)<3/4`$. This is shown by the left hand side of the dotted curve in Fig. 3. The stability in the $`m_2`$ direction requires $`3(b_2+1)+4a_1b_2(b_21)+\sqrt{3}\sqrt{34a_1(b_2+1)}>0`$, which is shown by the right hand side of the solid curve in Fig. 3. As we mentioned, $`b_2>3`$ is necessary for the stability of the potential at large $`f`$. As a result, the WFH phase is allowed at least as a local-minimum in the gray region in Fig. 3. It can be, however, shown that it is not a global minimum but a meta-stable phase. In fact, by comparing Eq. (16) with $`f=w`$ and Eq. (7), one finds $`a=a_1`$, $`b=3/2`$, and $`c=1+b_2>4`$. This parameter set is too restrictive and can not satisfy the condition Eq. (8) which is necessary for the absolute stability.

Even though it is only a meta-stable state, an important message here is that the 3d $`O(2)`$ $`\varphi ^6`$ theory in the HF approximation embodies a driving force (the Fock term) leading toward the WFH phase. Whether this new phase is indeed realized as a true ground state or not should be examined in an approach beyond the HF approximation or by numerical simulations. In $`O(N>2)`$ theories, the Fock term is suppressed relative to the Hartree term for large $`N`$. Therefore the WFH phase driven by the Fock term is expected in small $`N`$ instead of in large $`N`$.

In summary, we have studied higher dimensionful order parameters associated with the dynamical symmetry breaking which leaves manifest discrete symmetry. By taking a complex scalar field $`\varphi `$ and its Ginzburg-Landau potential as an example, we have given a general criterion to have a novel symmetry breaking pattern $`O(2)Z_2`$ characterized by the order parameter $`\mathrm{Re}\varphi ^20`$. We have demonstrated that the Hartree-Fock approximation to the 3d $`O(2)`$ $`\varphi ^6`$ theory gives such a novel phase at least as a meta-stable state. Whether the new phase exists as a true ground state beyond the HF approximation is an open question to be studied. We have assumed the stability around $`\varphi =0`$ in this Letter; further studies of the interplay among the three condensates, $`\varphi `$, $`|\varphi |^2`$, and $`\varphi ^2`$ should be also pursued. The lesson we can learn from our results in this Letter is that SSB through higher dimensionful order parameters is not so peculiar and may possibly be realized in more complex systems. Applications of the idea to e.g. gauge field theories and condensed matter systems are interesting future problems to be studied.

Acknowledgments: Y. W. is grateful to Ayumu Sugita for his crucial suggestions. K. F. is supported by Japan Society for the Promotion of Science for Young Scientists. This work is partially supported by the Grants-in-Aid of the Japanese Ministry of Education, Science and Culture (No. 15540254) and the U.S. Department of Energy (D.O.E.) under cooperative research agreement #DF-FC02-94ER40818.