format: text

text:
# X-wave mediated instability of plane waves in Kerr media

## Abstract

Plane waves in Kerr media spontaneously generate paraxial X-waves (i.e. non-dispersive and non-diffractive pulsed beams) that get amplified along propagation. This effect can be considered a form of conical emission (i.e. spatio-temporal modulational instability), and can be used as a key for the interpretation of the out of axis energy emission in the splitting process of focused pulses in normally dispersive materials. A new class of spatio-temporal localized wave patterns is identified. X-waves instability, and nonlinear X-waves, are also expected in periodical Bose condensed gases.

The dynamics of focused femtoseconds pulses (FFP) in optically nonlinear media has a fundamental importance and it is relevant in all of the applications of ultrafast optics. The basic mechanism, when dealing with the propagation in normally dispersive materials, is the splitting of the pulse, which has been originally predicted more than ten years ago Chernev and Petrov (1992); Rothenberg (1992) and recently reconsidered, due to the development of the physics of FFP. Ranka et al. (1996); Berge et al. (1996); Brabec and Krausz (1997); Trippenbach and Band (1997); Ranka and Gaeta (1998); Diddams et al. (1998); Zozulya and Diddams (1999); Litvak et al. (2000); Eisenberg et al. (2001); Berge et al. (2002); Ward and Berge (2003)

If $`z`$ is the direction of propagation, $`t`$ is the time in the reference frame where the pulse is still, and $`r=\sqrt{x^2+y^2}`$ the radial cylindrical coordinate, this process can be roughly divided into a series of steps. They describe the propagation of a gaussian (in space and time) pulse, which travels in a focusing medium, and undergoes relevant reshaping, due to the interplay of diffraction, dispersion and the Kerr effect (when the peak power is sufficiently greater than the critical power $`P_c`$ for self-focusing Chernev and Petrov (1992); Rothenberg (1992)): 1) the out of axis energy is focused towards $`r=0`$ around $`t=0`$; 2) the pulse at $`r=0`$ is compressed; 3) lobes appear in the on-axis temporal spectrum; 4) the pulse splits in the time domain.
Before the breakup, relevant out of axis energy emission and re-distribution occurs, as originally described by Rothenberg in Ref. Rothenberg, 1992. Looking to the spatio-temporal profile an X-shape (or hyperbola) is observed before the splitting (see, e.g., figures in Litvak et al. (2000); Berge et al. (2002)). This process has been theoretically described by different approaches Luther et al. (1994); Fibich et al. (1995); Litvak et al. (2000); Berge et al. (2002), and the onset of an X-shape can be ultimately related to the hyperbolic characteristics over which small perturbations evolve. Levi-Civita (1931) This can be checked, for example, by the hydrodynamical approach to the nonlinear Schrödinger equation.Staliunas (1993); Stringari (1996)

Recently, attention has been devoted to the existence of X-waves in optically nonlinear media. The latter are self-trapped (i.e. non-diffracting and non-dispersive) waves that are well known in the field of linear propagation in acustics Ziolkowski et al. (1989); Lu and Greenleaf (1992); Stepanishen and Sun (1997); Salo et al. (1999) and in electromagnetism Durnin et al. (1987); Recami (1998); Saari and Reivelt (1997); Sonajalg et al. (1997); Reivelt and Saari (2000); Mugnai et al. (2000); Sheppard (2001); Zamboni-Rached et al. (2002). The simplest X-wave has the shape of a double-cone, or clepsydra, that appears as an X when, e.g., a section is plotted in the plane $`(x,t)`$, and as V in the plane $`(r,t)`$, as shown in figure 1. Optical X-waves in nonlinear media have been theoretically predicted in Conti et al. (2002a) and experimental results, with a direct observation of the conical spatio-temporal shape, have confirmed their existence and generation in crystals for second-harmonic generation (SH-G) Valiulis et al. (2001); Conti et al. (2002b); Trapani et al. (2003). In Ref. Conti and Trillo, 2002 it has been theoretically shown how, during non-depleted-pump SHG, the phase matched spatio-temporal harmonics let the SH beam become an X-wave. 111For the sake of concreteness, I use the term “X-wave”, without making distinction among the many families of localized-wave patterns solutions to the linear Helmotz equation. See,e.g., Ref. Saari, 2001.

Conical emission (CE), or spatio-temporal modulational instability (MI) (Liou et al. (1992); Luther et al. (1994); Fonseca et al. (1999); Picozzi and Haelterman (2002)), has been addressed in Ref. Trillo et al., 2002 as a basic mechanism underlaying the spontaneous formation of an X-wave and, as a foundation for the understaing of the splitting, in Ref. Litvak et al., 2000. In this Letter I consider the formation of X-waves in Kerr media (or in quadratic media, in the regime where they mimic cubic nonlinearity, see, e.g., the chapter after Torruellas, Kivhsar and Stegeman in Trillo and Torruellas (2001) and Ref. Buriak et al., 2002). A new form of instability of plane waves can be introduced by directly involving self-localized spatio-temporal wave patterns. The process strongly resembles CE, i.e. the amplification of plane waves from noise, but in this case X-waves, instead of periodical patterns, emerge from the breakup of an unstable pump beam. This mechanism is responsible of the first stage of the pulse-splitting, i.e. the out-of-axis energy redistribution, such as MI breaks a continuos wave signal into a periodical pattern of solitons.Hasegawa (1984)

The wave equation describing the propagation in nonlinear Kerr media, whose refractive $`n`$ index obeys the law $`n=n_0+n_2I`$, with $`I`$ the optical intensity, in the framework of the paraxial and of the slowly varying envelope approximation, is written as

$$i_ZA+ik^{}_TA+\frac{1}{2k}_{XY}^2A\frac{k^{\prime \prime }}{2}_{TT}A+\frac{kn_2}{n_0}|A|^2A=0\text{,}$$
(1)
where $`A`$ is normalized such that $`|A|^2=I`$, and I have taken $`k^{\prime \prime }>0`$ (i.e. the medium is normally dispersive). $`X,Y,Z,T`$ are the real world variables, $`\lambda `$ the wavelength, $`n(\lambda )`$ the refractive index and $`k(\lambda )=2\pi n(\lambda )/\lambda `$. Eq. (1) can be casted in the adimensional form:

$$i_zu+\mathrm{\Delta }_{}u_{tt}u+\chi |u|^2u=0\text{,}$$
(2)
by defining $`z=Z/Z_{df}`$ with $`Z_{df}=2kW_0^2`$ and $`W_0^2`$ a reference beam waist; $`\mathrm{\Delta }_{}_{xx}+_{yy}`$ the tranverse Laplacian with $`(X,Y)=W_0(x,y)`$; $`t=(TZ/V_g)/T_0`$ the retarded time in the frame travelling at the group velocity $`V_g=1/k^{}`$ in units of $`T_0=(k^{\prime \prime }Z_{df}/2)^{1/2}`$. The optical field envelope $`A`$ is given by $`A=A_0u`$, with $`A_0=(n_0/k|n_2|Z_{df})^{1/2}`$. $`\chi >0`$ ($`\chi <0`$) identifies a focusing (de-focusing) medium being $`n_2>0`$ ($`n_2<0`$).
I start considering paraxial linear (i.e. $`\chi =0`$) X-wave solutions of (2) (paraxial X-waves have been considered in detail in Ref. Porras et al., 2002) , defined by

$$\psi (\mathrm{\Delta }_{}_{tt})\psi =0\text{.}$$
(3)
Introducing the complex variable $`v=(\mathrm{\Delta }it)^2+r^2`$, with $`\mathrm{\Delta }`$ a real valued arbitrary coefficient, I have from (3) the equation

$$6_v\psi +4v_{vv}\psi =0\text{,}$$
(4)
from which $`\psi =C_1/\sqrt{v}+C_2`$. $`C_1`$ and $`C_2`$ are arbitrary complex coefficients. Note that both the real and the imaginary parts of this solution, are real-valued X-wave profiles. The former being the simplest X-wave, given by (the branch cut for the square root is along the negative real axis)

$$\psi _X\mathrm{}(\frac{1}{\sqrt{v}})=\mathrm{}(\frac{1}{\sqrt{(\mathrm{\Delta }it)^2+r^2}})\text{.}$$
(5)
Remarkably $`\psi _X`$ still holds when referring to the Helmotz equation, instead of the paraxial wave equation (see e.g. Ref. Saari, 2001). Its plot is given in figure 1.

X-wave instability can be introduced in the same way as MI, i.e. by perturbing the plane-wave solution of (2): $`a=a_Pa_0\mathrm{exp}(ia_0^2z)`$, with $`a_0`$ a real-valued constant. Letting $`a=(a_0+ϵ(x,y,t,z))\mathrm{exp}(ia_0^2z)`$ I have, at first order in $`ϵ`$:

$$i_zϵ+ϵ+\chi a_0^2(ϵ^{}+ϵ)\text{.}$$
(6)
Writing $`ϵ=ϵ(r,t,z)=\psi _X(r,t)\mu (z)`$ gives

$$i_z\mu +\chi a_0^2(\mu +\mu ^{})=0\text{.}$$
(7)
The perturbation is thus

$$ϵ=[\alpha +i(\beta +2\chi a_0^2\alpha z)]\psi _X(r,t)$$
(8)
with $`\alpha `$ and $`\beta `$ arbitrary real-valued constants. Equation (8) represents an X-wave which grows linearly along propagation (independently on the sign of $`\chi `$), with amplification given by the intensity of the plane wave $`a_P`$. This situation strongly resembles MI, where plane waves are exponentially amplified at the expense of the pump beam, with a gain determined by the pump intensity. For this reason, it is natural to refer to this process as X-wave instability. As for MI, the amplified wave can be artificially externally feeded, or it can be generated by noise. 222 It is remarkable that both X-waves and MI plane waves, $`\mathrm{exp}(i\omega t+ik_xx+ik_yy)`$, fill all the space-time and have infinite energy. Furthermore X-waves may be enriched by orthogonality properties as described in Ref. Salo and Salomaa, 2001  Note also that X-wave instability can be triggered by conical emission. Indeed, as shown in Trillo et al. (2002), the latter generates the required spectrum to form an X-wave, which eventually gets amplified, as discussed above.

The previous treatment can be generalized in several ways. In particular it is possible to show that exponential amplification of X-wave-like beams can be attained. It is necessary to enlarge the definition of X-waves, i.e. Eq. (3), by introducing the equation:

$$\psi =\kappa \psi \text{,}$$
(9)
with $`\kappa `$ a real constant ($`\kappa 0`$ in the following). Letting $`ϵ=\mu (z)\psi (r,t)+\nu (z)^{}\psi (r,t)^{}`$ in (6), and setting to zero the coefficients of $`\psi `$ and $`\psi ^{}`$ the following linear system is obtained:

$$\begin{array}{c}i_z\mu +\kappa \mu +a_0^2\chi (\mu +\nu )=0\\ i_z\nu +\kappa \mu +a_0^2\chi (\mu +\nu )=0\text{.}\end{array}$$
(10)
If $`(\mu ,\nu )=(\widehat{\mu },\widehat{\mu })\mathrm{exp}(\gamma z)`$ I have

$$\left[\begin{array}{cc}i\gamma +\kappa +\chi a_0^2& \chi a_0^2\\ \chi a_0^2& i\gamma +\kappa +\chi a_0^2\end{array}\right]\left[\begin{array}{c}\widehat{\mu }\\ \widehat{\nu }\end{array}\right]=0\text{.}$$
(11)
The solvability condition yelds the allowed values for the gain $`\gamma =\pm \sqrt{\kappa (\kappa +2\chi a_0^2)}`$. The perturbation grows along $`z`$ if the following inequality is satisfied: $`\kappa (\kappa +2\chi a_0^2)>0`$. Thus, for a focusing (de-focusing) medium, generalized X-waves with $`2a_0^2<\kappa <0`$ ($`0<\kappa <2a_0^2`$) are exponentially amplified. The gain Vs $`\kappa `$, $`\gamma =\sqrt{\kappa (\kappa +2\chi a_0^2)}`$, plotted in figure 2, has a maximum value which corresponds to the most exploding self-localized packet.

Now the question arises as to whether or not Eq. (9) admits solutions resembling X-waves. Closed forms can be found proceeding as before: In term of $`v`$ Eq.(9) becomes:

$$6_v\psi +4v_{vv}\psi =\kappa \psi \text{,}$$
(12)
whose general solution is

$$\psi =C_1\frac{\mathrm{exp}(\sqrt{\kappa v})}{\sqrt{v}}+C_2\frac{\mathrm{exp}(\sqrt{\kappa v})}{\sqrt{v}}\text{.}$$
(13)
A real valued, localised, generalised X-wave is given by

