text:
# New empirical fits to the proton electromagnetic form factors

## Abstract

Recent measurements of the ratio of the elastic electromagnetic form factors of the proton, GEp/GMp, using the polarization transfer technique at Jefferson Lab show that this ratio decreases dramatically with increasing Q2, in contradiction to previous measurements using the Rosenbluth separation technique. Using this new high quality data as a constraint, we have reanalyzed most of the world $`ep`$ elastic cross section data. In this paper, we present a new empirical fit to the reanalyzed data for the proton elastic magnetic form factor in the region 0 $`<`$ $`Q^2`$ $`<`$ 30 GeV2. As well, we present an empirical fit to the proton electromagnetic form factor ratio, GEp/GMp, which is valid in the region 0.1 $`<`$ $`Q^2`$ $`<`$ 6 GeV2.

The elastic electromagnetic form factors are crucial to our understanding of the proton’s internal structure. Indeed, the differential cross section for elastic $`epep`$ scattering is described completely in terms of the Dirac and Pauli form factors, $`F_1`$ and $`F_2`$, respectively, based solely on fundamental symmetry arguments. Further, the Sachs form factors, $`G_{Ep}`$ and $`G_{Mp}`$, which are simply derived from $`F_1`$ and $`F_2`$, reflect the distributions of charge and magnetization current within the proton.

Until recently, the form factors of the proton have been determined experimentally using the Rosenbluth separation method , in which one measures elastic $`ep`$ cross sections at constant $`Q^2`$, and varies both the beam energy and scattering angle to separate the electric and magnetic contributions. In terms of the Sachs form factors, the differential cross section for elastic $`ep`$ scattering has traditionally been written as

$$\frac{d\sigma }{d\mathrm{\Omega }}=\sigma _{ns}\left(\frac{G_{Ep}^2+\tau G_{Mp}^2}{1+\tau }+2\tau G_{Mp}^2\mathrm{tan}^2(\theta _e/2)\right),$$
(1)
where $`\tau =Q^2/4M_p^2`$, $`\theta _e`$ is the in-plane electron scattering angle. For elastic $`ep`$ scattering, the so-called nonstructure cross section, $`\sigma _{ns}`$ is given by

$$\sigma _{ns}=\frac{\alpha _{em}^2\mathrm{cos}^2(\theta _e/2)}{4E^2\mathrm{sin}^4(\theta _e/2)}\frac{E^{}}{E}$$
(2)
where $`\alpha _{em}`$ is the electromagnetic coupling constant, and $`E^{}`$($`E`$) is the energy of the scattered (incident) electron.

From the measured differential cross section, one typically derives a “reduced cross section”, defined according to

$$\sigma _R=\frac{d\sigma }{d\mathrm{\Omega }}\frac{(1+\tau )ϵ}{\sigma _{ns}\tau }=\frac{ϵ}{\tau }G_{Ep}^2(Q^2)+G_{Mp}^2(Q^2),$$
(3)
where $`ϵ=\{1+2(1+\tau )\mathrm{tan}^2(\theta _e/2)\}^1`$ is a measure of the virtual photon polarization. Equation 3 is known as the Rosenbluth formula, and shows that fits to reduced cross section measurements made at constant $`Q^2`$ but varying $`ϵ`$ values may be used to extract both form factors independently.

With increasing $`Q^2`$, the reduced cross sections are increasingly dominated by the magnetic term $`G_{Mp}`$; at $`Q^2`$ $``$ 3 GeV2, the electric term contributes only a few percent of the cross section. Furthermore, referring to the open data points in the left panel of Fig. 1, one can see that the various Rosenbluth separation data sets  for the ratio $`\mu _pG_{Ep}/G_{Mp}`$, where $`\mu _p=2.79`$ is the magnetic moment of the proton, are not consistent with one another for $`Q^2`$ $`>`$ 1 GeV2. It is clear that a tremendous effort has gone into the analysis of these difficult experiments, however, one is forced to speculate that some of the experiments have underestimated the systematic errors. For example, the Rosenbluth experiments apply radiative corrections to their data at leading order, i.e. one hard photon emitted, which can vary significantly with $`ϵ`$ depending on the specific experiment. However, higher order radiative corrections, i.e. involving more than one photon, could in fact change the slope of the reduced cross section versus $`ϵ`$ plot, and thus have an impact on the extracted form factor ratio.

Due to the fundamental nature of the quantities at hand, a more robust method for measuring the proton electromagnetic form factors is certainly desirable. Over the last few years, focal plane polarimeters have been installed in hadron spectrometers in experimental facilities at Bates, Mainz, and Jefferson Lab. Specifically, one makes use of the polarization transfer method , in which one measures, using a focal plane polarimeter, the transverse ($`P_t`$) and longitudinal ($`P_{\mathrm{}}`$) components of the recoil proton polarization in $`{}_{}{}^{1}H(\stackrel{}{e},e^{}\stackrel{}{p})`$ scattering, using a longitudinally polarized electron beam.

The proton form factor ratio is given simply by

$$r=\mu _p\frac{G_{Ep}}{G_{Mp}}=\frac{P_t}{P_{\mathrm{}}}\frac{\mu _p\left(E_e+E_e^{}\right)}{2M_p}\mathrm{tan}\left(\theta _e/2\right).$$
(4)
Here, $`E_e`$ ($`E_e^{}`$) is the incident (scattered) electron energy. The polarization transfer method offers a number of advantages over the traditional Rosenbluth separation technique. Using the ratio of the two simultaneously measured polarization components greatly reduces systematic uncertainties. For example, a detailed knowledge of the spectrometer acceptances, something which plagues the cross section measurements, is in general not needed. Moreover, it is not necessary to know either the beam polarization or the polarimeter analyzing power, since both of these quantities cancel in measuring the ratio of the form factors. The dominant systematic uncertainty is the knowledge of spin transport, although in comparison to the size of the systematic uncertainties in cross section measurements, this too is small. Finally, the most recent theoretical calculations  of the effects of radiative processes verify the assertion that the form factor ratio is unaffected at the momentum transfers involved here.

As mentioned above, the proton form factor ratio has been measured at several facilities using the polarization transfer technique. It was first used by Milbrath et al, who determined $`r`$ at $`Q^2`$ = 0.38 and 0.50 GeV2, and the result was in good agreement with both the Rosenbluth separation results and a subsequent polarization measurement at $`Q^2`$ = 0.4 by Dieterich et al at Mainz. However, polarization transfer measurements up to $`Q^2`$=3.5 GeV2 in Hall A at Jefferson Lab  have revealed the somewhat surprising result shown in Fig. 1, that the form factor ratio decreases with increasing $`Q^2`$. Most recently, this trend has been confirmed in Jefferson Lab Experiment E99007 , which extends the form factor ratio measurement to $`Q^2`$=5.6 GeV2; these new data are also shown in Fig. 1. We have fit the Jefferson Lab data using a simple linear parameterization, i.e.,

$$r=1.0(0.130\pm 0.005)\left[Q^2(0.04\pm 0.09)\right],$$
(5)
for 0.04 $`<`$ $`Q^2`$ $`<`$ 5.6 GeV2. This empirical description, which gives an acceptable fit to the Jefferson Lab data with the smallest number of free parameters, is shown in the right panel of Fig. 1 using two dashed curves to represent the range of uncertainty in the fit parameters. Note that this description of the ratio is singularly inconsistent with the global fit from Ref.  shown in the left panel of Fig. 1. However, a similar decreasing ratio has been reported in a number of recent theoretical models. In addition, we also show in the right panel of Fig. 1 ratios calculated from recent fits to the electromagnetic form factors by Lomon  within the framework of the Gari-Krümpelmann model. The two curves shown differ somewhat from one another in their specific choice of the form of the hadronic form factors as well as in the parametrization of the behaviour at large $`Q^2`$, and as such represent upper and lower bounds on the extracted form factor ratio. However, these curves both lie significantly above the new data.

In the remainder of this paper, for the purposes of reanalyzing the cross section data, we have used an empirical prescription for the form factor ratio. For $`Q^2`$ $`<`$ 0.04 GeV2, we have used $`r=1`$; for 0.04 $`<`$ $`Q^2`$ $`<`$ 7.7 GeV2, we have employed Eq. 5; for $`Q^2`$ $`>`$ 7.7 GeV2, we have used $`r=0`$. The boundary of 0.04 GeV2 (7.7 GeV2) corresponds to the value of $`Q^2`$ where Eq. 5 predicts a ratio of 1 (0). The choice of setting $`r=0`$ for the largest $`Q^2`$ region is somewhat arbitrary, since no experimental data exists in this kinematic regime. However, since the electric contribution to the total cross section is minimal in this $`Q^2`$ region, our choice has in fact little impact on the extracted value of the magnetic form factor. In addition, our linear fit almost certainly has the wrong asymptotic behavior at very large $`Q^2`$, based on theoretical expectations, and therefore we have extended the empirical fit only to $`Q^2`$=7.7 GeV2 to match the higher $`Q^2`$ assumption of $`r=0`$.

