format: text

text:
# The distribution of shock waves in driven supersonic turbulence

## 1 Introduction

Many structures we come across in the Universe have been shaped by fluid turbulence. In astronomy, we often observe high speed turbulence driven by supersonic ordered motions such as jets, supernova shocks, galactic rotation and stellar winds (e.g. Franco & Carramiñana 1999). Although it has long been thought that astrophysical turbulence provides the best opportunity to investigate supersonic turbulence, the lack of a theory has stunted our attempts to understand the behaviour (von Hoerner 1962). Three dimensional high resolution numerical simulations now provide a method to make real progress. Driven turbulence is explored in this paper and the results compared to a sister study of decaying turbulence (Smith et al. 2000, hereafter Paper 1). Our aim here is to relate the type of turbulence to the properties of the shock waves.

We study here compressible turbulence without thermal conduction. No physical viscosity is modelled, but numerical viscosity remains present, and an artificial viscosity determines the dissipation in regions of strong convergence. Periodic boundary conditions were chosen for the finite difference ZEUS code simulations, fully described by Mac Low et al. (1998). Uniform three dimensional turbulence with an isothermal equation of state, an initially uniform magnetic field and periodic boundary conditions are investigated. The influence of self-gravity is also examined.

Our motivation here is to provide the observer with a means of recognising the type of turbulence from the properties of the generated shock waves. The general energetics of decaying and driven hydrodynamic turbulence have already been computed by Mac Low et al (1998) and Mac Low (1999), respectively, using ZEUS, a second-order, Eulerian hydrocode. Mac Low et al (1999) concluded that, since turbulence which is left to decay dissipates rapidly under all conditions, the motions we observe in molecular clouds must be continuously driven. Klessen et al. (2000) extended the hydrodynamic results by calculating self-gravitating models. A smoothed particle hydrodynamics (SPH) code was also employed to confirm the results for both the decaying and driven cases (Mac Low et al 1998, Klessen et al. 2000). Klessen et al. (2000) have provided the parameter scaling for applications to molecular clouds. Heitsch et al. (2000) present the magnetohydrodynamic extension to the self-gravitating case, and discuss the criteria for gravitational collapse.

Our immediate target is to derive the spectrum of shocks (the Shock Probability Distribution Function) generated by driven turbulence. With this knowledge, we will proceed to predict the spectroscopic properties in a following work. Observed individual bright, sharp features, such as arcs, filaments and sheets, are often interpreted as shock layers within which particles are highly excited (e.g. Eislöffel et al. 2000). Where unresolved, the excitation can still be explored quantitatively by employing spectroscopic methods. The gas excitation then depends on both the physics of shocks as well as the distribution of shock strengths.

Previous studies of compressible turbulence have concentrated on the density and velocity structure of the cold gas rather than the shocks (e.g. Porter et al. 1994; Falgarone et al. 1994, Vázquez-Semadeni et al. 1996, Padoan et al. 1998). This may be appropriate for the interpretation of clouds since, although the Mach number is still high, the shock speeds are too low to produce bright features. The simulations analysed here are also being interpreted by Mac Low & Ossenkopf (2000) in terms of density structure.

The method used to count shocks from grid-based simulations was developed in Paper 1. The one-dimensional counting procedure was verified through a comparison with full three-dimensional integrations of the dissipated energy in Paper 1. This method is appropriate for a ZEUS-type finite difference code for which shock transitions are spread out over a few zones. Here, we also study the shock transitions in all three directions and display the spatial distributions for the energy dissipated through the artificial viscosity in the shocks. We first present the shock jump PDFs and provide analytical fits (Sect. 2). Magnetohydrodynamic simulations (Sect. 3) with Alfv́en numbers of $`A=\mathrm{\hspace{0.17em}5}`$ and $`A=\mathrm{\hspace{0.17em}1}`$ are then explored. Note the definition of the Alfvén Mach number $`A=v_{\mathrm{rms}}/v_A`$, where the Alfvén speed $`v_A=(B^2/4\pi \rho )^{1/2}`$, and $`v_{\mathrm{rms}}`$ is the initial root mean square (rms) velocity.

Which shocks actually dissipate most of the energy? The shock power rather than number density will determine the spectral region in which the system can be best observed. Hence, we determine the energy dissipated as a function of the shock speed in Sect. 4. A large-scale simulation which included self-gravity is then analysed in Sect. 5. We next interpret the results in terms of the dynamical models (Sect. 6).

## 2 Hydrodynamic Turbulence

### 2.1 Model description

We first explore driven hydrodynamic turbulence. The three-dimensional numerical simulations on grids with $`128^3`$ and $`256^3`$ zones and periodic boundary conditions were initialized with an rms Mach number of $`M=5`$. The initial density is uniform and the initial velocity perturbations were drawn from a Gaussian random field, as described by Mac Low (1999). The power spectrum of the perturbations is flat and limited to small wavenumber ranges, $`k_{\mathrm{min}}<k<k_{\mathrm{max}}`$, with $`k_{\mathrm{min}}=7`$ and $`k_{\mathrm{max}}=8`$ in Fig. 1. The uniform driver is a simple constant rate of energy input with the distribution pattern fixed.

The simulations employ a box of length $`2L`$ and a unit of speed, $`u`$. The gas is isothermal with a sound speed of $`c_s=0.1u`$. Hence the sound crossing time is 20 L/u. We take the time unit as $`\tau =L/u`$. Assuming a mass m to be contained initially in a cube of side $`L`$, the dissipated power is given in units of $`\dot{E}_o=mu^2/\tau `$.

### 2.2 Shock number distribution

We calculate the one dimensional shock jump function, as discussed and justified in Paper 1. This is the number distribution of the total jump in speed across each converging region along a specific direction. This is written as $`dN/dv_j`$ where $`v_j`$ is the sum of the (negative) velocity gradients (i.e. $`\mathrm{\Sigma }\left[\delta v_x\right]`$ across a region being compressed in the x-direction). We employ the jump Mach number in the x-direction $`M_j=v_j/c_s`$ rather than $`v_j`$ since this is the parameter relevant to the dynamics. Thus, each bounded region of convergence in the x-direction counts as a single shock and the total jump in $`M_j`$ across this region is related to its strength.

Numerically, we scan over the whole simulation grid (x,y,z), recording each shock jump through

$$M_j=\underset{x=x_i}{\overset{x=x_f}{}}(\mathrm{\Delta }v_x/c_s)$$
(1)
with the condition that $`\mathrm{\Delta }v_x<0`$ in the range $`x_ix<x_f`$. This is then binned as a single shock element. The shock number distribution $`dN/dM_j`$ is obviously dimensionless.

Note that the shock number depends on the grid size. We find that on doubling the grid size, the shock number increases by a factor of four. This implies that the total shock front area in units of $`L^2`$ remains roughly constant. This holds for both driven and decaying turbulence and is an indication that the numerical resolution is sufficient to capture the vast majority of shocks.

### 2.3 Steady state description

The random Gaussian field at t = 0 rapidly transforms into a shock field (Fig 1). As can be seen, the shock distribution approaches a steady state,

reached by t = 1. The driving energy determines the break in the distribution and the maximum shock speed, but does not influence the distribution of shock speeds.

The driving wavenumber influences the steady state as shown in Fig. 2.

There is a moderate dependence of the shock number on the wavenumber, especially at high Mach numbers. This may be partly due to the time involved for the regenerated low wavenumber modes to steepen. Some of the longer wave modes may be damped by interactions with the turbulent field.

A power law distribution of shock velocities is uncovered. Note that for high Mach numbers, statistically, we can equate the shock Mach number to the jump Mach number to a good approximation. We find, as shown by the indicated line in the lower box of Fig. 1, an inverse square-root law $`dN/dM_jM_j^{0.5}`$. In detail, we find a fit of the form

$$\frac{dN}{dM_j}\mathrm{1.4\hspace{0.17em}10}^4k^{0.5}\dot{E}^{0.2}M_j^{0.5}$$
(2)
over the power law sections. The break in the power law is found to occur at $`M_j^{max}`$ given by

$$M_j^{max}7.6k^{0.2}\dot{E}^{0.4}.$$
(3)
The k-dependence for $`M_j^{max}`$ has been estimated by inspecting the energy dissipation diagrams below. The indicated value is consistent with that discernable in Fig. 2.

Integrating Eq. (2), using the limit Eq. (3), yields the remarkable result that the number of shocks is a constant: $`N178,000`$ for $`k=\mathrm{\hspace{0.17em}8}`$. That is, when the energy input is low, there are many more weak shocks. As the energy is increased, the number of shocks does not increase. Rather, the shock strengths increase. Hence, the number of shocks reaches a saturation level of about 180,000 on the 1283 grids.

Hence, the total shock surface area is $`Nk^{0.4}`$. This is confirmed directly from the shock counts. We find that the number of jumps does depend weakly on the chosen value for the minimum convergence. We take the case $`dE/dt=10`$ for illustrative numbers. Taking all converging flow regions, yields 181,000 shocks with 5.17 zones per shock. This gives 46% of the volume occupied by converging regions. Taking instead a minimum convergence of $`0.5c_s`$ as a criteria which relates to a steepened wave balanced by artificial viscosity, yields 161,000 regions with an average of 3.95 zones. That is, almost 90% of the converging regions are indeed associated with steepened waves, which occupy a constant 30% of the volume. This number density was also found in decaying turbulence: even though the shocks may weaken and interact, the total number is conserved, and the total shock area is independent of the grid size (Paper 1).

Many predictions and analytical fits to simulations have been published for the velocity gradients within turbulent flows. Here, for supersonic turbulence, we find a Gaussian high-speed tail to the shock distribution provides a rough fit. Figure 3

displays a combined power-law and Gaussian fit of the form

$$\frac{dN}{dM_j}\mathrm{1.4\hspace{0.17em}10}^4k^{0.5}\dot{E}^{0.2}M_j^{0.45}\mathrm{exp}\left[(M_j/15.1)^2\right]$$
(4)
for one case. To test if this is due to the initial Gaussian field we have also run a simulation with an exponential driver and found that a similar tail is still present. Hence the near-Gaussian could be due to the superposition of the initially randomly-located forcing waves, as expected from the Central Limit Theorem. The best fit we find , however, involves a tail of the form $`\mathrm{exp}(M^3)`$,

$$\frac{dN}{dM_j}\mathrm{1.4\hspace{0.17em}10}^4k^{0.5}\dot{E}^{0.2}M_j^{0.45}\mathrm{exp}\left[(M_j/16.6)^3\right].$$
(5)
Such steep tails have indeed been commonly found despite Gaussian driving, thought to be due to the more rapid decay of the faster shocks (e.g. Gotoh & Kraichnan 1998).

## 3 Magnetohydrodynamic Turbulence

We introduce a uniform magnetic field into the initial configuration. A weak field, in which the Alfvén and sound speeds are equal, has little overall influence (Fig. 4). The shock distribution is roughly isotropic and the shock number is not significantly altered from the equivalent hydrodynamic simulation. The power law section is not so well defined for the shock distribution parallel to the field.

A strong field introduces a strong anisotropy (Fig. 5). With a field such that the Alfvén speed equals the initial rms speed, the waves transverse to the field dominate. There are $`23`$ times more waves in the transverse direction for a given jump speed in each direction. The inverse square root power-law rule is again closely obeyed. We call these waves, rather than shocks, since the high Alfvén speed implies that a high fraction may be fast magnetosonic waves. The average number of zones, however, measured for each jump is only 4.7 (parallel) and 7.4 (transverse). This compares to an average of 5.4 for the hydrodynamic flow (and 23.7 zones for the initial Gaussian with k=4). Note these refer to the complete shock, not just the 2-3 zones across which the jump is highly non-linear and across which numerical viscosity is strong); there usually exists one or two zones on each side of the main jump across which the velocity joins smoothly onto the surrounding flow without oscillations). Hence, the zone measurements favour the interpretation that the dissipation is being carried out in short-wavelength non-linear magnetosonic waves.