$$\psi _\kappa =\mathrm{}(\frac{\mathrm{exp}(\sqrt{\kappa v})}{\sqrt{v}})\text{.}$$
(14)
(14) depends on two parameters $`\mathrm{\Delta }`$ and $`\kappa `$; while the first determines the decay constant as going far from the origin in the $`(r,t)`$ plane, the latter completely changes the shape of the wave. Examples for $`\mathrm{\Delta }=1`$ and $`\kappa =1`$ are shown in figure 3 and in figure 4 for $`\mathrm{\Delta }=1`$ and $`\kappa =1`$. Note that they are similar to the Bessel pulse beams described in Sheppard (2001).

To clarify the differences, I observe that the spatio-temporal spectrum develops around the curve $`\omega ^2=k_{}^2+\kappa `$, with $`\omega `$ the angular frequency corresponding to $`t`$ and $`k_{}`$ the transversal wavenumber. In figure 5 I show the spectra for different $`\kappa `$, with $`\kappa =0`$ corresponding to the simplest X-wave $`\psi _X`$. The appearance of this lines in the spatio-temporal far field is a clear signature of the X-waves instability, and can be directly observed in the experiments. Note that the spectra resemble the shape of the wave in the physical space as a consequence of their propagation invariance.

To show that X-wave instability can actually be observed in the experiments I numerically solved Eq. (2). I considered a gaussian input beam, $`A=A_0\mathrm{exp}(R^2/(2W_0^2)T^2/(2T_p^2))`$, whose intensity profile FWHM spot and duration are $`70\mu m`$ and $`100fs`$. Eq. (2) is integrated with reference to fused silica, with $`\lambda =800nm`$ ($`n_0=1.5`$,$`n_2=2.5\times 10^{20}`$ $`k^{\prime \prime }=360^{28}`$, SI units ), peak power $`P=1.5P_c`$, being $`P_c=(0.61\lambda )^2\pi /(8n_0n_2)2.6MW`$ the critical power for self-focusing. In figure 6 I show the spatio-temporal profile and the spectrum (in log scale) after $`3`$ diffraction lengths $`L_{df}=\pi n_0W_0^2/\lambda `$. Clearly an X-like profile is formed and the spectrum shows the features in Fig. 5.

In conclusion I have shown that a plane wave in Kerr media gives rise to linear and exponential amplification of X-waves, thus leading to a significat beam reshaping. A new class of X-waves is involved in this nonlinear process. The reported analysis provides insights for the interpretation of the pulse-splitting of focused femtosecond beams, and related phenomena, in the same way as MI is relevant for solitons generation. Indeed, the spostaneous formation of an X-wave can be another explanation of the out of axis energy redistribution typically observed. X-wave instability can also be an effective approach for the controlled generation of non-dispersive and non-diffractive pulses.

These results appear to be susceptible of several generalizations, as considering quadratric nonlinearity or vectorial effects, and have implications in all the fields encompassing nonlinear wave propagation, as acustics, hydrodynamics and plasma physics. For example, X-wave instability (as well as nonlinear X-waves) is also expected in periodical Bose-Einstein Condensates where Eq. (2) holds, being $`t`$ the direction of periodicity, in the presence of negative effective mass. Konotop and Salerno (2002); Kramer et al. (2002)

###### Acknowledgements.

I thanks S. Trillo and E. Recami for fruitful discussions, and the Fondazione Tronchetti Provera for the financial support.


text:
# 1 Introduction

## 1 Introduction

The ‘knee’, a rather sharp steepening in the primary cosmic ray ( CR ) energy spectrum at about 3 PeV, was inferred from the observation of a similar feature in the measured size spectrum of extensive air showers by Kulikov and Khristiansen (1958). The knee is commonly asserted to be due to an increasing failure of ‘Galactic containment’ of the CR generated by sources within the Galaxy, the containment being caused by the magnetic fields in the interstellar medium ( ISM ), However, it is the firmly-held view of the present authors that the knee is too sharp for this explanation and we have advanced what we claim to be a more realistic model. This is our ‘single source ( SS ) model’ (see Erlykin and Wolfendale, 1997, 2001b for recent details) which comprises cosmic ray acceleration up to the knee energy by supernova remnants, the knee itself being due to the truncation that occurs at 3 PeV for oxygen nuclei from a single, recent, nearby SNR. The other main accelerated nucleus at these energies is iron and its termination occurs at about 12 PeV where, it is claimed by us that there is a small second knee (when the spectrum is plotted as $`E^3I(E)`$ vs. $`E`$, the knees appear as small peaks). The remainder of the CR spectrum ( at least to some 109 GeV, or so ) is presumed due to ‘super’-SNR and other sources and their spatial distribution is such as to give a comparatively smooth spectrum in the PeV region.

Erlykin and Wolfendale (to be referred to henceforth as EW) have examined a variety of other cosmic ray data and concluded that there is either support for the model or that the data are neutral. Very recently, low energy gamma ray data have also been studied (EW, 2002) and the well-known ‘gamma ray excess’ in the Inner Galaxy, and deficit in the Outer Galaxy, have been explained in terms of propagation differences dependent on the conditions in the ISM from which the gamma rays come. The results relate to Galaxy-wide properties and, although the SNR acceleration hypothesis has been invoked, there is no significant information about the single source.

It is at higher gamma ray energies where potential problems exist (e.g. Drury et al., 1994). Most recently, Bhadra (2002) has argued that the single source should be visible in TeV gamma rays, and it is not. This is the topic to be addressed here. We use the results found in a very recent paper (EW, 2003a to be referred to as I), where we made predictions of the fluxes and the angular distribution of gamma ray intensity from SNR of different ages and at different distances from the Sun. The threshold energies were taken as 0.1 GeV and 1 TeV.

A critical feature of the Bhadra estimate was the ‘normalisation’ of the SNR ‘conditions’ so as to give the CR energy density created by the single source at earth. We regard this as a legitimate procedure and we follow this path, although other features of our model differ considerably from those adopted by Bhadra. Our calculations are thus not simply a ‘re-run of the Bhadra calculations with different values for the parameters’ but, rather, for what is certainly a more appropriate model of SNR acceleration and ( less certainly, perhaps ) a significantly different model of cosmic ray propagation.

## 2 The Bhadra estimate

Bhadra’s model is rather straightforward, in principle, at least: particles are accelerated by the SNR shock to give a differential spectrum $`AE^\gamma `$, with $`\gamma =2.0,`$ and these particles interact with the ambient (or swept-up) gas of density $`n`$, with $`n=1`$ cm-3. At this stage it is necessary to make critical remarks, however. Implicit in the Bhadra calculations is the assumption that the SNR shock accelerates the CR instantaneously, at ‘$`t=0`$’. Although this can be used to give a viable mathematical model, such a situation is certainly not appropriate to a real SNR where the acceleration occurs over an extended period: $`810^4`$ years in our model, and little different in other SNR acceleration models.

Nevertheless, we continue to describe the Bhadra calculations. It is assumed that the particles diffuse from the source a distance $`r`$ from the Sun in a normal, Gaussian fashion with diffusion coefficient $`D=10^{29}`$ cm2s-1. The parameters are chosen to give the required CR energy density at the Earth created by the single source. For instance, if the source is at the distance of 300 pc and it is 104 years old the total energy transferred from SN to CR is required to be 1.9$``$1050 erg. The Table summarizes the most important parameters.

The expected minimum gamma ray flux above 0.1 GeV rises with the age of the SN from $`0.210^7\mathrm{cm}^2\mathrm{s}^1`$ at 103 years to $`210^7\mathrm{cm}^2\mathrm{s}^1`$ at 105 years and from $`0.410^{11}\mathrm{cm}^2\mathrm{s}^1`$ to $`410^{11}\mathrm{cm}^2\mathrm{s}^1`$ above 1 TeV, respectively. For the real candidates for the single source discussed in EW, (1997) ( eg. Loop I, Clayton SNR ) the expected fluxes are substantially higher. Comparing these fluxes with the diffuse gamma ray background Bhadra found that for the present gamma ray telescopes it should have been possible to observe the single source. Since there has been no claim for an observation Bhadra concluded that the single source cannot be such a SNR.

Turning to our remarks on the validity of Bhadra’s model, in addition to the basic problem with the assumption about the instantaneous acceleration there are two further reasons why we cannot allow this conclusion to stand and, in fact, Bhadra made the neeeded reservations, viz. ‘the detection could be crucial, depending on the angular size of the object’ and ‘unless the source is in a lower density environment’. These are the points addressed in the following sections. We require, first, the likely whereabouts of the single source such that it can give the particle spectrum needed in the knee region.

## 3 The present treatment

### 3.1 The predicted particle spectrum

We have calculated cosmic ray energy spectra originating from SNR of different ages and distances. The model of acceleration was described in EW (2001a) and briefly in I. The propagation of accelerated cosmic rays through the ISM was calculated using two alternative assumptions about the mode of propagation from the source: ‘anomalous’ diffusion, viz. making allowance for the fractal-like nature of the ISM, and normal, Gaussian, diffusion. Following the work of Lagutin et al.(2001a,b) we distinguish these two modes by a parameter $`\alpha `$ determined by the fractal nature of the medium; $`\alpha `$ is equal to 1 or 2 for anomalous or normal diffusion, respectively.

The difference between the two modes can be seen most clearly in the shape of the lateral distribution function for the cosmic ray intensity: $`\frac{1}{(1+x^2)^2}`$ for $`\alpha =1`$ and $`exp(\frac{x^2}{4})`$ for $`\alpha =2`$, with $`x=\frac{r}{R_d}`$, $`r`$ being the distance from the radius $`R_s=100`$ pc where the particles start to diffuse and $`R_d`$ being the diffusion radius which is defined as $`R_d=H_z(\frac{t}{\tau (E)})^{\frac{1}{\alpha }}`$, i.e., there is a different time dependence for the two modes. $`H_z=1`$ kpc for the vertical scale of the galactic halo and $`t`$ and $`\tau (E)`$ are the diffusion time and time against escape, respectively, for the protons. Details of the propagation model are given in EW (2002) and in I.

The difference in the two lateral distributions is quite dramatic. Thus, for $`x=`$ 2, 4, 6 and 8, the ratio of the intensity for anomalous diffusion to that for normal diffusion changes as 0.11, 0.19, 5.9 and 2.1$``$103, respectively. The long tail for anomalous diffusion - the occasional considerable ‘penetration’ ( ’Levy flights’ ) - can have important implications. An example is that for the secondary to primary ratio; to our knowledge the implications have not been worked out.

We have determined the proton energy spectra expected for the two values of $`\alpha `$ and these are given in Figures 1 and 2. Also indicated in the Figures is the spectrum ‘required’ by the SS model. We have argued that the ‘needed’ particles are oxygen nuclei for the first ‘peak’ at 3 PeV and iron nuclei for the second ‘peak’ at $`12`$ PeV, and the requirement has been converted to rigidity before plotting, in Figures 1 and 2. An alternative association of the ‘peaks’ with helium and oxygen changes the indicated SS spectrum no more than by 17%. ( EW, 2003b ).

Dotted lines above and below the ‘SS’ curve indicate its uncertainty limits; the least uncertainty and therefore the most important constraint of the SS model is in the knee region. The limits come from the fit to the experimental data on the primary energy spectrum measured by means of the Cherenkov light emitted by extensive air showers ( EW, 2001b ). These data determine the magnitude of the uncertainty in the range of about a decade below the knee. At lower energies, - $`10^310^4`$ GeV - the upper limit is determined by the uncertainty of the direct measurements of the primary CR energy spectrum ( Biermann and Wiebel-Sooth, 1999 ). Although the shape of the SS spectrum has been adopted from the theoretical model of Berezhko et al.( 1996 ) its lower experimental limit at low energies is completely uncertain because the contribution of the Single Source to the total CR intensity at these energies is negligibly small and even consistent with zero. However, the actual value of the uncertainty at TeV energies is not important for this analysis because we have made a quantitative comparison of the calculated proton spectra with the SS spectrum only in the most important range of energies, viz. that covering a decade below the knee. The energy range of $`logE=4.65.6`$ used for this comparison is indicated in the Figures by ‘min’ and ‘max’.

Some comments are needed about the Figures. It is evident that $`\alpha =2`$ always gives a bad fit (Figure 2) whereas $`\alpha =1`$ can give a reasonable one, at least for energies above 104 GeV, up to the cut-off at 4$`10^5`$ GeV, for a range of age (T) and distance (R) values 111In spite of the fact that such an integral characteristic as $`\chi ^2/ndf`$ has a formally acceptable value for $`\alpha =2`$ at T$`85`$ kyear the slopes of the proton spectra are too steep to give the sharp knee and we do not consider it as a good contender.. Calculations made for a wide range of T and R allow us to estimate the range over which there is satisfactory agreement between calculations and the SS model; it is $`85<T<115`$ kyear and $`250<R<400`$ pc for the adopted set of input parameters.

It is relevant to point out that larger distances, too, would give a good fit to the spectral shape. The necessary upward movement in intensity could be effected by increasing the fraction of the shock energy going into CR. As was pointed out in I, Berezhko et al.,(1996) used much bigger values than the 10% used here \- their highest being 80%. Specifically, for 100 kyear an increase in CR yield by a factor 8 would allow the source to be at $`650`$ pc.