As stated above, the new Jefferson Lab data for the form factor ratio is in general disagreement with the higher $`Q^2`$ Rosenbluth separation measurements. Since the Rosenbluth separation measurements systematically attribute more strength in the cross section to the electric part (larger ratio, $`r`$), this means that their extracted values for the magnetic form factor, $`G_{Mp}`$, are potentially systematically too small.

We have reanalyzed the available cross section data using the following procedure:

* In the cases where experiments have extracted reduced cross section data at multiple $`ϵ`$ values at each $`Q^2`$ (Refs. ), we have reanalyzed these data using the aforementioned form factor ratio prescription. The net effect of this procedure is that the form factor ratio constraint fixes the ratio of the slope to the intercept of the graph of reduced cross section vs. $`ϵ`$. Therefore, in practice, one extracts only a single parameter from the new fit to the data.
* For the data of Sill et al. , which presents reduced cross section data at a single $`ϵ`$ value for each $`Q^2`$, we extract $`G_{Mp}`$ using the above form factor ratio prescription directly. The authors had assumed, quite reasonably at the time, $`r=1`$ in their extraction of $`G_{Mp}`$.
* It is important to recognize that Eqs. 1 and 2 describe the elementary single photon exchange electron scattering cross section, i.e., the Born cross section. However, other processes, such as vacuum polarization and radiative effects, contribute to the total cross section that one measures in an experiment. The technique that has been used to date is to correct the measured $`ep`$ cross section to account for the contribution from these extra processes, which are fully calculable within the framework of quantum electrodynamics, and thus extract the Born cross section . This being stated, the calculation of the vacuum polarization processes has varied between the different experiments. Specifically, only the most recent analyses  have included the $`\mu ^+\mu ^{}`$ and $`q\overline{q}`$ vacuum polarization contributions. These processes were not accounted for in the original work of Mo and Tsai . As a result, in our reanalysis of the cross section data we have included, using the same formalism as presented in , a correction to the extracted Born cross sections in the older experiments . The size of this correction is almost negligible for $`Q^2`$ $`<`$ 0.3 GeV2, but increases to about 1.3% at Q2=4 GeV2, which is the largest momentum transfer probed in the earlier measurments. We also point out that in the hadronic and particle physics communities, the vacuum polarization processes are often modelled in terms of a strengthening of the electromagnetic coupling, or alternatively in terms of a “screening” of the bare electron charge at low $`Q^2`$.

As an example of the effect of form factor ratio constraint, in Fig. 2, we show reduced cross section data which has been recalculated using the original data from Andivahis et al. for the four lowest $`Q^2`$ values of 1.75, 2.50, 3.25, and 4.00 GeV2. In each case, the dashed line is the best fit line using the direct Rosenbluth method. The solid line is the best fit using our form factor ratio constraint. The error bars shown are statistical and point-to-point systematic errors, as reported in Ref. , added in quadrature. The open and closed symbols in the figure are from two different spectrometers, known as the “1.6 GeV” and “8 GeV” spectrometers, respectively. In Ref. , the authors report that the data from these two spectrometers were cross-normalized at the two lowest $`Q^2`$ points, which resulted in a renormalization factor of $`0.956`$ being applied to the 1.6 GeV spectrometer data. In this paper, the authors also report an overall normalization uncertainty of $`\delta _{norm}=1.77\%`$. We have not included the renormalization uncertainty in each data point, but it has been included in the final uncertainty in the intercept of each graph, as

$$(\delta b_{final})^2=(\delta b_{raw})^2+(b_{raw}\delta _{norm})^2,$$
(6)
where $`b_{raw}`$ is the intercept of the straight-line fit to the data.

We note that our new fits would correspond to increasing the normalization factor of the 1.6 GeV spectrometer, or alternatively introducing a momentum dependence to the 1.6 GeV spectrometer acceptance. In an attempt to address this issue, we have used instead a renormalization factor of $`0.980`$ for the 1.6 GeV data, which represents an upper limit derived from the original normalization factor and uncertainty from Ref. .

In addition, in the final extraction of $`G_{Mp}`$ using the form factor ratio constraint, we have incorporated (in quadrature) the uncertainty in the form factor ratio itself, as expressed in Eq. 5. The results of the two extraction methods, including the final uncertainties, are summarized in Table I.

In Fig. 3, we show data for the proton magnetic form factor, expressed as $`G_{Mp}/\mu _pG_D`$, where $`G_D=(1+Q^2/0.71)^2`$ is the dipole form factor parameterization. In the left panel, we show the magnetic form factor as extracted using direct Rosenbluth separation (or using the assumption of $`r=1`$ in the case of the data of Sill et al.). In the right panel, we show the newly extracted data using the above constraint on the form factor ratio. In both panels, the dashed curve is the parameterization of Bosted , while in the right panel, the solid line is our new empirical fit, and is given by

$$\frac{G_{Mp}}{\mu _p}=\frac{1}{1+(0.116\pm 0.040)Q+(2.874\pm 0.098)Q^2+(0.241\pm 0.107)Q^3+(1.006\pm 0.069)Q^4+(0.345\pm 0.017)Q^5},$$
(7)
the form of which is consistent with Ref..

Imposing the form factor ratio constraint has the expected effect that the extracted magnetic form factor is systematically larger than when one uses the direct Rosenbluth method. Indeed, the data in the right panel of Fig. 3 lie approximately 1.5-3% above the Bosted parameterization. As mentioned earlier, the decreased electric strength implicit in a decreasing form factor ratio results in increased magnetic strength. However, perhaps the most striking feature of the reanalyzed data is that imposing the constraint on the form factor ratio results in uncertainties that are much reduced compared to using the direct Rosenbluth separation. This is due simply to the fact that we are extracting only a single parameter (proportional to $`G_{Mp}`$) from the cross section data, as opposed to extracting two parameters, as is done with the Rosenbluth method. Also, the Bosted parameterization and our new fit converge at large $`Q^2`$. As stated previously, the electric strength at large $`Q^2`$ decreases rapidly, and so indeed our choice of $`r=0`$, compared to the previous choice of $`r=1`$, has little effect on the magnetic form factor extracted from the data of Sill et al.

One also sees immediately that in using the new form factor constraint, the extracted magnetic form factor data from the various experiments are more consistent with one another, as well as with the new parameterization. Comparing the data extracted using direct Rosenbluth separation to the Bosted parameterization, we calculate $`\chi ^2/N_{d.f.}`$=0.97. This is somewhat larger than the value quoted in Ref. , since we have used a different data sample. Using the form factor constraint, and comparing to our new parameterization, $`\chi ^2/N_{d.f.}`$=0.82. Based on our new parameterizations of both $`G_{Mp}`$ and the form factor ratio, we have calculated the elastic hydrogen differential cross section directly using Eq. 1. We find that the deviation of the cross section calculated using the new fits from the measured world cross section data  form an approximately Gaussian distribution around -1%, with a standard deviation of $``$2.5%. As well, there is no significant $`Q^2`$ dependence of the deviation in the region 0.1 $`<`$ $`Q^2`$ $`<`$ 30 GeV2.

In conclusion, we have used recent measurements of the ratio of the elastic electromagnetic form factors of the proton, GEp/GMp, using the polarization transfer technique as a constraint in reanalyzing most of the world $`ep`$ elastic cross section data. We have presented a new empirical fit to the reanalyzed data for the proton elastic magnetic form factor in the region 0 $`<`$ $`Q^2`$ $`<`$ 30 GeV2, and find that over most of this kinematic region, the magnetic form factor is systematically 1.5-3% larger than had been extracted in previous analyses.

We thank R.L. Lewis for useful discussions regarding the vacuum polarization and radiative correction effects. This work was supported by the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# 1 Introduction

## 1 Introduction

In recent years a steadily growing interest in orbital angular momentum states of light can be observed. These light fields which are solutions of the scalar wave equation are mathematically described by Laguerre-Gaussian modes possessing a helical phase structure. As a result they have one or more phase singularities. The orbital angular momentum carried by these light fields is distinct from the angular momentum associated with polarization, it is quantized in units of $`\mathrm{}`$ and and can be converted into mechanical angular momentum . The possibility of using such light fields for driving micromachines, and their application as optical tweezers and optical traps make them possibly useful .

