text:
# Influence of Carbon Concentration on the Superconductivity in MgCxNi3

## Abstract

The influence of carbon concentration on the superconductivity ( SC ) in MgCxNi3 has been investigated by measuring the low temperature specific heat combined with first principles electronic structure calculation. It is found that the specific heat coefficient $`\gamma _n=C_{en}/T`$ of the superconducting sample ($`x1`$) in normal state is twice that of the non-superconducting one ($`x0.85`$). The comparison of measured $`\gamma _n`$ and the calculated electronic density of states ( DOS ) shows that the effective mass renormalization changes remarkably as the carbon concentration changes. The large mass renormalization for the superconducting sample and the low $`T_c`$( 7K ) indicate that more than one kind of boson mediated electron-electron interactions exist in MgCxNi3.

The recently discovered superconductor MgCNi3 HeT2001  with an anti-perovskite structure has attracted renewed attention, because it is a new intermetalic compound following the remarkable discovery of MgB2 NagamatsuJ2001  and is thought to be very close to ferromagnetism (FM) due to a large proportion of Ni in each unit cell RosnerH2002 . Experiments have shown that the superconductivity (SC) in MgCxNi3 is very sensitive to the content of carbon and disappears below about $`x=0.88`$ AmosTG2002 . Moreover, Ni-site doping with Cu (electron doping) and Co (hole doping) have very different effectHaywardMA2001 .

Many efforts have been focused on the origin of SC in MgCNi3. Some researchers explored the possibility of conventional phonon-coupled pairing ShimJH2001 ; DugdaleSB2001 . However, various experiments and theoretical calculations have not given a consistent coupling constant $`\lambda `$ in the frame of electron-phonon coupling (EPC) mechanism. For example, $`\lambda 11.6`$ was estimated from the calculated plasma frequency ($`\omega _p`$) and measured resistivity SinghDJ2001 , while the measured upper critical field $`H_{c2}`$ LiSY2001 ; WalteA2002 ; MaoZQ2002  suggested a much larger value of $`\lambda 2.53`$. Furthermore, two discrepant values of $`\lambda =0.7`$ HeT2001  and 2.2 WalteA2002  were determined in the similar specific heat measurements by taking different approximations. Besides these confusing results, it has been proposed that EPC is possibly not the only contribution to the coupling mechanism in MgCNi3 system SinghDJ2001 ; WalteA2002 . Opposed to the conventional pairing mechanism, it was suggested that FM and SC may coexist due to the ferromagnetic instability caused by the high density of electronic states (DOS) at the Fermi energy (EF) RosnerH2002 ; SheinIR2001 . NMR experiment confirmed the presence of substantial ferromagnetic spin fluctuations in MgCNi3 SingerPM2001 , while it also exhibited a clear coherence peak of the nuclear spin-lattice relaxation rate just below $`T_c`$, which is typical for an isotropic s-wave superconductor. More recently, Voelker et al. VoelkerK2002  proposed that MgCNi3 could be a multiband superconductor with a conventional phonon mechanism AgterbergDF1999  similar to MgB2 ChoiHJ2002 . Thus, the origin of SC in MgCNi3 is still controversial. One of the most illuminative methods is to study the influence of element-substitution on SC.

In this Letter, we compare the specific heat data of a superconducting MgCxNi3 with $`x1`$ and that of a non-superconducting one with $`x0.85`$. Using the tight-binding linear muffin-tin orbital (TB-LMTO) band method and coherent potential approximation (CPA) AndersenOK1985 ; TurekI1997 , we investigate the carbon concentration dependence of the electronic structure of MgCxNi3. Combining the results of the specific heat and band calculations, we propose a self-consistent picture to understand the reduced SC in non-stoichiometric MgCxNi3.

Poly-crystalline samples of MgCxNi3 were prepared by powder metallurgy method. Details of the preparation were published previously RenZA2002 . The superconducting transition of the stoichiometric sample( SC1 ) with $`T_c=7K`$ is shown in Fig. 1. A non-superconducting sample (NSC08) was synthesized by continually reducing the carbon component until the diamagnetism was completely suppressed. X-ray diffraction (XRD) patterns show nearly single phase in both samples. The lattice parameters $`a`$ determined from XRD are 3.812$`\AA `$ and 3.790$`\AA `$ for SC1 and NSC08, respectively. By comparing these lattice parameters with the reported $`ax`$ relation AmosTG2002 , we could estimate the carbon content as $`x=0.977`$ and 0.850 for SC1 and NSC08, respectively. According to the $`T_cx`$ relation presented in Ref AmosTG2002 , $`x=0.977`$ corresponds to $`T_c7K`$ and $`x=0.850`$ is very close to the critical stoichiometry of SC in MgCxNi3, which is in good agreement with our experimental results. The temperature dependence of the normalized resistivity is shown in Fig. 1. The heat capacity data presented here were taken with the relaxation method Bachmann1972  based on an Oxford cryogenic system Maglab, in which the heat capacity is determined by a direct measurement of the thermal time constant, $`\tau =(C_s+C_{add})/\kappa _w`$, here $`C_s`$ and $`C_{add}`$ are the heat capacity of the sample and addenda, respectively, while $`\kappa _w`$ is the thermal conductance between the chip and a thermal link. During the measurement the sample was cooled to the lowest temperature under a magnetic field ( field-cooling ) followed by data acquisition in the warming up process. The addenda includes a small sapphire chip, a tiny Cernox temperature sensor, small amount of Wakefield thermal conducting grease and gold leads ($`\varphi 25\mu `$m). All contributions to the heat capacity from the addenda have been measured separately and subtracted from the total specific heat before further analysis.

The typical low temperature specific heat at various magnetic fields up to 12 Tesla is presented in Fig. 2. For sample NSC08, all the data fall into an universal curve, indicating the independence of its specific heat on the magnetic field. For sample SC1, all the data above $`T_c`$ are independent on the field, and the normal state extends to the whole temperature region investigated here when $`H10T`$. The upper critical field $`H_{c2}(T)`$ can be determined from the position of the specific-heat jump, and $`H_{c2}(0)11T`$ is then derived by extrapolating $`H_{c2}(T)`$ to 0K. In order to describe the specific heat data of NSC08 and the normal state data of SC1, we use the following expression:

$$C_n(T)=\gamma _nT+\beta T^3+\delta T^5$$
(1)
where the linear-T term is the electronic contribution with $`\gamma _n`$ as Sommerfeld parameter, the second term represents the phonon contribution according to the Debye approximation, the last term is required to include deviations from the linear dispersion of the acoustic modes in extended temperature range. The fitting results are presented in Fig. 2, and the obtained parameters are shown in Table 1. Besides the obvious difference between the Debye temperatures ($`\mathrm{\Theta }_D`$) of SC1 and NSC08 (due to the larger elastic modulus of NSC08 than that of SC1), a more prominent difference is found between the Sommerfeld parameters $`\gamma _n`$ of these two samples, this is one of the main results in this work.

Before further analysis, we at first investigate the specific heat data of the superconducting sample SC1 in zero field, as shown in the lower inset of Fig. 2. The solid curve represents the theoretical fits, the part above $`T_c`$ is the fitting result to the normal state data of SC1 as mentioned above, while the part below $`T_c`$ is an attempt to describe the data using the specific heat expression of conventional superconductor: $`C_{es}=Aexp(\mathrm{\Delta }(0)/k_BT)/T^{1.5}`$, in which $`A`$ is a function of zero-temperature energy gap $`\mathrm{\Delta }(0)`$ and the normal-state DOS $`g(0)`$. This equation is a low temperature approximation and is not adequate above $`0.8T_c`$ where $`\mathrm{\Delta }(T)`$ has substantially deviated from $`\mathrm{\Delta }(0)`$. This conventional expression can fit our data below $`0.8T_c5K`$ quite well, though the returned $`\mathrm{\Delta }(0)`$ (1.46meV) leads to a large value of $`2\mathrm{\Delta }(0)/k_BT_c5`$ comparing with the weak coupling value 3.5. This indicates that MgCNi3 is a strong coupling superconductor, being consistent with the large value of $`\mathrm{\Delta }C/\gamma _nT_c1.7`$ ( the expected value for weak coupling is 1.43 ), where $`\mathrm{\Delta }C`$ is the specific heat jump at $`T_c`$.

As shown above, the specific heat coefficient $`\gamma _n`$ of SC1 is almost twice that of NSC08. In the frame of strong coupling mechanism, given the value of $`\gamma _n`$ and $`N(E_F)`$, an effective mass renormalization factor ($`\lambda `$) can be estimated from the following relation MigdalA1958 ; RickayzenG1980

$$\gamma _n=\frac{\pi ^2}{3}k_B^2(1+\lambda )N(E_F)$$
(2)
In order to obtain $`N(E_F)`$, we have calculated the electronic structures using a TB-LMTO method based on Green function formalism TurekI1997  within the atomic sphere approximation (ASA). The parameterized exchange-correlation potential of Ref VoskoSH1980  and the scalar-relativistic Dirac equation KoellingDD1977  were employed. The coherent potential approximation(CPA) used in our calculation allowed us to study the electronic structure of MgCxN3 for any Carbon concentration. In the crystal structure of cubic antiperovskite-type MgCN3, the atoms occupy the positions Mg (0,0,0), C (0.5,0.5,0.5), and Ni (0.5,0.5,0), (0.5,0,0.5), (0,0.5,0.5). Self-consistency was reached using 1000 k-points within the irreducible wedge of the simple cubic Brillouin zone. In the calculations, we adopted the experimental determined AmosTG2002  relation between the lattice parameter $`a`$ and the carbon concentration $`x`$. A fixed radii proportion as $`2.078:1:1.617`$ (Mg:C:Ni) was taken for the stoichiometric situation ($`x=1`$), and only the radii of carbon-site was reduced with decreasing $`x`$.