The shape and the absolute intensity of the CR energy spectrum give the most stringent constraints on the age and distance. The T-R region of SNR which could give an acceptable spectrum ( for our standard ‘10%’ ) is shown in Figure 3 by the 95% confidence level contour. Following Bhadra’s approach we have also used such an integrated characteristic of the spectrum as the energy density contained in it. Again, we used for the comparison just the last decade of the spectrum below the knee, because intensities at low energies are so poorly determined by the Single Source Model. The energy density contained in the spectrum of our single source between $`410^4`$ GeV and $`410^5`$ GeV is 1.84$`10^4\mathrm{eVcm}^3`$. Comparison of this value, allowing for uncertainty, with those expected for SNR of different ages and distances gives an acceptable T-R region indicated by the dashed line in Figure 3. It overlaps with the region deduced from the comparison of the spectral shapes. This proves the consistency of these results, although we must admit that the two methods are not completely independent, because consistent values for the spectra should inevitably give consistent energy densities contained in them.

In any case this analysis indicates that our single source should be located at about 300 - 350 pc from the Sun and should be about 90 - 100 kyears old. We know of no objection to such parameters ( see §4 for its likely location ).

### 3.2 The predicted gamma ray flux

#### 3.2.1 Calculations for the standard gas density

Fluxes of gamma rays from SNR of different ages and at different distances from the Sun have been calculated in I. In Figure 3 we show contours of the fluxes calculated for the interactions of protons, accelerated by the SNR, which propagate through the ISM by way of anomalous diffusion ( $`\alpha =1`$ ) and collide with nuclei of the ISM, The density of the ISM for these calculations has been taken to be our standard $`n=1\mathrm{c}\mathrm{m}^3`$. It is seen that the expected fluxes are about 3$`10^7\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >0.1`$ GeV and about 2$`10^{11}\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >1`$ TeV.

It is necessary now to study the implications of the fact that the cosmic rays from our Single Sourse are not only protons, but that the composition is mixed. At the same rigidity it consists of 21% P, 48% O, 13% ‘heavy’ nuclei ($`10<Z<23`$) and 18% Fe. Because nuclei of the same rigidity are more efficient in the production of gamma quanta ( see I ), the expected fluxes should be increased by a factor of 10.5. This gives expected fluxes of $`310^6\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >0.1`$ GeV and $`210^{10}\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >1`$ TeV. It is necessary to point out that the CR energy injection is higher by the same amount.

#### 3.2.2 The appropriate gas density to adopt

Since this is a first order effect, a detailed examination of the problem is required. As pointed out in I, there are several aspects that need attention, specifically:
(i) the effect of the pre-SNR winds from the progenitor star in excavating a ‘hole’ into which the SNR shock expands, this is essentially the origin of the HISM ( note that for vey young SNR, however, the progenitor star’s wind may have enhanced the density );
(ii) the general ISM density in the region where the progenitor star was situated and
(iii) the pressure, or otherwise, of clumpy gas ( often molecular ) in the vicinity of the SNR.

In the present case, it is almost certain that the progenitor star - at 300-350 pc from the sun - was in the Hot Interstellar Medium, where the density is often quoted as $`n310^3`$ cm-3 ( Berezhko et al., 1996; Cox and Reynolds, 1987, for example quote a density of $`410^3`$cm-3 ). Indeed, since all our model calculations assume this to be the case, a low density target material is a prerequisite. The assumption of the HISM comes from two factors:
(i) the absolute maximum particle energy needed to explain the knee ( $`310^6`$ GeV for oxygen nuclei ) appears only for this density, in the Berezhko et al. model, and
(ii) the HISM is eminently reasonable for a nearby source.

Factor (ii) can be examined in more detail. Insofar as the sun is located on the edge of a spiral arm, in one hemisphere ( the South ), at least, the gas density will be low. Frish ( 1997 ) has considered the situation in detail. Beyond the very local region, where there is the ‘Local Fluff’ ( of extent $``$3 pc ) in certain directions the column density of atomic ( and molecular ) hydrogen is very low. The ‘Local Bubble’ ( Loop I ) is quickly reached, this Bubble being caused by several SN over the past Myear. Here, one expects the HISM, with its density $`310^3`$ cm-3. Concerning the interarm region, Frish quotes an undisturbed part as having a density of $`410^4`$ cm-3, i.e. even smaller. However, much of the Local Bubble is here and the density will be higher because of material brought in from elsewhere.

Interestingly, Frish (1981) suggested, earlier, that the sun is embedded in one of the super-bubble shells associated with the formation of the Scorpius-Centaurus Association. Another possibility for the formation of Loop I (‘but very similar to that given above is that it was caused by activity in the ‘Upper Centaurus Lupus’ subgroup some 14-15 Myear ago ( De Geus, 1991 )).

A distinction must be made between a ‘recent’ SN exploding in Loop I, i.e. in the low density HISM caused by previous SNR and stellar activity, and the low level gamma ray flux from Loop I as a whole. This latter was considered in I and it was argued there that the mean density overall, allowing for molecular clouds inside the Loop, and the piled-up gas in the edges of the Loop, is $`n`$ 0.1 cm-3. There is no conflict with our $`n=310^3`$ cm-3 if the single source is well into the interior.

#### 3.2.3 The predicted gamma ray flux

With $`n=310^3`$, the predicted flux is reduced to $`10^8`$ cm-2s-1 for $`E_\gamma >`$0.1 GeV and $`0.710^{12}`$cm-2s-1 for $`E_\gamma >`$ 1 TeV.

If, surprisingly, the source is not in the HISM in the Local Bubble, but is isolated, then as was discussed in I, $`n0.1\mathrm{cm}^3`$ and the predicted fluxes are $`310^7\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >0.1`$ GeV and $`210^{11}\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >1`$ TeV.

The predicted fluxes are shown in Figure 4, for the various possibilities of n and Z.

### 3.3 The limiting sensitivity in practice

The sources must be detected against a background due to the Galactic diffuse emission arising from CR - ISM nucleus interactions ( and, in the case of TeV gamma rays, protons interacting in the atmosphere ). This background can be allowed for in a straightforward manner for ‘point’ sources, where a subtraction can be made of the signal nearby to the source, but for extended sources the problem is much more severe. The angular radius of a SNR, which is 90-100 kyear old, seen from a distance of 300-350 pc is about 20 ( see also Figure 5 in I ). For such a large size the determination of the background by a linear interpolation of the intensity between the adjacent regions is not at all accurate. Estimates of the background based on the known column density of the target gas, although more appropriate, are again not sufficiently accurate because the intensity of the initiating cosmic rays along the line of sight cannot be assumed to be strictly constant.

Estimates of the limiting fluxes were made from the available data and they are given in I. For a source of radius 20 the average limits are $`310^6\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >0.1`$ GeV and $`510^{11}\mathrm{cm}^2\mathrm{s}^1`$ for $`E_\gamma >1`$ TeV. These limiting sensitivities, which relate to observations made so far, are indicated in Figure 4.

The lower energy limit comes from an extrapolation of the EGRET results ( Hartman et al, 1999 ), which refer to small angular sizes, to 20. That for the upper energy band comes largely from the Tibet and HEGRA arrays ( Amenomori et al., 2001 and Lampeitl et al., 2001, respectively ).

It is instructive to examine these sensitivities in terms of the background fluxes. For $`E_\gamma >0.1`$ GeV, the intensity at low latitudes towards the Galactic Anticenter is about $`10^4\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1`$ and for our 20-radius source ( $`\mathrm{\Omega }`$ = 0.4 sr ) the flux would be $`410^5\mathrm{cm}^2\mathrm{s}^1`$. The plotted limit of $`310^6\mathrm{cm}^2\mathrm{s}^1`$ is thus 8% of the background. For the direction to the Galactic Center, where the background is higher by a factor of 3, this percentage is reduced to (2-3)%. When allowance is made for other sources of uncertainty our adopted limit seems reasonable.

For $`E_\gamma >1`$ TeV, the corresponding background in the Galactic Anticenter direction has been estimated to be $`1.610^9\mathrm{cm}^2\mathrm{s}^1`$ ( Porter and Protheroe, 1999 ), viz a flux of $`(67)10^{10}\mathrm{cm}^2\mathrm{s}^1`$ over a 20-radius source. Our adopted limit of $`510^{11}\mathrm{cm}^2\mathrm{s}^1`$ is thus the same: $``$8% of the background. In the direction of the Galactic Center the background is expected to be twice as large and the fraction of the limiting flux falls to $``$4%. This figure, too, is reasonable in view of misidentified proton contributions and a variety of technical problems associated with making absolute measurements at different zenith angles.

## 4 Conclusions

The conclusions are indicated in Figure 4. Only if the single source is in a ‘high’ density region ($`n0.1\mathrm{cm}^3`$) and the primaries are ‘heavy’ nuclei ( Z = 10 ) ( with consequent high CR energy injection ) will it be possible to detect the Single Source at high energies. Such a situation is not impossible but there would be problems for our Single Source model with such a high ambient density, specifically that it would not be possible with the model in use to reach the required 3 PeV energy at the knee. Insofar as we consider it very likely that the source is in the HISM and is extended, the chance of detecting it with contemporary instruments is considered, by us, to be very low. There are, however, hopes for the future. It is germane to consider in which direction future, improved, gamma ray detectors should be pointed in order to see it. The best that can be done here is to suggest the general direction of Loop I. However, since this structure occupies about 25% of the sky the ‘advice’ is not very helpful.

## Acknowledgments

The authors are grateful to The Royal Society and The University of Durham for financial support. Dr. Paula Chadwick, of the HESS collaboration, is thanked for providing useful information. Two unknown referees are also thanked for useful comments and suggestions.

## References

Aharonian, F.A. et al., 1997, Astroparticle Physics, 6, 369
Amenomori, M. et al., 2001, 27th Int. Cosm. Ray Conf. ( Hamburg ), 6, 2344
Berezhko, E.G., Elshin, V.K. and Ksenofontov, L.T., 1996, J. Exp. Theor. Phys., 82, 1
Bhadra, A., 2002, J. Phys. G: Nucl. Part. Phys., 28, 397
Biermann P.L. and Wiebel-Sooth B., 1999, Astronomy & Astrophysics - Interstellar Matter, Galaxy, Universe, Landolt-Bórnstein, Springer Verlag, Berlin/Heidelberg, 3, 37
Cox D.P. and Reynolds R.J., 1987, Annual Rev. Astron. Astrophys. 25, 303
De Geus, E.J., 1991, Astron. Astrophys., 262, 258
Drury, L.O’C., Aharonian, F.A. and Völk, 1994, Astron. Astrophys., 287, 959
Erlykin, A.D. and Wolfendale, A.W., 1997, J. Phys. G: Nucl. Part. Phys., 23, 979
Erlykin, A.D. and Wolfendale, A.W., 2001a, J. Phys. G: Nucl. Part. Phys., 27, 941
Erlykin, A.D. and Wolfendale, A.W., 2001b, J. Phys. G: Nucl. Part. Phys., 27, 1005
Erlykin, A.D. and Wolfendale, A.W., 2002, J. Phys. G: Nucl. Part. Phys., 28, 2329
Erlykin, A.D. and Wolfendale, A.W., 2003a ( submitted to J. Phys. G )
Erlykin, A.D. and Wolfendale, A.W., 2003b, Nucl. Phys. B ( Proc. Suppl. ) ( in press )
Frisch, P.C., 1981, Nature, 293, 377
Frisch, P.C., 1997, astro-ph/9705231
Hartman, R.C. et al., 1999, Astron. Astrophys. Suppl. Ser., 123, 179
Kulikov, G.V. and Khristiansen, G.B., 1958, J. Exp. Theor. Phys., 35, 635
Lagutin, A.A. et al. 2001a, Nucl. Phys. B ( Proc. Suppl. ), 97, 267
Lagutin, A.A. et al. 2001b, Proc. 27th Int. Cosm. Ray Conf. ( Hamburg ), 5, 1900
Lampeitl, H. et al., 2001, 27th Int. Cosm. Ray Conf. ( Hamburg ), 6, 2348
Porter T.A. and Protheroe R.J., 1999, 26th Int. Cosm. Ray Conf., Salt Lake City, 4. 306


text:
# Discrete Fourier Transform in Nanostructures using Scattering

## I Introduction