One may be led to believe that the orbital angular momentum of photon pairs created by spontaneous parametric down-conversion is not conserved if using classical detection methods. However from the quantum physics perspective most interesting are those applications which exploit the quantum properties of photons with orbital angular momentum. As already shown  spontaneous parametric down-conversion creating pairs of photons conserves orbital angular momentum on the individual photon level. Also, the two down-converted photons in a pair have been demonstrated to be in an entangled state with respect to their orbital angular momentum. As the Laguerre-Gaussian modes can be used to define an infinitely dimensional Hilbert space, orbital angular momentum entangled photons provide access to multi-dimensional entanglement which involves many orthogonal quantum states. Multi-dimensional entangled states are of considerable importance in the field of quantum information and quantum communication enabling, for example, quantum cryptography with higher alphabets.

For these quantum applications of orbital angular momentum it is necessary to be able to analyze orbital angular momentum entangled states of individual photons and one has to have quantitative measures for the multi-dimensional entanglement. Since analyzing entanglement locally always implies analyzing superpositions it is essential to have experimental techniques for preparing and analyzing superpositions of different orbital angular momentum eigenstates. One quantitative measure for multi-dimensional entanglement is a Bell inequality experiment generalized to more than two dimensions.

In the following we will present techniques, already realized in experiment, for preparing superpositions of Laguerre-Gaussian modes with arbitrary amplitude and phase ratios. This was done employing two different techniques, using computer generated holograms on the one hand and an interferometric method on the other hand. We will also present how these techniques will be applied in a Bell-inequality experiment as mentioned above, which is in progress in our laboratory.

## 2 Mathematical description of the Laguerre-Gaussian modes

The well-known Gaussian beam is not the only solution of the scalar paraxial wave equation. Also, the Hermite Gaussian (HG)  and the Laguerre-Gaussian (LG) modes for which the electromagnetic field amplitude is given by

$$u_{p,l}(r,\theta ,z)=\sqrt{\frac{2p!}{\pi (p+|l|)!}}\frac{1}{w(z)}\left(\frac{r\sqrt{2}}{w(z)}\right)^{|l|}L_p^{|l|}\left(\frac{2r^2}{w(z)^2}\right)e^{\frac{r^2}{w(z)^2}}e^{\frac{ikr^2}{2R(z)}}e^{i(2p+|l|+1)\mathrm{arctan}(\frac{z}{z_R})}e^{il\theta }$$
(1)
are solutions of this equation. With their two-fold infinite number of the indices both the LG and the HG modes build an orthogonal basis set for describing any paraxial transversal mode of the free propagation. For our further considerations we will only focus on the LG modes.

An LG mode (1) is characterized by its two indices $`p`$ and $`l`$ and by the standard Gaussian beam parameter definitions for the spot size $`w(z)`$, the radius of wavefront curvature $`R(z)`$ and the Rayleigh length $`z_R`$. The $`L_p^l(x)`$ term in (1) is a generalized Laguerre polynomial. The indices $`p`$ and $`l`$ are referred to as the radial and azimuthal mode index respectively, $`p+1`$ is the number of radial nodes and $`l2\pi `$ is the phase variation along a closed path around the beam center (Figure 1). This phase variation which is due to the $`e^{il\theta }`$ term in (1) results in a helical structure of the wave front. In consequence there is a phase singularity in the beam center for $`l0`$ and in order to satisfy the wave equation the intensity has to vanish there. Therefore such states are also called doughnut modes. Since the Laguerre-Gaussian modes (1) are angular momentum eigenstates they carry an orbital angular momentum of $`\mathrm{}l`$ per photon. The fact Laguerre Gaussian modes carry an orbital angular momentum was predicted and it was experimentally verified .

As usual in quantum mechanics, external variables describe the quantum state in real space while internal variables refer to additional variables. In that sense it is important to stress that the angular momentum carried by Laguerre Gaussian modes with $`l0`$ is an external angular momentum distinct from the internal angular momentum of the photons associated with their spin .

## 3 Production of the Laguerre-Gaussian modes

There are several experimental methods like cavity induced production, astigmatic mode conversion  and the use of computer generated holograms  for creating LG modes. In this article we restrict ourselves to describing the use of computer generated holograms.

A hologram is a recording of the interference pattern between the desired field and some reference field. The simplest possible reference field is the plane wave.

$$R=R_0e^{ik_xx+ik_zz}$$
(2)
The interference pattern produced by such a beam propagating at an angle $`\xi =\mathrm{arctan}(\frac{k_x}{k_z})`$ and e.g. an $`LG_{01}`$ mode propagating along the z-direction can be calculated numerically (Figure 2a). It is a line grating with one dislocation in the center. Now, if this grating is illuminated by the reference beam, which is sufficiently well realized if this hologram is placed at the waist of a Gaussian beam, the $`LG_{01}`$ mode is reproduced. Intuitively speaking the phase dislocation exerts a “torque” onto the diffracted beam because of the difference of the local grating vectors in the upper and lower parts of the grating. This “torque” depends on the diffraction order n and on the number of dislocations $`\mathrm{\Delta }m`$ of the hologram. Consequently the right and left diffraction orders gain different handedness and the associated orbital angular momentum values differ in their sign. The definition of the sign can be chosen by convention. The modulation of the incoming beam can be either done in the phase using transmission phase gratings or reflection gratings or in the amplitude using absorption gratings the latter being rather inefficient.

In our experiments we used transmission phase gratings with a diffraction efficiency of 70$`\%`$ after blazing. The binary structure of the hologram in (Figure 2a) is modified by blazing such that it results in a pattern as shown in (Figure 2b). The transmission function of such a hologram in polar coordinates is

$$T(r,\varphi )=e^{i\delta \frac{1}{2\pi }\mathrm{mod}(\mathrm{\Delta }m\varphi \frac{2\pi }{\mathrm{\Lambda }}r\mathrm{cos}\varphi ,2\pi )}$$
(3)
where $`\delta `$ is the amplitude of the phase modulation and the second factor in the superscript is the actual pattern of the blazed hologram. $`\mathrm{\Lambda }`$ is the spacing period of the grating and $`\mathrm{mod}(a,b)=ab\mathrm{Int}\left(\frac{a}{b}\right)`$. As mentioned above the fraction of intensity diffracted into higher orders of the hologram that is for a grating with one dislocation to the higher order LG modes with $`l1`$ can be decreased by blazing. But the blazing has no influence on the composition of the output beam’s radial index p-terms. The relative amplitudes of the p-terms depend on the choice of the beam parameters of the input beam to output beam . However for our experimental applications using holograms with one dislocation the relative amplitudes of the $`p0`$-terms become negligible.

## 4 Superpositions of the Gaussian and the Laguerre-Gaussian modes

Central to many, if not all, quantum experiments is the concept of superposition. It is therefore important to be able to both produce and analyze superpositions of the various states of a chosen basis. We therefore discuss now superpositions of the LG modes presented above. There are several experimental methods for producing such superpositions. One simple way is to use a Mach-Zehnder interferometer as sketched in (Figure 3a). After the input mode is split by the first beam splitter the beam of each arm is sent through a hologram causing the desired mode transformations which are in general different transformations in the two arms. For producing superpositions of the Gaussian and the $`LG_{01}`$ mode it is sufficient to place a hologram with one dislocation only in one of the arms of the interferometer. The two beams are brought together on a second beam splitter where they are superposed. An advantage of this method is that by attenuating each of the arms and using a phase plate superpositions with arbitrary amplitudes and relative phases can be produced without changing the setup. The resulting interference pattern for a superposition of an $`LG_{00}`$ mode (=Gaussian mode) and an $`LG_{01}`$ mode is shown in Figure 3b.

However the interferometric preparation of superposition modes has also some disadvantages. The experimental setup becomes too complex and too difficult to control when one has to create and analyze superpositions several times. It would be necessary to keep the relative phase in the interferometer arms stable and one also has to take the different divergences of the interfering modes into account. A more convenient but less general method for creating superposition modes is to use a displaced hologram . This method is particularly suitable for producing superpositions of an $`LG_{0l}`$ with the Gaussian mode which may also be seen as an $`LG`$ mode with $`l=0`$. The transmission function of a hologram which is designed to transform a Gaussian mode into an $`LG_{01}`$ mode is given by (3) with $`\mathrm{\Delta }m=1`$. In order to tranform an incoming Gaussian beam into an $`LG_{01}`$ beam the hologram should be placed at the waist of the Gaussian beam and the beam should be sent through the center of the hologram where the dislocation is located. The intensity pattern of such an $`LG_{01}`$ mode posseses a centrally located singularity. By shifting the dislocation out of the beam center step by step one can experimentally observe that the singularity becomes eccentric resulting in the same pattern achieved by the interferometric setup.