The DOS of MgCxN3 for $`x=1`$ are presented in the upper inset of Fig. 3, which is consistent with the reports of other groups SinghDJ2001 ; ShimJH2001 ; SheinIR2002 ; Kim2002 . The $`\pi ^{}`$ antibonding states of Ni-$`3d`$ and C-$`2p`$ are located just below $`E_F`$, yielding a high DOS peak whose height is proved to be sensitive to the exact method. Fortunately, the obtained value of $`N(E_F)=4.53`$ eV-1cell-1 is not far away from the result of general potential linearized augmented planewave method (4.99 eV-1cell-1) SinghDJ2001  and that of full-potential LMTO method (4.57 eV-1cell-1) SheinIR2001 , which is a good starting point for the following consideration of the doping effect.

Fig. 3 shows the evolution of the DOS around Fermi level with the change of carbon concentration $`x`$ in MgCxNi3. It is found that the Fermi level slightly shifts towards low energy with the decreasing $`x`$, while the Ni$`3d`$-C$`2p`$ antibonding peak is depressed remarkably. In the nearest-neighbor approximation, the degeneration of this DOS peak indicate that the carbon-site vacancies locally break the hybridization between Ni-$`3d`$ and C-$`2p`$ orbitals and hence lead to a redistribution of the electronic states. Atom-resolved DOS shows that the near-$`E_F`$ DOS peak of Ni-$`3d`$ component and that of C-$`2p`$ component do not depart from each other with decreasing $`x`$, indicating that carbon-site vacancies does not destroy the Ni$`3d`$-C$`2p`$ bond for all carbon concentrations investigated here. Therefore, the decrease of the peak possibly means the enhanced itinerancy of electrons. As shown in Fig. 3, the depressed DOS peak just below $`E_F`$ results in a notable reduction of $`N(E_F)`$. The dependence of $`N(E_F)`$ on the carbon concentration is presented in Fig. 4. We can see a obvious decrease of $`N(E_F)`$ from 4.26eV-1 at $`x=0.977`$ to 3.14eV-1 at $`x=0.85`$. Moreover, the doping dependence of $`N(E_F)`$ is linear above a critical doping $`x_c0.88`$ which is very close to the carbon content in NSC08 and also the lower limit of the carbon concentration for SC in MgCxNi3 AmosTG2002 . Below $`x_c`$, the doping dependence of $`N(E_F)`$ becomes much weaker than that above $`x_c`$.

By inserting the experimental value of $`\gamma _n`$ and the calculated $`N(E_F)`$ into Eq( 2), we can determine the coupling constant as $`\lambda =2.37`$ and 1.35 for the samples SC1 and NSC08, respectively. The former is close to that reported by Wälte et al. WalteA2002 , while the remarkable difference in the coupling strength between a superconducting sample and a carbon-deficient non-superconducting one is for the first time found in the MgCxNi3 system. However, applying $`\lambda =2.37`$ to the McMillan’s $`T_c`$ formula McMillanWL1968  $`T_c=\frac{\mathrm{\Theta }_D}{1.45}e^{\frac{1+\lambda }{\lambda }}`$ gives rise to $`T_c>30K`$, which seems to be too large as compared to experimental value $`T_c=7K`$. In other words, $`T_c=7K`$ corresponds to $`\lambda =0.42`$, which is much smaller than the value determined here. This contradiction can be resolved if we adopt the picture that there are two kinds of boson-mediated electron-electron interactions, for example, electron-phonon coupling (EPC) and spin fluctuation (SF) SinghDJ2001 ; WalteA2002 . Due to the different effect of the spin fluctuations on mass renormalization and superconducting properties, it should add in the mass renormalization term and subtract in the pairing term for SC. So in the frame of strong coupling, the specific heat determined $`\lambda `$ has two contributions: $`\lambda =\lambda _{ph}+\lambda _{spin}`$. Similarly, the simplified McMillan’s $`T_c`$ formula should be modified as SinghDJ2001

$$T_c=\frac{\mathrm{\Theta }_D}{1.45}e^{\frac{1+\lambda _{ph}+\lambda _{spin}}{\lambda _{ph}\lambda _{spin}}}$$
(3)
Using the determined parameters for sample SC1 (namely, $`\lambda =2.2`$, $`T_c=7K`$ and $`\mathrm{\Theta }=351K`$), one obtains $`\lambda _{ph}\lambda _{spin}=0.95`$ from Eq.( 3). Combining this result with $`\lambda _{ph}+\lambda _{spin}=2.37`$, we determine the strength of both EPC and electron-paramagnon coupling as $`\lambda _{ph}=1.66`$ and $`\lambda _{spin}=0.71`$. Firstly, by assuming that $`\lambda _{spin}`$ is invariable with the decreasing carbon concentration, we could obtain $`\lambda _{ph}=0.64`$ for sample NSC08. This suggests a strong dependence of the EPC strength on the carbon concentration in MgCxNi3. For NSC08, $`\lambda _{ph}=0.64`$ has become smaller than $`\lambda _{spin}=0.71`$, which is closely associated with the disappearance of SC. Secondly, on the assumption that $`x=0.85`$ is the critical stoichiometry for SC, i.e., $`\lambda _{ph}=\lambda _{spin}`$, we obtain $`\lambda _{ph}=\lambda _{spin}=0.665`$. This indicates a $`6\%`$ decrease of $`\lambda _{spin}`$ from $`x=0.977`$ to $`x=0.85`$. However, $`\lambda _{ph}`$ is lost near $`60\%`$ in the same system, which seems to be responsible for the suppression of SC in MgCxNi3.

We have been aware of some recent calculations SheinIR2002 ; SheinIR2001  suggesting that both Co doping and Cu doping in MgCNi3 lead to a reduction of $`N(E_F)`$, which seems to be critical in depressing SC in this system (MgCNi3-xYx \[Y=Co or Cu\]). However, these investigations were based on the ordered substitutions, i.e., x=1, 2 or 3, which is somewhat away from the key problem because the solubility of Cu is limited to $`3\%`$ ($`x=0.1`$) in technology and the bulk SC disappears for only $`1\%`$ Co doping ($`x=0.03`$). Our further calculations did not show obvious reduction of $`N(E_F)`$ in this doping range. Therefore, some necessary investigations on the variation of the EPC in this system must be supplemented.

In summary, we have investigated the specific heat data in MgCxNi3 system and found a remarkable difference between the Sommerfeld parameters of a superconducting sample ($`x1`$) and a non-superconducting one ($`x0.85`$). Band calculations reveal a distinct decrease of the density of the electronic states at Fermi level with the decreasing carbon concentration. Further analysis indicate that there is more than one kind of boson-mediated electron-electron interactions existing in MgCxNi3 system. In the frame of strong coupling theory, a substantial depression of the electron-phonon coupling caused by the decrease of carbon concentration is for the first time found in this system and seems to be responsible for the impairment of its superconductivity.

###### Acknowledgements.

This work is supported by the National Science Foundation of China (NSFC 19825111, 10074078), the Ministry of Science and Technology of China ( project: NKBRSF-G1999064602 ), and Chinese Academy of Sciences with the Knowledge Innovation Project.


format: text

text:
# Research in Econophysics

## I Introduction

During the last several years, in addition to his research in condensed matter theory, Victor Yakovenko has been working in a recently emerged field of studies often called “econophysics”. Econophysics applies statistical physics methods to economical, financial, and social problems. Detailed references to the econophysics research in Victor Yakovenko’s group are given in his Web page http://www2.physics.umd.edu/~yakovenk/econophysics.html. The results have been published in refereed journals DY-money ; DY-income ; DY-wealth ; DY-Heston ; SY-Heston ; DY-statistical  and presented at international conferences and seminars. His research has been reviewed in a popular article in the American Scientist magazine Hayes , and Australian Financial Review, the leading Australian business newspaper, has published an op-ed column about his studies Price . His student Adrian Dragulescu received Ph.D. in 2002 and now works as a risk analyst at the Constellation Energy Group in Baltimore, which is the owner of Baltimore Electric and Gas Company. Currently Victor Yakovenko works with another graduate student A. Christian Silva.

## II Statistical Mechanics of Money, Income, and Wealth

In this Section, we overview Refs. DY-money ; DY-income ; DY-wealth ; DY-statistical , which use an analogy with statistical physics to describe probability distributions of money, income, and wealth.

The equilibrium statistical mechanics is based on the Boltzmann-Gibbs law, which states that the probability distribution function (PDF) of energy $`E`$ is $`P(E)=Ce^{E/T}`$, where $`T`$ is the temperature, and $`C`$ is a normalizing constant. The main ingredient in the textbook derivation of the Boltzmann-Gibbs law is conservation of energy. When two economic agents make a transaction, some amount of money is transferred from one agent to another, but the sum of their money before and after transaction is the same: $`m_1+m_2=m_1^{}+m_2^{}`$. Then, by analogy with statistical physics, the equilibrium PDF of money $`m`$ in a closed system of agents should have the Boltzmann-Gibbs form $`P(m)=Ce^{m/T}`$, where $`T`$ is the effective “money temperature” equal to the average amount of money per agent. This exponential distribution is indeed observed in computer simulations DY-money , as shown in Fig. 1.

It is interesting to compare this result with the actual PDF of money in the society. Unfortunately, it is very difficult to find the data on distribution of money $`m`$. On the other hand, a lot of statistical data is available for distribution of income $`r`$ (for revenue). Fig. 2 shows that the PDF of individual income in USA is very well fitted by the exponential function $`P(r)=Ce^{r/T}`$ DY-income .