A very similar difference between parallel and transverse shock numbers is found in the case of decaying turbulence (Paper 1). It is clear that the magnetic pressure is not damping the shock waves. In contrast, the extra magnetic energy which becomes tied up in the waves also helps maintain them.

## 4 Energy dissipation

The increasing number of strong shocks with increasing wavenumber is qualitatively consistent with the finding of Mac Low (1999) that the total energy dissipation rate increases with the driving wavenumber. To discuss the energy dissipation rate, however, we must first also consider the density $`\rho `$ within each shock. We have then calculated the power dissipated by artificial viscosity within each shock front, as described in Paper 1. This yields the power dissipated per unit shock speed. We actually calculate in this section the component of the dissipation along the x-axis and the jump Mach number along the x-axis which, given the statistics, represents the true three dimensional ‘power distribution function’. Numerically, over the whole simulation grid (x,y,z), we identify each compact range $`x_ix<x_f`$ along the x-axis for which $`\mathrm{\Delta }v_x=v_{x+1}v_x<0`$. For each jump we then find the energy dissipated by artificial viscosity as

$$P_j=\frac{C}{L^2}\underset{x=x_i}{\overset{x=x_f}{}}\rho _x(\mathrm{\Delta }v_x)^3,$$
(6)
where C measures the number of zones over which artificial viscosity will spread a shock. This is then binned as a single shock element. The shock number distribution $`dP/dM_j`$ is multiplied by three to account for the energy dissipated in each of the three dimensions. Further details of the method can be found in Paper 1, where its reliability was also verified.

The surface brightness of the shocks, as well as the column density of all the gas within the cube, is displayed in Fig. 6 for a pure

hydrodynamic driven example. Note that these are not slices, but we have integrated through the z-direction. Much of the gas has been swept up into a flattened cloud. Many shocks appear sharper when only the x-direction is accounted for since this emphasizes the shocks transverse to the line of sight.

Dissipating shocks and the density distribution do not correlate very well, although the main cloud contains the main elongated region of dissipation. The most powerful shocks are mainly within the cloud.

We find that the energy dissipation rate also takes on a power law $`M_j`$ dependence (Fig. 7). In contrast to the number

distribution, it is independent of the driving power input. There is also a wavenumber dependence, as shown in Fig. 8.

The total energy dissipated in the power law section (summing the losses for each direction) can be written

$$\frac{dP}{dM_j}\mathrm{1.05\hspace{0.17em}10}^2k^{0.5}M_j^{1.5}.$$
(7)
The error in the power law index of 1.5 is $`0.12`$. The break in the power law is found to occur at $`M_j^{max}`$ as given by Eq. (3). Integrating Eq. (7), we obtain Eq. (3). Thus we have acquired a self-consistent mathematical description.

## 5 Gravitational Collapse

Can we distinguish a collapsing self-gravitating cloud from a turbulent region within which the Jeans mass is not reached? We here take a simulation of driven hydrodynamic turbulence in which self-gravity is switched on after a steady driven state has been reached (Model D2 from Klessen et al. 2000). The particular parameters chosen are shown in the caption to Fig. 9.

The number and distribution of shock speeds is not particularly different from the equivalent non-gravity numerical experiment. The power-law section is limited, consistent with other low wavenumber simulations (see Fig. 2).

In stark contrast, the dissipated energy transforms from the typical non-gravity case (upper box of Fig. 10) to one dominated by a few accretion shocks (lower box of Fig. 10). A

relatively small number of well-defined shocks generate strong individual signatures. These shocks are of typical cloud speeds but are strongly dissipative because they propagate into very dense regions.

The distribution of the shocks at time t = 6 is displayed in Fig. 11. The elongated cloud at time t = 2 (Fig. 6) has now collapsed into narrow filaments and cores. There is high shock dissipation around the collapsing cores although other strong shocks are still quite widespread. The code does not distinguish between a shock and a collapsing flow. The artificial viscosity prescription treats all converging flow regions as dissipative. Nevertheless, a collapsing flow does physically dissipate its energy as it shocks onto a growing core. Hence the dynamical evolution is correctly modelled. The surface brightness of the shocks is, of course, limited by the grid resolution. The high SPH resolution, however, verifies that the general properties are realistic (Klessen et al. 2000). A new problem arising in the shock analysis is that a single converging flow onto a core may actually be a double shock layer. We would here record this as a single shock in the counting procedure and in Fig. 10. We have not tried to correct for this.

## 6 Interpretation

Turbulence has proven easy to interpret, but the interpretations have been hard to prove. Our aim here is limited to answer: why a universal inverse square root power law? Supersonic turbulence possesses simplifying characteristics which we here employ to understand the shock distribution.

First, it is significant that the flow is dominated by the integrated number and power of the high-speed shocks at the knee in the shock jump distributions. These shocks are generated by the driver and propagate at high speed into the fluid. In contrast, a shock from the power-law section plays a passive rôle: when overtaken by a strong shock, its strength decays according to the speed and density of the oncoming material. Hence the mean lifetime, $`\tau `$, of a shock layer, and hence of the shock momentum, is a constant, independent of its own absolute speed, $`v`$. Expressed mathematically, for a steady state, this yields

$$\frac{}{v}\left[\frac{N}{v}\frac{1}{\tau }\right]=0,$$
(8)
from which we find the distribution of absolute shock speeds (i.e. the speed of the layers within the box) is $`N/v=constant`$. This is confirmed from the simulations - Fig. 12 shows that, below

the high-speed knee, the distribution of absolute velocities within the shocks is indeed independent of the absolute speed.

The relationship between the absolute speed, $`v`$, and the jump speed, $`v_j`$, for the ensemble of shocks is expected to resemble that of Burger’s turbulence (Paper 1, see also Gotoh & Kraichnan 1993). The steepening of these independently-moving shocks is proportional to the absolute speed simply because higher amplitude waves produce stronger shocks, whereas low amplitude waves can only lead to weak jumps. Secondly, the higher amplitudes imply that the velocity gradients are also steeper. Hence the material swept up into a shock is from a wider region which possesses a larger velocity difference. These two linear effects lead to a relationship for the ensemble shock jump and absolute speeds of the form $`v_jv^2`$. This is the essence of mapping closure theory, which then uses this relationship to transform initial Gaussian distributions into exponentials (Chen et al 1989, Kraichnan 1990). Here, we transform the distribution of absolute shock speeds generated by the uniform driving, to obtain

$$\frac{N}{v_j}v_j^{1/2}.$$
(9)
This agrees with the numerical experiments. Note that the result depends on the nature of the driving: uniform driving produces strong shocks locally which then propagate through the ambient medium. In contrast, a white-noise driver would continually create weak small-scale turbulence throughout the space and might thus produce, for example, some other power-law behaviour.

We found here that the high-velocity tail is somewhat steeper than Gaussian. A similar result has been obtained for other types of turbulence, with the same form occurring, $`\mathrm{exp}(v_j^3)`$, as uncovered here (see Gotoh & Kraichnen 1998). This is due to the Gaussian forcing being modified by the high dissipation rates at these speeds.

The total energy ‘radiated’ in the shocks can be found on integrating Eq. (7) over the jump Mach Number. This yields $`P=0.669\dot{E}`$. The remaining injected energy is also lost in the shocks where numerical diffusion is, of course, strong. As shown in Table 1 of Paper 1, we expect approximately 0.68 of the energy to be radiated in the shocks. Hence, the description of driven turbulence is fully consistent.

## 7 Conclusions

We have analysed numerical simulations of uniformly driven supersonic, magnetohydrodynamic and self-gravitating turbulence. Below a critical jump speed, we find a power law distribution of fast shocks with the number of shocks inversely proportional to the square root of the shock jump speed. Hence, unlike the decaying case, the driven case possesses an ‘inertial range’ of shock strengths. This range is, however, dynamically passive, being mediated by strong shocks injected at a higher knee in the distribution. The knee is mainly responsible for the dissipation of energy, and thus the power-law range may only have a weak observational signature. This will be explored in the following paper of this series.

These results contrast with the exponential distribution and slow shock dissipation associated with decaying turbulence described in Paper 1.

A strong magnetic field does not alter the shape of the shock number distributions. It enhances the shock number transverse to the field direction at the expense of parallel shocks. The distribution of parallel shocks is, however, extended to higher jump speeds. Transverse waves thus appear to find support in the magnetic pressure.

A simulation with self-gravity demonstrates the development of a number of highly dissipative accretion shocks. These shocks are rare and are not apparent in the number distributions. The shock number distribution demonstrates a gradual curvature rather than a power-law. This implies that high-speed shocks can be maintained by gravitational acceleration.

Finally, we have shown how the power-law behaviour may arise. No substantial theory has been developed to predict this inverse square root law. The power law is closely related to the distribution in absolute shock speeds, which is shown to be flat. The particular power law is predicted to depend on the type of driver.

The energy is dissipated in shocks mainly associated with the denser regions of the clouds formed by the large scale wave motions. Numerous other shocked regions are scattered throughout. Very strong shocks are associated with the collapsing cores in the self-gravitating case.

Our results must be considered within the context of the imposed models. While turbulence has been pinned down in many studies, further studies have simply extended the number of pins required. Extensions to this work include: (1) the shock distributions corresponding to non-uniform drivers, (2) the spatial propagation of turbulent energy and (3) the observable signatures: the spectral energy distributions derived from the shock distributions,

###### Acknowledgements.

MDS benefitted greatly from the hospitality of the Max-Planck-Institut für Astronomie. We thank R. Klessen and E. Zweibel for advice and discussions. Computations were performed at the MPG Rechenzentrum Garching. M-MML holds a CAREER fellowship from the US National Science Foundation, grant number AST 99-85392.


text:
# Metal-insulator transitions in cyclotron resonance of periodic nanostructures due to avoided band crossings

## I Introduction

Since the pioneering work of Bloch  and Landau  the investigation of spectral properties and transport features of Bloch electrons in magnetic fields  have attracted much attention. While infinitely degenerate Landau levels are found for free electrons in a static magnetic field, a spatially periodic potential without field yields a band structure in the spectrum. The general situation of a periodic potential and a magnetic field is still a matter of intensive research. A particular simple realization for a two-dimensional (2D) periodic potential is, e.g., given by

$$V(𝐫)=V_x\mathrm{cos}\left(\frac{2\pi x}{a}\right)+V_y\mathrm{cos}\left(\frac{2\pi y}{b}\right),$$
(1)
where $`a,b`$ are the periods in the corresponding directions. In a one-band approximation this potential leads to the well-known Harper model , where a metal–insulator transition (MIT) is observed upon changing the modulation amplitudes from $`V_y>V_x`$ to $`V_x>V_y`$ . At the critical point, i.e., $`V_x=V_y`$, the spectrum and the eigenfunctions are multifractals  and give rise to anomalous diffusion . The one-band Harper model, however, has an integrable classical limit, whereas Bloch electrons in a magnetic field show chaotic dynamics in the classical limit. In quantum mechanics this is reflected in the coupling of Landau bands. This coupling leads to avoided band crossings, which induce an unusual transport phenomenon: For $`V_xV_y`$ electrons show ballistic transport in the direction of strong potential modulation and localization in the direction of weaker modulation .

Analogous phenomena can be found in driven systems with a chaotic classical limit even within the one-band approximation. A well studied example is the kicked Harper model (KHM) . For small kicking strength its behavior is analogous to the Harper model , while for increasing kicking strength avoided band crossings due to the classical non-integrability induce a variety of MITs . Recently, it was demonstrated that a variety of kicked models, including the KHM, can be realized for electrons on a lattice in the presence of a magnetic field driven by a smooth electric field . The effective kicking is a result of a resonant interaction between the electronic motion and the driving field. These models correspond to cyclotron resonance experiments in antidot arrays and in organic metals . By investigating the dynamics, a localization-delocalization transition induced by the amplitude of the electric field was observed  and remained unexplained so far.