Numerical simulations show that there are also higher order $`LG_{0l}`$ components present in the case when the superposition is produced by a displaced hologram. The traces in the bottom of Figure 5c are corresponding to the amplitudes of the $`LG_{01}`$, $`LG_{02}`$ and $`LG_{02}`$ modes. This is also required by the unitary of the procedure. Nevertheless the relation (5) between the relative amplitudes and the position of the singularity is still a good approximation because the amplitudes of these higher orders are small. In the case where the beam is sent through a border region of the hologram far away from the dislocation the output beam is again a Gaussian beam because there the hologram acts as an ordinary grating.

An important question is whether in the experiment one actually observe coherent superpositions rather than incoherent mixtures. The distinction between coherent superposition and incoherent mixture of Gaussian and LG modes is that the latter posses no phase singularity. This is because adding the spatial intensity distributions of these two modes will yield a finite intensity everywhere in the resulting pattern. In contrast, in a coherent superposition the amplitudes are added and therefore the phase singularity must remain and is displaced to an eccentric location (Figure 3b). It will appear at that point where the amplitudes of the two modes are equal with opposite signs. Therefore the radial distance of the singularity from the beam center is a measure of the amplitude ratio of the Gaussian to the $`LG`$ components whereas the angular position of the singularity is determined by their relative phase.

To obtain quantitative results we calculated the intensity distribution of a normalized superposition mode of an $`LG_{00}`$ and an $`LG_{01}`$ mode described by

$$\frac{1}{\sqrt{1+\gamma ^2}}\left[|u_{00}+\gamma e^{i\phi }|u_{01}\right]$$
(4)
and looked for the position of its phase singularity. Here $`|u_{00}`$ and $`|u_{01}`$ denote the amplitudes of an $`LG_{00}`$ and an $`LG_{01}`$ respectively as given by (1), $`\gamma `$ is the relative amplitude of the $`LG_{01}`$ mode and $`\phi `$ is the relative phase of the interfering modes. After inserting the corresponding $`LG`$-amplitudes (1) into (4) we found for the position of the singularity the cylindrical coordinates

$$r=\frac{w_0}{\gamma \sqrt{2}}\text{,}\theta =\phi ,$$
(5)
where $`w_0`$ denotes the waist size of the Gaussian beam.(Figure 3b).

Although this expression only holds for superpositions of $`LG_{00}`$ and $`LG_{0\pm 1}`$ modes it can easily be generalized to superpositions of higher order LG mode containing more than just one phase singularity.

In order to prove that the eccentric mode is indeed a superposition mode one has to project the superposition onto the orthogonal basis states. Experimentally this was done by sending a Gaussian laser beam (HeNe, 632nm) through a displaced hologram producing a superposition of the Gaussian and the $`LG_{01}`$ mode. In the next step the output mode was projected onto the Gaussian mode by coupling into a mono-mode optical fiber. Since all other modes have a larger spatial extension, only the Gaussian mode can propagate in the mono-mode fiber and it therefore acts as a filter for the $`LG_{0l}`$ modes with $`l0`$.

Having only a filter transmitting the Gaussian mode the $`LG_{01}`$ mode had to be identified via an additional step. This was done by introducing a second hologram making a mode transformation on the output beam reducing the azimuthal index $`l`$ by one before coupling into mono-mode optical fibers(Figure 4). For each position of the displaced hologram the transmitted intensity to the Gaussian and the $`LG_{01}`$ detector was measured.

The results are shown in Figure 5a and 5b. When the hologram is centered the incoming Gaussian mode is transformed into an $`LG_{01}`$ mode. As a result the intensity of light coupled into the Gauss-detector is a minimum and the one coupled into the $`LG_{01}`$-detector a maximum. As the hologram is shifted the intensity at Gauss-detector increases and the intensity at the $`LG_{01}`$-detector drops. The asymmetry in Figure 5a) is a result of the imperfection of the holograms . However the extinction ratio $`e`$ for all measurements was always far better than $`1:300`$.

These results are in agreement with our numerical calculations of the superposition modes(Figure 5c). The action of the hologram on the incoming field is characterized by (3). The transmitted beam directly after the hologram is given by

$$u_{out}(r,\varphi )=T(r,\varphi )u_{in}(r,\varphi )$$
(6)
Denoting the relative position of the hologram with respect to the beam center with $`(r_0,\varphi _0)`$ one finds the projection of a transmitted $`LG_{0l}`$ mode onto the $`LG_{0L}`$-mode to be

$$a_L^l(r_0,\varphi _0)=_{\mathrm{}}^+\mathrm{}_0^{2\pi }r𝑑r𝑑\varphi \left(u_{0L}(r,\varphi ,0)e^{i\mathrm{\Delta }m\frac{2\pi }{\mathrm{\Lambda }}r\mathrm{cos}\varphi }\right)^{}T(rr_0,\varphi \varphi _0)u_l(r,\varphi ,0)$$
(7)
This expression represents the amplitude of an $`LG_{0L}`$-mode when a hologram with one dislocation placed at $`(r_0,\varphi _0)`$ is illuminated by an $`LG_{0l}`$ mode. The numerical simulation of a setup which corresponds to our experiment is given in Figure 5c.

As mentioned above, the superposition does not only consist of two modes. The output mode will also contain higher order modes (Figure 5c bottom). This is a necessary consequence of the fact that the action of the hologram is unitary. In general, the relative amplitudes of the higher order modes will increase when the hologram has many dislocations or when the single dislocation hologram is illuminated by an $`LG_{0l}`$ beam with $`|l|>1`$. However they are negligible (Figure 5c bottom) for the case that the hologram is illuminated by an $`LG_{01}`$ mode.

## 5 Conclusion and Outlook

In this work we showed techniques for producing and analyzing superpositions of Laguerre-Gaussian modes. Our experimental results show that it is possible to achieve mode detection of high distinction ratio ($`1:300`$) with the technique described above. As recently shown the orbital angular momentum of photons is conserved in parametric down-conversion  and the down-converted photons are found to be in an entangled state with respect to the orbital angular momentum. Since orbital angular momentum entangled photons can be used as qu-nits they open a practical approach to multi-dimensional entanglement where the entangled states do not only consist of two orthogonal states but of many of them. We expect such states to be of importance for the current efforts in the field of quantum computation and quantum communication. For example, quantum cryptography with higher alphabets could enable one to increase the information flux through the communication channels.

However the ultimate confirmation of the entanglement of the orbital angular momentum will be a Bell inequality experiment generalized to more states . Such an experiment will also give a quantitative measure of the multi-dimensional entanglement produced by parametric down-conversion. Employing the techniques presented here for preparing and analyzing superpositions of orbital angular momentum states such an experiment is presently in progress in our laboratory.

Acknowledgment

The Authors would like to thank Marek Zukowski for discussion and Jochen Arlt for correspondence. This work was supported by the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF).


text:
# NUMERICAL SIMULATIONS OF INCOMPRESSIBLE MHD TURBULENCE

## I INTRODUCTION

Astrophysical flows are complicated and dynamic. For example, molecular clouds are clumpy over a vast range of scales and observations of velocity spectral line broadening indicate that dynamic pressure ($`\rho 𝐯^2`$/2) usually dominates the thermal pressure ($`nkT`$). Furthermore, velocity spectral lines show multiple components. These facts suggest that cloud internal velocity possesses disordered components. In molecular clouds, the characteristic Reynolds number is estimated to be larger than $`10^{10}`$, which is consistent with the turbulent state in the clouds. In general, when the Reynolds number ($`LV/\nu `$; L=characteristic size of the system, V=velocity dispersion, and $`\nu `$=viscosity) is larger than an order of 10 - 100, the system becomes unstable and turbulent.

In laboratory systems, fluids are usually incompressible and unmagnetized. Kolmogorov phenomenology provides an excellent insight for such systems. Suppose that we ‘disturb’ the fluid at a scale $`L`$. We call this scale as energy injection scale or largest energy containing eddy scale. Then the energy injected to the scale $`L`$ cascades to progressively smaller and smaller scales. Ultimately, the energy will reach the molecular dissipation scale $`l_d`$ and it will be lost there. The scales between $`L`$ and $`l_d`$ is called the inertial range. Suppose a scale $`l`$ lies in the inertial range. Let the characteristic velocity associated with the scale be $`v_l`$. Kolmogorov theory states that the kinetic energy ($`v_l^2`$) is transferred to one-level smaller scale within one eddy turnover time ($`l/v_l`$). It is natural that the cascade rate ($`v_l^2/(l/v_l)`$) be scale-independent. Therefore, we have

$$v_ll^{1/3}.$$
(1)
(One-dimensional) Energy spectrum $`E(k)`$ is one of the most important quantities in turbulence theories. Note that $`E(k)dk`$ is the amount of energy between the wavenumber $`k`$ and $`k+dk`$. When $`E(k)`$ follows a power law, $`kE(k)`$ is the energy near the wavenumber $`k1/l`$. Since $`v_l^2`$ represents a similar energy, we have $`v_l^2kE(k)`$. Therefore, equation (1) becomes