The standard plot of PDF inevitably puts an upper limit on the horizontal axis (120 k$/year in Fig. 2). A standard way of representing the whole income distribution without any truncation is the so-called Lorenz curve shown in Fig. 3. The horizontal axis of the Lorenz curve, $`x(r)`$, represents the fraction of population with incomes below $`r`$, and the vertical axis $`y(r)`$ represents the fraction of the total income this population accounts for. As $`r`$ changes from 0 to $`\mathrm{}`$, $`x(r)`$ and $`y(r)`$ change from 0 to 1 and parametrically define the Lorenz curve in the $`(x,y)`$ space. The diagonal line $`y=x`$ represents the Lorenz curve in the case where all population has equal income. The inequality of the actual income distribution is characterized by the Gini coefficient $`0G1`$, which is the area between the diagonal and the Lorenz curve, normalized to the area of the triangle beneath the diagonal. For the exponential PDF, the Lorenz curve and the Gini coefficient can be easily calculated DY-income :

$$y=x+(1x)\mathrm{ln}(1x),G=1/2.$$
(1)
The solid line in Fig. 3 shows the theoretical Lorenz curve given by Eq. (1), and the points show the income data for 1979–1997. The agreement is quite good, in the first approximation, given that the curve (1) has no fitting parameters. The inset shows that the Gini coefficient is close to the theoretical value 1/2, although the inequality does increase during the last 20 years.

One may notice that discrepancy between the theory and the data occurs at the upper end of Fig. 3. The origin of this discrepancy becomes clear when we look at the cumulative distribution of income up to 1 M$/year shown in Fig. 4.

It is clear from Fig. 4 that income distribution for the great majority of population (more than 97%) is described by the exponential Boltzmann-Gibbs law. However, for a small fraction of population (less than 3%) with income above 100 k$/year, the PDF changes to the Pareto power law. The extra income in the upper tail of the distribution can be considered as a “Bose condensate”, and the Lorenz curve should be modified as DY-statistical

$$y=(1f)[x+(1x)\mathrm{ln}(1x)]+f\delta (1x),$$
(2)
where the last term is the delta-function, and $`f`$ is the fraction of income in the “Bose condensate”. As shown in Fig. 5, Eq. (2) gives an excellent fit of the data, and $`f=16\%`$ in 1997.

Thus far we discussed the distribution of individual income. By taking a convolution of two exponential distributions, it is easy to show that the PDF of family income is given by the modified exponential formula $`P(r)=Cre^{r/T}`$ DY-income . As Fig. 6 shows, this formula is in excellent agreement with the data. The corresponding Lorenz curve for family income is shown in Fig. 7 and compared with the data from the Bureau of Census for 1947–1994. It is amazing that the shape of the income distribution remains the same for half a century, and it is in agreement with the theoretical formula.

The theoretically calculated Gini coefficient for family income is 3/8=37.5% DY-income . The inset in Fig. 7 shows that the Gini coefficient in USA reported by the Bureau of Census for 1947–1994 is very close to the theoretical value. The average Gini coefficients for different regions of the World in 1988 and 1993 are shown in Fig. 8. For the well-developed market economies of West Europe and North America, the Gini coefficient is very close to the calculated value 37.5% and does not change in time. In other regions of the World, income inequality is higher. The special case is the former Soviet Union and East Europe, where inequality was lower before the fall of communism and has greatly increased afterwards.

In statistical physics, the exponential Boltzmann-Gibbs distribution is the equilibrium one, because it maximizes entropy. The data shown above demonstrate that probability distribution of income is also described by the Boltzmann-Gibbs law, and the equilibrium state of maximal entropy has been achieved in developed market economies.

## III Probability Distribution of Stock-Market Fluctuations

The first theory of stock-market fluctuations was proposed in 1900 in the Ph.D. thesis of the French mathematical physicist Louis Bachelier Taqqu . (Henri Poincar was on his Ph.D. committee.) His thesis developed the concept of Brownian motion (before the famous Einstein’s paper of 1905) for stock-market prices. A modern version of this theory is routinely used in financial literature. The theory predicts a Gaussian probability distribution for stock-price fluctuations. On the other hand, it is well known that the tails of the distribution are not Gaussian (the so-called “fat tails”). To improve agreement, it was proposed that the diffusion coefficient of the Brownian motion is not a constant, but itself is a stochastic variable. One popular model was proposed by Steve Heston, who is a faculty of the Department of Finance, University of Maryland.

In Ref. DY-Heston , Dragulescu and Yakovenko derived the probability distribution of price changes for the Heston model and compared it with the data. Fig. 9 shows the probability distribution of log-return $`x`$ for the Dow-Jones index during the 20-years period 1982–2001. (In the main panel of Fig. 9, the curves are offset vertically for clarity; the inset shows the same curves without offset.) The log-return $`x=\mathrm{ln}(S_2/S_1)`$ is the logarithm of the ratio of the stock prices $`S_2`$ and $`S_1`$ for two moments of time $`t_2`$ and $`t_1`$ with the average market growth subtracted. The probability distribution $`P_t(x)`$ depends on the time lag $`t=t_2t_1`$, which is indicated near each curve in Fig. 9. The solid curves show the analytically derived distribution, and the points show the Dow-Jones data. We see that the Dragulescu-Yakovenko formula DY-Heston  very well describes $`P_t(x)`$ for a broad range of time lags from one day to one year (252 trading days).

Dragulescu and Yakovenko also found that for times $`t`$ longer than the relaxation time of the model, $`P_t(x)`$ becomes a function of a single combination $`z`$ of the two variables $`x`$ and $`t`$. Thus, when plotted vs. $`z`$, the points for different time lags should collapse on a single scaling curve. That indeed happens, as show in Fig. 10. The solid line is the theoretically calculated scaling curve expressed in terms of a modified Bessel function. Notice that the agreement extends over seven orders of magnitude on the vertical axis.

In the recent paper SY-Heston , Silva and Yakovenko found that the same results hold for Nasdaq and S&P 500 in 1980s and 1990s. By analyzing the statistics of fluctuations, they concluded that the decline of stock market after 2000 is a long-term change of regime, not a temporary fluctuation, unlike the crash of 1987.

## IV Conclusions

We have demonstrated that methods and techniques of statistical physics can be successfully applied to economical and financial problems. The great experience of physicists in working with experimental data gives them a unique advantage to uncover quantitative laws in the statistical data available in economics and finance. The interdisciplinary field of econophysics is bringing new insights and new perspectives, which are likely to revolutionize the old social disciplines.


text:
# Unusual Subpulse Modulation in PSR B0320+39

## 1 Introduction

For a number of pulsars it is known that each pulse is composed of a number of subpulses, and that for each successive pulse the subpulses appear to “drift” by a given amount across the profile (Drake & Craft 1968). The rate at which subpulses drift across the profile is not constant but rather depends on their pulse longitude. This can be understood in the context of the widely adopted “carousel” model of drifting subpulses (Ruderman 1972), which postulates the existence of a ring of equally spaced “sparks” above the polar cap that give rise to “tubes” of plasma streaming upward from the surface, ducted along magnetic field lines. At a certain altitude these particles emit microwave radiation that is beamed along tangents to the local field lines, giving rise to a beam pattern that consists of a ring of “subbeams” (bright spots), reflecting the polar cap spark configuration. As the pulsar rotates the distant observer samples emission along a certain path in the beam pattern. Under the right viewing geometry, the sight-line makes a tangential pass along the ring of subbeams, giving rise to the reception of a sequence of one to a few subpulses every time the star rotates. Due to the non-linear mapping between (spin) longitude and magnetic azimuth, the subpulses should appear more closely spaced around the point at which the sight-line makes its closest approach to the magnetic pole. As the carousel slowly rotates about the magnetic pole, the longitude at which the emission from a given spark is seen exhibits a monotonic drift with time, which by extension of the above argument is greater in magnitude further from the magnetic pole.

By stacking the intensity time series modulo the pulse period to produce a two-dimensional array in pulse number and pulse longitude, and taking Fourier transforms along constant-longitude columns (forming a longitude-resolved fluctuation spectrum; LRFS), Backer (1970a, b) was able to measure not only the characteristic periodicity in pulse number ($`P_3`$, the time taken for one subpulse to drift to the former position of its neighbour), but also its amplitude and phase as a function of pulse longitude. Many subsequent studies employed this method but ignored the phase information, however those that examined the longitude-phase relation (Backer 1970a; Wright 1981; Davies et al. 1984; Biggs et al. 1987; Ashworth 1988) found behaviour qualitatively consistent with that expected under the carousel model.

The longitude-phase relationship expected under the carousel model depends on the number of sparks present, the misalignment angle between the spin and magnetic axes, and their orientation with respect to the line of sight. Hence studies of the relationship have the potential to provide valuable information about the pulsar, particularly in comparison to the observed longitudinal dependence of the polarisation position angle which, under the magnetic pole model of Radhakrishnan & Cooke (1969), depends on the same geometric parameters.

Motivated by these goals, we embarked on a program of analysis of several pulsars with drifting subpulses using the Westerbork Synthesis Radio Telescope (WSRT). One of the target sources was PSR B0320+39. This pulsar was discovered independently by Damashek et al. (1978) and Izvekova et al. (1982) in surveys of the northern sky at 400 and 102.5 MHz respectively. Further observations revealed the emission of very regular drifting subpulses (Izvekova et al. 1982) which occur in two distinct pulse longitude intervals, with a non-drifting component present in the average profile (Izvekoza et al. 1993; Kouwenhoven 2000). In this paper we report a most unexpected result from our study of the longitude-phase relationship of the drifting subpulses of PSR B0320+39.

## 2 Observations and Analysis