The fast developing fields of quantum computing and quantum information processing have attracted much interest recently . Most of the papers in the field of quantum computing, such as Shor’s and Grover’s about factoring numbers and searching a database, use many-particle physics (i.e. qubits) for their implementation. However, it is also possible to implement quantum information processes in the so-called unary representation of dimension $`N`$, which consists of a statistical ensemble of non-interacting $`N`$-level systems that can be described by a single-particle Hamiltonian. So only interference, but no entanglement is needed for the implementation of quantum algorithms in the unary representation. It was already shown that the Grover algorithm can be performed in a single shot using the unary representation. While in Ref.  the quantum database search has been done experimentally with atoms in a beam that serve as non-interacting $`N`$-level systems, Refs.  and  show theoretically that the large spin of molecular magnets and nuclear spins in semiconductors can be used as non-interacting $`N`$-level systems, respectively. Here we propose a method to perform the classical discrete Fourier transform (DFT) in the time ensemble of quantum or classical scattering events, where the scattering potential provides the classical $`N`$-level system for the read-in of information and the $`N`$ detectors provide the classical $`N`$-level system for the read-out of the information. In the case of quantum scattering a single-particle Hamiltonian describes the information processing. In contrast to Shor’s algorithm our approach requires an exponential overhead of the hardware. However, once our device is built the DFT can be computed very efficiently, i.e. in a single shot. In our method we propose to scatter a coherent matter or laser beam off a two-dimensional potential consisting of $`N`$ concentric rings or $`N`$ peaks, the amplitudes of which represent the elements of a vector in an $`N`$-dimensional vector space where the information is encoded. The Fourier transformation of this vector can then be read out immediately by detectors surrounding the 2D potential. Our method is closely related to the general field of Fourier optics, where a diffraction pattern is Fourier transformed into a desired image pattern by means of laser scattering. We go a step further in the sense that we describe how to use a 2D scattering potential to perform the DFT efficiently. The concentric rings and peaks can be produced by charge or spin distributions. The calculation time $`\tau _c`$ is roughly given by the diameter of the detector ring (e.g. 1 mm) divided by the light velocity $`c=3\times 10^8`$ m/s, i.e. $`\tau _c10^{11}`$ s. Using for example a laser with a wavelength of 500 nm and a spot radius of 100 $`\mu `$m enables the Fourier transformation of a function with 100 sampling points. The larger the spot size and the smaller the wavelength of the laser, the more sampling points can be achieved. The advantage of this device would be that it can be used already for example as DFT coprocessor to speed up present computers operating at room temperature. Rough estimates show a possible speedup of $`10^6`$ or more (see below).

In Sec. II we give an overview on the DFT and the Fast Fourier Transform (FFT). Sec. III and Sec. IV describe the quantum-mechanical coherent matter scattering and the classical laser beam scattering off a specially designed 2D potential that is used for the implementation of the DFT.

## II Discrete Fourier Transformation

The discrete Fourier transformation (DFT) has a wide range of applications in the field of digital signal processing, such as spectral analysis and filtering. The DFT enables the possibility to extract the period of a function, which can be used to factor numbers. Thus finding efficient ways to compute the DFT is the key to the efficient factorization of numbers. The calculation time of a DFT can be greatly reduced by the classical algorithm called Fast Fourier Transform (FFT), which gets rid of the redundancies found in the DFT. In a similar way Shor’s quantum algorithm reduces the calculation time of the discrete quantum Fourier transformation (DQFT) by means of an efficient quantum circuit that uses qubits as information carriers. In our proposal the superposition principle of quantum or classical wavefunctions is used to perform an $`N`$-dimensional DFT with $`N`$ shots or even a single shot in one scattering device. The FFT provides then the possiblity to link several scattering devices in parallel. This means that the speedup factor $`N`$ can be multiplied by the number of the parallel scattering devices.

First it is instructive to give the definition of a DFT. When we calculate the DFT of a function $`f(z)`$ that maps complex numbers on complex numbers, we can first discretize the function $`f(z)`$ in $`M`$ total sampling points $`z_1,z_2,\mathrm{},z_M`$, leading to $`M`$ values $`f_1=f(z_1),f_2=f(z_2),\mathrm{},f_M=f(z_M)`$. Then the DFT of $`f`$ is given by

$$F_k=\underset{j=0}{\overset{M1}{}}e^{2\pi ijk/M}f_j.$$
(1)
The matrix elements of a DFT consist of the irreducible representation of the cyclic group produced by a rotation R2π/M about the angle $`2\pi /M`$, where $`M`$ is the dimension of the vector space, which is equal to the number of sampling points. For example a DFT of dimension $`M=4`$ is represented by

$$\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& e^{\frac{1}{4}2\pi i}& e^{\frac{1}{2}2\pi i}& e^{\frac{3}{4}2\pi i}\\ 1& e^{\frac{1}{2}2\pi i}& e^{2\pi i}& e^{\frac{1}{2}2\pi i}\\ 1& e^{\frac{3}{4}2\pi i}& e^{\frac{1}{2}2\pi i}& e^{\frac{1}{4}2\pi i}\end{array}\right]$$
(2)
It will turn out that the exponential overhead of the hardware imposes a physical limit to the number of sampling points $`N`$ that can be processed by our DFT device. However, the Danielson-Lanczos Lemma, which is widely used for the FFT algorithm, states that the DFT for $`M`$ sampling points can be separated into the sum over the DFT for the $`M/2`$ even sampling points and the DFT for the $`M/2`$ odd sampling points, i.e.

$$F_k=\underset{j=0}{\overset{M/21}{}}e^{2\pi ik(2j)/M}f_{2j}+\underset{j=0}{\overset{M/21}{}}e^{2\pi ik(2j+1)/M}f_{2j+1}.$$
(3)
Usually, the FFT uses this Danielson-Lanczos Lemma recursively down to Fourier transforms of length 20. Therefore it is essential that $`M=2^m`$, where $`m`$ is a positive integer. As we want to use 2D scattering potentials with real amplitudes, we can separate the real and imaginary part of the DFT:

$$F_k=\underset{j=0}{\overset{M1}{}}e^{2\pi ijk/M}g_j+i\underset{j=0}{\overset{M1}{}}e^{2\pi ijk/M}h_j,$$
(4)
where $`f(z_j)=g(z_j)+ih(z_j)`$ with real functions $`g(z_j)`$ and $`h(z_j)`$.

Since our device can perform the DFT efficiently for $`N`$ sampling points, we can stop the recursion at $`N=2^n`$ sampling points, where $`n`$ is a positive integer that is smaller than $`m`$. Using our device as coprocessor of an existing computer, the remaining work of the main processor consists in combining the $`M/N`$ partial Fourier transforms, which takes of order $`\mathrm{log}_2(M/N)M/N=(mn)M/N`$ operations. Thus, compared with initially necessary $`M\mathrm{log}_2(M)`$ operations a speedup of at least a factor $`N`$ is achieved. In this estimate of the speedup, the necessary time to calculate the DFT for the $`N`$ sampling points is assumed to be of the order of 10 ps, i.e. negligible. Alternatively, one could say that an existing computer equipped with our device is capable of factoring numbers that are $`N`$ times larger than the largest nowadays factorizable number. If $`M/N`$ DFT devices are linked in parallel, a speedup of factor $`M`$ is achieved.

## III Coherent matter scattering

Now we describe the various possible implementations of our DFT device, in which a coherent matter or laser beam is scattered off a 2D potential. This potential consists of radial $`\delta `$-functions or $`\delta `$-peaks on a straight line, the amplitudes of which both sample the points of the real or imaginary part of our function $`f(z)`$ and represent the amplitudes of our basis vectors. Let $`𝐯=_{j=1}^Na_j𝐰_j`$ be the input vector, where $`𝐰_j`$ are the basis vectors with amplitudes $`a_j=g(z_j)`$ or $`a_j=h(z_j)`$. Then our 2D potential is of the form $`V(𝐱)_{j=1}^Na_j\delta (𝐱𝐱_j)`$. The read-out is done by detectors surrounding this potential. While the coherent matter beam needs a high precision of the potential, the laser beam provides a high fault-tolerance of the potential since faults that are smaller than the wavelength of the laser beam remain undetected. We give the derivation for both the coherent matter and the laser beam scattering. We first start with the coherent matter beam. We use the nomenclature of Ref. . First we calculate the 2D scattering amplitude, from which we can infer the precise shape of the 2D potential required for the implementation of the DFT. For this we have to start from the well-known Lippman-Schwinger equation

$$𝐱|\psi ^{(\pm )}=𝐱|\varphi +\frac{2m}{\mathrm{}^2}d^2x^{}G_\pm (𝐱,𝐱^{})𝐱^{}\left|V\right|\psi ^{(\pm )},$$
(5)
where $`|\varphi `$ is the incoming wave, $`V`$ is the 2D scattering potential, and

$$G_\pm (𝐱,𝐱^{})=\frac{\mathrm{}^2}{2m}𝐱\left|\frac{1}{EH_0\pm iϵ}\right|𝐱^{}$$
(6)
is the Green’s function, which also solves the Helmholtz equation $`(\mathbf{}^2+k^2)G_\pm (𝐱,𝐱^{})=\delta ^{(2)}(𝐱𝐱^{})`$. For the evaluation of the Green’s function we transform first to the momentum representation, i.e.

$$G_\pm (𝐱,𝐱^{})=\frac{\mathrm{}^2}{2m}\frac{d^2p^{}}{(2\pi \mathrm{})^2}\frac{e^{\frac{i}{\mathrm{}}𝐩^{}(𝐱𝐱^{})}}{E\frac{p^2}{2m}\pm iϵ},$$
(7)
where we have used the 2D overlap $`𝐱|𝐩=e^{i𝐩𝐱/\mathrm{}}/2\pi \mathrm{}`$. Inserting $`E=\mathrm{}^2k^2/2m`$, $`𝐩^{}=\mathrm{}𝐐`$, and going to polar coordinates $`(Q,\varphi )`$, we get

$`G_\pm (𝐱,𝐱^{})`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle _0^{\mathrm{}}}Q𝑑Q{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle \frac{e^{iQ|𝐱𝐱^{}|\mathrm{cos}\varphi }}{k^2Q^2\pm iϵ}}`$ (8)
$``$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle _0^{\mathrm{}}}Q𝑑Q{\displaystyle \frac{2^{3/2}\sqrt{\pi }}{\sqrt{Q|𝐱𝐱^{}|}}}`$ (10)
$`\times {\displaystyle \frac{\mathrm{cos}(Q|𝐱𝐱^{}|\frac{\pi }{4})}{k^2Q^2\pm iϵ}}`$
$`=`$ $`{\displaystyle \frac{1}{4\sqrt{\pi }}}{\displaystyle \frac{e^{\pm i(k\left|𝐱𝐱^{}\right|\frac{\pi }{4})}}{\sqrt{2k\left|𝐱𝐱^{}\right|}}}.`$ (11)
For the approximation we have used the integral representation of the Bessel function $`J_0(Q\left|𝐱𝐱^{}\right|)=\frac{1}{2\pi }_0^{2\pi }𝑑\varphi e^{iQ\left|𝐱𝐱^{}\right|\mathrm{cos}\varphi }\sqrt{\frac{2}{\pi Q\left|𝐱𝐱^{}\right|}}\mathrm{cos}(Q\left|𝐱𝐱^{}\right|\frac{\pi }{4})`$, which holds in the asymptotic regime $`Q\left|𝐱𝐱^{}\right|>>1`$, i.e. the local 2D potential $`𝐱^{}\left|V\right|𝐱^{\prime \prime }=V(𝐱^{})\delta ^{(2)}(𝐱^{}𝐱^{\prime \prime })`$ and the detectors at $`𝐱`$ are far apart from each other. Then we obtain with Eq. (5)

$$𝐱|\psi ^{(+)}=\frac{1}{2\pi }\left[e^{i𝐤𝐱}+\frac{e^{i(kr\frac{\pi }{4})}}{\sqrt{kr}}f(𝐤^{},𝐤)\right],$$
(12)
where $`r:=\left|𝐱\right|`$, $`r^{}:=\left|𝐱^{}\right|`$, $`𝐤^{}:=k\frac{𝐱}{\left|𝐱\right|}`$, and

$$f(𝐤^{},𝐤)=\frac{\sqrt{2\pi }}{4}\frac{2m}{\mathrm{}^2}d^2x^{}e^{i𝐤^{}𝐱^{}}V(𝐱^{})𝐱^{}|\psi ^{(+)}$$
(13)
is the 2D scattering amplitude. The amplitude of the 2D spherical wave decays as $`1/\sqrt{r}`$. This is in agreement with the topological argument that the circumference of a circle is $`2\pi r`$ and that the probability has to be conserved on the circle. In first-order Born approximation, where $`|f(𝐤^{},𝐤)|\sqrt{kr}`$, we can insert $`𝐱^{}|\psi ^{(+)}=\frac{e^{i𝐤𝐱^{}}}{2\pi }`$ into Eq. (13), and we define $`𝐪:=𝐤^{}𝐤`$ as the scattering wavevector. Eq. (13) becomes then

$$f^{(1)}(𝐤^{},𝐤)=\frac{1}{4\sqrt{2\pi }}\frac{2m}{\mathrm{}^2}d^2x^{}e^{iqx^{}\mathrm{cos}\phi }V(𝐱^{}),$$
(14)
where $`\phi `$ is the angle between $`𝐪`$ and $`𝐱^{}`$.