$$E(k)k^{5/3}.$$
(2)
This is the well-known Kolmogorov spectrum. Astrophysical fluids are compressible and magnetized. Nevertheless, many astrophysical quantities follow this Kolmogorov spectrum. It is important to note that Kolmogorov theory predicts: 1. the velocity dispersion ($`v_l`$) decreases with scales as $`v_ll^{1/3}`$; 2. the slope of the energy spectrum is $`5/3`$.

Earlier works on ISM focused on the velocity dispersions and derived a qualitative agreement with the Kolmogorov theory. Larson (1981) and Myers (1983) found that velocity dispersion decreases as the size of the clumps decreases, which may support turbulent model of ISM. Scalo (1984) showed that two-point velocity correlation increases as the separation decreases, which is also qualitatively consistent with turbulent model. All the above mentioned studies are for scales larger than or similar to $``$ 1 pc.

On the other hand, other studies address the issue of power spectra. Interstellar scintillation observations indicate electron density spectrum follows a power law over 7 decades of length scales (see Armstrong et al. 1995). The slope of the spectrum is very close to $`5/3`$ for $`10^6m`$ \- $`10^{14}m`$. The spectrum is somewhat uncertain for scales larger than $`10^{14}m`$ ($`0.01pc`$). However, recently, Lazarian & Pogosyan (2000) and Starnimirovic & Lazarian (2001) showed the existence of the $`5/3`$ velocity power spectrum over pc-scales in HI. Solar-wind observations provide in-situ measurements of the power spectrum of magnetic fluctuations. Leamon et al. (1998) obtained a slope of $`1.7`$.

Therefore it is evident that turbulence is prevalent in astrophysical fluids. Understanding the nature of such turbulence is important because it affects many astrophysical processes. For example, star formation rate depends on the turbulence decay time-scale. Heat and cosmic ray transport also depend on the statistics of turbulence. Recently anisotropic MHD Turbulence has also been evoked to provide fast magnetic reconnection rate (Lazarian & Vishniac 1999).

One approach to studying the ISM is to perform time-dependent numerical simulations to model the ISM, including as many of the interacting phenomena as practical. Of course, the physics included in such models must necessarily be highly simplified, and it is difficult to determine which features of the final model result from which physical assumptions (or initial conditions). Our approach is to use simplified numerical simulations to study the influences of various physical phenomena in isolation. We want to obtain a physical feeling for the general effects that each phenomenon has on the nature of the ISM. Here we review our previous works on incompressible MHD turbulence. It is obvious that the real ISM is compressible, but we want to separate the effects of magnetic turbulence from those involving compression. In later works we will include compression for comparison with the present models, thereby isolating its importance directly.

As we said earlier, we will numerically study incompressible MHD turbulence. We believe Fourier spectral method is one of the best numerical methods for such calculations. In II, we show how the spectral method works. In III, we introduce our recent works done with the spectral code: anisotropy, decay speed of MHD turbulence, and magnetic structures below the viscous cut-off. In IV, we give conclusion.

## II NUMERICAL METHOD

### (a) PSEUDO-SPECTRAL METHOD

For simplicity, let us first consider incompressible unmagnetized fluid. Then, the governing equation is the Navier-Stocks equation:

$$\frac{𝐯}{t}=(\mathrm{𝐯𝐯})+\nu ^2𝐯+P,$$
(3)
where $`𝐯`$ is velocity, $`\nu `$ is viscosity, and $`P`$ is gas pressure. Since we are considering an incompressible fluid, we can assume that $`\rho =1`$.

The Fourier spectral method is highly suitable for incompressible flows. Simply put, the idea of the method is simple: we transform all physical space variables into Fourier (or wave-vector) space and solve the differential equation there. Hereafter, we assume that we are solving the equation on a uniform rectangular grid. Some advantage of the method is:

1. It shows fast convergence for smooth functions.
2. It is easy to implement divergence-free conditions.
3. It is easy to implement periodic boundary condition.
4. Fourier space analyses (e.g. calculation of energy
spectrum) are natural consequences.

We note that the gradient operator becomes $`i𝐤`$ and that multiplication becomes convolution sum in Fourier space. Here $`i=\sqrt{1}`$. Using these facts, we obtain the Navier-Stocks equation in Fourier space:

$$\frac{\widehat{𝐯}(𝐤)}{t}=i𝐤(\widehat{\mathrm{𝐯𝐯}})(𝐤)i\nu k^2\widehat{𝐯}(𝐤)+i𝐤\widehat{P}(𝐤),$$
(4)
where

$$(\widehat{\mathrm{𝐯𝐯}})(𝐤)_{𝐩+𝐪=𝐤}𝑑𝐤\widehat{𝐯}(𝐩)\widehat{𝐯}(𝐪).$$
(5)
The incompressibility condition becomes

$$i𝐤\widehat{𝐯}(𝐤)=0,$$
(6)
which means that the velocity vector stays perpendicular to the wave vector $`𝐤`$. Since the pressure term in equation (4) is parallel to $`𝐤`$, we can drop the term. On the other hand, the nonlinear term (i.e the term containing $`\widehat{\mathrm{𝐯𝐯}}`$) may have both parallel and perpendicular components. Therefore, we need to remove the component parallel to $`𝐤`$ by applying an appropriate projection operator.

The calculation of the nonlinear term is the most time consuming part of the Fourier spectral methods. For simplicity, let us consider the following one dimensional convolution sum:

$$_{p+q=k}𝑑k\widehat{f}(p)\widehat{g}(q).$$
(7)
Then the discrete representation of the convolution sum is

$$\underset{p+q=k}{}\widehat{f}_p\widehat{g}_q.$$
(8)
The pseudo-spectral method approximates the convolution sum as

$$\underset{p+q=k}{}\widehat{f}_p\widehat{g}_q\left[F^1\left[F(\widehat{f})F(\widehat{g})\right]\right]_k,$$
(9)
where $`F`$ and $`F^1`$ denote forward/inverse Fourier transform. The cost of such calculation is $`O(N\mathrm{log}N)`$, where $`N`$ is the number of grid points.

Once we obtain the convolution sum, we can use any techniques for ordinary differential equations to solve the time evolution of the Navier-Stocks equation.

### (b) ACTUAL SETUP

In actual calculations, we have adopted a pseudospectral code and calculated the time evolution of incompressible magnetic turbulence subject to a random driving force per unit mass:

$$\frac{𝐯}{t}=(\times 𝐯)\times 𝐯(\times 𝐁)\times 𝐁+\nu ^2𝐯+𝐟+P^{},$$
(10)
$$\frac{𝐁}{t}=\times (𝐯\times 𝐁)+\eta ^2𝐁,$$
(11)
$$𝐯=𝐁=0,$$
(12)
where $`𝐟`$ is a random driving force, $`P^{}P/\rho +𝐯𝐯/2`$, $`𝐯`$ is the velocity, and $`𝐁`$ is magnetic field divided by $`(4\pi \rho )^{1/2}`$. In this representation, $`𝐯`$ can be viewed as the velocity measured in units of the r.m.s. velocity, v, of the system and $`𝐁`$ as the Alfven speed in the same units. The time $`t`$ is in units of the large eddy turnover time ($`L/v`$) and the length in units of $`L`$, the inverse wavenumber of the fundamental box mode. In this system of units, the viscosity $`\nu `$ and magnetic diffusivity $`\eta `$ are the inverse of the kinetic and magnetic Reynolds numbers respectively. The turbulence is driven by 21 random forcing components in Fourier space. The peak of energy injection occurs at $`k2.5`$. The amplitudes of the forcing components are tuned to ensure $`v1`$ We use exactly the same forcing terms for all simulations. We use up to $`256^3`$ collocation points. We use an integration factor technique for kinetic and magnetic dissipation terms and a leap-frog method for nonlinear terms. At $`t=0`$, the magnetic field has only its uniform component and the velocity field is restricted to the range $`2k4`$ in wavevector space.

For some runs, we use hyperviscosity and hyperdiffusivity for the dissipation terms. The power of hyperviscosity is set to either 2 or 8, so that the dissipation term in the above equation is replaced with

$$\nu _h(^2)^h𝐯,$$
(13)
where $`h=2`$ or $`8`$. A similar expression is used for the magnetic dissipation term. In subsections III(a) and III(b), we use unit magnetic Prandtl number ($`\nu =\eta `$) and the notation NY-$`B_0`$Z. In subsection III(c), we use physical viscosity and physical- or hyper-diffusion and the notation NXY-$`B_0`$Z. Here N = 256, 144 refers to the number of grid points in each spatial direction; X = P refers to physical viscosity; Y = P, H2, H8 refers to the form of diffusion (physical-, the 2nd power hyper-, or the 8th power hyper-diffusion, respectively); Z=0.5 or 1 refers to the strength of the external magnetic field in Alfven velocity unit.