In this paper we investigate these dynamical properties in terms of spectral properties in a wide range of driving strengths. We show that the previously found localization-delocalization transition is a consequence of avoided band crossings due to the change of the driving field strength. Moreover, the spectral analysis predicts a cross-like ballistic transport, which is numerically confirmed. We start in Section II by briefly introducing the model and explain the numerical approach in Section III. In Section IV we present our main results, which give a clear picture for the MIT in terms of the spectrum. We finally summarize our findings in Section V.

## II The model

The one particle Hamiltonian, describing the cyclotron resonance of an electron in the 2D periodic potential $`V(𝐫)`$ of Eq. (1) can be written in the following form:

$$=\frac{1}{2m^{}}\left(𝐩e𝐀(𝐫,t)\right)^2+V(𝐫),$$
(2)
where $`𝐩=(p_x,p_y)`$ is the momentum of an electron with effective mass $`m^{}`$ and charge $`e`$ moving in the $`xy`$-plane perpendicular to the constant magnetic field $`B`$ pointing in $`z`$-direction. We choose the Landau gauge, which yields a vector potential $`𝐀(𝐫,t)=(\frac{E_x}{\nu }\mathrm{sin}\nu t,xB\frac{E_y}{\nu }\mathrm{cos}\nu t)`$ containing also the alternating electric field $`𝐄=(E_x\mathrm{cos}\nu t,E_y\mathrm{sin}\nu t)`$ with frequency $`\nu `$.

Applying a gauge transformation to the Hamiltonian (2), the time dependence can be shifted to the periodic potential $`V(𝐫)V(𝐫,t)`$ : If the fields $`B,E_x,E_y`$ satisfy $`\omega ^{}E_y=\nu E_x`$ then one can rewrite the Hamiltonian as $`=_0(𝐩,x)+V(𝐫,t)`$, where $`_0=_0(𝐩,x)`$ describes the cyclotron motion with frequency $`\omega ^{}=eB/m^{}`$. The perturbation $`V(𝐫,t)`$ splits the infinitely degenerate Landau levels of $`_0`$ into a mini-band structure. For strong electric fields and $`\nu =1`$ this Hamiltonian $``$ reduces to a kicked system

$`\stackrel{~}{}=L\mathrm{cos}p+K{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}`$ $`\{\mathrm{cos}(q\kappa _0)\delta (t2n+1)`$ (3)
$`+`$ $`\mathrm{cos}(q+\kappa _0)\delta (t2n)\},`$ (4)
where the canonical operators $`q`$ and $`p`$ corresponding to the spatial coordinates $`x,y`$ satisfy $`[p,q]=i2\pi \mathrm{}_{\mathrm{eff}}`$ and $`\kappa _0=(2\pi eE_y/(bm^{})\pi /4)\mathrm{mod}2\pi `$; $`L,K`$ denote the rescaled amplitudes of the potential (1). The effective Planck constant $`\mathrm{}_{\mathrm{eff}}=\mathrm{}/(abeB)`$, corresponds to the inverse number of magnetic flux quanta through a unit cell $`a\times b`$. The system described by Eq. (3) has a chaotic classical limit  and for $`\kappa _0=0`$ reduces to the KHM.

## III Numerical Approach

The time evolution operator for one period corresponding to (3) reads