For the implementation of the DFT we start first with a 2D potential of the form $`V(𝐱^{})=V(r^{})(1+\mathrm{cos}\phi )`$. If the scattering is elastic, we choose $`\phi =\pi /2\theta /2`$. If the material allows inelastic scattering, we choose $`\phi =\pi /2\theta `$. It will turn out that for this potential $`𝐪`$ has to point always into the same direction. For example $`𝐪`$ can be always perpendicular to the x-axis. Then we can define $`\alpha `$ as the angle between the x-axis and $`𝐱^{}`$, which turns our potential into $`V(𝐱^{})=V(r^{})(1+\mathrm{sin}\alpha )`$. This 2D potential is shown in Fig. 1. Inserting this potential into Eq. (14) leads to

$$f^{(1)}(𝐤^{},𝐤)=\frac{\sqrt{2\pi }}{4}\frac{2m}{\mathrm{}^2}r^{}𝑑r^{}\left[J_0(qr^{})+iJ_1(qr^{})\right]V(r^{}),$$
(15)
which in the asymptotic limit $`r^{}1/q`$ becomes

$$f^{(1)}(𝐤^{},𝐤)=\frac{1}{2\sqrt{q}}\frac{2m}{\mathrm{}^2}\sqrt{r^{}}𝑑r^{}e^{i(qr^{}\frac{\pi }{4})}V(r^{}).$$
(16)
By choosing $`V(r^{})=_{j=1}^Na_j\delta (r^{}r_j)/\sqrt{r^{}}`$ we can encode the information into the $`N`$ amplitudes $`a_j`$ and perform a DFT

$$f^{(1)}(𝐤^{},𝐤)=\frac{1}{2\sqrt{q_\nu }}\frac{2m}{\mathrm{}^2}e^{i\frac{\pi }{4}}\underset{j=1}{\overset{N}{}}a_je^{iq_\nu r_j}$$
(17)
In order for it to be a discrete Fourier transformation (DFT), we have to discretize the variables. Let us make a first attempt: $`r_j=j\frac{l}{N},j=1,2,\mathrm{},N`$, and $`q_\nu =\nu \frac{2\pi }{l},\nu =1,2,\mathrm{},N`$, with $`\nu `$ referring to the $`\nu `$th detector. The smallest optically resolvable distance between the radii $`r_j`$ is $`\lambda /2`$. So it is safe to choose $`l=N\lambda `$, where $`\lambda =2\pi /k`$ is the deBroglie wavelength of the coherent beam. Then we obtain $`\frac{\nu }{N}=2\mathrm{sin}\frac{\theta }{2}`$ for the elastic scattering and $`\frac{\nu }{N}=\mathrm{tan}\theta `$ for the inelastic scattering (see below). The radii of the rings are multiples of the wavelength $`\lambda `$, i.e. $`r_j=j\lambda `$. There is still the requirement $`r_j1/q_\nu `$ that has to be satisfied. Therefore we have to choose $`j=N+1,N+2,\mathrm{},2N`$, which still produces a DFT, since the DFT is $`2\pi `$-periodic.

Our 2D potential $`V(𝐱^{})`$ has the $`1+\mathrm{cos}\phi `$ dependence, so that we obtain $`J_0`$ and $`J_1`$ Besselfunctions in Eq. (15). Since we do not want to adjust the orientation of $`V(𝐱^{})`$ according to the scattering vector $`𝐪`$, we have to keep the direction of $`𝐪`$ constant in space. This requirement can be satisfied in two ways, which can be identified by elastic or inelastic scattering. In the case of elastic scattering $`𝐤`$ and $`𝐤^{}`$ have the same length. So $`q:=\left|𝐤𝐤^{}\right|=2k\mathrm{sin}\frac{\theta }{2}`$, where $`\theta `$ is the angle between $`𝐤`$ and $`𝐤^{}`$. For every incoming beam of angle $`\vartheta _{\mathrm{in}}=\theta /2`$ the resulting vector element after the DFT can be read out in the detector at the outgoing angle $`\vartheta _{\mathrm{out}}=\vartheta _{\mathrm{in}}=\theta /2`$. The number of incoming and outgoing angles is exactly $`N`$. Thus the DFT can be performed with $`N`$ shots, which is shown in Fig. 2. In the case of inelastic scattering, $`𝐤`$ and $`𝐤^{}`$ do not have the same length anymore. It is then possible to perform the DFT in a single shot. For this $`q:=\left|𝐤𝐤^{}\right|=k\mathrm{tan}\theta `$. Since $`k^{}=n_rk`$ and $`k=k^{}\mathrm{cos}\theta `$, we obtain an angle-dependent refractive index $`n_r=1/\mathrm{cos}\theta `$, which is probably difficult to realize experimentally. At least in the case of a laser beam (see below) it may be possible to use photonic crystals that mimick $`n_r=1/\mathrm{cos}\theta `$. The inelastic scattering is shown in Fig. 3.

Instead of $`V(𝐱^{})=(1+\mathrm{sin}\alpha )_{j=1}^Na_j\delta (r^{}r_j)/\sqrt{r^{}}`$ consisting of concentric rings, we can alternatively use peaks of the form $`V(𝐱^{})=_{j=1}^Na_j\delta (r^{}r_j)\delta (\alpha \alpha _j)/r^{}`$, which is shown in Fig. 4. After putting the peaks on a straight line ($`\alpha _1=\alpha _2=\mathrm{}=\alpha _N=\pi /2`$), we obtain from Eq. (14)

$$f^{(1)}(𝐤^{},𝐤)=\frac{1}{4\sqrt{2\pi }}\frac{2m}{\mathrm{}^2}\underset{j=1}{\overset{N}{}}a_je^{iq_\nu r_j\mathrm{cos}\phi _\nu },$$
(18)
where $`\alpha _\nu =\phi _\nu +\pi /2`$ is the angle between the scattering vector $`𝐪_\nu `$ and the x-axis. With this potential we can discretize the variables similarly as before: $`r_j=j\frac{l}{N},j=1,2,\mathrm{},N`$, and $`q_\nu \mathrm{cos}\phi _\nu =\nu \frac{2\pi }{l},\nu =1,2,\mathrm{},N`$. The advantage of this method is that the peaks can be much easier produced experimentally than the rings. Also, we can start from the peak with $`j=1`$, because we do not require $`r_j1/q_\nu `$ anymore. There is also no scaling factor $`1/\sqrt{q_\nu }`$ needed for the detectors (see Eq. (17)) and the scattering can be done elastically with a single shot (see Fig. 5) because the scattering vectors $`𝐪_\nu `$ do not need to point into the same direction. However, the potential made of peaks decreases as $`1/r^{}`$, and not as $`1/\sqrt{r^{}}`$ as for the rings.

For the experimental realization of the 2D scattering off a 2D potential one could use a 2D electron gas that is scattered off an electrostatic potential. Like in the case of the laser beam scattering (see next section), it would be also possible to use 3D scattering off a 2D potential, such as the scattering of spin-polarized neutrons off a spatial electronic spin distribution.

## IV Laser beam scattering

Since coherent matter waves have a very short deBroglie wavelength, the resolution of the 2D potential must be very high ($`<1`$ Å), which is beyond today’s technology. In order to overcome this difficulty, one can use 3D scattering of a laser beam off a 2D potential, where the wavelength of the laser can be chosen according to the accuracy of the 2D potential. For this we have to start from Maxwell’s equations. Following Ref.  we need to solve the Helmholtz equation

$`(\mathbf{}^2+k^2)𝐃`$ $`=`$ $`\mathbf{}\times \mathbf{}\times (𝐃ϵ_0𝐄)`$ (20)
$`{\displaystyle \frac{iϵ_0\omega }{c}}\mathbf{}\times (𝐁\mu _0𝐇),`$
where $`𝐃`$ and $`𝐁`$ are the electrostatic and magnetic fields within the sample, $`𝐄`$ and $`𝐇`$ are the respective fields of the laser beam, $`ϵ_0`$ is the dielectric constant, $`\mu _0`$ is the permeability, and $`k=\sqrt{ϵ_0\mu _0}\omega /c`$ is the wavevector of our laser beam. The result reads

$$𝐃=𝐃^{(0)}+𝐀_se^{ikr}/r,$$
(21)
where the scattering amplitude in Born approximation is given by

$`{\displaystyle \frac{\mathit{ϵ}^{}𝐀_s^{(1)}}{D_0}}`$ $`=`$ $`{\displaystyle \frac{k^2}{4\pi }}{\displaystyle }d^3x^{}e^{i𝐪𝐱^{}}\{\mathit{ϵ}^{}\mathit{ϵ}_0{\displaystyle \frac{\delta ϵ(𝐱^{})}{ϵ_0}}`$ (23)
$`+(𝐧\times \mathit{ϵ}^{})(𝐧_0\times \mathit{ϵ}_0){\displaystyle \frac{\delta \mu (𝐱^{})}{\mu _0}}\},`$
where $`\delta ϵ(𝐱^{})ϵ_0`$ and $`\delta \mu (𝐱^{})\mu _0`$. Here $`𝐧_0`$ and $`\mathit{ϵ}_0`$ are the direction and the polarization of the incoming beam, $`𝐧`$ and $`\mathit{ϵ}`$ are the direction and the polarization of the outgoing beam, respectively. Let us restrict ourselves to variations in the dielectric constant $`\delta ϵ(x^{})`$, i.e. $`\delta \mu (x^{})=0`$. In analogy to the 2D potential for the coherent matter beam, we can use a dielectric variation of the ring form $`\delta ϵ(𝐱^{})=(1+\mathrm{sin}\alpha )_{j=1}^Na_j\delta (r^{}r_j)/\sqrt{r^{}}`$ or of the peak form $`\delta ϵ(𝐱^{})=_{j=1}^Na_j\delta (r^{}r_j)\delta (\alpha \alpha _j)/r^{}`$. Performing the integration in cylindrical coordinates, $`d^3x^{}=r^{}𝑑r^{}𝑑\phi 𝑑z`$, for the peak form leads to

$$\frac{\mathit{ϵ}^{}𝐀_s^{(1)}}{D_0}=\frac{k^2d}{4\pi ϵ_0}\mathit{ϵ}^{}\mathit{ϵ}_0\underset{j=1}{\overset{N}{}}a_je^{iq_\nu r_j\mathrm{cos}\phi _\nu },$$
(24)
where $`d`$ is the thickness of the sample. So 3D scattering off a 2D potential provides another means to perform the DFT.

For the experimental implementation the spatial dielectric variation $`\delta ϵ(𝐱^{})`$ can be produced for example electrostatically by means of the Kerr effect. As another example one could use a crystal with spin-orbit interaction where an imbalance of spin-up and spin-down electrons leads to different dielectric constants for left and right circularly polarized light. In this case a circularly polarized laser beam is scattered off a spatial dielectric variation $`\delta ϵ(𝐱^{})`$ that is created by a spatial electronic spin distribution. It is also worth to mention that the laser scattering can be done at room temperature.

Now we estimate the amount of sampling points (=$`N`$) that can be used with today’s technology. The state-of-the-art lasers of e.g. $`\lambda =500`$ nm wavelength produce an intensity of J=1 GW/cm2 at the sample. This implies a photon number density of $`n=\frac{J}{c\mathrm{}\omega }10^{10}`$ photons per (100 $`\mu `$m)3. If only 10-2 of the laser beam is scattered, about 104 photons remain at the detectors surrounding the sample at $`r=1`$ mm. Within a radius of 100 $`\mu `$m$`=2N\lambda `$ we can insert $`N=100`$ rings or peaks, i.e. use a basis of $`N=100`$ sampling points.

Until now we have taken a 2D variation of the dielectric constant $`\delta ϵ(𝐱^{})=_{j=1}^Na_j\delta (r^{}r_j)\delta (\alpha \alpha _j)/r^{}`$. For experiments it is more realistic to assume a 2D variation of the dielectric constant $`\delta ϵ(𝐱^{})=_{j=1}^Na_j\delta ^{(w)}(r^{}r_j)\delta ^{(w)}(\alpha \alpha _j)/wr^{}`$ made of delta functions of width $`w`$, such as $`\delta ^{(w)}(r^{}r_j)=w/\pi [(r^{}r_j)^2+w^2]`$. Then Eq. (24) becomes approximately

$`{\displaystyle \frac{\mathit{ϵ}^{}𝐀_s^{(1)}}{D_0}}`$ $`=`$ $`{\displaystyle \frac{k^2d}{4\pi ϵ_0}}\mathit{ϵ}^{}\mathit{ϵ}_0{\displaystyle \underset{j=1}{\overset{N}{}}}a_j{\displaystyle 𝑑r^{\prime \prime }\frac{we^{iq_\nu (r^{\prime \prime }+r_j)\mathrm{cos}\phi _\nu }}{\pi (r^{\prime \prime 2}+w^2)}}`$ (25)
$`=`$ $`{\displaystyle \frac{k^2de^{q_\nu |\mathrm{cos}\phi _\nu |w}}{4\pi ϵ_0}}\mathit{ϵ}^{}\mathit{ϵ}_0{\displaystyle \underset{j=1}{\overset{N}{}}}a_je^{iq_\nu r_j\mathrm{cos}\phi _\nu },`$ (26)
where we have made the substitution $`r^{\prime \prime }=r^{}r_j`$. As long as $`w1/2k=\lambda /4\pi 1/q_\nu |\mathrm{cos}\phi _\nu |`$, this exponential damping factor is negligible. Using a laser of wavelength $`\lambda =500`$ nm implies that $`w50`$ nm. It is also possible to scale up the proposed setup to longer wavelengths without changing the results. One could say that the required accuracy of the 2D potential is given by the wavelength of the laser. Therefore our scheme is very flexible and very fault-tolerant.
Since $`\frac{\nu }{N}=2\mathrm{sin}\frac{\theta }{2}`$ for elastic scattering, the detectors are to be arranged within a scattering angle of $`\pi /3`$. The smallest angle between two detectors is $`\theta _{\mathrm{min}}=\pi /3N=\pi /300`$, which leads to a minimum distance of $`r\theta _{\mathrm{min}}=10`$ $`\mu `$m between the detectors.