Diagnostics for our code can be found in Cho and Vishniac (2000b).

## III PROPERTIES OF MHD TURBULENCE

In this section, we discuss the properties of MHD turbulence in the presence of a strong uniform background field.

### (a) ANISOTROPY

Historically hydrodynamic turbulence in an incompressible fluid was successfully described by the eddy cascade (Kolmogorov 1941), but MHD turbulence was first modeled by wave turbulence (Iroshnikov 1963, Kraichnan 1965). This theory assumes isotropy of the energy cascade in Fourier space. However, the magnetic field defines a local symmetry axis since it is easy to mix field lines in directions perpendicular to the local $`𝐁`$ and much more difficult to bend them. In a turbulent medium, the kinetic energy associated with large scale motions is greater than that of small scales. However, the strength of the local mean magnetic field is almost the same on all scales. Therefore, it becomes relatively difficult to bend magnetic field lines as we consider smaller scales, leading to more pronounced anisotropy. A self-consistent model of MHD turbulence which incorporates this concept of scale dependent anisotropy was introduced by Goldreich & Sridhar (1995; hereinafter GS95).

The major predictions of the GS95 model are as follows. First, anisotropy is scale dependent. Roughly speaking, the semi-major axis ($`1/k_{}`$) and the semi-minor axis ($`1/k_{}`$) of an eddy satisfy the relation $`k_{}k_{}^{2/3}`$, which means smaller eddies are relatively more elongated than larger ones. Here $`k_{}`$ and $`k_{}`$ are wave numbers perpendicular and parallel to the background field. Second, the model predicts that the one-dimensional energy spectrum is Kolmogorov-type if expressed in terms of perpendicular wavenumbers, i.e. $`E(k_{})k_{}^{5/3}`$.

Numerical simulations by Cho & Vishniac (2000a) and Maron & Goldreich (2001) have mostly supported the GS95 model and helped to extend it. Both analyses stressed the point that scale dependent anisotropy can be measured only in local coordinate frames which are aligned with the locally averaged magnetic field direction. Cho & Vishniac (2000a) calculated the structure functions of velocity and magnetic field in the local frames, and found that the contours of the structure functions do show scale dependent anisotropy, consistent with the predictions of the GS95 model. Fig 1. shows that the semi-major axis ($`1/k_{}`$) is proportional to the 2/3 power of the semi-minor axis ($`1/k_{}`$), implying that $`k_{}k_{}^{2/3}`$. One dimensional energy spectrum follows Kolmogorov spectrum: $`E(k)k^{5/3}`$ (see Cho & Vishniac 2000a).

### (b) DECAY OF MHD TURBULENCE

Turbulence plays a critical role in molecular cloud support and star formation and the issue of the time scale of turbulent decay is vital for understanding these processes.

If MHD turbulence decays quickly then serious problems face researchers attempting to explain important observational facts, i.e. turbulent motions seen within molecular clouds without star formation (see Myers 1999) and rates of star formation (Mckee 2000). Earlier studies attributed the rapid decay of turbulence to compressibility effects (Mac Low 1999). Cho et al. (2001), as well as earlier ones (Cho & Vishniac 2000a, Maron & Goldreich 2001), shows that turbulence decays rapidly even in the incompressible limit. This can be understood in the framework of GS95 model where mixing motions perpendicular to magnetic field lines form hydrodynamic-type eddies. Such eddies, as in hydrodynamic turbulence, decay in one eddy turnover time.

Below we consider the effect of imbalance on the turbulence decay time scale. MHD turbulence can be described by opposite-traveling wave packets. ‘Imbalance’ means that wave packets traveling in one direction have significantly larger amplitudes than the other. In astronomy, many energy sources are localized. For example, SN explosions and OB winds are typical point energy sources. Furthermore, astrophysical jets from YSOs are believed to be highly collimated. With these localized energy sources, it is natural to think that interstellar turbulence is typically imbalanced.

In this subsection, we explicitly relate the degree of imbalance and the decay time scale of turbulence in the presence of a strong uniform background field.

In Figure 2 we demonstrate that an imbalanced cascade does extend the lifetime of MHD turbulence. We use the run 144H-$`B_0`$1 to investigate the decay time scale. We ran the simulation up to t=75 with non-zero driving forces. Then at t=75, we turned off the driving forces and let the turbulence decay. At t=75, there is a slight imbalance between upward and downward moving components ($`E_+=0.499`$ and $`E_{}=0.40`$). This results from a natural fluctuation in the simulation. The case of $`(E_+)_{t_0}=80\%(E_{})_{t_0}`$ corresponds to the simulation that starts off from this initial imbalance. In other cases, we either increase or decrease the energy of $`𝐳^{}`$ components and, by turning off the forcing terms, let the turbulence decay. We can clearly observe that imbalanced turbulence extends the decay time scale substantially. Note that we normalized the initial energy to 1. The y-axis is the total (=up $`+`$ down) energy.

In this section, we found that turbulence decay time can be slow. If we consider the interstellar medium at larger scales, it is definitely compressible, has a whole range of energy injection/dissipation scales (see Scalo 1987), and the relative role of vortical versus compressible motions being unclear. Nevertheless, we believe that our simplified treatment may still elucidate some of the basic processes. To what extend this claim can be justified will be clear when we compare compressible and incompressible results. However, if we accept that fast and slow MHD modes are subjected to fast collisionless damping (see Minter & Spangler 1997) the remaining modes are incompressible Alfven modes. Those should be well described by our model when turbulence is supersonic but sub-Alfvenic. Our preliminary results (Cho & Lazarian 2002) show that the coupling of the modes is marginal even in compressible regime.

### (c) STRUCTURES BELOW VISCOUS CUT-OFF

In hydrodynamic turbulence viscosity sets a minimal scale for motion, with an exponential suppression of motion on smaller scales. Below the viscous cutoff the kinetic energy contained in a wavenumber band is dissipated at that scale, instead of being transferred to smaller scales. This means the end of the hydrodynamic cascade, but in MHD turbulence this is not the end of magnetic structure evolution. For viscosity much larger than resistivity, $`\nu \eta `$, there will be a broad range of scales where viscosity is important but resistivity is not. On these scales magnetic field structures will be created through a combination of large scale shear and the small scale motions generated by magnetic tension. As a result, we expect a power-law tail in the energy distribution, rather than an exponential cutoff. To our best knowledge, this is a completely new regime of MHD turbulence. Ambipolar diffusion damps turbulent motions in ISM. Therefore, the study of such a viscous damped MHD turbulence is particularly important for many astrophysical processes below the ambipolar diffusion scale in ISM.

In Cho, Lazarian, & Vishniac (2002), we numerically demonstrate the existence of the power-law magnetic energy spectrum below the viscous damping scale. We use the same pseudo-spectral incompressible MHD code. We use physical viscosity for velocity. The kinetic Reynolds number is around 100. We use three different values of magnetic diffusion: physical diffusion, hyper-diffusion with order 2, and hyper-diffusion with order 8. We will present a theoretical model for this regime in an upcoming paper (Lazarian, Vishniac, & Cho 2002).

In Figure 3, we plot energy spectra. The spectra consist of several parts. First, the peak of the spectra corresponds to the energy injection scale. Second, for $`2<k<6`$, kinetic and magnetic spectra follow a similar slope. This part may be the short inertial range in general sense. Third, magnetic and kinetic spectra decouple at $`k6`$. Forth, for $`20<k`$, it appears that a new damped-scale inertial range emerges. The emergence of the new inertial range is conspicuous in the order 8 hyperdiffusion simulation.

## IV SUMMARY

We described a pseudo-spectral MHD code and demonstrated its usefulness in astrophysics through our recent works. We showed that the rate at which MHD turbulence decays depends on the degree of energy imbalance between wave packets traveling in opposite directions. A substantial degree of imbalance can substantially extend the decay time scale of the MHD turbulence, which might be useful for the study of star formation. We also showed that magnetic fields can have rich structures below the viscous cut-off scale, which implies that magnetic fields have stochastic structures below the ambipolar diffusion scale in ISM and, therefore, affect many astrophysical processes, such as cosmic ray transports, magnetic reconnection, and thermal diffusion.

ACKNOWLEDGEMENTS I thank Alex Lazarian, Ethan T. Vishniac, and Peter Goldreich for valuable discussions. The introduction of this proceeding is partly based on email communication with John Scalo. This work was partially supported by National Computational Science Alliance under CTS980010N and AST000010N and utilized the NCSA SGI/CRAY Origin2000.