We observed PSR B0320+39 on 2000 November 3 with the Westerbork Synthesis Radio Telescope (WSRT). Signals from fourteen telescopes were added (after appropriate delays) in each sense of linear polarisation in a 10 MHz band centred at $`328`$ MHz. These were processed using the Dutch Pulsar Machine (PuMa), operating as a digital filterbank with 64 channels and $`819.2`$ $`\mu `$s sampling in total intensity (Stokes I). For details of the system see Voûte et al. (2002).

In offline analysis we de-dispersed the data and binned the resultant time series into an array in pulse longitude and pulse number using the ephemeris of Arzoumanian et al. (1994). A total of 3400 pulses in 3701 longitude bins was recorded, with a region of 512 bins including the on-pulse used in further analysis. The samples were normalised to produce an average profile with a zero baseline and a peak value of unity. The first 100 pulses are shown in a greyscale form in Fig. 1.

The techniques employed in the further analysis are based on those presented in detail by Edwards & Stappers (2002). We computed the two-dimensional Discrete Fourier Transform (DFT) of the first 2048 pulses of the pulse longitude – pulse number array to yield what we call the two-dimensional fluctuation spectrum (2DFS). This is mathematically equivalent to the harmonic resolved fluctuation spectrum (HRFS; Deshpande & Rankin 1999; Edwards & Stappers 2002)111Without the use of “twiddle factors” (Cooley & Tukey 1965), the 2DFS will differ from the HRFS. However, it is easily shown that this merely amounts to the provision of slightly different sampling lattices over the same continuous function, due to the assignment of integer pulse numbers (2DFS) versus integer harmonic numbers (HRFS)., and the form of the power spectrum found here for PSR B0320+39 matches that found by Kouwenhoven (2000) for the HRFS. For the purposes of clear display we also computed a spectrum using data from which we subtracted the average profile from each pulse, and padded with zeros in the longitude axis, in order to improve the resolution. The resolution in the time-associated axis was unnecessarily high (2048 elements), so to improve the signal-to-noise ratio we smoothed the power spectrum by convolution with a $`7\times 1`$ -bin ($`\nu _t\times \nu _l`$) boxcar before plotting the relevant portion of the spectrum222The sampling interval used allows longitude-associated frequencies up to 1850 cycles/period, however all significant power is present in the portion of the spectrum shown. in Fig. 2. The spectrum is dominated by a pair of components associated with the drift modulation; two are present because the spectrum derives from the Fourier transform of real-valued data (see Edwards & Stappers 2002). Also present with high significance are a pair of responses corresponding to the second harmonic of the subpulse response. The features near $`\nu _t=0.02`$ and $`0.04`$ cycles/period are significant, however they also appear in an analysis of off-pulse data, indicating that they are caused by periodic interference.

After performing the 2D DFT, the complex spectrum was shifted in both axes by an amount which placed the primary response to the drifting subpulses at zero frequency. We then masked other components present in the spectrum (namely the complex conjugate “mirror” of the primary response, the pair of second harmonics and the DC component) by multiplication of the coefficients with a transfer function:

$$f(\nu _t)=\{\begin{array}{cc}0\hfill & |\nu _t\nu _{tc}|/w<0.5\hfill \\ 2|\nu _t\nu _{tc}|/w1\hfill & 0.5|\nu _t\nu _{tc}|/w1\hfill \\ 1\hfill & |\nu _t\nu _{tc}|/w>1\hfill \end{array},$$
(1)
where $`\nu _t`$ is the time-associated axis, $`\nu _{tc}`$ is the centre frequency of the component to mask, and $`w`$ is the frequency extent of the region of zero transmission. The results that follow derive from filters with $`w=0.05`$ cycles/period, although as expected they were found to be insensitive to $`w`$ unless it becomes too small to mask adequately or so large that it masks some of the desired signal. We then computed the inverse DFT of the result.

We represent the subpulse signal as the real part of the product of a pure two-dimensional sinusoid and a two-dimensional complex modulation envelope. By the convolution theorem, this envelope is the result of the inverse transform performed above. It describes the deviations from a pure periodicity of constant amplitude that arise due to amplitude modulation at the pulse period, scintillation or nulling, non-uniform spacing of subpulses in pulse longitude (e.g. due to sight-line curvature) and slow variations in the drift rate. Using the iterative scheme of Edwards & Stappers (2002) we decomposed the two-dimensional envelope into the product of two modulation envelopes that vary in pulse longitude and pulse number respectively. Such a decomposition is able to model the behaviour of a rotating carousel, where the geometry defines the longitude-phase dependence and the carousel rotation determines the time-phase dependence, and its basic assumptions – that the longitudinal spacing of the subpulses is time invariant, and that the time spacing of the subpulses is longitude independent – are in agreement with previous studies; see Edwards & Stappers (2002). The validity of this approach in the case of PSR B0320+39 was confirmed by computation of the 2DFS of the difference between the observed signal and that predicted using the real part of product of the two one-dimensional envelopes. No significant power remained in the portion of the spectrum previously occupied by the primary drifting subpulse response.

The spectral shifting performed earlier has the effect of removing a constant phase slope in each of the envelopes. Since the slope varies over the envelope, the choice of value to use in the shifting is somewhat arbitrary and was refined after viewing the resultant envelopes. The periodicity with pulse number was well defined and a value of $`\widehat{P_3}=P_1/0.118`$ was used. The response to the drifting component in the other axis is very broad, and in fact extends well into negative frequencies (Fig. 2). This complicates the choice of value for $`1/\widehat{P_2}`$ to be used in shifting the spectrum in the longitude axis. We chose a value of $`\widehat{P_2}=P_1/60`$ to give the resultant envelope a relatively small phase slope over most of the pulse333Since this value does not lie on the natural grid of the DFT, the shift was actually performed in the time-longitude domain by multiplication with a complex exponential..

The inferred longitude-dependent envelope (Fig. 3) shows that nearly linear subpulse modulation occurs in two distinct longitude intervals, with a striking phase offset of $``$180° between them. In terms of the longitude-time diagram, the effect of the offset is to cause the drift band peaks in one component to extrapolate to troughs in the other and vice versa. To check this result we computed the (complex) LRFS and examined the amplitudes and phases in frequency bins near $`\widehat{P_3}=P_1/0.118`$. Since the modulation is not a pure periodicity, the power is spread over several bins in the spectrum. Although of a much lower signal-to-noise ratio, the shape of the phase envelope in each bin was consistent with that inferred with our preferred technique, which makes optimal use of all modulation power. Similarly, division of the observation into segments in time (which results in wider frequency bins and potentially also in reduced frequency drift in each segment) and examination of the phase relation in the peak bin of each segment also gave consistent results (regardless of the segment size). In order to improve the signal-to-noise ratio of these other methods, we also used an iterative algorithm based on multiplication with the complex conjugate of a template envelope to determine the correct phases with which to add envelopes from different frequency bins or time segments (producing in effect coarse-grained versions of the algorithm of Edwards & Stappers 2002). The resultant summed envelopes were again consistent with the result of our preferred algorithm (Fig. 3). As a final check, we folded the longitude-time data about the pulse number axis (modulo $`1/0.118`$ pulses) to produce an average profile as a function of drift phase and pulse longitude (Fig. 4). As with the LRFS methods, much of the modulation was lost due to the unstable drift motion, however some residual modulation remained and is sufficient to confirm the unusual phase behaviour.

Knowing that there is a $``$180° offset between the two halves of the subpulse phase envelope, we may explain two otherwise curious facts about PSR B0320+39. The first concerns the values of $`P_2`$ previously reported, which are around $`P_1/120`$ (Izvekoza et al. 1993 report $`3.1\mathrm{°}`$ of longitude, i.e. $`P_2P_1/116`$ at 406 MHz; Kouwenhoven 2000 reports $`24.5`$ ms, i.e. $`P_2123`$, at 328 MHz). These strongly differ from our nominal $`\widehat{P_2}=P_1/60`$ because they are based on measurements of actual subpulse separations, which will nearly always be taken between subpulses on opposite sides of the phase shift. Since the subpulse in the trailing component leads (or lags) the position it would take in a continuous linear phase envelope by $`180\mathrm{°}`$ of subpulse phase, the derived spacing differs from that expected from the typical phase slope by a factor of 1/2 (and probably sometimes 3/2). The second curiosity explained by the phase offset is the local minimum in the subpulse HRFS response at around 60 cycles/period (Fig. 2). The complex longitude-dependent modulation envelope (once the constant phase slope is subtracted) is quite close to an odd real-valued function, which causes its Fourier transform to have local minimum in its power at DC. This is the cause of the local minimum around $`60`$ cycles/period, the typical longitudinal drift frequency. Without the phase offset the envelope would be close to a real-valued even function and its Fourier transform would not have a strong minimum at its centre.

Furthermore, we note that the equivalence of the 2DFS and HRFS provides an explanation for the apparent “secondary” component at $`0.118`$ cycles/period, suggested by Kouwenhoven (2000) to be due to the ambiguity between which drift band joins to which between the two components. In fact, the subpulse response is very broad in the longitude-associated frequency axis, and extends significantly into negative frequencies. The symmetry of the two dimensional Fourier transform of real functions is such that the power at any frequency pair $`\nu _l,\nu _t`$ is equal to that at $`\nu _l,\nu _t`$ (where $`\nu _l`$ and $`\nu _t`$ are the frequencies in the longitude-associated and time-associated axes respectively). If only the positive frequencies in $`\nu _l`$ are plotted, as is the case with the HRFS, no information is lost but individual components that span the $`\nu _l=0`$ axis will appear as two spatially separated components at $`\nu _t`$ values of opposite sign. When the full spectrum is plotted as in Fig. 2 the proper interpretation becomes apparent.

## 3 Discussion & Interpretation

### 3.1 Polar Cap Models

The phenomenon of drifting subpulses has been used repeatedly in investigations of pulsar emission as a potentially strong clue about the geometry of the system. The fact that those pulsars that show nearly linear drifting subpulses tend to have a double, “unresolved double” or “single” average pulse profile morphologies, with many showing a transition from resolved to unresolved double as the observing frequency is increased, is usually seen as strong support for some form of “hollow-cone” polar cap model (Radhakrishnan & Cooke 1969; Ruderman & Sutherland 1975; Backer 1976; Rankin 1983, 1986; Lyne & Manchester 1988). PSR B0320+39 shares these characteristics of regular drifting subpulses and a transition from a resolved double to an unresolved double average profile morphology between 102 MHz (Izvekova et al. 1982)444Due to an improved system response the profile of Kuz’min & Losovskii (1999) shows much stronger component separation than that of Izvekova et al. (1982); see fig. 6. and 328 MHz (this work; also Damashek et al. 1978; Izvekoza et al. 1993).

The local steepening of the phase envelope of PSR B0320+39 in the transition region is inconsistent with the carousel model in its most basic form. For a single ring of “sparks”, as long as the sparks are elongated in either magnetic azimuth or opening angle (or neither, but not both), the subpulse modulation seen in a ring of constant opening angle has the same phase as that in any other ring, if measured from the same fiducial azimuth. Therefore, the observed longitude-dependent subpulse phase envelope is directly tied to to the magnetic azimuth sampling of the sight-line (via multiplication by the number of sparks, $`N`$) (Edwards & Stappers 2002). The only way to produce a sudden swing in magnetic azimuth (and hence subpulse phase) is to pass close to the magnetic pole. Such a configuration is inconsistent with the usual interpretation of unresolved double profile morphologies, and with the gradual linear polarisation position angle swing observed in this pulsar (Suleimanova & Pugachev 2002). To prove this assertion we attempted to fit the observed subpulse phase swing with the geometric formula of Edwards & Stappers (2002), using a $`\chi ^2`$ minimisation algorithm. From the figure provided by Suleimanova & Pugachev (2002) we conservatively estimate the range of polarization position angle gradients consistent with the data as $`\mathrm{d}\chi /\mathrm{d}\varphi =12\pm 3\mathrm{°}/\mathrm{°}`$ (where $`\chi `$ and $`\varphi `$ are the position angle and pulse longitude respectively). The parameters of the fit were required to be consistent with this gradient via the geometric model of Radhakrishnan & Cooke (1969). The result (Fig. 5) clearly fails to describe the observed behaviour.

One way to produce a jump in subpulse phase is to adopt a nested ring polar cap spark configuration (Gil & Sendyk 2000) and to assign the two halves of the profile to rings of different opening angles and numbers of sparks. However, this would result in a measurably different $`P_3`$ unless the circulation rates of the two carousels were different by a fortuitous ratio (i.e. the inverse of the ratio of the numbers of sparks in the rings). In addition, the transition from resolved to unresolved double profile with increasing frequency could not be explained by standard radius-to-frequency mapping (RFM; Cordes 1978). Some other mechanism such as the finite beam-width of the elementary emission mechanism ($`2\gamma ^1`$), the finite band-width of emission at a given radius (with attendant RFM in the effective range of radii visible), or refractive broadening would be required to produce the profile evolution.

### 3.2 Double Imaging of Magnetospheric Origin

#### 3.2.1 Basic Considerations

Given the failure of the standard model to account for the observed subpulse drift behaviour of PSR B0320+39, we develop in this section a pair of explanations based on a simple and plausible extension of the model. Namely, we argue that the most important features of the longitude-dependent subpulse modulation envelope – a local reduction of amplitude accompanied by a rapid swing in the phase angle – are the result of destructive interference between two superposed drifting subpulse signals of nearly opposite phase. The observed behaviour and this explanation of it can be understood by analogy with the phenomenon of orthogonal polarisation mode transitions, where the pulsar signal is seen to consist of the incoherent sum of two signals with nearly orthogonal polarisation states (e.g. Manchester et al. 1975; Stinebring et al. 1984). In both cases the total intensity of the two signals simply add in the average profile, while the near anti-parallel nature of the vectors representing the complex subpulse modulation signal (in the Argand plane) or the polarised part of the signal (in the Poincaré sphere) explains the reduction in magnitude (subpulse amplitude, polarised intensity) and the swing in orientation (subpulse phase, polarisation position angle) during a transition of dominance from one mode to the other. Other than the obvious features of the shift in phase and the attenuation in amplitude, the continuity of both the average profile and the absolute value of the gradient of the subpulse amplitude envelope serve as supporting evidence for this hypothesis.

Since the polar cap cannot be physically populated with overlapping spark systems, we argue that the superposed drifting components of PSR B0320+39 are associated with a form of “double imaging” of the polar cap excitation pattern. Whereas the polar cap may consist of a single, drifting ring of excitation points (e.g. sparks), the physics that map this to a radiated beam pattern must do so under two different transformations that are received in incoherent addition by the observer. This gives rise to a subpulse pattern that is the sum of two overlapping components with differing peak pulse longitudes, and phases that differ in their contribution at a given pulse longitude by a constant offset of $`180\mathrm{°}`$. This also implies that the average pulse profile is the simple sum of two overlapping components approximately equal to the corresponding subpulse amplitude envelopes.

In this section we discuss two potential mechanisms for such an effect. The first associates the images with different ray paths in a refractive magnetosphere (Petrova & Lyubarskii 2000), while the second attributes them to emission originating at two discrete altitudes in the magnetosphere (Rankin 1993). Before going into the details of each scenario, given that they both produce a beam consisting of a pair of nested cones, we are are in a position to examine the consistency between these models and the observed pulse profile morphology. This is important because, in assigning the two components of PSR B0320+39 to cones of different opening angles, we are abandoning the standard “single cone” explanation favoured for most pulsars with regular drifting subpulses. We suggest instead that a cut through an outer cone (from the outside to the inside) produces a leading component, while a trailing component is produced by a pass along the edge of the inner cone. At most radio frequencies the difference in the opening angles of the cones is small and the components overlap, but radius-to-frequency mapping causes the outer cone to expand at lower frequencies, giving the components sufficient separation to be resolved at 100 MHz. Clearly one should also expect a third component to be present on the trailing edge of the profile as the sight-line passes across the outer cone for the second time. This is actually observed at frequencies at and above 600 MHz (Gould & Lyne 1998), giving further support to our assertion that more than one cone of emission is visible from this pulsar. As shown in Fig.6, the pulse profiles observed for PSR B0320+39 match this picture well. We now turn to the specifics of the two mechanisms.

#### 3.2.2 Refraction

Refractive effects have long been considered potentially important in pulsar magnetospheres (e.g. Melrose 1979; Barnard & Arons 1986), particularly in the production of circular polarisation (Allen & Melrose 1982) and superposed orthogonally polarised radiation (McKinnon 1992). Only recently, however, has the potential for refraction to explain pulse profile morphologies seen detailed study (Lyubarskii & Petrova 1998; Petrova & Lyubarskii 2000; Petrova 2000). Based on the simple argument that the magnetospheric plasma density has minima both at the magnetic pole and at the edges of the cone defined by the last open magnetic field lines, Petrova & Lyubarskii (2000) showed that rays emitted poleward of the peak of the plasma distribution could be refracted across the pole to exit with an angle to the pole of the opposite sign to that with which they were emitted. Rays emitted outward of the plasma peak are simply refracted to somewhat greater opening angles. This is illustrated schematically in Fig. 7. Clearly, at any pulse longitude the observer may potentially receive two rays, one from either side of the magnetic pole.

In an axisymmetric plasma distribution, refraction occurs purely in the plane of the magnetic field lines. Petrova (2000) shows that the magnetosphere tends to concentrate rays from within a given set of co-planar field lines into two components of different peak opening angle, corresponding to the different ray paths described above. These are associated with inner and outer components of the conal beam, which in terms of explaining the morphology of pulsar average profiles serves the same purpose as the double cone of Rankin (1993) (see below). This grouping of rays occurs within a given plane of field lines; all possible origins for visible rays must be co-planar, along with the line of sight and the magnetic axis. Therefore, to the extent that the plasma distribution is axisymmetric, the geometric conveniences of the basic magnetic pole model (Radhakrishnan & Cooke 1969) are preserved, including the polarisation behaviour, and the mapping of azimuthal intensity modulations (e.g. from a plasma flow initiated by a polar cap spark carousel) to observable drifting subpulses.

Moreover, if the polar cap pattern consists of a ring of an odd number of points of excitation, the subpulse pattern seen in the component associated with pole-crossing rays should be close to $`180\mathrm{°}`$ out of phase with that in the other component. The slight deviation from purely out of phase signals required to explain the finite longitude interval over which the phase transition occurs might be explained as a retardation effect or a refractive delay. Alternatively, inclusion of magnetospheric rotation effects or refraction under a moderately non-axissymmetric plasma distribution (Petrova & Lyubarskii 2000 model both but only in the context of polarisation wave mode coupling; see also Blaskiewicz et al. 1991) could potentially provide the kind of symmetry-breaking needed to produce subpulse phase offsets other than 0 or 180°. The form of the beam pattern and the manner in which a rotating carousel polar cap system maps to drifting subbeams under such conditions is yet to be investigated, to our knowledge.

We note that this model can also explain a recent result concerning the drifting subpulses of PSR B0809+74. Ramachandran et al. (2002) found that over the course of each subpulse, the polarisation state jumps in a regular manner from one polarisation state to another, almost orthogonal state. This is similar to the behaviour seen in PSR B0320+39, in that the observed total intensity drifting subpulse pattern can be understood as the sum of two components with a phase offset between them. However, in PSR B0320+39 the drift components cannot be associated with orthogonal polarisation modes since the linearly polarised average profile at 408 MHz peaks in intensity at the centre of the profile (Gould & Lyne 1998), where maximal depolarisation ought to take place, were the components orthogonally polarised. The single-pulse polarimetric work of Suleimanova & Pugachev (2002) confirms that the situation is similar at 103 MHz: significant power is present in two modes only in the profile wings. The model of Petrova (2001) for the origin of orthogonal polarisation modes offers a simple explanation for this disparity. Under that model the secondary polarisation mode is produced via the conversion of power from the ordinary mode in regions of propagation that are quasi-longitudinal to the magnetic field. If radiation in one of the conal components experiences significant conversion while the other does not, as could well occur since the ray trajectories involved are very different between the two, the situation is as seen in PSR B0809+74. If the configuration is not conducive to conversion then both components have similar polarisation, as in PSR B0320+39. Since the conversion process also alters the orientation (i.e. position angle, ellipticity) of the modes, the longitude-dependent non-orthogonality of the polarisation modes reported by Ramachandran et al. (2002) might also receive an explanation via this effect.

#### 3.2.3 Dual Emission Altitudes

A second potential explanation for how we receive the sum of two related subpulse signals is that the two signals arise at different altitudes in the magnetosphere. This is the essence of the “double conal” model of Rankin (1993), so-called because the production of radiation cones at two distinct altitudes in the magnetosphere gives rise to a beam pattern that consists of two nested cones. The radial separation of the emitting regions and the differential tangent angle between the active field lines and the magnetic axis that arises as a result causes the polar cap plasma flux cone to produce a wider radiation beam at the upper emission region, compared to the lower one. If the height difference is large, up to four well-resolved conal profile components are produced, while reduced vertical separation causes the components to overlap and produce a “boxy” profile shape. Fig. 8 illustrates this model as it might apply for PSR B0320+39 at 102 MHz, with the conal beams overlapping to a lesser degree than at higher frequencies.

Clearly this model is able to explain how we see superposed drifting patterns, but we have not yet offered an explanation for why there should be a subpulse phase offset between them. They arise from the same pattern of sparks (or plasma flux production) on the polar cap, yet we observe a phase difference of $`180\mathrm{°}`$ near the profile centre. For the subpulse pattern to be delayed by $`P_2/230`$ ms due to simple retardation (Cordes 1978) requires that the emission regions are separated by $`9000`$ km, 1–2 orders of magnitude greater than the values typically inferred from profile morphologies (e.g. Gangadhara & Gupta 2001; Mitra & Rankin 2002; Gupta & Gangadhara 2003) . Including differential aberration reduces this by no more than a factor of 2 (Cordes 1978). Hence we consider this explanation unlikely. An alternative explanation worth pursuing in the future, but outside the scope of this work, is that a significant electric field exists throughout the magnetosphere, causing particles to accelerate upward mostly along magnetic field lines, but with some azimuthal $`E\times B`$ drift (e.g. Wright 2003). The drift would cause the beam pattern radiated at a given height to be rotated in magnetic azimuth by a certain amount, which may account for the subpulse phase offset seen here.

In any case, under this model the phase difference must ultimately be due to the different altitudes at which the components originate. Since RFM must apply to the outer (i.e. upper) component to explain the 102 MHz pulse profile, the difference in subpulse phase between the components ought to be a strong function of frequency. Although this is the first study to examine the longitude dependence of subpulse phase in this pulsar, the fact that the longitude-dependent subpulse amplitude envelopes measured at 102 MHz and 408 MHz by Izvekoza et al. (1993) drop to close to zero in the profile centre, just as at 328 MHz, suggests that at all three frequencies the subpulses add destructively in this interval, and therefore that the phase difference between them is close to $`180\mathrm{°}`$, regardless of frequency. The origin of the phase offset therefore remains a mystery under this scenario.

## 4 Conclusions

The drifting subpulses of PSR B0320+39 occur in two distinct intervals of pulse longitude. Using a sensitive new technique we have shown that the extrapolation of subpulse drift from a given longitude interval is $`180\mathrm{°}`$ out of phase with the subpulses of the other. This was confirmed by analysis of the longitude-resolved fluctuation spectrum and by folding the data into profiles corresponding to different phases of the subpulse drift.

This suggests the presence of two components of nearly opposite drift phase, and based on the forms of the average pulse profile and the longitudinal dependence of the amplitude of subpulses, we argued that the modes are superposed at all longitudes. The adding of intensity contributions from the two modes results in attenuation of the modulations in the longitude interval where the modes are of similar intensity, and produces a steep transition of $`180\mathrm{°}`$ in the longitudinal subpulse phase envelope at the point where the component intensities are equal.

Two specific models are suggested to explain this behaviour. Both produce a nested pair of conal beams, one via magnetospheric refraction, the other via emission at two distinct altitudes. The leading and trailing components are associated with outer and inner cones respectively. The second cut of the outer cone does not give rise to a component at $`102`$ MHz or $`328`$ MHz, but does produce a third, trailing component at frequencies above $`600`$ MHz. The component separation increases rapidly at frequencies below $`328`$ MHz, in accordance with the expectations of radius to frequency mapping. The refractive scenario can explain the observed phase offset as a simple consequence of having an odd number of subbeams in the carousel, whilst under the dual-altitude emission model the phase offset is more difficult to understand. Planned observations of the time and frequency dependence of the phase offset, and the behaviour of subpulses in the third component present at and above $`600`$ MHz, should prove most helpful in the further evaluation of these models.

###### Acknowledgements.

The authors wish to thank P. Weltevrede for enlightening discussions, and the referee whose persistent skepticism led to the correction of an error in our analysis. RTE is supported by a NOVA fellowship. BWS is supported by NWO Spinoza grant 08-0 to E. P. J. van den Heuvel. The Westerbork Synthesis Radio Telescope is administered by ASTRON with support from the Netherlands Foundation For Radio Astronomy. Part of this research has made use of the data base of published pulse profiles maintained by the European Pulsar Network, available at http://www.mpifr-bonn.mpg.de/pulsar/data/.


text:
# The extended empirical process test for non-Gaussianity in the CMB, with an application to non-Gaussian inflationary models

## I Introduction

Much information is expected from ongoing and future experiments aimed at measuring the CMB radiation. In this framework, a considerable amount of attention has been drawn in the recent literature by Gaussianity tests on CMB data. Such tests have been motivated by the prediction of the so-called theory of inflation, currently the most popular theory of the very early universe. Standard inflation predicts that the temperature fluctuations of the CMB, should be (very close to) Gaussian distributed (see the reviews Narlikar and Padbanabahn (1991); Watson ). On the other hand more recent, non-standard models of inflation predict small non-Gaussianities; detecting such effects would be important for understanding the physics of the primordial epochs Martin et al. (2000); Contaldi et al. (2000); Linde and Mukhanov (1997); Bartolo et al. (2002); Gupta et al. (2002); Gangui et al. (2002). More radical deviations from Gaussianity are predicted by models where topological defects are induced by phase transitions (See Ref.Gangui et al. (2001) for a review). Detecting non-Gaussianities is also relevant as a tool for controlling systematic effects in CMB data. Among the many procedures advocated in the literature, many are based upon topological properties of spherical Gaussian fields Novikov et al. (2000); Gott et al. (1990); Eriksen et al. (2002); Heavens and Gupta (2000); Scaramella and Vittorio (1991); Doré and S. Colombi ; Barreiro et al. (2001); others are based on harmonic space properties Phillips and Kogut (2000); Komatsu et al. (2002); Winitzki and Wu ; Hu (2001); Kunz et al. (2001); Barreiro and Hobson (2001) For empirical results, see Wu et al. (2001); Polenta et al. (2002); Ferreira et al. (1998); Banday et al. (2000).

In a recent paper Hansen et al. (2002) we proposed a new procedure to detect non-Gaussianity in harmonic space using so-called empirical processes Shorack and Wellner (1986); Marinucci and Piccioni (2002). More precisely, let $`T(\theta ,\phi )`$ denote the CMB fluctuations field, which we assume, as always, to be homogeneous and isotropic, for $`0<\theta \pi ,`$ $`0<\phi 2\pi `$. Assuming that $`T(\theta ,\phi )`$ has zero mean and finite second moments, it is well-known that the following spectral representation holds:

$$T(\theta ,\phi )=\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}mY_\mathrm{}m(\theta ,\phi )\text{ ,}$$