## V Conclusion

In conclusion, we have shown that by means of scattering a coherent matter or laser beam off a 2D potential it is possible to perform the DFT on $`N`$ sampling points. So the scattering is done quantum-mechanically or classically. The 2D potential can be made of concentric rings or peaks on a straight line. The scattering off the peaks should be much easier to implement experimentally since it can be done elastically with a single shot. Note that the DFT is performed classically in both cases since we use a classical 2D potential to encode classical information. This provides the great advantage that our device can already be used as DFT coprocessor for existing classical computers, without having to worry about the decoherence of quantum states, like in the case of quantum computing or spintronics devices. In addition, $`M/N`$ devices can be linked together to perform the DFT in parallel on $`M`$ sampling points, thereby making use of the FFT. It should be possible with today’s technology to build a DFT device with the capability to Fourier transform an $`N=1024`$-dimensional vector. If an existing classical computer was equipped with $`M/N=1024`$ such DFT devices that are linked together, a speedup of a factor $`M=1024N10^6`$ could be achieved.

## VI Acknowledgments

This work has been supported in part by the Swiss NSF, NCCR Nanoscience, U.S. NSF, and DARPA/ARO. We would like to thank Ryan Epstein and Martino Poggio for useful discussions.


text:
# I Introduction

## I Introduction

There are abundant data from solar, atmospheric, laboratory and long baseline neutrino experiments on neutrino mass and mixing. The present experimental data, including recent results from KamLAND and K2K, on neutrino masses and mixing angles can be summarized as follow. The 99.3% C.L. level allowed ranges for the mass-squared differences are constrained to be: $`1.5\times 10^3`$ eV2 $`\mathrm{\Delta }m_{atm}^25.0\times 10^3`$ eV2, and $`2.2\times 10^5`$ eV2 $`\mathrm{\Delta }m_{solar}^22.0\times 10^4`$ eV2, with the best fit values given by $`\mathrm{\Delta }m_{atm}^2=3.0\times 10^3`$ eV2, and $`\mathrm{\Delta }m_{solar}^2=7.0\times 10^5`$ eV2. The mixing angles are in the ranges of $`\mathrm{sin}^22\theta _{atm}>0.85`$ and $`0.18\mathrm{sin}^2\theta _{solar}0.37`$. Also the CHOOZ experiment gives a bound of $`0.22`$ on the $`\nu _e\nu _x`$ (where $`\nu _x`$ can be either $`\nu _\mu `$ or $`\nu _\tau `$ or a linear combination) oscillation parameter.

Present data can be explained by oscillations between three active neutrinos ***There are additional evidences for oscillation between electron and muon neutrinos from LSND experiment. If confirmed more neutrinos are needed to explain all the data. with the atmospheric neutrino and K2K data explained by oscillation between the muon and the tauon neutrinos, and the solar neutrino and KamLAND data explained by oscillation between the electron and muon neutrinos. Neutrino oscillations provide direct evidence of non-zero neutrino masses and mixing between different species of neutrinos.

We will assume that the neutrinos are Majorana as favored by some theoretical considerations . The relevant term in the Lagrangian describing Majorana neutrino masses is $``$ $`=\nu _L^TCM\nu _L+H.C.`$, where $`\nu _L=(\nu _{eL},\nu _{\mu L},\nu _{\tau L})`$ are the left-handed neutrinos and $`C`$ is the charge conjugation operator. The mass matrix $`M`$ is symmetric due to Fermi statistics. A convenient basis to study neutrino masses and mixing is the weak basis where all charged leptons are already diagonalized. The unitary mixing matrix $`V`$ and the mass eigenvalues in this basis are related by

$`D=V^TMV,`$ (1)
where $`D`$ is a diagonal matrix. The diagonal entries $`m_i`$ of $`D`$ are the mass eigenvalues which can always be made real by an appropriated choice of phase convention. For three generations of neutrinos, the mixing matrix can be described by three rotation angles and three phases. Note that in general $`TrMTrD`$.

Although there is a lot of data on neutrinos, more data are needed to determine detailed properties of neutrinos. Oscillation experiments, as indicated above, cannot measure the values $`m_i`$, but only the mass-squared differences $`\mathrm{\Delta }m_{21}^2=m_{\nu _2}^2m_{\nu _1}^2`$ and $`\mathrm{\Delta }m_{32}^2=m_{\nu _3}^2m_{\nu _1}^2`$ and some mixing angles. In analyzing the present data, it is generally assumed that $`\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{solar}^2`$, and $`\mathrm{\Delta }m_{32}^2\mathrm{\Delta }m_{atm}^2`$ and thus $`|\mathrm{\Delta }m_{21}^2||\mathrm{\Delta }m_{32}^2|`$. The sign of $`\mathrm{\Delta }m_{21}^2`$ has now been determined to be positive from matter effect in the solar neutrino data, but the sign of $`\mathrm{\Delta }m_{32}^2`$ is not determined yet.

There is at present certainly no information on any of the three CP violating phases, and in the foreseeable future no set of experiments can fully determine all the parameters in the neutrino mass matrix. Certain theoretical inputs have to be employed to reconstruct the neutrino mass matrix. Several proposals have been made to reduce the parameters, such as texture zero and determinant zero requirement for the mass matrix.

Here we propose another way of reducing the number of unknown parameters, by imposing a condition on the mass eigenvalues,

$`m_{\nu _1}+m_{\nu _2}+m_{\nu _3}=0.`$ (2)
This “zero sum” condition that the neutrino masses sum to zero implies that some of the eigen-masses must have opposite signs. These signs of course can always be arranged by making chiral phase redefinition of the neutrino fields.

If there is no CP violation in the neutrino mass matrix $`M`$, the mass matrix can always be made real and it can be diagonalized by an orthogonal transformation, namely $`V^TV=I`$. In this case the “zero sum” condition (2) is equivalent to the traceless condition, $`TrM=Tr(V^{}V^{}D)=TrD=0`$. If CP is not conserved, the “zero sum” and traceless conditions are different. One needs to be careful about the phase definitions. In this paper we simply explore the phenomenological consequences of requiring the neutrino masses to satisfy the “zero sum” condition with CP conservation without speculating on its theoretical origin. We note, however, that it holds if $`M=[A,B]`$, that is, the mass matrix can be expressed as a commutatorIn the simplest versions of models proposed in ref., $`M`$ results from radiative correction and comes out to be the commutator of a coupling matrix and the mass-squared matrix of the charged leptons. However there are a number of variations of this model and the condition $`m_{\nu _1}+m_{\nu _2}+m_{\nu _3}=0`$ fails to hold in most of them. of two matrices $`A`$ and $`B`$.

## II “Zero Sum” Condition and Neutrino Masses

Direct measurement of neutrino masses is a very difficult experimental task. If the “zero sum” condition $`m_{\nu _1}+m_{\nu _2}+m_{\nu _3}=0`$ or equivalently $`TrM=0`$, is applied, all the neutrino masses are determined in terms of the mass-squared differences. We have

$`m_{\nu _1}^2={\displaystyle \frac{1}{3}}\left[2\mathrm{\Delta }m_{21}^2+\mathrm{\Delta }m_{32}^2\pm 2\sqrt{(\mathrm{\Delta }m_{32}^2)^2+\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{32}^2+(\mathrm{\Delta }m_{21}^2)^2}\right],`$ (3)
$`m_{\nu _2}^2={\displaystyle \frac{1}{3}}\left[\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{32}^2\pm 2\sqrt{(\mathrm{\Delta }m_{32}^2)^2+\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{32}^2+(\mathrm{\Delta }m_{21}^2)^2}\right],`$ (4)
$`m_{\nu _3}^2={\displaystyle \frac{1}{3}}\left[\mathrm{\Delta }m_{21}^2+2\mathrm{\Delta }m_{32}^2\pm 2\sqrt{(\mathrm{\Delta }m_{32}^2)^2+\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{32}^2+(\mathrm{\Delta }m_{21}^2)^2}\right].`$ (5)
The choice of the signs in front of the square root is decided by the requirement that all $`m_{\nu _i}^2`$ must be larger or equal to zero. From the expression for $`m_{\nu _3}^2`$, one can determine that the “+” has to be chosen for the expressions for $`m_{\nu _{2,3}}`$, and “$``$” for the expression for $`m_{\nu _1}`$. The relative signs of the eigen-masses $`m_i`$ are determined by the condition (2). We will use a convention such that $`m_{\nu _3}0`$ in our later discussions.

For $`|\mathrm{\Delta }m_{21}^2||\mathrm{\Delta }m_{32}^2|`$, the above simplify to, to leading order in $`\epsilon \mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{32}^2`$, we have

$`m_{\nu _1}^2={\displaystyle \frac{2}{3}}[|\mathrm{\Delta }m_{32}^2|(1+{\displaystyle \frac{\epsilon }{2}}){\displaystyle \frac{1}{2}}\mathrm{\Delta }m_{32}^2(1+2\epsilon )],`$ (6)
$`m_{\nu _2}^2={\displaystyle \frac{2}{3}}[|\mathrm{\Delta }m_{32}^2|(1+{\displaystyle \frac{\epsilon }{2}}){\displaystyle \frac{1}{2}}\mathrm{\Delta }m_{32}^2(1\epsilon )],`$ (7)
$`m_{\nu _3}^2={\displaystyle \frac{2}{3}}(|\mathrm{\Delta }m_{32}^2|+\mathrm{\Delta }m_{32}^2)(1+{\displaystyle \frac{\epsilon }{2}}).`$ (8)
At present the sign of the measured $`\mathrm{\Delta }m_{21}^2`$ is determined to be positive, but the sign of $`\mathrm{\Delta }m_{32}^2`$ is not determined, there are two possible solutions corresponding to the sign of $`\mathrm{\Delta }m_{32}^2`$. In Table I we list all solutions for the best fit values of the mass-squared differences.

We see that the mass eigenvalues exhibit two types of hierarchies,

$`I)m_{\nu _3}2m_{\nu _1}2m_{\nu _2}0.064\text{eV}`$ (9)
$`II)m_{\nu _1}m_{\nu _2}0.054\text{eV},\text{and}m_{\nu _3}0.00064\text{eV}.`$ (10)
The sign of $`\mathrm{\Delta }m_{32}^2`$ decides which mass hierarchy the solutions belong to. Note that the “natural” sign $`\mathrm{\Delta }m_{32}^2>0`$ corresponds to scenario I), in which the masses are of the same order of magnitude, in contrast to scenario II), in which $`m_{\nu _3}`$ is two order of magnitude smaller than $`m_{\nu _1}`$ and $`m_{\nu _2}`$. We would like to suggest that I) is more favored than II).

In contrast to oscillation experiments, the contribution of the neutrinos to the energy density of the universe, $`\mathrm{\Omega }_\nu _i|m_i|/(90\text{eV})`$ depends on the absolute values of $`|m_i|`$ of course. With the neutrino masses predicted with the central values of $`\mathrm{\Delta }m_{21}^2`$ and $`\mathrm{\Delta }m_{32}^2`$, we would have, $`\mathrm{\Omega }_\nu 1`$. Even if one uses the upper bound of $`\mathrm{\Delta }m_{32}^2=5.0\times 10^3`$, the contribution to the energy density is still much smaller than one. Other astrophysical and cosmological data can give information on the neutrino masses, such as the CMB power spectrum and large scale structure survey data. The absolute neutrino mass sum $`|m_{\nu _1}|+|m_{\nu _2}|+|m_{\nu _3}|`$ predicted by the “zero sum” condition is only a factor of two smaller than the present bound of 0.23 eV obtained from combining WMAP and Galaxy sky survey data, and is close to the sensitivity of 0.12 eV of the combined PLANCK CMB data with the SDSS sky survey. A future sky survey with an order of magnitude larger survey volume would allow the sensitivity to reach 0.03 eV. The mass ranges predicted by the condition (2) may be tested in the future.

## III Mixing and Masses