### REFERENCES

Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995, ApJ, 443, 209

Cho, J. & Vishniac, E. T. 2000a, ApJ, 539, 273

Cho, J. & Vishniac, E. T. 2000b, ApJ, 538, 217

Cho, J., Lazarian, A., & Vishniac, E. T. 2001, ApJ, in press

Cho, J. & Lazarian, A. 2002, in preparation

Cho, J., Lazarian, A., & Vishniac, E. T. 2002, in preparation

Goldreich, P. & Sridhar, H. 1995, ApJ, 438, 763 (GS95)

Iroshnikov, P. 1963, Astron. Zh., 40, 742 (English version: 1964, Sov. Astron., 7, 566)

Kolmogorov, A. 1941, Dokl. Akad. Nauk SSSR, 31, 538

Kraichnan, R. 1965, Phys. Fluids, 8, 1385

Larson, R. B. 1981, MNRAS, 194, 809

Lazarian, A. & Pogosyan, D. 2000, ApJ, 537, 720L

Lazarian, A.& Vishniac, E. T. 1999, ApJ, 517, 700

Lazarian, A.& Vishniac, E. T., & Cho, J. 2002, in preparation

Leamon, R. J., Smith, C. W., Ness, N. F., & Matthaeus, W. H. 1998, J. Geophys. Res., 103, 4775

Mac Low, M. 1999, ApJ, 524, 169

Maron, J. & Goldreich, P. 2001, ApJ, 554, 1175

McKee, C.F. 1999, astro-ph/9901370

Minter, A. & Spangler, S. 1997, ApJ, 485, 182

Myers, P.C. 1983, ApJ, 270, 105

Myers, P.C. 1999, in The Origin of Stars and Planetary Systems, ed. by Charles J. L. and Nikolaos D.K., Kluwer Academic Publishers, 1999, p.67

Scalo, J. M. 1984, ApJ, 277, 556

Scalo, J. M. 1987, in Interstellar Processes, ed. D. J. Hollenbach & H. A. Thronson (Dordrecht: Reidel), 349

Stanimirovic, S., & Lazarian, A. 2001, ApJ, in press


text:
# Self-consistent quantum effects in the quark meson coupling model

## I Introduction

Usually the frame of Quantum Hadrodynamics (QHD)  is the departure to the study of the nuclear many-body problem describing nucleons interacting with scalar and vector mesons. This meson field theory has quite successfully described the properties of the nuclear matter and finite nuclei using the mean field approximations for the meson fields. The vacuum polarization corrections arising from the nucleon fields as well as the meson fields have also been considered to study the nuclear matter. This is one way of obtaining a softer equation of state yielding a lower compressibility than would be reached without quantum effects.

While descriptions of the nuclear phenomena have been efficiently formulated using some hadronic degrees of freedom as in QHD, there have been interesting observations which reveals the medium modification of the internal structure of the nucleon. For this, degrees of freedom from the fundamental theory of strong interacting systems, QCD, are expected to be considered. Due to the complex structure of this theory we are lead to formulate effective models which have the main properties and symmetries of QCD as chiral symmetry and its spontaneously symmetry breaking. One of the first models put forward along these lines was the quark-meson coupling (QMC)  model. It is an extension of QHD which includes the explicit quark structure of the baryons. This model describes nuclear matter with nucleons as non-overlapping MIT bags containing quarks inside them which interact with the scalar and the vector mesons. Recently, the model has been extended to finite temperatures , neutron stars  as well as to study beyond the mean field level by the incorporation of the exchange terms . The density dependent bag constant has also been investigated within this frame

In the present work, the vacuum polarization corrections, described in ref. , are included in the quark meson coupling model taking into account the quantum fluctuations of the scalar field (as done in ) self-consistently as well. A similar approach is being investigated in the frame of QHD . The equation of state for dense matter is derived. We organize the paper as follows: In section II, we derive the EOS for dense matter including the quantum effects arising from the nucleons and the $`\sigma `$-mesons within QMC model. In section III, we discuss in detail the numerical results obtained in the present work and discuss possible outlook.

## II Theory

The details of the QMC model have been given in Ref. . Since we now include the vacuum polarization effects, we give here a few important steps for completeness.

In this model, the nucleon in nuclear matter is assumed to be described by a static MIT bag in which quarks interact with the scalar ($`\sigma `$) and the vector ($`\omega `$) mesons. The quark field $`\psi _q(\stackrel{}{r},t)`$ inside the bag then satisfies the equation

$$\left[i\gamma ^\mu _\mu (m_q^0g_\sigma ^q\sigma )g_\omega ^q\omega \gamma ^0\right]\psi _q(\stackrel{}{r},t)=0,$$
(1)
where $`m_q^0`$ is the current quark mass and $`g_\sigma ^q`$ and $`g_\omega ^q`$ are the quark couplings with the $`\sigma `$ and $`\omega `$ mesons.

The normalized ground state for a quark (in an s-state) in the bag, which has a radius $`R`$, is given as

$$\psi _q(\stackrel{}{r},t)=N\mathrm{exp}(i\frac{ϵ_qt}{R})\left(\begin{array}{c}j_0(xr/R)\\ i\beta _q\stackrel{}{\sigma }\widehat{r}j_1(xr/R)\end{array}\right)\frac{\chi _q}{\sqrt{4\pi }},$$
(2)
where $`x`$ is the dimensionless quark momentum and the single particle quark energy, in units of $`R^1`$, is

$$ϵ_q=\mathrm{\Omega }_q+g_\omega ^q\omega R,\beta _q=\sqrt{\frac{\mathrm{\Omega }_qRm_q^{}}{\mathrm{\Omega }_q+Rm_q^{}}},$$
(3)
with $`\mathrm{\Omega }_q=(x^2+R^2m_{q}^{}{}_{}{}^{2})^{1/2}`$; $`m_q^{}=m_q^0g_\sigma ^q\sigma `$ is the effective quark mass, $`\chi _q`$ is the quark spinor and $`N`$ is the normalization factor.

The boundary condition at the bag surface is given by:

$$i\gamma n\psi _q=\psi _q.$$
(4)
This, for the ground state, reduces to

$$j_0(x)=\beta _qj_1(x)$$
(5)
which determines the eigen frequency, $`x`$, of this lowest mode in the medium. The form of the quark wavefunction in equation (2) is almost identical to that of the solution in free space. However the parameters in the expression have been substantially modified by the surrounding nuclear medium. Thus the quarks in the nucleon embedded in the nuclear medium are more relativistic than those in a free nucleon.

The energy of the nucleon bag is

$$M^{}=3\frac{\mathrm{\Omega }_q}{R}\frac{Z}{R}+\frac{4}{3}\pi R^3B,$$
(6)
where $`B`$ is the bag constant and $`Z`$ parametrizes the sum of the center-of-mass (c.m.) motion and the gluonic corrections. Note that this center-of-mass treatment is different from that of Jin and Jennings . The bag radius $`R`$ is that which minimizes the nucleon bag energy through

$$\frac{M^{}}{R}=0.$$
(7)
We now proceed to study the equation of state (EOS) for nuclear matter including the vacuum polarisation effects from nucleon and sigma mesons at zero temperature. The details of the theory have already been discussed in ref. . Only a few important steps are given here. The energy density after subtracting out the pure vacuum contribution then becomes

$$ϵ_0=ϵ_{MFT}+\mathrm{\Delta }ϵ,$$
(8)
with

$$ϵ_{MFT}=\frac{g}{(2\pi )^3}\underset{|𝐤|<k_F}{}𝑑𝐤(k^2+M_{}^{}{}_{}{}^{2})^{1/2}+\frac{1}{2}m_\sigma ^2\sigma _0^2+\frac{1}{2}m_\omega ^2\omega _0^2,$$
(9)
and

$$\mathrm{\Delta }ϵ=\frac{g}{(2\pi )^3}𝑑𝐤\left[(k^2+M_{}^{}{}_{}{}^{2})^{1/2}(k^2+M^2)^{1/2}\frac{g_\sigma \sigma _0M}{(k^2+M^2)^{1/2}}\right].$$
(10)
The above expression for the energy density is divergent. After renormalisation by adding the counter terms , we have the expression for the finite renormalised energy density,

$$ϵ_{ren}=ϵ_{MFT}+\mathrm{\Delta }ϵ_{ren},$$
(11)
where

$$\mathrm{\Delta }ϵ_{ren}=\frac{g}{16\pi ^2}\left[M_{}^{}{}_{}{}^{4}\mathrm{ln}\left(\frac{M^{}}{M}\right)+M^3(MM^{})\frac{7}{2}M^2(MM^{})^2+\frac{13}{3}M(MM^{})^3\frac{25}{12}M(MM^{})^4\right].$$
(12)
The baryon density is given by