where $`Y_\mathrm{}m(\theta ,\phi )`$ denotes the spherical harmonics. The random coefficients (amplitudes) $`\left\{a_\mathrm{}m\right\}`$ have zero-mean with variance $`|a_\mathrm{}m|^2=C_{\mathrm{}}`$. They are uncorrelated over $`\mathrm{}`$ and $`|m|`$: $`a_\mathrm{}ma_\mathrm{}^{}m^{}^{}=C_{\mathrm{}}\delta _{\mathrm{}}^{\mathrm{}^{}}\delta _m^m^{}`$, and $`a_{\mathrm{},m}=a_{\mathrm{},m}^{}`$. The sequence $`\left\{C_{\mathrm{}}\right\}`$ denotes the angular power spectrum of the random field and the asterisk complex conjugation. Furthermore, if $`T(\theta ,\phi )`$ is Gaussian, the $`\left\{a_\mathrm{}m\right\}`$ have a complex Gaussian distribution. Upon observing $`T(\theta ,\phi )`$ on the full sky, the random coefficients can be obtained through the inversion formula

$$a_\mathrm{}m=_\pi ^\pi _0^\pi T(\theta ,\phi )Y_\mathrm{}m^{}(\theta ,\phi )\mathrm{sin}\theta d\theta d\phi \text{ , }m=0,\pm 1,\mathrm{},\pm l\text{ },\text{ }l=1,2,\mathrm{}\text{ .}$$
(1)
Our idea in that paper was to study the empirical distribution function for the $`\left\{a_\mathrm{}m\right\},`$ and to use these results to implement tests for non-Gaussianity in harmonic space, in the presence of an unknown angular power spectrum. Although the Monte Carlo evidence produced in Hansen et al. (2002) was rather encouraging, several important issues were left for further research. Indeed, it must be noted that although the joint distribution of the spherical harmonic coefficients is invariant under rotations (under Gaussianity), their numerical value is not; as a result, the sample outcome of the test in Hansen et al. (2002) is not rotationally invariant. In this paper this remark is exploited to implement an extended, more powerful version of the empirical process procedure, based upon the combination of results from several different choices of coordinate systems. Furthermore, empirical process ideas are also implemented in four different setups: whereas in Hansen et al. (2002) we focused only on the squared spherical harmonic coefficients, here we consider also the real and imaginary parts as separate random variables, which allows us to gain further efficiency improvements. More important, we consider the empirical process on the random phases of the Gaussian field. This approach is especially convenient, because under Gaussianity not only the random phases are model independent (compare Chiang and Naselsky ) but the Smirnov uniformization (see below) is exact, i.e. no bias arises. Moreover the joint information on the random phases and the random amplitudes uniquely identifies the distribution of the field, under Gaussianity. We address also some experimentally meaningful situations, such as the presence of instrumental noise and a galactic cut in the map; concerning the latter, we advocate a solution which may have some independent interest for other harmonic space methods of statistical inference. Our tests are implemented on models of non-Gaussianity based upon ongoing research on non-standard inflationary models.