To obtain more information about neutrino properties, one needs to have information from mixing. Before the SNO and KamLAND data, assuming three active neutrino oscillations, one of the solutions which can account for known data is the bi-maximal mixing matrix with $`|V_{e2}|=|V_{\mu 3}|=1/\sqrt{2}`$. Experimental data from SNO and the recent data from KamLAND however disfavor the maximal mixing for the $`V_{e2}`$ entry. We were thus led to propose the following mixing matrix,

$`V=\left(\begin{array}{ccc}\hfill \frac{2}{\sqrt{6}}& \hfill \frac{1}{\sqrt{3}}& \hfill 0\\ \hfill \frac{1}{\sqrt{6}}& \hfill \frac{1}{\sqrt{3}}& \hfill \frac{1}{\sqrt{2}}\\ \hfill \frac{1}{\sqrt{6}}& \hfill \frac{1}{\sqrt{3}}& \hfill \frac{1}{\sqrt{2}}\end{array}\right).`$ (14)
This mixing matrix (but with the first and second column interchanged) was first suggested by Wolfenstein more than 20 years ago. It has subsequently been studied extensively by Harrison, Perkins and Scott, and Xing.

As mentioned before, oscillation experiments can not determine the relative signs of the mass eigenvalues which implies that one can multiply a phase matrix $`P=Diag(e^{i\rho },e^{i\sigma },1)`$ from the right on $`V`$. With CP invariance, $`\sigma `$ and $`\rho `$ can take the values of zero or $`\pm \pi /2`$.

With the above information on the mixing matrix, let us estimate two observables related to neutrino mass measurements, the effective mass electron neutrino mass $`m_e^2`$ measured in tritium beta decay end point spectrum experiments, and the effective Majorana electron neutrino mass $`m_{ee}`$ in neutrinoless double beta decays.

The effective mass measured in tritium end point spectrum experiments is given by, $`m_e^2=(m_{\nu _1}^2|V_{e1}|^2+m_{\nu _2}^2|V_{e2}|^2+m_{\nu _3}^2|V_{e3}|^2)`$. Using the values for the neutrino masses in Table I, we find that $`m_e`$ is below the sensitivity of $`0.12`$ eV for the proposed experiment KATRIN. However, neutrinoless double decay experiments may be sensitive to the predicted ranges. The amplitude of neutrinoless double beta decay is proportional to the effective electron-neutrino Majorana mass $`m_{ee}`$, namely $`M_{11}`$, given by

$`|m_{ee}|=|m_{\nu _1}V_{e1}^2e^{2i\rho }+m_{\nu _2}V_{e2}^2e^{2i\sigma }+m_{\nu _3}V_{e3}^2|.`$ (15)
From the above expression we see that the value $`m_{ee}`$ depends on the Majorana phases $`\rho `$ and $`\sigma `$. If the phases $`\rho `$ and $`\sigma `$ are all zero, the neutrinoless double beta decays would have the smallest rate because the cancellation imposed by the “zero sum” condition. If the phases take values such that both $`m_{\nu _1}e^{2i\rho }`$ and $`m_{\nu _2}e^{2i\sigma }`$ are positive, the neutrinoless double beta decays would have the largest rate possible. To have an idea of what possible effects of the Majorana phases can have on neutrinoless double beta decays, we calculated the range for $`|m_{ee}|`$ for the best fit values of the mass squared differences with $`V`$ given in eq.(14) for arbitrary $`\rho `$ and $`\sigma `$. The allowed ranges are also listed in Table I at the last column. The ranges obtained are well below the present upper bounds of 0.4 eV, but can almost be fully covered by future experiments, such as GENIUS, MAJORANA, EXO, MOON or COURE, where sensitivity as low as 0.01 eV seems possible.

## IV Neutrino Mass Matrix

We now study implications of the “zero sum” neutrino mass condition on the neutrino mass matrix $`M`$. If the neutrino masses and mixing matrix are known to a good precision, one can invert the eigen-masses according to eq. (1) to obtain the mass matrix $`M`$. To have some feeling how this may provide important information about the mass matrix, we present some details of the mass matrices which produce the mixing matrix $`V`$ in eq. (14).

The most general mass matrix which can produce the mixing matrix $`V`$ in eq. (14) can be specified in the following form by the mass eigenvalues,

$`M={\displaystyle \frac{m_{\nu _1}}{6}}\left(\begin{array}{ccc}4& 2& 2\\ 2& 1& 1\\ 2& 1& 1\end{array}\right)+{\displaystyle \frac{m_{\nu _2}}{3}}\left(\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right)+{\displaystyle \frac{m_{\nu _3}}{2}}\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 1\\ 0& 1& 1\end{array}\right).`$ (25)
Being real symmetric (and so a fortiori Hermitian) the above three matrices generate a $`U(1)U(1)U(1)`$ subgroup of $`U(3)`$.

With the “zero sum” condition $`m_{\nu _1}+m_{\nu _2}+m_{\nu _3}=0`$, the above matrix can be written as

$`M={\displaystyle \frac{m_{\nu _1}}{3}}\left(\begin{array}{ccc}2& 1& 1\\ 1& 1& 2\\ 1& 2& 1\end{array}\right)+{\displaystyle \frac{m_{\nu _2}}{3}}\left(\begin{array}{ccc}1& 1& 1\\ 1& 1/2& 5/2\\ 1& 5/2& 1/2\end{array}\right).`$ (32)
Since $`|\mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{32}^2|1`$, it is instructive to work with the case $`\mathrm{\Delta }m_{21}^2=0`$ as the first approximation and to see how the obtained mass matrix can be perturbed to produce the desired $`\mathrm{\Delta }m_{21}^2`$. With this approximation the two mass hierarchies I) and II) are then $`(m_{\nu _1}=m_{\nu _2}=m_{\nu _3/}2=2a)`$, and $`(m_{\nu _1}=m_{\nu _2}=2a,m_{\nu _3}=0)`$, respectively. The corresponding mass matrices are given by

$`I)M_0=a\left(\begin{array}{ccc}2& 0& 0\\ 0& 1& 3\\ 0& 3& 1\end{array}\right);II)M_0={\displaystyle \frac{1}{3}}a(\begin{array}{ccc}2& 4& 4\\ 4& 1& 1\\ 4& 1& 1\end{array}).`$ (39)
Note that the unperturbed mass matrix $`M_0`$ looks “simpler” in the “natural” hierarchy I) than in the “inverted” hierarchy II).

The mass matrix in case I) has been studied in a previous paper by us. The desired mass-squared difference $`\mathrm{\Delta }m_{21}^2`$ can be obtained by a small perturbation of the form

$`\delta M_T=\epsilon a\left(\begin{array}{ccc}0& 1& 1\\ 1& 0& 1\\ 1& 1& 0\end{array}\right),`$ (43)
with the perturbed eigenvalues given by $`m_{\nu _1}=2a(1\epsilon /2)`$, $`m_{\nu _2}=2a(1+\epsilon )`$, and $`m_{\nu _3}=4a(1+\epsilon /4)`$. The parameter $`\epsilon `$ to the lowest order is given by $`\epsilon =\mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{32}^2`$.

For case II), adding $`\delta M_T`$ to $`M_0`$, the eigen-masses are given by $`m_{\nu _1}=2a(1\epsilon /2)`$, $`m_{\nu _2}=2a(1\epsilon )`$, and $`m_{\nu _3}=a\epsilon `$ with $`\epsilon \mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{32}^2`$.

The perturbation $`\delta M_T`$ preserves the “zero sum” condition in eq. (2). One can also consider situations where the perturbations break this “zero sum” condition, but still produce the mixing matrix $`V`$ in eq. (14). For example, adding a “democratic” perturbation,

$`\delta M_D=\epsilon a\left(\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right),`$ (47)
produces a mixing matrix of the form given by $`V`$ in eq. (14), but different mass eigenvalues, ($`m_{\nu _1}=2a`$, $`m_{\nu _2}=2a(1+3\epsilon /2)`$, $`m_{\nu _3}=4a`$, $`\epsilon =\mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{32}^2`$), and ($`m_{\nu _1}=2a`$, $`m_{\nu _2}=2a(13\epsilon /2)`$, $`m_{\nu _3}=0`$, $`\epsilon =\mathrm{\Delta }m_{21}^2/3\mathrm{\Delta }m_{32}^2`$) for case I) and case II), respectively.

We would like to comment that the reconstruction of mass matrix depends crucially on knowing the mixing matrix even if all the eigenvalues for the masses are known through experimental measurements or theoretical considerations such as $`m_{\nu _1}+m_{\nu _2}+m_{\nu _3}=0`$ or the traceless condition for real Majorana neutrino mass matrix proposed here. Had the mixing matrix been given by the bi-maximal mixing matrix which was allowed by pre-SNO and KamLAND data, the mass matrix would be very different. We mention that the simplest version of the model proposed in ref. with all diagonal entries zero in the mass matrix can easily produce a bi-maximal mixing, but it is difficult to produce the mixing matrix given by $`V`$ in eq. (14).

## V Conclusions

We have studied the consequences of neutrino masses with the “zero sum” condition $`m_{\nu _1}+m_{\nu _2}+m_{\nu _3}=0`$. With this condition the neutrino masses can be determined from measured mass-squared differences from oscillation experiments. We find that this condition predicts only two types of neutrino mass hierarchies with one of them characterized by $`m_{\nu _3}2m_{\nu _1}2m_{\nu _2}0.063`$ eV, and another by $`m_{\nu _1}m_{\nu _2}0.054`$ eV and $`m_{\nu _3}0.0064`$ eV. These masses although small, can be probed by experiments from CMB measurements and large scale structure survey, and can also be probed by neutrinoless double beta decay experiments. In conjunction with information on neutrino mixing, the “zero sum” condition also predicts simple mass matrices for neutrinos.

###### Acknowledgements.

We thank H. Murayama, S. Pakvasa, R. Volkas and L. Wolfenstein for helpful discussions and comments on the manuscript. This work was supported in part by NSC under grant number NSC 91-2112-M-002-42, and by the MOE Academic Excellence Project 89-N-FA01-1-4-3 of Taiwan, and by NSF under grant number PHY 99-07949 of USA. AZ thanks Professor Pauchy Hwang and the Department of Physics of the National Taiwan University, where this work was initiated, for warm hospitality.
Note Added The Ansatz discussed in this paper has been proposed earlier by D. Black, A. Fariborz, S. Nasri and J. Schechter. Their theoretical motivation and the values of the masses they obtained are however rather different.


text:
# Avoided Antiferromagnetic Order and Quantum Critical Point in CeCoIn5

## Abstract

We measured specific heat and resistivity of heavy fermion CeCoIn5 between the superconducting critical field $`H_{c2}=5T`$ and 9 T, with field in the  direction, and at temperatures down to 50mK. At 5T the data show Non Fermi Liquid behavior down to the lowest temperatures. At field above 8T the data exhibit crossover from the Fermi liquid to a Non Fermi Liquid behavior. We analyzed the scaling properties of the specific heat, and compared both resistivity and the specific heat with the predictions of a spin-fluctuation theory. Our analysis leads us to suggest that the NFL behavior is due to incipient antiferromagnetism (AF) in CeCoIn5 with the quantum critical point in the vicinity of the $`H_{c2}`$. Below $`H_{c2}`$ the AF phase which competes with the paramagnetic ground state is superseded by the superconducting transition.

preprint: submitted to Physical Review Letters LA-UR-03-0364

When the symmetry of the ground state of a system changes as a function of an external or internal parameter, the system is said to undergo a quantum phase transition. If, in addition, this transition is second order, the system has a Quantum Critical Point (QCP) at the critical value of the parameter. The competition between the nearly degenerate ground states determines the behavior of the system over a range of temperatures and tuning parameter values in the vicinity of QCP. In this region of the phase diagram the properties of the system differ from those on either side of the transition, and often exhibit unusual dependence on the temperature and the tuning parameter. This has made quantum critical phenomena a subject of intense current interest.

Study of quantum critical points in heavy fermion systems has been a focus of particular attention (for a recent review see Ref. Stewart, 2001). In these materials the competition typically takes place between a paramagnetic and a magnetically ordered ground states. The unconventional behavior near QCP is manifested in the deviation of the temperature dependence of measured properties from those of metals described by the Landau Fermi Liquid (FL) theory. In that theory the electronic specific heat is linear in temperature, $`C(T)=\gamma T`$, and the resistivity increases quadratically from a residual value, $`\rho =\rho _0+AT^2`$. In systems tuned to QCP the Sommerfeld coefficient, $`\gamma (T)=C/T`$, commonly diverges as the temperature goes to zero, and has been variously argued to behave as either $`\mathrm{log}T`$ or $`T^\alpha `$, with $`\alpha <0`$. Resistivity with an exponent less than two is also ubiquitous in these compounds.

Tuning the system through a QCP can be accomplished experimentally by varying sample’s composition v. Löhneysen et al. (1994); Maple et al. (1995), applying pressure Mathur et al. (1998), or applying magnetic field Grigera et al. (2001). In non-stoichiometric compounds the Kondo disorder, where a range of Kondo temperatures $`T_K`$ appears due to different environments of the f-electron ions, is an important mechanism leading to a Non Fermi Liquid (NFL) behavior Bernal et al. (1995); MacLaughlin et al. (1996). In these compounds it is not easy to separate this origin of NFL behavior from the consequences of the proximity to a QCP. Hence the stoichiometric compounds receive great deal of attention in the field of quantum criticality.