$$\rho _B=\frac{gk_F^3}{6\pi ^2}$$
(13)
In the above, $`g`$ is the spin-isospin degeneracy factor which is equal to $`4`$ for nuclear matter and to $`2`$ for neutron matter.

Next, we consider the quantum corrections due to the scalar mesons. Including a quartic scalar self-interaction, the Hamiltonian density for the scalar mesons becomes

$$_\sigma =\frac{1}{2}_\mu \sigma ^\mu \sigma +\frac{1}{2}m_\sigma ^2\sigma ^2+\lambda \sigma ^4,$$
(14)
with $`m_\sigma `$ and $`\lambda `$ being the bare mass and coupling constant respectively. We calculate the expectation value of the Hamiltonian density and perform the renormalisation according to prescription of ref. . The resulting gap equation for $`M_\sigma ^2`$, which minimizes the energy, in terms of the renormalised parameters $`m_R^2`$ and $`\lambda _R`$ can be written as

$$M_\sigma ^2=m_R^2+12\lambda _R\sigma _0^2+12\lambda _RI_f(M_\sigma ),$$
(15)
where

$$I_f(M_\sigma )=\frac{M_\sigma ^2}{16\pi ^2}\mathrm{ln}\left(\frac{M_\sigma ^2}{\mu ^2}\right).$$
(16)
Using the above equations we obtain the energy density for the $`\sigma `$ in terms of $`\sigma _0`$ which is given by

$$ϵ_\sigma =3\lambda _R\left(\sigma _0^2+\frac{m_R^2}{12\lambda _R}\right)^2+\frac{M_\sigma ^4}{64\pi ^2}\left(\mathrm{ln}\left(\frac{M_\sigma ^2}{\mu ^2}\right)\frac{1}{2}\right)3\lambda _RI_f^22\lambda \sigma _0^4.$$
(17)
Where $`\mu `$ is a mass scale. The above expression is given in terms of $`\sigma `$ mass $`m_R`$ and $`\lambda _R`$ except for the last term which is still in terms of the bare coupling constant $`\lambda `$ and did not get renormalised because of the structure of the gap equation . However, from the renormalisation procedure it is easy to see in this last work that when $`\lambda _R`$ is kept fixed, the bare coupling $`\lambda 0_{}`$. Therefore the last term in eq. (17) will be neglected in the numerical calculations.

After subtracting the vacuum contribution we obtain:

$`\mathrm{\Delta }ϵ_\sigma `$ $`=`$ $`ϵ_\sigma ϵ_\sigma (\sigma _0=0)`$ (18)
$`=`$ $`{\displaystyle \frac{1}{2}}m_R^2\sigma _0^2+3\lambda _R\sigma _0^4+{\displaystyle \frac{M_\sigma ^4}{64\pi ^2}}\left(\mathrm{ln}\left({\displaystyle \frac{M_\sigma ^2}{\mu ^2}}\right){\displaystyle \frac{1}{2}}\right)3\lambda _RI_f^2{\displaystyle \frac{M_{\sigma ,0}^4}{64\pi ^2}}\left(\mathrm{ln}\left({\displaystyle \frac{M_{\sigma ,0}^2}{\mu ^2}}\right){\displaystyle \frac{1}{2}}\right)+3\lambda _RI_{f0}^2,`$ (19)
where $`M_{\sigma ,0}`$ and $`I_{f0}`$ are the expressions given by eqs. (15) and (16) with $`\sigma _0=0`$.

The energy density and pressure with baryon and the sigma condensate $`\sigma _0`$ are respectively done by

$$ϵ_{ren}=ϵ_0^{finite}+\mathrm{\Delta }ϵ_{ren},$$
(20)
and

$$P=\frac{g}{3(2\pi )^3}\underset{|𝐤|<k_F}{}𝑑𝐤\frac{k^2}{(k^2+M_{}^{}{}_{}{}^{2})^{1/2}}+\frac{1}{2}m_\omega ^2\omega _0^2\mathrm{\Delta }ϵ_\sigma \mathrm{\Delta }ϵ_{ren},$$
(21)
where

$$ϵ_0^{finite}=\frac{g}{(2\pi )^3}\underset{|𝐤|<k_F}{}𝑑𝐤(k^2+M_{}^{}{}_{}{}^{2})^{1/2}+\frac{1}{2}m_\omega ^2\omega _0^2+\mathrm{\Delta }ϵ_\sigma $$
(22)
with $`\mathrm{\Delta }ϵ_{ren}`$ given by eq. (12) and $`\mathrm{\Delta }ϵ_\sigma `$ by eq. (19).

The energy density from the $`\sigma `$ field as given by eq. (19) is still in terms of the renormalisation scale $`\mu `$ which is arbitrary. We choose this to be equal to the renormalised sigma mass $`m_R`$ in doing the numerical calculations. This is because changing $`\mu `$ would mean changing the quartic coupling $`\lambda _R`$, and $`\lambda _R`$ here enters as a parameter to be chosen to give the incompressibility for nuclear matter in the correct range. The parameters $`g_\sigma ^q`$ and $`g_\omega `$ are fitted so as to describe the ground-state properties of nuclear matter correctly. For a given baryon density $`\rho _B`$, the energy density, the density dependent radius of the nucleon and the nucleon effective mass are calculated at zero temperature.

## III Results and discussion

We now proceed with the numerical calculations for the nuclear matter. We start fixing the bag properties in the vacuum. They are given in Table 1. We next calculate the ground state properties of the nuclear matter and fit the scalar and vector coupling constants $`g_\sigma ^q`$ and $`g_\omega (=3g_\omega ^q)`$ to get the correct saturation properties for a given renormalised $`\sigma `$ mass and coupling, $`m_R`$ and $`\lambda _R`$. The omega and sigma couplings for given $`\lambda _R`$ are tabulated in Table II. Results from Relativistic Hartree approximation (RHA) are also shown. Little change is noted in the values of the parameters for the range of $`\lambda _R`$ which is considered here.

Using these values, we plot the binding energy $`(E_B=ϵ/\rho _BM_N)`$ for nuclear matter as a function of density in Figure 1. In the same figure we also plot the results for the relativistic Hartree approximation (RHA). Clearly, including baryon and $`\sigma `$-meson quantum corrections leads to a softer equation of state which is further softer for a higher value of $`\lambda _R`$. The equation of state, pressure, $`P`$ versus as a function of energy density $`ϵ`$ is displayed in Fig. 2 for the different cases. For comparison, the causal limit $`P=ϵ`$ is also shown in the figure. All the cases studied here respect the causal condition $`P/ϵ1`$, so that the speed of sound remains lower than the speed of light.

In figure 3 we plot the effective nucleon mass as a function of density. At the saturation density we get $`M^{}=0.817M`$ and $`0.83M`$ for $`\lambda _R=3.0`$ and $`\lambda _R=4.5`$ respectively. These values may be compared with the results of $`M^{}=0.775M`$ for normal QMC model and of $`0.793M`$ with the relativistic Hartree approximation in QMC model. This influence is much higher at high nuclear densities. We can conclude that quantum effects, at the level we consider, increase the effective mass $`M^{}`$.

We plot the in medium effective radius of the nucleon ($`R^{}`$) as a function of density in figure 4. $`R^{}`$ is also increased with relation to the mean field approximation mainly for higher densities. It is possible to understand this result as an effect which prevents QMC from deconfinement at high densities.

In figure 5, we plot the in-medium $`\sigma `$-meson mass, $`M_\sigma `$, as a function of density. $`M_\sigma `$ increases with density as $`\lambda _R`$ is positive. The higher is the coupling, higher is the increase of the sigma mass with density. This would be rather an indication of a further chiral symmetry breaking instead of its restoration .

To summarize, we have used a non-perturbative approach to include the quantum effects in nuclear matter in the framework of QMC. The calculations of the scalar meson quantum corrections was done here in a self-consistent manner including multi-loop effects. This leads to a softening of the equation of state. We have also calculated the effective mass of the $`\sigma `$-meson as modified by the quantum corrections. The effective sigma mass is seen to increase with density. The results may be suggesting that, at high densities, quantum fluctuations prevents the model from the chiral symmetry restoration as well as from deconfinement. These results will be extensively studied in another work where the coupling constant $`\lambda _R`$ and its influence on the chiral symmetry behavior -whose order parameter can be considered to be $`\sigma _0`$\- at high densities and temperatures will be shown in the frame of the model worked out here.

## IV Acknowledgments

P.K.P. would like to acknowledge FAPESP (Processo- 99/08544-0) for financial support and the IFT, São Paulo, for kind hospitality. F.L.B. acknowledges support from FAPESP (Processo- 97/01317-3).