The plan of the paper is as follows: in section (II) we review the empirical process method; the extensions provided in this paper are presented in (III). Section (IV) presents the results from Monte-Carlo simulations on 100 realizations of a non-standard inflationary model, with varying levels of non-Gaussianity. Comments and directions for further research are collected in the final section (V).

## II The empirical process method

The details of the empirical process approach to detect non-Gaussianity in the CMB were given in Hansen et al. (2002). In short, the method consists of a family of tests which focus on the total distribution of $`a_\mathrm{}m`$ and checks for dependencies between $`k`$ $`\mathrm{}`$-rows. The first step is to transform the spherical harmonic coefficients into variables $`u_\mathrm{}m`$ which have an approximate uniform distribution between $`0`$ and $`1`$, given that the $`a_\mathrm{}m`$ were initially Gaussian distributed. This is done using the Smirnov transformation, defined as

$$u_\mathrm{}0=\mathrm{\Phi }_1\left(\frac{|a_\mathrm{}0|^2}{\widehat{C}_{\mathrm{}}}\right)\text{ , }u_\mathrm{}m=\mathrm{\Phi }_2\left(\frac{2|a_\mathrm{}m|^2}{\widehat{C}_{\mathrm{}}}\right)\text{ , }m=1,2,\mathrm{},l,\text{ }l=1,2,\mathrm{}L\text{ ,}$$

where $`\mathrm{\Phi }_n`$ is the cumulative distribution function of a $`\chi ^2`$ with $`n`$ degrees of freedom and $`\widehat{C}_{\mathrm{}}`$ are the power spectrum coefficients estimated from the data. The error introduced by using estimated $`\widehat{C}_{\mathrm{}}`$ instead of the real underlying $`C_{\mathrm{}}`$ is dealt with using a bias-subtraction, as described in Hansen et al. (2002).

Then the joint empirical distribution function for row $`\mathrm{}`$ is formed,

$$\widehat{F}_{\mathrm{}\mathrm{}\mathrm{}+\mathrm{\Delta }_{\mathrm{},k1}}(\alpha _1,\mathrm{},\alpha _k)=\frac{1}{(\mathrm{}+1)}\underset{m=0}{\overset{\mathrm{}}{}}\left\{\underset{¯}{\mathrm{𝟏}}(\widehat{u}_\mathrm{}m\alpha _1)\underset{i=2}{\overset{k}{}}\underset{¯}{\mathrm{𝟏}}(\widehat{u}_{\mathrm{}+\mathrm{\Delta }_{\mathrm{},i1},m+\mathrm{\Delta }_{mi}}\alpha _i)\right\}\text{ , }\mathrm{\Delta }_{mi}0\text{ ,}$$

where $`\mathrm{\Delta }_{\mathrm{},i}`$ determines the spacing between the rows for which the dependencies are tested and $`\mathrm{\Delta }_{mi}`$ denotes the difference in $`m`$ for row $`i`$. The parameters $`\alpha _i`$ run over the interval $`[0,1]`$. The empirical process is expressed using the centered and rescaled $`\widehat{F}_{\mathrm{}\mathrm{}\mathrm{}+\mathrm{\Delta }_{\mathrm{},k1}}`$ given as

$$\widehat{G}_{\mathrm{}\mathrm{}\mathrm{}+\mathrm{\Delta }_{\mathrm{},k1}}(\alpha _1,\mathrm{},\alpha _k)=\sqrt{(\mathrm{}+1)}\left\{\widehat{F}_{\mathrm{}\mathrm{}\mathrm{}+\mathrm{\Delta }_{\mathrm{},k1}}(\alpha _1,\mathrm{},\alpha _k)\underset{i=1}{\overset{k}{}}\alpha _i\right\}\text{ .}$$

The intuition behind this procedure is as follows: if the $`a_\mathrm{}m`$s are Gaussian, $`\widehat{G}`$ converges to a well-defined limiting process, whose distribution can be readily tabulated. On the other hand, for non-Gaussian $`a_\mathrm{}m`$s $`\{\widehat{F}(\alpha _1,\mathrm{},\alpha _k)_{i=1}^k\alpha _i\}`$ and thereby $`\widehat{G}`$ will take ‘high’ values over some parts of $`\alpha `$-space. Thus, the analysis of some appropriate functional of $`\widehat{G}`$ can be used to detect non-Gaussianity. To combine the information over all multipoles into one statistic, we define

$$\widehat{K}_L(\alpha _1,\mathrm{},\alpha _k,r)=\frac{1}{\sqrt{L\mathrm{\Delta }_{\mathrm{},k1}}}\underset{\mathrm{}=1}{\overset{[(L\mathrm{\Delta }_{k1})r]}{}}\widehat{G}_{\mathrm{},\mathrm{},\mathrm{}+\mathrm{\Delta }_{\mathrm{},k1}}(\alpha _1,\mathrm{},\alpha _k)\text{ ,}$$

where $`L`$ is the highest multipole where the data is signal dominated.

The method can then be summarized as follows: the distribution of $`\mathrm{sup}|\widehat{K}_L|`$ is found using Monte-Carlo simulations of Gaussian distributed $`a_\mathrm{}m`$. Then, for a given observed set of $`a_\mathrm{}m`$, the value $`k_{\mathrm{max}}=\mathrm{sup}|\widehat{K}_L|`$ is found and compared to the distribution obtained from Monte-Carlo. The consistency of the data with a Gaussian distribution can then be estimated to any suitable $`\sigma `$-level.

## III Improvements of the empirical process method

The process described in the previous section can be improved in several ways. One important thing to notice is that if the $`a_\mathrm{}m`$ are non-Gaussian distributed, their joint distribution will change under a rotation of the sky, and therefore it will not be invariant with respect to the choice of coordinate system. In spherical harmonic space, the rotation of the sky with the Euler angles $`(\mathrm{\Phi }_1,\mathrm{\Theta },\mathrm{\Phi }_2)`$ can be written as,

$$a_\mathrm{}m^R=\underset{m^{}}{}D_{mm^{}}^{\mathrm{}}(\mathrm{\Phi }_1,\mathrm{\Theta },\mathrm{\Phi }_2)a_\mathrm{}m^{},$$
(2)
where the rotation matrices $`D_{mm^{}}^{\mathrm{}}(\mathrm{\Phi }_1,\mathrm{\Theta },\mathrm{\Phi }_2)`$ are described in Appendix (A). This new sets of $`a_\mathrm{}m`$ will yield a different (although not independent) value of $`k_{\mathrm{max}}`$; averaging the $`k_{\mathrm{max}}`$ obtained from the empirical process over several rotations may thus increase the probability for detecting non-Gaussianity in the data. In Table 1, we indeed prove that this effect is noticeable, considering a non-Gaussian inflationary model with $`f_{NL}=300`$. To obtain this data, we used 100 Gaussian Monte-Carlo simulations to find the distribution of $`k_{\mathrm{max}}`$. Then we used 100 realizations of the non-Gaussian model and for each realization we averaged the $`k_{\mathrm{max}}`$ value obtained over 45 rotation in the $`\mathrm{\Theta }`$-direction. The results are shown using the process to test Gaussianity and dependency between 1, 2 and 3 $`\mathrm{}`$-rows. We define the $`1\sigma `$, $`2\sigma `$ and $`3\sigma `$ detection levels as the $`k_{\mathrm{max}}`$ values over which we had $`32\%`$, $`5\%`$ and $`1\%`$ of the hits in the Gaussian simulations respectively.

As seen from equation (2), a rotation in the $`\mathrm{\Phi }`$-direction will only change the phase of the $`a_\mathrm{}m`$ for which the empirical process as defined above is insensitive. The reason is that we define the $`u_\mathrm{}m`$ using the modulus of the $`a_\mathrm{}m`$. But by defining a complex $`u_\mathrm{}m`$ for which the real and imaginary parts are just the uniformized real and imaginary parts of the $`a_\mathrm{}m`$, the phase information can be kept, and the test will be sensitive also to rotations in the $`\mathrm{\Phi }`$ direction. The new $`u_\mathrm{}m`$ can be obtained by a similar Smirnov transformation,