$`𝒰`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{iL}{\mathrm{}_{\mathrm{eff}}}}\mathrm{cos}p\right)\mathrm{exp}\left({\displaystyle \frac{iK}{\mathrm{}_{\mathrm{eff}}}}\mathrm{cos}(q\kappa _0)\right)`$ (6)
$`\times \mathrm{exp}\left({\displaystyle \frac{iL}{\mathrm{}_{\mathrm{eff}}}}\mathrm{cos}p\right)\mathrm{exp}\left({\displaystyle \frac{iK}{\mathrm{}_{\mathrm{eff}}}}\mathrm{cos}(q+\kappa _0)\right).`$
In order to obtain its spectrum, one has to approximate the irrational effective Planck constant by a rational approximant, $`\mathrm{}_{\mathrm{eff}}=2\pi M/N`$. Following standard techniques for kicked systems , this leads to a $`N\times N`$-matrix for the evolution operator $`𝒰`$, which depends on two Bloch phases $`\theta _q,\theta _p[0,2\pi ]`$. These Bloch phases arise from periodic boundary conditions in the $`q`$\- and $`p`$-direction, respectively. Diagonalization of $`𝒰`$ for all pairs $`(\theta _q,\theta _p)`$ gives its spectrum. Keeping $`\theta _q`$ fixed and varying $`\theta _p`$ yields a part of this spectrum, from now on referred to as the “$`p`$-spectrum”. Eigenfunctions, which are localized in $`p`$-direction will be associated with very narrow bands in the $`p`$-spectrum, whereas extended states give rise to broad bands. Since for higher approximants of $`\mathrm{}_{\mathrm{eff}}`$ the width of the narrow bands decreases exponentially, we will call them levels. In an analogous way one can define the “$`q`$-spectrum”, which reflects the properties of the eigenfunctions in $`q`$-direction. Figure 1 shows a $`q`$\- and $`p`$-spectrum, where one observes broad bands and levels.

In order to investigate dynamics, we will study the spreading of wave packets, initially $`\delta `$-localized in either $`q`$\- or $`p`$-direction by determining the corresponding variances $`M_q(t)`$ and $`M_p(t)`$. For the spreading in $`p`$-direction we use periodic boundary conditions in $`q`$-direction after one unit cell, giving a discrete lattice in momentum space with spacing $`\mathrm{}_{\mathrm{eff}}`$. An analogous construction is used for the spreading in $`q`$-direction. When considering transport in phase space in both directions, we use periodic boundary conditions after $`r`$ unit cells in $`q`$-direction, which gives $`r`$ different Bloch phases $`\theta _q=2\pi j/r`$ ($`j=0,\mathrm{},r1`$). The complete wave packet is then obtained by combining the wave packets for each Bloch phase and it is defined on a grid with spacing $`\mathrm{}_{\mathrm{eff}}/r`$ in momentum space.

## IV Results

We will study the spectrum of Eq. (6) as a function of $`\kappa _0`$ and its consequences for dynamical properties. In particular, we will focus on the parameters used in Ref.  and show, that they reflect the generic behavior of the model (3).

In Ref.  for $`K=L=5and\mathrm{}_{\mathrm{eff}}=2\pi /(7+\sigma _g)`$, where $`\sigma _g=(\sqrt{5}1)/2`$ denotes the golden mean, a transition from diffusive to ballistic dynamics was found, when $`\kappa _0`$ was varied from $`0`$ to $`\pi \sigma _g`$ (cf. Fig. 1 of Ref. ). In Fig. 1 we present the $`q`$\- and $`p`$-spectrum for one period of $`\kappa _0`$. Apparently, the duality of $`q`$\- and $`p`$-spectrum is conserved for all $`\kappa _0`$, i.e., the $`q`$-spectrum shows wide bands, where the $`p`$-spectrum is level-like and vice versa.

For $`\kappa _0=0`$ the variances increase almost linearly in both $`q`$\- and $`p`$-direction, corresponding to diffusive spreading in either direction (top of Fig. 2a). This well-known behavior of the KHM is a direct consequence of the multifractal nature of spectrum and eigenfunctions.

For $`\kappa _0=\pi \sigma _g`$ bands occur in the $`p`$-spectrum explaining the ballistic transport in $`p`$-direction observed in Ref. . There are, however, also levels, which are related to bands in the $`q`$-spectrum, predicting a ballistic spreading in $`q`$-direction. This is confirmed in Fig. 2b. All eigenstates are extended in either of the directions and localized in the other one. A typical initial wave packet will excite both types of eigenfunctions. The part of the wave packet consiting of eigenfunctions extended in the $`q`$-direction gives ballistic transport in the $`q`$-direction. The other part, which is a superposition of eigenfunctions extended in the $`p`$-direction, yields ballistic transport in the $`p`$-direction. Therefore one finds a superposition of ballistic transport in $`q`$\- and $`p`$-direction, namely a cross-like transport (bottom of Fig. 2b). The cross-like form of the wave packet is clearly different from the isotropic form for $`\kappa _0=0`$ (Fig. 2a). There are several other transport features one can observe from the Husimi plot in Fig. 2b. The weight of the wave packet spreading ballistically along the $`p`$-direction is much bigger than the corresponding weight for the $`q`$-direction. This is a signature of the spectrum, namely that there are more bands in the $`p`$-spectrum than there are bands in the $`q`$-spectrum at $`\kappa _0=\pi \sigma _g`$, as can be seen in Fig. 1. It should be noted, that the velocities in $`p`$\- and $`q`$-direction (which can be inferred from the variances $`M_p(t)`$ and $`M_q(t)`$) are in general not related to the number of bands but are determined by the band widths.

Now we will turn away from the symmetry line $`K=L`$ and focus on the parameters $`K=2,L=4,\mathrm{}_{\mathrm{eff}}=2\pi /(7+\sigma _g)`$. In Fig. 3 we show the $`q`$\- and $`p`$-spectrum as a function of $`\kappa _0`$, which again appears to be dual. For the KHM ($`\kappa _0=0`$) the $`q`$-spectrum shows broad bands, whereas the $`p`$-spectrum is level-like, indicating transport in $`q`$-direction, but localization in $`p`$-direction. This is confirmed by the dynamics (Fig. 4a). Variation of $`\kappa _0`$ induces avoided band crossings leading to the occurrence of bands in the $`p`$-spectrum (see magnification in Fig. 3). This is associated with a localization-delocalization transition in $`p`$-direction as observed in Ref. . As there are still bands in the $`q`$-spectrum one again finds the generic cross-like transport (Fig. 4b). Here, the weight of the wave packet spreading along the $`q`$-direction is bigger than the corresponding weight in $`p`$-direction, as is implied by the much larger number of bands in the $`q`$-spectrum (Fig. 3). Another characteristic of the wave packet is the increased localization length in $`p`$-direction of the part of the wave packet spreading in $`q`$-direction compared to Fig. 4a. This is due to the increased localization length of the eigenstates corresponding to the levels in the $`p`$-spectrum. This change in the localization length can be explained by avoided band crossings.

## V Summary

In conclusion, we have shown here, that the previously found MIT in cyclotron resonance in a 2D-periodic potential in the presence of a magnetic field can be fully understood by investigating the spectrum of the corresponding time evolution operator. We have presented numerical evidence, that the duality of the spectra is conserved under changes of the driving strength. We were able to explain the previously observed dynamics , which has shown a transition from diffusive to ballistic transport as well as a transition from localization to ballistic transport in the $`p`$-direction. Furthermore, the spectra allowed us to predict the transport behavior in $`q`$-direction leading to a cross-like transport. These predictions were also confirmed numerically. We would like to mention, that cross-like transport is generic for systems with dual spectra and a chaotic classical limit.

The observed transport phenomena in phase space of the effective Hamiltonian (3) correspond to analogous transport phenomena in the $`xy`$-plane. The required conditions for realization of the observed MIT may well be achieved in experiments on cyclotron resonance  of a 2D electron gas embedded in lateral superlattices fabricated on GaAs heterostructures .

This research was supported by the Niedersachsen Ministry of Science. A.I. and S.F. were also supported by the Israel Science Foundation founded by the Israel Academy of Sciences and Humanities. A.I. thanks Prof. T. Geisel for his hospitality at the MPI für Strömungsforschung.


text:
# Regularizing the Divergent Structure of the Light-Front Currents

## I Introduction

With the recent advances in the Hamiltonian renormalization program, Light-Front Dynamics (LFD) appears to be a promising technique to impose the relativistic treatment of hadrons. In LFD a Fock-space expansion of bound states is made to handle the relativistic many-body effects in a consistent way . The wave function $`\psi _n(x_i,\stackrel{}{k}_i,\lambda _i)`$ describes the component with $`n`$ constituents, with longitudinal momentum fraction $`x_i`$, perpendicular momentum $`\stackrel{}{k}_i`$ and helicity $`\lambda _i`$, $`i=1,\mathrm{},n`$. It is the aim of LFD to determine those wave functions and use them in conjunction with hard scattering amplitudes to describe the properties of hadrons and their response to electroweak probes. Important steps were taken towards a realization of this goal. However, at present there are no realistic results available for wave functions of hadrons based on QCD alone. In order to calculate the response of hadrons to external probes, one might resort to the use of model wave functions. The variational principle enabled the solution of a QCD-motivated effective Hamiltonian and the constructed LF quark-model provided a good description of the available experimental data spanning various meson properties . The same reasons that make LFD so attractive to solve bound-state problems in field theory make it also useful for a relativistic description of nuclear systems. Presently, it is realized that a parametrization of nuclear reactions in terms of non-relativistic wave functions must fail. LF methods have the advantage that they are formally similar to time-ordered many-body theories, yet provide relativistically invariant observables. Furthermore, as far as the concerned amplitude is unconditionally (or absolutely) convergent, the LF Hamiltonian approach must yield the same result as the one obtained by a covariant Feynman approach.

However, not all is well. As we have recently shown, the amplitudes that are only conditionally convergent must be treated with care. A case in point discussed in our work was the calculation of a current matrix element in quantum field theory. A typical amplitude is given by the triangle diagram. One encounters this diagram e.g. when computing the meson form factor (see Fig. 1(a)).

The vertices denoted by $`\mathrm{\Phi }`$ are coupling constants in covariant perturbation theory. The hard scattering process is the absorption of a photon of momentum $`q`$ by a (anti-)quark. In the LFD approach the covariant amplitude is replaced by a series of LF time-ordered diagrams. In the case of the triangle diagram they are depicted in Figs. 1(b) and 1(c). The first (Fig. 1(b)) of these two diagrams is easily interpreted in terms of the LF wave functions $`\mathrm{\Psi }`$. However, the other diagram (Fig. 1(c)) has one vertex that can again be written in the same way as before, but it contains also another vertex, denoted by $`\mathrm{\Psi }^{}`$, that cannot be written as a LF wave function. The necessity of this new element $`\mathrm{\Psi }^{}`$ in LFD has also been discussed in the cases of semileptonic meson decays and deeply virtual Compton scattering. One may call $`\mathrm{\Psi }`$ and $`\mathrm{\Psi }^{}`$ the vertices of leading order and non-leading order, respectively, in the sense that the leading order vertex corresponds to the lowest Fock-state whereas the non-leading order vertex takes into account the higher Fock-states. The diagrams with $`\mathrm{\Psi }^{}`$ are designated as nonvalence diagrams while those with vertices of type $`\mathrm{\Psi }`$ only are designated as valence diagrams. In order to obtain the invariant form factor, the two LF form factors must be added. It depends on the situation whether one can limit oneself to a single component of the current $`J^\mu `$ to extract the invariant objects. For the electromagnetic current of a spin-0 particle any single component would suffice to extract the unique form factor. On the other hand, in situations like semileptonic pseudoscalar meson decay, which involves two independent form factors, or the electromagnetic current of a (axial-)vector particle that is described by three independent form factors, one must use information from several current components to determine the invariant amplitudes.

Earlier, we presented an analysis of contributions from the nonvalence diagrams. We constructed both leading and non-leading order vertices using pointlike covariant ones. The model that we used was essentially an extension of Mankiewicz and Sawicki’s $`(1+1)`$-dimensional quantum field theory model , which was later reinvestigated by several others . While their model was a simple 1+1 dim. scalar field theory, it included a binding effect of the two-body bound-state. Indeed, in Ref. , the relativistic two-body bound-state form factor was discussed in the full range of the binding energy. The starting model wave function was the solution of the covariant Bethe-Salpeter equation in the ladder approximation with a relativistic version of the contact interaction . The covariant model wave function was a product of two free single particle propagators, the overall momentum-conserving Dirac delta function, and a constant vertex function. Consequently, all our form factor calculations were various ways of evaluating the Feynman triangle diagram in quantum field theory. As pointed out in Ref. , however, the elastic electromagnetic form factors of a bound-state computed from the triangle diagram and from the Hamiltonian front-form dynamics are the same. Since our aim was to analyze the taxonomy of the triangle diagram, we didn’t choose any particular gauge for the electromagnetic gauge field but presented the equivalence of the physical form factor $`F(q^2)`$ in any choice of the electromagnetic gauge ***For a recent advocacy of the anti-light-cone gauge, $`A^{}=0`$, see M. Morara, R. Soldati, and G. McCartor, AIP Conference Proceedings $`\mathrm{𝟒𝟗𝟒}`$, 284-290 (1999). or in any choice of the current component.

Our conditionally convergent example of LFD with a fermion-loop showed that the bad component of the current, $`J^{}`$, with spin-1/2 constituents exhibits a persistent end-point singularity in the contribution from the nonvalence diagram. However, the calculation carried out so far was semi-realistic as the model was 1+1 dimensional and only a point-vertex was considered.

In the present work, we extend our analysis of the fermion-loop to the case of 3+1 dimensions. In 3+1 dimensions both the covariant and the LF calculations are divergent and the model without any smeared vertex for the fermion loop is not well defined. In the recent literature, the fermion-loop was regulated by smearing the $`q\overline{q}`$ bound-state vertex in a covariant manner. However, the vertex function was not symmetric in the four momenta of the constituent quarks and could hardly be considered a realistic approximation of a $`q\overline{q}`$ bound state. It was regarded only as a convenient cutoff prescription which makes the one-loop integrals finite. Furthermore, the calculation of the meson decay constant reveals that the end-point singularity is not completely cancelled by such an asymmetric choice of the vertex function. (See the next section for more details.) We show in this work that the fermion-loop can also be regulated by taking a non-local gauge-boson vertex. With this method satisfying the Ward-Takahashi (W-T) identity, we found the complete cancellation of the end-point singularity not only in the electromagnetic form factor but also in the decay constant. The non-local gauge-boson vertex remedies also the conceptual difficulty associated with the asymmetric way of treating $`q`$ and $`\overline{q}`$ in the previous calculations.

Nevertheless, one should distinguish the bound-state from the confined-state. In this work, we are treating the mesons as bound-states rather than confined-states, because we do not yet know how to make a covariant regularization for the confined-state. Thus, our emphasis here is the inclusion of the non-leading order vertex rather than the model-building of a realistic meson wavefunction.

We have performed the LF calculation in parallel to the covariant Feynman calculation. Our light-front results entail the similarity between the vertex-smearing-technique (either for the bound-state vertex or the gauge-boson vertex) and the Pauli-Villars regularization. For the bosonic loop calculation, the two methods turn out to be identical. However, for the fermionic loop calculation, the vertex-smearing-technique shares only the same structure of the denominators with the Pauli-Villars regularization. Using the gauge-boson vertex-smearing technique, we found that the persistent end-point singularity is removed even if the smeared-vertex is taken to the limit of the point-vertex.

A significant entity of our work is the taxonomy of valence and non-valence contributions substantiating the finiteness of each contribution when the gauge-boson vertex is regulated. Our results satisfy current conservation. However, we note that each contribution individually depends on the reference-frame even though the sum is always frame-independent. We thus elaborate the frame-dependence of individual contributions. Also, the zero-mode contribution should be distinguished from the instantaneous contribution. For the numerical estimates of physical observables, we present the electromagnetic form factors of the $`\pi ,K`$, and $`D`$ mesons.

In the next Section (Section II), we present both the covariant Feynman calculations and the LF calculations using the LF energy integration for the electromagnetic form factors of a pseudoscalar meson with spin-1/2 constituents. Section III contains the numerical taxonomy of both the valence and nonvalence diagrams to the electromagnetic form factors of $`\pi ,K`$ and $`D`$-mesons. The conclusion and discussion follow in Section IV. The general formula including unequal mass cases are summarized in Appendix A and the analytic behavior of the valence and nonvalence contributions at $`Q^20`$ limit in various frames is summarized in Appendix B.

## II Calculations

The electromagnetic form factors can be extracted from the matrix elements of the current $`J^\mu `$

$$p^{}|J^\mu |p=ie_m(p^\mu +p^\mu )F(q^2),$$
(1)
where $`e_m`$ is the charge of the meson and $`q^2=(p^{}p)^2`$ is the square of the four momentum transfer. If one uses the plus-component, $`J^+=J^0+J^3`$, the LF calculation gives two finite contributions, the LF valence form factor and the LF nonvalence form factor, that add up to the covariant result, as expected. The importance of the nonvalence contribution varies strongly with the momentum transfer and depends sensitively on the binding energy of the meson. For small values of $`q^2`$ and small binding energy, the valence part is dominant, but elsewhere the nonvalence diagram is essential for agreement between the LF calculation and the covariant results.

The form factor can also be extracted from the minus-component of the current, $`J^{}=J^0J^3`$. Covariance guarantees that it makes no difference whether the form factor is determined using the plus or the minus current matrix element. As LFD is not manifestly covariant, it may happen that $`J^{}`$ leads to a form factor different from the one determined using $`J^+`$. As we showed in Ref. , the matrix element of $`J^{}`$ diverges even in 1+1 dimensional LFD. Unless one regulates $`J^{}`$, the current cannot be conserved. To assure current conservation, it was crucial to identify the term that causes the divergence. We have identified this term exactly and found that it is an infinite function of the momentum transfer. If this infinite term is subtracted, the two LF contributions become finite as it must be in a conserved current. Moreover, their sum equals again the covariant result as expected. Still, the regularized LF contributions obtained from $`J^{}`$ are different from the ones extracted from the plus current. The differences grow with increasing binding energy.

The covariant fermion triangle-loop (Fig. 1(a)) in 3+1 dimension is divergent if all the vertices are point-like. In the recent literature, the fermion-loop was regulated by smearing the $`q\overline{q}`$ bound-state vertex in a covariant manner. However, the vertex function used in Refs. was not symmetric in the four momenta of the constituent quarks as we discussed in the introduction. We note that the conceptual difficulty associated with the asymmetry could be remedied if the cutoff prescription is used in the gauge-boson vertex rather than in the $`q\overline{q}`$ bound-state vertex. In Fig. 1(a), we thus replace the point photon-vertex $`\gamma ^\mu `$ by a non-local (or smeared) photon-vertex $`S_\mathrm{\Lambda }(kp)\gamma ^\mu S_\mathrm{\Lambda }(kp^{})`$, where $`S_\mathrm{\Lambda }(p)=\mathrm{\Lambda }^2/(p^2\mathrm{\Lambda }^2+iϵ)`$ and $`\mathrm{\Lambda }`$ plays the role of a momentum cut-off. Our method is gauge invariant, and satisfies the W-T identity. Even though we have computed the covariant amplitude (Fig. 1(a)) with unequal constituent masses, for the clarity of presentation we will focus in this Section on the equal mass case, i.e. $`m_q=m_{\overline{q}}=m`$. The basic formulas for the general case are given in the Appendix.

The covariant amplitude (Fig. 1(a)) for a pseudoscalar meson is given in the equal-mass case by

$$p^{}|J^\mu |p=4N_cg^2\mathrm{\Lambda }^4\frac{d^4k}{(2\pi )^4}\frac{(m^2k^2+pp^{})k^\mu +(k^2m^2kp^{})p^\mu +(k^2m^2kp)p^\mu }{D(k)D(kp)D(kp^{})D_\mathrm{\Lambda }(kp)D_\mathrm{\Lambda }(kp^{})},$$
(2)
where $`N_c`$ is the number of colors and $`g`$, modulo the obvious charge factor $`e_m`$ is the normalization constant fixed by the unity of the form factor at zero momentum transfer and the denominator factor $`D(k)`$ from the quark propagator with momentum $`k`$ is given by $`D(k)=k^2m^2+iϵ`$. While the bound-state vertex is still pointlike, it satisfies corresponding relativistic bound-state equation and a binding effect of the two-body bound-state is thus included in our analysis. We note that one can split the denominators in Eq. (2) into four terms (see below) to show the similarity between the methods of photon-vertex-smearing and the Pauli-Villars regularization, namely the vertex-smearing-technique shares the same structure of the denominators with the Pauli-Villars regularization. To investigate the issue of the end-point singularity, we present the $`J^{}`$ calculation in the following.

First, the valence contribution shown in Fig. 1(b) is obtained in the range $`0<k^+<p^+`$. By picking up the pole of the spectator quark, i.e. $`k^{}=\frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}iϵ}{k^+}`$, we obtain

$`p^{}|J^{}|p`$ $`=`$ $`{\displaystyle \frac{4\pi iN_cg^2\mathrm{\Lambda }^4}{(\mathrm{\Lambda }^2m^2)^2}}{\displaystyle \frac{dk^+d^2\stackrel{}{k}_{}}{(2\pi )^4}\frac{\{k^+p^{}p_{}^{}{}_{}{}^{}+\stackrel{}{k}_{}\stackrel{}{p^{}}_{}p^{}+\stackrel{}{k}_{}\stackrel{}{p}_{}p_{}^{}{}_{}{}^{}\stackrel{}{p}_{}\stackrel{}{p^{}}_{}k^{}\}}{k^+(k^+p^+)(k^+p_{}^{}{}_{}{}^{+})}}`$ (4)
$`\times \left[{\displaystyle \frac{1}{E(p,\mathrm{\Lambda })E(p^{},\mathrm{\Lambda })}}{\displaystyle \frac{1}{E(p,\mathrm{\Lambda })E(p^{},m)}}{\displaystyle \frac{1}{E(p,m)E(p^{},\mathrm{\Lambda })}}+{\displaystyle \frac{1}{E(p,m)E(p^{},m)}}\right],`$
where the energy denominator $`E(p,\mathrm{\Lambda })`$ is defined as

$$E(p,\mathrm{\Lambda })=\frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{k^+}p^{}\frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{p}_{})^2}{k^+p^+}.$$
(6)
As we see from Eq. (LABEL:eq2.3), the result depends on the reference frame. However, as we will see later the sum of the valence and nonvalence contributions will of course be identical in any frame. In the present Section, we choose the frame where $`q^+=\alpha p^+`$ and $`\stackrel{}{p}_{}=0`$ and the momentum fraction $`x`$ is defined by $`k^+=(1x)p^+`$. Then the LF valence form factor (see Eq. (1)) is given by

$`F_{\mathrm{val}}^{}(q^2)`$ $`=`$ $`{\displaystyle \frac{2N_cg^2\mathrm{\Lambda }^4}{(2+\alpha )+\stackrel{}{q}_{}^{}{}_{}{}^{2}/M^2}}{\displaystyle _0^1}𝑑x{\displaystyle \frac{d^2\stackrel{}{k}_{}}{(2\pi )^3}\frac{(1x)(M^2+\stackrel{}{q}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}})(1+\alpha )\stackrel{}{k}_{}\stackrel{}{q}_{}}{(1x)x^2(x+\alpha )^2}}`$ (7)
$`\times `$ $`{\displaystyle \frac{1}{\left(\frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{1x}+\frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{x}M^2\right)\left(\frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{1x}+\frac{\mathrm{\Lambda }^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{x}M^2\right)}}`$ (8)
$`\times `$ $`{\displaystyle \frac{1}{\left(\frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{1x}+\frac{m^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{x+\alpha }\frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}}{1+\alpha }\right)\left(\frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{1x}+\frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{x+\alpha }\frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}}{1+\alpha }\right)}},`$ (9)
where $`M`$ is the meson mass. The limit to the point vertex can be taken by letting $`\mathrm{\Lambda }\mathrm{}`$ and we find that the result is finite in this frame.

Next, the nonvalence contribution shown in Fig. 1(c) is obtained in the range $`p^+<k^+<p_{}^{}{}_{}{}^{+}`$. Here, the pole is not taken as the spectator but as $`k^{}=p_{}^{}{}_{}{}^{}+\frac{m^2+(\stackrel{}{k}_{}\stackrel{}{p^{}}_{})^2iϵ}{k^+p_{}^{}{}_{}{}^{+}}`$ (for the term corresponding to the Pauli-Villars particle, $`k^{}=p_{}^{}{}_{}{}^{}+\frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{p^{}}_{})^2iϵ}{k^+p_{}^{}{}_{}{}^{+}}`$). Following a procedure similar to the one described in the valence case, we find the LF nonvalence form factor to be given by

$`F_{\mathrm{nv}}^{}(q^2)`$ $`=`$ $`{\displaystyle \frac{2N_cg^2\mathrm{\Lambda }^4(1+\alpha )}{(\mathrm{\Lambda }^2m^2)[M^2(2+\alpha )+\stackrel{}{q}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}]}}`$ (11)
$`\times {\displaystyle _0^\alpha }dx{\displaystyle }{\displaystyle \frac{d^2\stackrel{}{k}_{}}{(2\pi )^3}}{\displaystyle \frac{1}{x^2(1+x)(\alpha x)}}[{\displaystyle \frac{N(m)}{D0D1D2}}{\displaystyle \frac{N(\mathrm{\Lambda })}{D5D3D4}}],`$
where

$`N(m)`$ $`=`$ $`\left[m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}(1+x)\left({\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2}}{1+\alpha }}{\displaystyle \frac{m^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{\alpha x}}\right)\right]\left[{\displaystyle \frac{m^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{\alpha x}}+M^2\right]`$ (13)
$`+{\displaystyle \frac{(1+x)M^2(M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2})}{1+\alpha }}\stackrel{}{k}_{}\stackrel{}{q}_{}M^2,`$
$`N(\mathrm{\Lambda })`$ $`=`$ $`\left[m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}(1+x)\left({\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2}}{1+\alpha }}{\displaystyle \frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^{\mathrm{\hspace{0.17em}2}}}{\alpha x}}\right)\right]\left[{\displaystyle \frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{\alpha x}}+M^2\right]`$ (15)
$`+{\displaystyle \frac{(1+x)M^2(M^2+\stackrel{}{q}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}})}{1+\alpha }}\stackrel{}{k}_{}\stackrel{}{q}_{}M^2,`$
$`D0`$ $`=`$ $`{\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2}}{1+\alpha }}{\displaystyle \frac{m^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^{\mathrm{\hspace{0.17em}2}}}{\alpha x}}{\displaystyle \frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{1+x}},`$ (16)
$`D1`$ $`=`$ $`M^2{\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2}}{1+\alpha }}+{\displaystyle \frac{m^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^{\mathrm{\hspace{0.17em}2}}}{\alpha x}}+{\displaystyle \frac{\mathrm{\Lambda }^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{x}},`$ (17)
$`D2`$ $`=`$ $`M^2{\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2}}{1+\alpha }}+{\displaystyle \frac{m^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{\alpha x}}+{\displaystyle \frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{x}},`$ (18)
$`D5`$ $`=`$ $`{\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}}{1+\alpha }}{\displaystyle \frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{\alpha x}}{\displaystyle \frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{1+x}},`$ (19)
$`D3`$ $`=`$ $`M^2{\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2}}{1+\alpha }}+{\displaystyle \frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{\alpha x}}+{\displaystyle \frac{\mathrm{\Lambda }^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{x}},`$ (20)
$`D4`$ $`=`$ $`M^2{\displaystyle \frac{M^2+\stackrel{}{q}_{}^{}{}_{}{}^{2}}{1+\alpha }}+{\displaystyle \frac{\mathrm{\Lambda }^2+(\stackrel{}{k}_{}\stackrel{}{q}_{})^2}{\alpha x}}+{\displaystyle \frac{m^2+\stackrel{}{k}_{}^{}{}_{}{}^{2}}{x}}.`$ (21)
We find from Eq. (11) that the two terms in the integrand are individually finite as long as the parameter $`\mathrm{\Lambda }`$ is finite. There is neither an ultra-violet divergence nor an end-point singularity. However, in the limit to the point vertex (i.e. $`\mathrm{\Lambda }\mathrm{}`$), we find that each term has not only a linear divergence in the $`\stackrel{}{k}_{}`$ integration, but also an end-point singularity in the $`x`$ integration, which cancel each other exactly. Thus, in the point vertex limit, we find that the end-point singularity is completely removed even though the result is logarithmically divergent as it must be in the 3+1 dimensional fermion-loop with point vertices. This shows a striking difference from the calculation without relying on the vertex regularization from the beginning. The critical reason for this is that the end-point singularity for the fermion-loop is a consequence of the bottomless nature of the Dirac sea and the vertex-smearing effectively provides the weighting in the Dirac sea de-emphasizing the lower part. The previously identified end-point singularity is exactly cancelled by the identical end-point singularity from the term generated by the vertex-smearing corresponding to the Pauli-Villars particle. Therefore, in the regularized case, all the physical degrees of freedom are taken into account.

In addition, the sum of the valence and nonvalence contributions $`F(q^2)=F_{\mathrm{val}}(q^2)+F_{\mathrm{nv}}(q^2)`$ is of course identical to the results obtained by other components of the current, i.e. either $`J^+`$ or $`J^{}`$.In the electromagnetic form factor, there is no zero-mode contribution in $`J_{}`$ current. Also, the net result $`F(q^2)`$ is independent of the choice of reference frame.

A similar calculation can be made for the pseudoscalar-meson decay constant $`f`$. Taking a non-local gauge-boson vertex, we verified again the exact cancellation of linear and logarithmic divergences in the (two-point) fermion loop. The result for the equal mass case such as the pion is given by

$$f=\frac{N_cgm\mathrm{\Lambda }^4}{4\sqrt{2}\pi ^2(\mathrm{\Lambda }^2m^2)^2}_0^1𝑑x\mathrm{log}\left[\frac{C(m,m)C(\mathrm{\Lambda },\mathrm{\Lambda })}{C(m,\mathrm{\Lambda })C(\mathrm{\Lambda },m)}\right],$$
(22)
where $`C(m_1,m_2)`$ is given by

$$C(m_1,m_2)=x(1x)M^2(1x)m_1^2xm_2^2.$$
(23)
If the asymmetric $`q\overline{q}`$ bound-state vertex-smearing is used as in Refs., the above result, Eq.(22), is replaced by

$$f=\frac{N_cgm\mathrm{\Lambda }^2}{4\sqrt{2}\pi ^2(\mathrm{\Lambda }^2m^2)}_0^1𝑑x\mathrm{log}\left[\frac{C(m,m)}{C(m,\mathrm{\Lambda })}\right].$$
(24)
We note that the logarithmic divergence is not completely cancelled in the unequal mass case if the minus component of the weak current is used with the asymmetric $`q\overline{q}`$ bound-state vertex-smearing. The non-local gauge-boson vertex used in this work does not suffer from such an incomplete cancellation of divergences no matter which component of the current is used.

In the next section, we present the numerical taxonomy of the $`\pi `$, $`K`$, and $`D`$ meson electromagnetic LF form factors choosing various reference-frames, as well as the values of the decay constants $`f_\pi ,f_K`$ and $`f_D`$.

## III Numerical Results

In Table I, we present our model parameters such as the constituent quark masses ($`m`$), and the cutoff values ($`\mathrm{\Lambda }`$), as well as the decay constants that we calculated here and compared to the experimental data . The meson masses ($`M`$) are taken as the experimental values . Our model parameters have been chosen to fit both the charge radii and the decay constants to the experimental data well, although the available data for the decay constant of the $`D`$-meson is only an upper limit. We also compare in Table I the decay constants from our symmetric non-local gauge-boson vertex ($`f^{\mathrm{sym}}`$) with those from the asymmetric $`q\overline{q}`$ bound-state vertex ($`f^{\mathrm{asym}}`$) smearing case. However, we note that our calculation here is limited in value because the zero-range approximation is used for the bound-state vertices.

Three different reference frames were considered: the Drell-Yan-West (DYW) frame, the target-rest frame (TRF) and the Breit frame. The corresponding kinematics is given in Table II. The DYW frame has gained some popularity in deep-inelastic scattering calculations because in that frame $`q^+=0`$ identically. This frame can be obtained by taking the $`\alpha 0`$ limit in the frame presented in section II. In the limit $`\alpha 1`$, the frame presented in section II coincides with the target-rest frame with $`\theta =0`$.

Now we comment on our results shown in the figures below. (In all these figures, we use thick solid line for the covariant form factor, thick dashed line for $`F_{\mathrm{val}}^+`$, thick dot-dashed line for $`F_{\mathrm{nv}}^+`$, thin dashed and dot-dashed lines for the corresponding minus LF form factors.) We show in Fig. 2 the results of our numerical calculations using the DYW frame and compare with the experimental data  Our total results, represented by the solid line in each figure, are also in very good agreement with the experimental data of the pion and kaon form factors, respectively.

In the DYW frame $`F_{\mathrm{nv}}^+`$ vanishes identically. Remarkably, we find that the nonvalence part of the minus current, which in this reference frame coincides with the zero-mode contribution, makes a very important contribution to the total form factor and may even dominate over the valence part in the whole $`Q^2`$-range considered. There are quantitative differences between the results obtained for the different mesons, which are to large extent due to the difference in binding: the tighter the binding, the more important $`F_{\mathrm{nv}}^{}`$ becomes. We checked the binding effect in the case of the pion. By varying the quark mass alone from the realistic value of $`0.25`$ GeV to $`0.07007`$ GeV, thus lowering the binding energy to 0.1% of the pion mass, we found that the value of $`F_{\mathrm{val}}^{}`$ at $`Q^2=0`$ was increased from 6% to 69%. Still, for $`Q^2>0.3`$ GeV2 $`F_{\mathrm{val}}^{}`$ is lower than $`F_{\mathrm{nv}}^{}`$. This indicates that for the larger values of $`Q^2`$, the relativistic effects can still be large.

The frame dependence of the different components of the current can be studied by comparing the results of calculations in different frames and at different values of the polar angle $`\theta `$. It is worth mentioning that the results must be independent of the azimuthal angle $`\varphi `$, because rotations about the $`z`$-axis are kinematical transformations. (We used this fact as a check on the correctness of our codes.)

In Fig. 3 we show the results of our numerical calculations in the Breit frame at $`\theta =0`$, which differs from the target-rest frame at $`\theta =0`$ by a boost in the $`z`$-direction only, so the results are identical in these two frames. The first thing we notice is the great difference with the DYW frame. Now there is a sizeable contribution from $`F_{\mathrm{nv}}^+`$, which dominates at higher values of the momentum transfer: at $`Q^2=0`$, $`F_{\mathrm{val}}^+`$ coincides with the covariant form factor. It crosses $`F_{\mathrm{nv}}^+`$ at some value of $`Q^2`$, the crossover point being smaller for larger binding energy. It is of special interest to separate the instantaneous part, i.e. the contribution to $`F^+`$ from diagrams with one internal instantaneous propagator. They are given in the figures by crosses ($`\times `$). It turns out that they give a very large contributions to the pion form factor but become negligible for the heavy $`D`$-meson case.

Turning to the minus current we see that in this reference frame the relative importance of $`F_{\mathrm{nv}}^{}`$ becomes more prominent than in the DYW frame for the more tightly bound mesons $`\pi `$ and $`K`$. In the case of the pion the dominance of $`F_{\mathrm{nv}}^{}`$ is so strong, that it can hardly be distinguished from the covariant form factor. The figures (in Fig. 3) suggest that at $`Q^2=0`$ the values obtained are frame independent. This is indeed the case, as we found in our calculations and can be understood as the equality of the form factors in the long-wave-length limit.

A first glance at the angle dependence is given in Fig. 4. In both the target-rest frame and the Breit frame for $`\theta =\pi /2`$ the transferred momentum is purely transverse. As the pure Lorentz boosts in the transverse direction, $`K_1`$ and $`K_2`$, are dynamical in LFD, we must expect a strong frame dependence and this is indeed what we find. As $`q^+=0`$ in the Breit frame for $`\theta =\pi /2`$, the part $`F_{\mathrm{nv}}^+`$ vanishes identically. Still, the DYW frame results differ for the minus current, which is clear from a comparison with Fig. 2. The reason is that they are not connected by a kinematical transformation. For the same reason the zero-mode contributions differ in the two frames. In the DYW frame it is very close to the total form factor, while in the Breit frame it even overshoots the covariant form factor at large values of $`Q^2`$ by a factor of almost 2. In this connection we want to mention the work of Frederico et al. , who performed a calculation similar in spirit to ours, but claimed that in the Breit frame at $`\theta =\pi /2`$ $`F_{\mathrm{val}}^{}`$ is always very close to the covariant form factor. This discrepancy was due to the difference of their $`F_{\mathrm{val}}^{}`$ definition in the Breit frame at $`\theta =\pi /2`$, where they removed the term that is odd under the transformation $`p^{}k^{}(p^{}k^{})`$ in their definition. However, our definition of $`F_{val}^{}`$ is general for any angle $`\theta `$ and our results on $`\theta `$-dependence are smooth as shown in Fig. 9.

We show the systematics of the angle dependence for the case of the pion only. Our results are depicted in Figs. 5-8. ( The thick solid line for $`\theta =0`$, thin dotted line for $`\theta =\pi /4`$, the thick dashed line for $`\theta =\pi /4`$, the thin dot-dashed for $`\theta =3\pi /4`$, and the the thick long-dashed line for $`\theta =\pi `$.) One sees immediately that the angle dependence is smooth but can be very strong, both for the valence and the nonvalence parts, calculated from either the plus or the minus current.

One might try to exploit the angle dependence to optimize the calculation of the form factor in a noncovariant framework. However, as the figures above clearly show, there is no value for the angle $`\theta `$ where both $`F_{\mathrm{nv}}^+`$ and $`F_{\mathrm{nv}}^{}`$ are negligible, or even suppressed, compared to the valence parts, for all values of $`Q^2`$. On the contrary, as the values of the form factor components at $`Q^2=0`$ are frame independent, we can be sure that $`F_{\mathrm{nv}}^{}`$ must be very important for an important part of the $`Q^2`$ range.

As a summary of what we found concerning the angle dependence we show in Fig. 9 the complete angle dependence for two values of the momentum transfer, $`Q^2=0`$ and 0.1 GeV2. All curves in this figure are clearly smooth and demonstrate the fact that there is no preferred value for the polar angle $`\theta `$.

## IV Conclusion and Discussion

In this paper, we have analyzed all the components of the current quantized on the light-front to compute the electromagnetic form factors of pseudoscalar mesons with spin-1/2 consituents. Since our aim in this work was to analyze the taxonomy of the triangle diagram, we did not choose any particular electromagnetic gauge but just presented the equivalence of the physical form factor $`F(q^2)`$ in any choice of the current component. The divergence appearing in the 3+1 dimensional fermion-loop calculations was regulated by the covariant vertex-smearing-technique. Performing the light-front calculation, we verified that the vertex-smearing-technique is similar to the Pauli-Villars regularization. In the $`J^{}`$ computation, we find that the critical smearing effect persists even in the limit to the point vertex because the end-point singularity existing otherwise is completely removed once the limit is taken at the end of the calculation. If the limit is taken at the beginning of the calculation, however, we have alreay shown that the end-point singularity in the nonvalence contribution leads to an infinitely different result from that obtained by the covariant Feynman calculation . Our taxonomical analysis demonstrated that each individual contribution, whether valence or nonvalence, is finite regardless of which component of the current is considered. However, we stress that each contribution depends on the reference-frame even though the sum does not. Of course, the invariance of the sum ensures the current conservation. Also, the zero-mode contribution should be distinguished from the instantaneous contribution as we have numerically estimated the differences.

¿From our numerical calculations we can conclude that to obtain agreement with the covariant form factor one needs both the valence and the nonvalence parts. For tightly bound states the nonvalence parts dominate in an important part of the range of $`Q^2`$-values that we studied. It is natural that this result runs counter to nonrelativistic intuition, which says that the valence parts should dominate, because the tightly bound states are not expected to be non-relativistic. In Ref.  it was demonstrated in the 1+1-dimensional case, that by weakening the binding $`F_{\mathrm{val}}^+`$ and $`F_{\mathrm{val}}^{}`$ approximate the covariant form factor more and more closely. Here we found that indeed for less tightly bound state the valence parts come closer to the covariant result, but even at $`Q^2=0`$ $`F_{\mathrm{val}}^{}`$ gives only 69% of the covariant form factor.

If two reference frames can be connected by a kinematical Lorentz transformation the LF form factors calculated in these two frames must be the same. Otherwise they must differ. We found that the angle dependence within a chosen reference frame is always smooth, although it may be very strong. For some values of the polar angle $`\theta `$ two reference frames may be connected by a kinematical Lorentz transformation, e.g., the target-rest frame and the Breit frame for $`\theta =0`$. On the other hand, the DYW frame and the Breit frame for $`\theta =0`$ are not connected by a kinematical transformation and the form factor components consequently do not coincide. More details of the frame-dependence can be found in Ref. .

While the calculations carried out in this work is more realistic than the case of 1+1 dim with the point vertex, they are still semi-realistic as the model uses a bound-state rather than a confined-state. Moreover, we used for the vertices those obtained from a Bethe-Salpeter equation using a contact interaction which are of zero range and may emphasize the importance of the nonvalence contributions rather differently from a more realistic model. Thus, there is still much room for extending our model towards a more realistic model. However, the essential conclusions about the similarity to the Pauli-Villars regularization and the cancellations of both ultra-violet divergence and the end-point singularity remain intact. Nevertheless, the numerical results may differ from a more realistic model calculation. This point is presently under investigation.

###### Acknowledgements.

This work was supported in part by a grant from the US Department of Energy and the National Science Foundation. This work was started when HMC and CRJ visited the Vrije Universiteit and they want to thank the staff of the department of physics at VU for their kind hospitality. BLGB wants to thank the staff of the department of physics at NCSU for their warm hospitality during a stay when this work was completed. The North Carolina Supercomputing Center and the National Energy Research Scientific Computer Center are also acknowledged for the grant of Cray time.

## A Unequal-mass case

If we number the momenta of the internal lines of the fermion triangle as $`k_1`$ (spectator), $`k_2=p+k_1`$, and $`k_3=p^{}+k_1`$ (struck quark) respectively, and the corresponding masses as $`m_1`$, $`m_2`$, and $`m_3`$ then we find for the trace appearing in the numerator of the covariant integral, the expression

$$T^\mu =4\left[(m_2m_3\stackrel{}{k}_2\stackrel{}{k}_3)k_1^\mu +(m_1m_3\stackrel{}{k}_2\stackrel{}{k}_3)k_2^\mu +(m_1m_2\stackrel{}{k}_2\stackrel{}{k}_3)k_3^\mu \right].$$
(A1)
The valence diagram is obtained if one calculates the integral over $`k_1^{}`$ by closing the contour around the pole corresponding to putting the spectator $`k_1`$ on shell. This corresponds to the following values for the minus-components of the momenta to be used in the expression for $`T^\mu `$

$$k_1^{}=\frac{\stackrel{}{k}_1^2+m_1^2}{k_1^+},k_2^{}=p^{}+k_1^{},k_3^{}=p^{}+k_1^{}.$$
(A2)
To obtain the nonvalence diagram one closes the contour around the pole corresponding to $`k_3`$. Then one gets for the part of the propagator with the mass $`m_3`$

$$k_3^{}=\frac{\stackrel{}{k}_3^2+m_3^2}{k_3^+},k_1^{}=p^{}+k_3^{},k_2^{}=p^{}p^{}+k_3^{}.$$
(A3)
For the part of the propagator with the cutoff $`\mathrm{\Lambda }^{}`$ the same formula can be used, but with $`m_3`$ replaced by $`\mathrm{\Lambda }^{}`$.

As the smearing we use affects only the denominator parts of the propagator of the struck quark, the replacement of $`m_3`$ by $`\mathrm{\Lambda }^{}`$ in $`T^\mu `$ occurs only in the minus components.

The energy denominators are easy to find. One obtains for the valence parts

$`D(kp)`$ $`=`$ $`p^{}+{\displaystyle \frac{\stackrel{}{k}_1^2+m_1^2}{k_1^+}}{\displaystyle \frac{\stackrel{}{k}_2^2+m_2^2}{k_2^+}},`$ (A4)
$`D(kp^{})`$ $`=`$ $`p^{}+{\displaystyle \frac{\stackrel{}{k}_1^2+m_1^2}{k_1^+}}{\displaystyle \frac{\stackrel{}{k}_3^2+m_3^2}{k_3^+}}.`$ (A5)
The nonvalence part has the same denominator $`D(kp^{})`$, but $`D(kp)`$ is changed to

$$D^{}(kp)=p^{}p^{}+\frac{\stackrel{}{k}_3^2+m_3^2}{k_3^+}\frac{\stackrel{}{k}_2^2+m_2^2}{k_2^+}.$$
(A6)
The energy denominators $`D_\mathrm{\Lambda }`$ and $`D_\mathrm{\Lambda }^{}^{}`$ are obtained by the substitutions $`m_2\mathrm{\Lambda }`$ and $`m_3\mathrm{\Lambda }^{}`$, respectively.

In the case where $`m_1=m_2=m_3`$ and $`\mathrm{\Lambda }=\mathrm{\Lambda }^{}`$ the final formulas are much simplified. The explicit formulas in the main text are valid for the equal-mass case. The general case can be easily constructed from the expressions (A1 \- A6).

In the explicit formulas we use the notation

$$k^\mu =k_1^\mu ,x=\frac{k^+}{p^+}(\mathrm{valence}),x=\frac{k^+p^+}{p^+}(\mathrm{nonvalence}),\alpha =\frac{p^+p^+}{p^+}.$$
(A7)
## B Analyticity

¿From the covariant expression for the amplitude one can prove that the form factors are analytic functions of $`Q^2`$. This proof is not valid for the LF-time-ordered amplitudes. One can, however, expand the expressions that one obtains for the different parts of the form factor in terms of $`Q`$ and determine the $`Q`$-dependence at small values of $`Q^2`$. It turns out that the results depend on the kinematics: it matters in which reference frame one does the calculations.

We use the fromulas for the momenta given in Table II and expand the trace and the energy denominators in powers of $`Q`$. From this expansion the analyticity properties of the amplitudes follow. We have numerically verified that a blow-up of Figs. 5-8 for small $`Q^2`$ illustrates the analyticity properties of the amplitudes discussed in this Appendix. The three different reference frames are discussed consecutively.

### 1 Target-rest frame

We discuss the valence part of the plus-current in detail; the other LF form factors can be treated in a similar way. First we expand the momenta.

$`q^\pm `$ $`=`$ $`q\mathrm{cos}\theta +𝒪(Q^2),\stackrel{}{q}_{}=Q\mathrm{sin}\theta \widehat{n}+𝒪(Q^3),`$ (B1)
$`p^\pm `$ $`=`$ $`M+Q\mathrm{cos}\theta +𝒪(Q^2),\stackrel{}{p}_{}^{}=Q\mathrm{sin}\theta \widehat{n}+𝒪(Q^3).`$ (B2)
The trace consists of a piece that is independent of $`Q`$ and a piece that consists of two parts, one that is proportional to $`Q\stackrel{}{k}_1\widehat{n}\mathrm{sin}\theta `$ and a part proportional to $`Q\mathrm{cos}\theta `$. The denominator $`D(kp^{})`$ has a similar behaviour, but $`D(kp)`$ is independent of $`Q`$.

Upon integration over $`\stackrel{}{k}_1`$ the terms proportional to $`\stackrel{}{k}_1\widehat{n}\mathrm{sin}\theta `$ vanish. Consequently, the valence part $`F_{\mathrm{val}}^+`$ has the small-$`Q`$ behaviour

$$F_{\mathrm{val}}^+F_{\mathrm{val}}^{+\mathrm{\hspace{0.17em}0}}+F_{\mathrm{val}}^{+\mathrm{\hspace{0.17em}1}}Q\mathrm{cos}\theta .$$
(B3)
The other cases, nonvalence plus-current, valence and nonvalence minus-current, show the same pattern. There are pieces independent of $`Q`$, parts with the $`Q`$-dependence $`Q\stackrel{}{k}_1\widehat{n}\mathrm{sin}\theta `$ and ones proportional to $`Q\mathrm{cos}\theta `$. As the pieces proportional to $`\stackrel{}{k}_1\widehat{n}`$ vanish upon integration, all the components of the form factor show a behaviour similar to Eq. (B3). So, only for $`\theta =\pi /2`$ do we find that the LF time-ordered amplitudes calculated in the target-rest frame are analytic in $`Q^2`$.

### 2 Breit frame

In the case of the Breit frame we can follow the same line as in the case of the target-rest frame. The only difference is that now no terms of the form $`Q\stackrel{}{k}_1\widehat{n}\mathrm{sin}\theta `$ appear. As those terms give no contribution to first order in $`Q`$ anyway, this does not alter the result: also in the Breit frame the components of the form factor have the same small-$`Q`$ behaviour as in Eq. (B3).

### 3 Drell-Yan-West frame

In the Drell-Yan-West frame there is no angle dependence. However, there is a term linear in $`Q`$. It is proportional to $`Q\widehat{n}\stackrel{}{k}_1`$. Of course, this term also vanishes upon integration over $`\stackrel{}{k}_1`$. Therefore, in the Drell-Yan-West frame we find no term linearly dependent on $`Q`$, so the amplitudes are analytic in $`Q^2`$.

### 4 Zero mode

The zero mode is defined as the contribution to the nonvalence amplitude that survives the limit $`q^+0`$. It is easy to see that only the minus-current can have a zero-mode part, because the integral defining $`F_{\mathrm{nv}}^+`$ has an integrand that remains finite when $`q^+`$ goes to 0. However, the integrand defining $`F_{\mathrm{nv}}^{}`$ diverges when the limit $`q^+0`$ is taken. In order to determine the limit, one may expand the integrand in powers of $`q^+`$. As $`q^+`$ is taken to be zero in the Drell-Yan-West frame, the analysis cannot be done in that frame, but it can most easily be carried out in the Breit frame.

The algebra being straightforward but tedious, we shall not give the details. Rather we quote the final result. If we consider the smeared case, the most divergent part of the integrand for $`F_{\mathrm{nv}}^{}`$ has the behaviour $`1/q^{+\mathrm{\hspace{0.17em}2}}`$. As the integration over $`k_1^+`$ ranges from $`p^+`$ to $`p^+`$, we can scale the integration variable as in Eq. (A7) and obtain an integral over $`x`$ from 0 to 1. The Jacobian being $`q^+`$, it cancels one factor $`q^+`$ in the denominator, so we only need to show that the leading term vanishes to prove that the zero-mode contribution is finite.

Upon carrying out the algebra we find the following result. The leading part is proportional to

$$x^2(\stackrel{}{k}_3^2+m_3^2)(\stackrel{}{k}_3^2+\mathrm{\Lambda }^{\mathrm{\hspace{0.17em}2}})(1x)^2(\stackrel{}{k}_2^2+m_2^2)(\stackrel{}{k}_2^2+\mathrm{\Lambda }^2).$$
(B4)
This part of the integrand does not vanish, because $`\stackrel{}{k}_2=\stackrel{}{p}_{}\stackrel{}{k}_{}\stackrel{}{k}_3=\stackrel{}{p}_{}^{}\stackrel{}{k}_{}`$. However, if we take $`\mathrm{\Lambda }^{}=\mathrm{\Lambda }`$ and $`m_2=m_3`$, being the mass of the struck quark, and as $`\stackrel{}{k}_2`$, $`\stackrel{}{k}_3`$, and $`\stackrel{}{k}_{}`$ differ only by a constant vector, this function vanishes after integration over $`x`$ from 0 to 1 and $`\stackrel{}{k}_{}`$ over the whole of $`𝖨𝖱_2`$. So we see that the coefficient of the contribution proportional to $`1/q^+`$ vanishes. The remaining part, the piece that survives the limit $`q^+0`$, is the zero mode contribution. In the Breit frame as well as the target-rest frame it is well-defined and finite.


text:
# Effective hadronic Lagrangian for charm mesons

## 1 Introduction

Suppression of charmonium production in high energy heavy ion collisions is a possible signature for the quark-gluon plasma (QGP) formed in ultra-relativistic heavy ion collisions . Recent data from the Pb+Pb collision at $`P_{\mathrm{lab}}=158`$ GeV$`/c`$ in the NA50 experiment at CERN  have shown an anomalously large $`J/\psi `$ suppression in high $`E_\mathrm{T}`$ events. While there are suggestions that this anomalous suppression may be due to the formation of the QGP , more conventional mechanisms based on $`J/\psi `$ absorption by comoving hadrons have also been proposed as a possible explanation . Since the latter scenario depends on the value of $`J/\psi `$ absorption cross sections by hadrons, it is important to have better knowledge of these cross sections in order to understand the observed anomalous charmonium suppression . Even in heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC), where the QGP is most likely to be formed, the effect of hadronic absorption of $`J/\psi `$ is still non-negligible . Furthermore, one needs to know if $`J/\psi `$ can also be produced from the hot hadronic matter in the later stage of heavy ion collisions .

Also, charm mesons such as $`D`$ and $`D^{}`$ are expected to be abundantly produced in heavy ion collisions at RHIC energies and beyond. It has been shown that the high mass ($`M>2`$ GeV) dilepton spectrum at RHIC may be dominated by decays from charm and bottom hadrons . Since charm quarks may lose appreciable energies in a quark-gluon plasma via gluon radiations, study of charm meson spectrum in heavy ion collisions could provide useful information on the properties of the QGP . However, charm mesons may interact strongly with hadrons during later stage of heavy ion collisions, leading to possible changes in their final spectra. As a result, the dilepton spectrum from the decay of charm meson pairs could also be modified . To study the energy loss of charm quarks in the QGP thus requires the understanding of charm meson interactions in hadronic matter.

Various approaches have been used in evaluating the charmonium absorption cross sections by hadrons. In the quark-exchange model, an earlier study  has shown that the $`J/\psi `$ absorption cross section by pion, $`\sigma _{\pi \psi }`$, has a peak value of about $`7`$ mb at $`E_{\mathrm{kin}}\sqrt{s}m_\pi m_\psi 0.8`$ GeV, but a recent study  gives a peak value of only $`\sigma _{\pi \psi }1`$ mb at the same $`E_{\mathrm{kin}}`$ region. On the other hand, the perturbative QCD approach  predicts that the $`J/\psi `$ dissociation cross section increases monotonously with $`E_{\mathrm{kin}}`$ and has a value of only about $`0.1`$ mb around $`E_{\mathrm{kin}}0.8`$ GeV. Charmonium absorption cross sections by hadrons have also been studied in meson-exchange models based on effective hadronic Lagrangians. Using pseudoscalar-pseudoscalar-vector-meson (PPV) couplings and without form factors at the interaction vertices, Matinyan and Müller  have found $`\sigma _{\pi \psi }0.3`$ mb at $`E_{\mathrm{kin}}=0.8`$ GeV. In a recent study, Haglin  has included also the three-vector-meson (VVV) and four-point couplings (or contact terms), and obtained much larger values for the $`J/\psi `$ absorption cross sections. Large discrepancies in the magnitude of $`\sigma _{\pi \psi }`$ (as well as $`\sigma _{\rho \psi }`$) thus exist among the predictions from these three approaches. In this study, we use a meson-exchange model as in  but treat differently the VVV and four-point couplings in the effective Lagrangian. For charm meson scattering cross sections with pion and rho meson, previous studies from the meson-exchange model include only the PPV interactions . We now also use the extended hadronic Lagrangian to study these cross sections.

## 2 The effective Lagrangian

The free Lagrangian for pseudoscalar and vector mesons in the limit of SU(4) invariance can be written as $`_0=\mathrm{Tr}\left(_\mu P^{}^\mu P\right)\mathrm{Tr}\left(F_{\mu \nu }^{}F^{\mu \nu }\right)/2`$, where $`F_{\mu \nu }=_\mu V_\nu _\nu V_\mu `$, and $`P`$ and $`V`$ denote, respectively, the $`4\times 4`$ pseudoscalar and vector meson matrices in SU(4) . To obtain the couplings between pseudoscalar and vector mesons, we introduce the minimal substitution

$$_\mu P𝒟_\mu P=_\mu P\frac{ig}{2}\left[V_\mu P\right],F_{\mu \nu }_\mu V_\nu _\nu V_\mu \frac{ig}{2}[V_\mu ,V_\nu ].$$
(1)
The effective Lagrangian is then given by

$``$ $`=`$ $`_0+ig\mathrm{Tr}\left(^\mu P[P,V_\mu ]\right){\displaystyle \frac{g^2}{4}}\mathrm{Tr}\left([P,V_\mu ]^2\right)`$ (2)
$`+`$ $`ig\mathrm{Tr}\left(^\mu V^\nu [V_\mu ,V_\nu ]\right)+{\displaystyle \frac{g^2}{8}}\mathrm{Tr}\left([V_\mu ,V_\nu ]^2\right).`$
Since the SU(4) symmetry is explicitly broken by hadron masses, terms involving hadron masses are added to Eq.(2) using the experimentally determined values.

The effective Lagrangian in Eq. (2) is generated by minimal substitution, which is equivalent to treating vector mesons as gauge particles. The gauge invariance leads to the current conservation; i.e., in the limit of zero vector meson masses, degenerate pseudoscalar meson masses, and SU(4) invariant coupling constants, one has $`^{\lambda _k\mathrm{}\lambda _l}p_{j\lambda _j}=0`$ for any given process, where the index $`\lambda _j`$ denotes the external vector meson $`j`$. We have checked that all the amplitudes without form factors satisfy the requirement of current conservation.

The above effective Lagrangian allows us to study various interactions of charm mesons and $`J/\psi `$ with hadrons. These include the charm meson scattering such as $`\pi D\rho D^{}`$, the charm meson production and annihilation such as $`\pi \pi D\overline{D}`$, and the charmonium absorption and annihilation such as $`\pi \psi D\overline{D}^{}`$. This effective hadronic Lagrangian has also been extended to SU(5) to study $`\mathrm{{\rm Y}}`$ absorption in hadronic matter . In the following, we show the results for $`J/\psi `$ absorption and charm meson scattering by pion and rho meson.

## 3 $`J/\psi `$ absorptions

For $`J/\psi `$ absorption by $`\pi `$ and $`\rho `$ mesons, we study the following processes:

$`\pi \psi D^{}\overline{D},\pi \psi D\overline{D^{}},\rho \psi D\overline{D},\rho \psi D^{}\overline{D^{}}.`$
The full amplitudes for these processes can be found in .

From the $`D^{}`$ decay width , the coupling constant $`g_{\pi DD^{}}`$ is found to be $`g_{\pi DD^{}}=4.4`$. Using the vector meson dominance (VMD) model, we determine other three-point coupling constants as

$`g_{\rho DD}=g_{\rho D^{}D^{}}=2.52,g_{\psi DD}=g_{\psi D^{}D^{}}=7.64.`$ (3)
For the four-point coupling constants, there is no empirical information, and we thus use the SU(4) relations to determine their values as

$`g_{\pi \psi DD^{}}=g_{\pi DD^{}}g_{\psi DD},g_{\rho \psi DD}=2g_{\rho DD}g_{\psi DD},g_{\rho \psi D^{}D^{}}=g_{\rho D^{}D^{}}g_{\psi D^{}D^{}}.`$ (4)
Form factors are introduced at interaction vertices to take into account the composite nature of hadrons. Unfortunately, there is no empirical information on form factors involving charmoniums and charm mesons. We thus take the form factors as the usual monopole form at the three-point $`t`$ and $`u`$ channel vertices, i.e., $`f_3(t\mathrm{or}u)=\mathrm{\Lambda }^2/(\mathrm{\Lambda }^2+𝐪^2)`$, where $`\mathrm{\Lambda }`$ is a cutoff parameter, and $`𝐪^2`$ is the squared three momentum transfer in the c.m. frame, given by $`(𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟑})_{\mathrm{c}.\mathrm{m}.}^2`$ and $`(𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟒})_{\mathrm{c}.\mathrm{m}.}^2`$ for $`t`$ and $`u`$ channel processes, respectively. For simplicity, we use the same value for all cutoff parameters, i.e., $`\mathrm{\Lambda }_{\pi DD^{}}=\mathrm{\Lambda }_{\rho DD}=\mathrm{\Lambda }_{\rho D^{}D^{}}=\mathrm{\Lambda }_{\psi DD}=\mathrm{\Lambda }_{\psi D^{}D^{}}\mathrm{\Lambda }`$, and choose $`\mathrm{\Lambda }`$ as either $`1`$ or $`2`$ GeV to study the uncertainties due to form factors. We also assume that the form factor at four-point vertices is given by $`f_4=\left[\mathrm{\Lambda }^2/(\mathrm{\Lambda }^2+\overline{𝐪^2})\right]^2`$, where $`\overline{𝐪^2}`$ is the average value of the squared three momentum transfers in $`t`$ and $`u`$ channels, $`\overline{𝐪^2}=p_{i,\mathrm{c}.\mathrm{m}.}^2+p_{f,\mathrm{c}.\mathrm{m}.}^2`$.

Figure 1a shows the cross section as a function of the c.m. energy of initial-state mesons. Although form factors strongly suppress the cross sections, the $`J/\psi `$ absorption cross sections remain appreciable after including form factors. The values for $`\sigma _{\pi \psi }`$ and $`\sigma _{\rho \psi }`$ are roughly $`7`$ and $`3`$ mb, respectively, and are comparable to those used in phenomenological studies of $`J/\psi `$ absorption by comoving hadrons in relativistic heavy ion collisions . The thermal average of $`J/\psi `$ absorption cross sections is shown in figure 1b. At the temperature of $`150`$ MeV, for example, $`\sigma _{\pi \psi }v`$ and $`\sigma _{\rho \psi }v`$ are about $`1`$ and $`2`$ mb, respectively.

## 4 Charmed meson scattering

We consider the following processes for charm meson scattering by $`\pi `$ and $`\rho `$ mesons:

$`\pi D\rho D^{},\pi D\pi D,\pi D^{}\pi D^{},\pi D^{}\rho D,\rho D\rho D,\rho D^{}\rho D^{}.`$
There are also similar processes for anticharm mesons. The full amplitudes for the above processes can be found in .

For the coupling constants, we use $`g_{\rho \pi \pi }=6.1`$ . The SU(4) symmetry then gives $`g_{\rho \rho \rho }=g_{\rho \pi \pi }`$ and the following relations for the four-point couplings:

$`g_{\pi \rho DD^{}}=g_{\pi DD^{}}g_{\rho DD},g_{\pi \pi D^{}D^{}}=2g_{\pi DD^{}}^2,g_{\rho \rho DD}=2g_{\rho DD}^2,g_{\rho \rho D^{}D^{}}=g_{\rho D^{}D^{}}^2.`$ (5)
Form factors at the $`s`$ channel vertices are taken as $`f_3(s)=\mathrm{\Lambda }^2/(\mathrm{\Lambda }^2+𝐤^2)`$ , with $`𝐤`$ denoting the three momentum of either the incoming or outgoing particles in the c.m. frame, i.e., $`𝐤^2=p_{i,\mathrm{c}.\mathrm{m}.}^2`$, or $`p_{f,\mathrm{c}.\mathrm{m}.}^2`$. For form factors at four-point vertices, we take them to be $`f_4^{}=\left[\mathrm{\Lambda }^2/(\mathrm{\Lambda }^2+\overline{𝐤^2})\right]^2`$, where $`\overline{𝐤^2}`$ is the average value of the squared three momenta in the form factors for the $`s`$, $`t`$, and $`u`$ channels, i.e., $`\overline{𝐤^2}=5(p_{i,\mathrm{c}.\mathrm{m}.}^2+p_{f,\mathrm{c}.\mathrm{m}.}^2)/6`$.

The charm meson scattering cross sections and their thermal averages are shown in figure 2a and 2b, respectively. As expected, the magnitude of the cross sections decreases with decreasing cutoff parameter. For the cutoff parameters used here, the cross sections for $`\pi D,\pi D^{},\rho D`$ and $`\rho D^{}`$ scattering are all roughly between 10 and 20 mb. We note that form factors only suppress modestly (by a factor of two and less) these total cross sections and their thermal averages. This is due to the fact that these cross sections are dominated by elastic processes, which involve small momentum transfer near the threshold. In contrast, the process $`\pi \psi D^{}\overline{D}`$ has a large threshold, and form factors suppress its cross section by as much as a factor of 8.

## 5 Summary

In summary, we have introduced an effective hadronic Lagrangian, that includes the charm mesons, to study the interactions of charmed mesons and $`J/\psi `$ in hadronic matter. In particular, we have calculated the absorption cross sections of $`J/\psi `$ and scattering cross sections of charmed mesons by $`\pi `$ and $`\rho `$ mesons. We find that the $`J/\psi `$ absorption cross sections are much larger than those in a previous study, where only pseudoscalar-pseudoscalar-vector-meson couplings were considered. Including form factors at the interaction vertices, the values for $`\sigma _{\pi \psi }`$ and $`\sigma _{\rho \psi }`$ are about $`7`$ mb and $`3`$ mb, respectively, and their thermal averages at the temperature of $`150`$ MeV are roughly $`1`$ mb and $`2`$ mb, respectively. These values suggest that the absorption of $`J/\psi `$ by comoving hadrons may play an important role in $`J/\psi `$ suppression in relativistic heavy ion collisions. We also find that the scattering cross sections of $`D`$ or $`D^{}`$ by $`\pi `$ or $`\rho `$ mesons are all about $`10`$ mb and thus expect these scatterings to significantly modify the charm meson spectra in heavy ion collisions.

## Acknowledgments

This work was supported in part by the National Science Foundation under Grant No. PHY-9870038, the Welch Foundation under Grant No. A-1358, and the Texas Advanced Research Program under Grant No. FY99-010366-0081.

## References


text:
# Universality limits on bulk fermions

## Acknowledgements

It is a pleasure to thank Ll. Ametller, A. Delgado, M. Masip, M. Pérez-Victoria, A. Pomarol, M. Quirós and F. Zwirner for discussions. JS thanks the Dipartimento di Fisica Galileo Galilei and INFN Sezione di Padova for their warm hospitality. This work has been supported by CICYT, by Junta de Andalucía and by the European Union under contract HPRN-CT-2000-00149. JS thanks MECD for financial support.