One class of such materials are Ce-based compounds, which have an antiferromagnetic (AF) ground state at ambient pressure. Hydrostatic pressure suppresses the magnetic ordering temperature $`T_N`$ to zero at a critical pressure of the QCP. Such approach was used successfully for CeCu2Ge2 Jaccard et al. (1992) CeRh2Si2 Movshovich et al. (1996), CePd2Si2, CeIn3 Mathur et al. (1998), and some other compounds. Alternatively, AF order can be suppressed by an applied magnetic field. When this was done in YbRh2Si2, an NFL behavior was revealed again on the paramagnetic side of the QCP Gegenwart et al. (2002).

In this Letter we present the results of the specific heat and resistivity measurements on stoichiometric CeCoIn5, which show that this compound has a QCP with magnetic field as a tuning parameter. We show that the quantum critical behavior is most likely due to the proximity of an antiferromagnetic state. This case is particularly interesting because CeCoIn5 is an ambient pressure superconductor. We argue that magnetic order which competes with the paramagnetic state is superseded by the superconducting state.

CeCoIn5 is a tetragonal, quasi-2D compound, with layers of CeIn3 separated by layers of CoIn2. It is an ambient pressure heavy fermion superconductor Petrovic et al. (2001) with $`T_c=2.3`$ K, the highest value for this class of compounds. Superconductivity in CeCoIn5 is unconventional, with lines of nodes in the energy gap, as demonstrated by specific heat and thermal conductivity Movshovich et al. (2001) and NQR Kohori et al. (2001) measurements. Pauli limiting analysis Movshovich et al. (2002) and thermal conductivity modulations in magnetic field Izawa et al. (2001) indicate that CeCoIn5 is a singlet, $`d_{x^2y^2}`$ superconductor. Recently it was shown that the superconducting transition in CeCoIn5 in magnetic fields close to the upper critical field for  direction, $`H_{c2}=4.95T`$, is first order below $`T=0.7`$Izawa et al. (2001); Bianchi et al. (2002). NFL behavior of CeCoIn5 at the field of 5 T persists over a large region of temperature Kim et al. (2001), with $`\gamma \mathrm{log}(T)`$ between 0.4 K and 8 K.

Here we concentrate on the low temperature region down to 50 mK and magnetic fields between $`H_{c2}=4.95`$ T and 9 T. Figure 1(a) shows Sommerfeld coefficient of CeCoIn5 versus temperature for several magnetic fields between 5 T and 9 T. Large low temperature Schottky anomaly tails due to the quadrupolar and magnetic spin splitting of In and Co nuclei were subtracted Movshovich et al. (2001). At high magnetic field and low temperature Schottky term accounts for most of the specific heat measured, and its subtraction results in an increased scatter of the data points at low temperature. The data for 5 T follow logarithmic temperature behavior below 1 K down to the lowest temperature studied. As the magnetic field is increased above 5 T, the low temperature data starts to deviate from the logarithmic behavior at ever higher temperature. Clear FL regime of $`\gamma =`$ const. at low temperature is recovered at 8 Tesla. Both 8 T and 9 T data exhibit sharp crossovers between the NFL (logarithmic) and the FL regimes at a temperature $`T_{FL}^C`$. For fields below 8 T $`\gamma `$ continues to rise with decreasing temperature and does not show saturation within the temperature range shown. This behavior is different from that of YbRh2Si2, where $`\gamma `$ saturates below some temperature $`T^{}(H)`$ to a constant value for all $`HH_c`$, with a rather sharp break from the $`\mathrm{log}(T)`$ behavior for $`T>T^{}`$ Trovarelli et al. (2000).

The logarithmic divergence of $`\gamma (T,H)`$ at 5T and $`\gamma `$’s increase followed by a crossover to the FL behavior at low $`T`$ for higher fields is indicative of a field-tuned QCP near 5T. This observation is strongly supported by the scaling properties of $`\gamma `$. For a phase transition with $`T_c=0`$, the temperature alone sets the energy scale of fluctuations. Scaling of dynamical properties with energy, $`E`$, as $`E/T`$ is taken as good evidence of quantum criticality Varma et al. (1989); Sachdev (1999).

In our case the tuning parameter is the applied magnetic field. Hence the energy scale of fluctuations away from the QCP depends on $`|HH_c|`$. Therefore the entropy, and hence the specific heat, should scale as $`\gamma (H)\gamma (H_c)f((HH_c)/T^\beta )`$. Experimentally, at fields sufficiently away from the QCP, we achieve the FL regime at low $`T`$. In that regime $`\gamma `$ depends on the field but is $`T`$-independent. Consequently, the scaling relation has to take the form $`\gamma (H)\gamma (H_c)(HH_c)^\alpha f((HH_c)/T^\beta )`$. The form of scaling is similar to that obtained for U0.2Y0.8Pd3 Andraka and Tsvelik (1991), CeCu5.8Ag0.8 Heuser et al. (1998), and for YbRh2Si2 Trovarelli et al. (2000).

We find that with the choice of $`H_c=5`$ T, the best scaling is achieved for $`\alpha =0.71`$ and $`\beta =2.5`$; it is shown in Fig. 1(b). It is a remarkably good scaling, spanning both FL and NFL regimes, with all four data sets for different fields overlapping each other in the entire experimental temperature range. This scaling is a strong indication that the behavior of CeCoIn5 in the part of the phase space explored in these experiments is governed by the QCP very close to $`H_c=5`$ T. The high power of temperature in the argument of the scaling function (the scaling dimension of the magnetic field) suggests that the field is very efficient in suppressing the critical fluctuations.

Importantly, scaling implies that the behavior observed here is associated with a second order transition. Therefore it is unlikely that first order transition from the superconducting to normal state above 4.7 T Bianchi et al. (2002) controls the properties in the parameter range investigated here, and we need to consider alternative competing orders.

The crossover from the NFL to FL regimes is also clear from resistivity measurements in the same magnetic field range. Fig. 2(a) shows resistivity, $`\rho `$, of CeCoIn5 below 5 K for fields between 5.6 T and 9 T, as well as the zero-field data. The data at 5.6 T and above is not affected by the proximity to the superconducting transition, as seen from the magnetoresistance measurements of CeCoIn5 at 100 mK. Movshovich  To focus on the low temperature behavior, we show the data below 1 K in Fig. 2(b). The data for 8 T below 300 mK and for 9 T below 500 mK display the FL-like $`\mathrm{\Delta }\rho \rho \rho (T=0)=AT^2`$ behavior, consistent with the linear in $`T`$ specific heat data for these fields. In contrast, for fields of 5.6 T and 6 T we found $`\mathrm{\Delta }\rho T`$ below 200 mK, and the data for 7 T show intermediate behavior, the best fit being the sum of linear and quadratic terms in this temperature range. At these fields the linear in $`T`$, rather than quadratic, dependence of the resistivity is consistent with the NFL behavior of the specific heat. Qualitative difference in the temperature behavior of $`\rho (T)`$ at 6T and at 9T in the low temperature range is emphasized in the upper inset in Fig. 2(b), where resistivity is plotted versus $`T^2`$ for two fields, 9 T and 6 T. The 9 T data are linear in $`T^2`$ over the entire temperature range shown, while the 6 T data has a pronounced negative curvature.

To further analyze the FL region of the phase diagram of CeCoIn5 in the inset (d) of Fig. 2(b) we plot $`\gamma ^2`$ versus $`A`$ for the 8 T and 9 T data. Within error bars the two values lie on the line $`\gamma ^2=.3\times A`$ $`\mathrm{J}^2/\mathrm{mol}^2\mathrm{K}^2\mu \mathrm{\Omega }\mathrm{cm}`$. The slope is about three times larger than the Kadowaki-Woods ratio Kadowaki and Woods (1986), approximately obeyed by many heavy fermion compounds. This discrepancy is within the typical range of the scatter of the data for other HF compounds, and confirms the FL ground state of CeCoIn5 for 8 T and 9 T.

Having established the existence of the QCP near $`H_c5`$ T, we proceed to identify the order competing with the paramagnetic FL ground state. Since a closely related compound, CeRhIn5, is an ambient pressure antiferromagnet with the Neel temperature $`T_N=3.8`$ K, AF order is a natural candidate. It is also consistent with the $`d_{x^2y^2}`$ superconductivity.

While the complete theory of quantum criticality in itinerant antiferromagnets is still lacking, the singular contribution of the critical AF spin fluctuations to the specific heat and scattering has been considered in Ref. Moriya and Takimoto, 1995 in the framework of a self-consistent renormalization (SCR) theory. This theory was used to analyze the specific heat of CeCu5.2Ag0.8 near the QCP Stewart (2001). The input parameter of the theory is the distance from a QCP, $`y_0`$, ($`y_0=0`$ at the QCP), related to the inverse square of the magnetic correlation length. The specific heat and the resistivity are functions of the reduced temperature, $`T/T_0`$, where $`2\pi ^2T_0`$ is of the order of the exchange interaction. We use this theory to fit our data.

Close to the QCP we obtained a good fit to both the specific heat data (shown in Fig. 1(a)), and resistivity (shown in the inset of Fig. 2(a)) for the same parameter values. To estimate $`T_0`$, we compared CeCoIn5 with CeRhIn5. Pressure of 16 kbars suppresses the AF state in CeRhIn5, at which point the superconducting state with $`T_c=2.1`$ K emergesHegger et al. (2000). This lead Sidorov et al. to suggest that physical properties of CeCoIn5 are very close to those of CeRhIn5 at 16 kbars Sidorov et al. (2002). The $`T_0=0.4`$ K we chose is consistent with the exchange energy of the order of a few Kelvin, expected from the $`T_N=3.8`$ K for CeRhIn5. In the fit we assumed an additive field-independent contribution to $`\gamma (T)`$ of 0.2 J/mol K2 from non-critical fermions. For the 5 T data the theoretical fits are identical for any $`y_00.01`$ over the temperature range studied. Consequently, CeCoIn5 is very close to a QCP at this field.

Since in our experiments the magnitude of the applied field (used to tune the system to the QCP) is comparable to the exchange interaction, we expect that the changes in the spectrum of the spin fluctuations and the quasiparticle scattering rate lead to deviations from the predictions of the SCR theory when the field is increased away from QCP. Indeed, at higher fields the discrepancy between the fits and the resistance data becomes more pronounced (especially at higher temperature). However, the overall S-shape of the resistance curves agrees with the data. We therefore conjecture that the QCP is due to the tendency towards the AF order.

Finally, in Fig. 3 we plot the phase diagram of CeCoIn5 in the $`TH`$ plane. The normal-superconducting phase boundary was determined in Refs. Petrovic et al., 2001; Bianchi et al., 2002. At fields of 8T and 9T we see a clear transition from the NFL to FL behavior in both specific heat ($`T_{FL}^C`$) and resistivity ($`T_{FL}^\rho `$). Below 8T, in the temperature range studied, the heat capacity and resistance deviate from the critical behavior (5 T data), although the FL regime is not reached. For these fields we define a temperature $`T_{cr}`$ below which $`\gamma `$ begins to deviate from $`\mathrm{log}(T)`$ behavior, leading to a broad crossover region. If the FL behavior exists at these fields, it occurs over a narrow range of very low temperatures.

In conclusion we find that, while there is no signature of the AF ordered phase, the behavior of CeCoIn5 above the superconducting upper critical field, $`H_{c2}=4.95`$T, at low temperatures is controlled by a field-tuned antiferromagnetic quantum critical point. The existence of the QCP is confirmed by scaling analysis of the specific data. Incipient AF state in CeCoIn5 is strongly suggested by the the fits to both the specific heat and the resistivity based on the spin-fluctuation theory Moriya and Takimoto (1995). The underlying physical picture is that the AF spin fluctuations promote anti-alignment of the electron spins. Whether time-reversal breaking (antiferromagnetism) or singlet formation with phase coherence (superconductivity) occurs first depends on the details of the system. In CeCoIn5 the superconductivity prevails, preventing the magnetic order from developing. However, there are no critical fluctuations associated with the (first order Bianchi et al. (2002)) superconducting transition at $`H_{c2}`$, so that the behavior above $`H_{c2}`$ is controlled by the proximity to the AF critical point. Our results imply that the transition temperature of the avoided antiferromagnetic order vanishes at fields very close to the superconducting upper critical field. This concomitant suppression of the AF and SC orders is unusual, and may reflect an important aspect of the underlying physics, presenting further challenges to theory and experiment.

We thank J. D. Thompson, C. Varma, Q. Si, P. Coleman, and G. R. Stewart for stimulating discussions. We are also grateful to J. D. Thompson for help with selecting Indium-free samples. Work at Los Alamos National Laboratory was performed under the auspices of the U.S. Department of Energy.