$$u_\mathrm{}m^r=G\left(\frac{a_\mathrm{}m^r}{\sqrt{\widehat{C}_{\mathrm{}}}}\right)\text{ , }u_\mathrm{}m^i=G\left(\frac{a_\mathrm{}m^i}{\sqrt{\widehat{C}_{\mathrm{}}}}\right)\text{ , }m=1,2,\mathrm{},l,\text{ }l=1,2,\mathrm{}L\text{ ,}$$

where $`G`$ is the cumulative distribution function of a Gaussian variable with zero mean and unit variance and the $`r`$ and $`i`$ indices are used to indicate real and imaginary parts.

With these new $`u_\mathrm{}m`$ the standard empirical process can be applied, now also allowing the check for dependencies between real and imaginary parts. We will in this paper arrange the $`u_\mathrm{}m`$ in two ways. First we will double the number of rows, putting the imaginary rows in between the real rows, then we apply the standard empirical process using $`\mathrm{\Delta }\mathrm{}=\mathrm{\Delta }m=250`$ for $`L=500`$ as in Hansen et al. (2002) (from now on we will label this the C test). Otherwise we will just double the length of the $`\mathrm{}`$-rows, putting the imaginary parts on the negative $`m`$ side. Then we apply the standard empirical process still with $`\mathrm{\Delta }\mathrm{}=\mathrm{\Delta }m=250`$ for $`L=500`$ (called the B test). In both ways we check dependencies between real and imaginary parts. In Table 2 we show results from these two new tests applied to 100 realizations of the $`f_{NL}=300`$ model used above (to be compared with the results of test A without rotations in Table 1). Even if the number of detections in each test vary, the specific realizations detected by each of the tests also vary. This suggests that by combining the results from all tests, the number of detections can be increased. We suggest to make the average $`k_{\mathrm{max}}`$ over all 3 tests as a new statistic. The average has to be weighted with the mean $`\overline{k}_{\mathrm{max}}`$ for each test, obtained from Monte-Carlo. Thus the new statistic would be:

$$k_{\mathrm{max}}=k_{\mathrm{max}}^A+\frac{\overline{k}_{\mathrm{max}}^A}{\overline{k}_{\mathrm{max}}^B}k_{\mathrm{max}}^B+\frac{\overline{k}_{\mathrm{max}}^A}{\overline{k}_{\mathrm{max}}^C}k_{\mathrm{max}}^C,$$

where the letters indicate test A, B and C.

As we showed above, the old test can be improved significantly using rotations in the $`\mathrm{\Theta }`$-directions. The two new tests will also be sensitive to rotations in the $`\mathrm{\Phi }`$-directions, as these rotations only change the phases of the $`a_\mathrm{}m`$. In Table 3 we show a similar increase in the number of detections using 30 rotations in the $`\mathrm{\Phi }`$-rotation. We also show the total number of detections, combining test A averaged over 30 $`\mathrm{\Theta }`$-rotations and test B and C averaged over 30 $`\mathrm{\Phi }`$-rotations. Clearly the highest number of detections can be obtained combining all 3 tests with as many rotations as possible for each, and in this way testing all possible dependencies in spherical harmonic space.

Finally the same sort of ideas can be implemented by working directly on the random phases by taking the uniformized coefficients to be $`u_\mathrm{}m=\mathrm{arctan}(a_\mathrm{}m^i/a_\mathrm{}m^r)/\pi `$. From a theoretical point of view, this choice is particularly appealing because the uniformization is exact. Random phases as a test for non-Gaussianity were considered also by Chiang and Naselsky ; these authors however do not consider the empirical process approach. We label this test D and we report the related results in Table 7.

In practice, a relevant issue is the coupling between spherical harmonic coefficients introduced by partial sky coverage from a galactic cut. It turns out that this can be well controlled if the Monte-Carlo calibration of the $`k_{\mathrm{max}}`$ distribution are done with the same galactic cut. The introduction of the galactic cut has the effect of reducing the number of detections, but this can be significantly improved if the gap is refilled, using the two borders between the gap and the observed area as a mirror. This creates a discontinuity at the galactic equator, but improves the power of the procedure. In our test, using this refilling technique, we found no noticeable reduction in the number of detections with and without a galactic cut of $`5^{}\times 2`$ in the 100 realizations used in the simulations.

## IV Detection of non-Gaussian inflationary models

In this section we will show the results of the 4 versions of the empirical process method applied to 100 realizations of a non-standard inflationary model. Indeed, as discussed in Bartolo et al. (2002), one way of parameterizing the possible presence of non-Gaussianity in the primordial gravitational potential $`\mathrm{\Phi }`$ is to represent it in the following way,

$$\mathrm{\Phi }=\varphi +f_{NL}(\varphi ^2<\varphi ^2>),$$
(3)
where $`\varphi `$ is a zero mean Gaussian random field and $`f_{NL}`$ is a non-linear parameter which can be observationally constrained. Our empirical examples are based on this model, using for simplicity a pure Sachs-Wolfe angular power spectrum ($`C_{\mathrm{}}\mathrm{}(\mathrm{}+1)`$). It should be noted that because $`\varphi `$ is of the order $`10^5`$, the quadratic term is five orders of magnitudes smaller than the linear one.

We will use 3 levels of non-Gaussianity, $`f_{NL}=100`$, $`f_{NL}=300`$ and $`f_{NL}=1000`$ (a map for each $`f_{NL}`$ with its corresponding histogram is shown in figures 1-6). On these realizations we put non-uniform noise corresponding to the noise-level as expected from the Planck LFI $`100\mathrm{G}\mathrm{H}\mathrm{z}`$ channel. We also include a $`5^{}\times 2`$ galactic cut. The Monte-Carlo for determination of the distribution of $`k_{\mathrm{max}}`$ will when not otherwise indicated be made using the same maps with $`f_{NL}=0`$.

The results of the tests are reported in Tables 4-7. We see that the 3-rows version of the test provides the best performance. Among the 4 tests, none appears to detect significantly more than others. Down to a level $`f_{NL}=300`$, more than half of the maps are detected at the $`2\sigma `$ level. However at $`f_{NL}=100`$, the detection level at $`2\sigma `$ has dropped to about $`30\%`$ at the best. Nevertheless, as the spectrum we used went only to multipole $`\mathrm{}=500`$ and that radiative transfer was not taken into account, we expect that these numbers can improve significantly.

## V Comments and conclusions

In Hansen et al. (2002) we argued that the empirical process method enjoyed a number of advantages over existing methods, such as being completely model independent, allowing for a rigorous statistical analysis, detecting the location of non-Gaussianity in harmonic space, being easily amendable to extensions and generalizations. In this paper the potentials of this approach are further investigated. Three further classes of tests in the same spirit as in the previous paper are proposed. Changes of coordinates are implemented to make the procedure rotationally invariant and more powerful against non-Gaussian alternatives. We have also considered experimentally meaningful situations such as the presence of gaps and noise. Finally, we have tested the power of these procedures against physically reasonable models such as those expected by non-standard inflation; we show that for the model considered, significant detections of non-Gaussianities are achieved for non linear coupling parameters down to about 100. Overall we believe that empirical process methods emerge as valuable tools for further studies on the non-Gaussianity of the CMB.

###### Acknowledgements.

We acknowledge the use of the Healpix package<sup>1</sup><sup>1</sup>1http://www.eso.org/science/healpix/. We are very grateful to Sabino Matarrese and Michele Liguori for supplying us with non-Gaussian maps. We are also very grateful to Amadeo Balbi, Paolo Natoli and Paolo Cabella for very useful discussions. FKH acknowledges financial support from the CMBNET Research Training Network.

## Appendix A Rotation Matrices

A spherical function $`T(\theta ,\varphi )`$ is rotated by the operator $`\widehat{D}(\alpha \beta \gamma )`$ where $`\alpha \beta \gamma `$ are the three Euler angles for rotations Risbo (1996) and the inverse rotation is $`\widehat{D}(\gamma \beta \alpha )`$. For the spherical harmonic functions, this operator takes the form,

$$Y_\mathrm{}m(\theta ^{},\varphi ^{})=\underset{m^{}=\mathrm{}}{\overset{\mathrm{}}{}}D_{m^{}m}^{\mathrm{}}(\alpha \beta \gamma )Y_\mathrm{}m^{}(\theta ,\varphi ),$$
(4)
where $`D_{m^{}m}^{\mathrm{}}`$ has the form

$$D_{m^{}m}^{\mathrm{}}(\alpha \beta \gamma )=e^{im^{}\alpha }d_{m^{}m}^{\mathrm{}}(\beta )e^{im\gamma }.$$
(5)
Here $`d_{m^{}m}^{\mathrm{}}(\beta )`$ is a real coefficient with the following property:

$$d_{m^{}m}^{\mathrm{}}(\beta )=d_{mm^{}}^{\mathrm{}}(\beta ).$$
(6)
The D-functions also have the following property:

$$D_{m^{}m}^{\mathrm{}}(\alpha \beta \gamma )=\underset{m^{\prime \prime }}{}D_{m^{}m^{\prime \prime }}^{\mathrm{}}(\alpha _2\beta _2\gamma _2)D_{m^{\prime \prime }m}^{\mathrm{}}(\alpha _1\beta _1\gamma _1),$$
(7)
where $`(\alpha \beta \gamma )`$ is the result of the two consecutive rotations $`(\alpha _1\beta _1\gamma _1)`$ and $`(\alpha _2\beta _2\gamma _2)`$.

The complex conjugate of the rotation matrices can be written as

$$D_{mm^{}}^{\mathrm{}}=(1)^{m+m^{}}D_{(m)(m^{})}^{\mathrm{}}.$$
(8)