text:
# Enhanced inverse bremsstrahlung heating rates in a strong laser field

## I Introduction

The description of electron-ion collisions in the presence of a strong oscillatory electric field has been an essential part of every laser plasma interaction model. Processes such as collisional laser energy absorption, bremsstrahlung radiation emission and transport include electron scattering on ions in the presence of a high frequency electromagnetic field. Advances in laser technology, particularly those related to the generation of ultra-short laser pulses and progress in inertial confinement fusion studies have challenged our understanding of scattering processes over a wide range of conditions and plasma parameters. Nonetheless, heating rate calculations do ; silin ; bunkin ; kw ; kull ; pert1 ; pert2 ; decker ; shvets ; mulser  in strong laser fields have involved simplifications that are equivalent to a Born approximation for electron trajectories in the Coulomb fields. This approximation, although well-accepted in the physics community, should be revisited in view of recently discovered nonlinear processes such as harmonic generation and above threshold ionization corkum0 ; brabec ; ivanov  that have been explained as being due to strong modifications of electron orbits in the laser field during electron-ion scattering events. Similar nonlinear modifications of particle trajectories have recently been discussed in numerical studies of inverse bremsstrahlung heating rates fraiman ; fraimannew . Our paper continues this analysis; we examine classical electron trajectories that contribute to enhanced plasma heating fraiman  and discuss their importance for calculations of effective collision frequencies with certain classes of electron distribution functions. We also summarize existing theories of inverse bremsstrahlung heating and compare them with results of our test particle calculations.

For the purpose of making comparisons, seminal works by Dawson and Oberman do  and Silin silin  and numerous subsequent publications on high field effects bunkin ; kw ; kull ; pert1 ; pert2 ; decker ; shvets ; mulser  are classified into several broad categories: classical-mechanical analysis of binary electron-ion collisions in the dielectric approximation do ; pert2 ; decker ; shvets ; kinetic theory which uses the Landau-type collision integral in the oscillating electron field silin ; ballistic model (or the model of instantaneous electron-ion collision) bunkin ; pert1 ; mulser ; ”first principle” quantum-mechanical calculations bunkin ; kw ; and, more recently, the quantum kinetic theory in the dielectric approximation kull . All theories consider the two body scattering problem as a starting point and their predictions lead to remarkably similar expressions which may differ only by a logarithmic factor.

The agreement between these approaches follows from two common characteristics inherent in all theories, i.e. the Born approximation is used for the electron orbits in a Coulomb field and the cross-section is assumed to have a weak dependence on the laser field phase. The first approximation, with notable exception of the low frequency approximation kw , has been employed in most quantum mechanical calculations bunkin  where it applies to fast electrons with the velocity $`v_0`$ satisfying the condition $`Ze^2/\mathrm{}v_01`$ fedorov , where $`Ze`$ is the ion charge. The Born approximation is equivalent to small angle scattering and small momentum exchange in classical kinetic theory studies do ; silin ; shvets . This corresponds to a straight line electron trajectory approximation, that is also a part of ideal plasma collision theories that lead to the Landau or Balescu-Guernsey-Lenard operators ichimaru . A simplified description of electron orbits removes the sensitivity of scattering processes to the initial phase of an electric field. However,it has been demonstrated fraiman  that groups of electrons entering an ion interaction sphere at specific phase can significantly modify the collisional cross-section and the heating rates.

Irregular electron trajectories in the combined, high frequency field of a laser and a Coulomb field of an ion have been found by Wiesenfeld num1 . His numerical studies have identified stochastic trajectories for electrons with quiver velocity larger than the initial drift velocity. This effect has been elaborated upon in detail in recent studies by Fraiman et al. fraiman ; fraimannew . In addition to irregular trajectories that involve quasi-capture of electrons in complicated orbits encircling ions, a phase space analysis has revealed the existence of initial conditions corresponding to large angle correlated collisions with anomalously large energy transfer from the field to particles. It was shown recently in Ref. jetp  that the effect of correlated collisions exists not only in the Coulomb potential but also in short-range potentials in the presence of a strong laser field.

There is a similarity between these processes and the interpretation of such atomic physics phenomena as multi-photon ionization and harmonic generation proposed by Corkum corkum0 . It was suggested that free electrons that are created in the process of tunnel ionization can be brought back to the original atom by the laser field and undergo multiple correlated collisions. This process is sensitive to the phase of an electron in the oscillatory field and is responsible for processes such as high-harmonic generation, two-electron ejection, double ionization, etc. From the point of view of quantum mechanics, the high probability of secondary scattering was explained by the phenomenon of Coulomb focusing brabec ; ivanov . This is a process by which the Coulomb attraction of the ion compensates for the natural dispersion of the electron wave packet. A classical analog involves the slowly moving electron that is brought repeatedly into close proximity of an ion as it performs large amplitude oscillations in the laser field. Each encounter results in a small deflection of an electron trajectory toward the ion. These small changes accumulate into a deflection of an electron trajectory leading to a head-on collision producing large momentum change and causes an electron to leave the ion interaction sphere. An important question in the context of laser produced plasmas is determining under what circumstances such trajectories are statistically significant and can alter macroscopic heating rates.

Our paper examines effects of correlated collisions and irregular scattering trajectories on the inverse Bremsstrahlung heating rate. By using mechanisms such as the Coulomb focusing corkum0 ; brabec , the parachute effect fraiman  or quasi-capture num1 , we analyze strongly modified electron orbits and classify them in terms of their initial phase, impact parameters, and drift velocities. By means of test particle simulations we examine the statistical importance of these trajectories. They indeed strongly modify the heating rates for all three different energy distribution functions used in the simulations, i.e. mono-energetic, Maxwellian and anisotropic electron distribution functions.

As we show in this paper, there are large discrepancies between existing analytical theories and test particle simulations of heating rates. They warrant further investigations that may include molecular dynamics simulations dufty ; batishchev . Classical models of electron-ion interactions in molecular dynamics models include an effective potential uhlenbeck ; deutsch  to prevent the collapse of such system due to electrostatic attraction and properly account for quantum diffraction effects at short distances. We will also perform test particle simulations with this effective potential and will compare them with Coulomb potential scattering studies.

The paper is organized as follows. In Sec. II we summarize existing theories of inverse bremsstrahlung heating and present three representative expressions. Section III presents a calculation model and gives examples of scattering orbits. Section IV describes the case of mono-energetic, beam-like electron distribution functions. We discuss different wave polarizations and field strengths. Heating rates for the case of a Maxwellian electron distribution are calculated in Section V. Section VI deals with anisotropic electron distributions that are encountered in photo-ionized plasmas. Section VII contains a summary and conclusions.

## II Description of the energy exchange in the laser field

In the introduction, we have identified several groups of analytical results describing electron-ion collisions in the presence of a homogeneous oscillating electric field. This division is a convenient characterization for various formalisms that are used to calculate the heating rate for electrons. All theories lead to remarkably similar heating rates, particularly for energetic electrons and small angle collisions, where electron scattering is well described by the Born approximation.

The Dawson-Oberman model do ; pert2 ; decker ; shvets  of the electron-ion correlation function is based on the classical Born approximation of binary collisions and on an explicit average over random ion positions. The Coulomb field of the immobile ion of the charge $`Ze`$ is considered as a weak perturbation to the electron motion which is comprised of a drift velocity $`𝐯_0`$ along a straight trajectory and a quiver velocity $`𝐯_E=e𝐄/m\omega `$ in a laser field of amplitude $`𝐄`$ abd the frequency $`\omega `$. Here $`e`$ and $`m`$ are the charge and the mass of the electron. The Dawson-Oberman collision operator describes electron-ion scattering as an instantaneous event resulting in a small angle deflection of the electron trajectory and accounts for the screening effect of other electrons and the high frequency field in the dielectric approximation. The electron energy gain has been found in the weak field limit do , and also for an arbitrary strength of the oscillating field decker  with a Maxwellian distribution of electrons do ; decker  or with an arbitrary distribution function pert2 ; shvets . Derivations by Pert pert2  and Shvets and Fisch shvets  arrive at similar expressions starting from the test particle, two body scattering problem. The energy gain for a given electron in a plasma with the ion density $`n_i`$ has the following form pert2 ; shvets :

$$\frac{dϵ}{dt}=\frac{2n_iZ^2e^4}{mv_0^2}\underset{l_0}{\overset{+\mathrm{}}{}}l\frac{d𝐤}{k^2}J_l^2\left(\frac{\mathrm{𝐤𝐯}_E}{\omega }\right)𝐯_0\frac{}{𝐤}\delta \left(\frac{𝐤𝐯_0}{\omega }l\right),$$
(1)
where $`J_l`$ is the Bessel function and the integration over the transferred momentum $`𝐤`$ involves the upper cut-off limit at $`k_{max}=mv_0^2/2Ze^2`$, which corresponds to the assumption of small angle scattering.

Silin silin  has calculated the high frequency nonlinear conductivity and the effective collision frequency for a fully ionized plasma by using a kinetic equation with a collision integral which accounts for small-angle scattering. In this approach, the high frequency electric field defines the cut-off limit of the impact parameter and the effective collision velocity, but does not affect the collision event itself. The Silin’s collision frequency is in close agreement with the Dawson-Oberman model decker  in both limits of weak and strong amplitude high frequency fields. A similar result for the electron energy gain has been obtained from the solution to the electron kinetic equation with the Landau collision operator after averaging the result with respect to the high frequency oscillations by Catto catto .

Several classical models of the effective collision frequency have been derived from explicit treatment of the dynamics of electron-ion scattering assuming instantaneous and elastic interactions bunkin ; pert1 ; mulser . We will identify them as the impact approximation pert1 ; bunkin , or the ballistic model mulser . From the conservation of energy for such collisions, one finds that the change in the average electron kinetic energy, which is related to the electron drift velocity, is proportional to the change in the electron momentum. The energy gain is averaged with respect to the phase of laser field at the time of the electron-ion collision and expressed in terms of the transport cross section:

$$\frac{dϵ}{dt}=\frac{2n_iZ^2e^4}{m}_0^{2\pi }𝑑\varphi \frac{\mathrm{\Lambda }(𝐯_E𝐯_0\mathrm{cos}\varphi +v_E^2\mathrm{cos}^2\varphi )}{(v_0^2+2𝐯_E𝐯_0\mathrm{cos}\varphi +v_E^2\mathrm{cos}^2\varphi )^{3/2}}.$$
(2)
Here, the Coulomb logarithm

$$\mathrm{\Lambda }=\mathrm{ln}(1+\rho _{max}^2/\rho _{min}^2)^{1/2}$$
(3)
depends on the full electron velocity $`𝐯=𝐯_0+𝐯_E\mathrm{cos}\varphi `$ mulser ; pert3  and $`\rho _{max}=v/\omega `$ and $`\rho _{min}=\mathrm{max}\{2Ze^2/mv^2,\mathrm{}/mv\}`$ are the classical cut-offs at large and small impact parameters. According to Eq. (2), there is a singularity in the integrand if the quiver velocity, $`v_E`$, approaches the drift velocity, $`v_0`$. This singularity corresponds to orbits passing close to the ion and first appears for electrons launched in the direction of the electric field. Correlated collisions in presence of high frequency electric field have been observed in numerical simulations of particle trajectories fraiman ; num1  for the parallel launch, $`𝐯_0𝐯_E`$. They remove this singularity and lead to an electron energy gain that is much larger than one that follows from Eq. (2), if $`v_0<v_E`$.

A quantum-mechanical description of electron-ion collisions in the Coulomb and laser fields employing the low frequency approximation has been developed by Bunkin and Fedorov bunkin  and Kroll and Watson kw . Agreement between the quantum mechanical and classical derivations of the electron heating rate has been demonstrated by Ferrante et al. ferrante . Within the theoretical framework of the Born approximation, this is a consequence of the well-known fact landau3  that the electron scattering cross section in the Coulomb field is the same in classical and quantum-mechanical calculations. The energy gain averaged with respect to the Maxwellian electron velocity distribution function schlessinger  also shows a good agreement between classical and quantum results. The recent quantum theory of the inverse Bremsstrahlung heating rate kull , based on the dielectric approximation, compares well with classical results, provided the quantum cut-off parameters are introduced into the classical kinetic theory.

Quantum mechanical calculations bunkin ; schlessinger  describe the energy exchange rate between an electron and a laser field by using the cross sections derived by Kroll and Watson kw  for electron-ion scattering. The result is given by

$`{\displaystyle \frac{dϵ}{dt}}`$ $`=`$ $`{\displaystyle \frac{n_iZ^2e^4}{2mv_0}}{\displaystyle \underset{l=l_{min}}{\overset{+\mathrm{}}{}}}{\displaystyle }d𝐧^{}{\displaystyle \frac{\xi l(1+\xi l)^{1/2}}{[1+\xi l/2(1+\xi l)^{1/2}𝐧𝐧^{}]^2}}\times `$ (4)
$`J_l^2\left({\displaystyle \frac{2𝐯_E}{\xi v_0}}[(1+\xi l)^{1/2}𝐧^{}𝐧]\right),`$
where $`𝐧`$ and $`𝐧^{}`$ refer to the propagation directions of the electron before and after collision, respectively, $`\xi =2\mathrm{}\omega /mv_0^2`$, $`v_0`$ is the initial drift velocity, and $`l_{min}=1/\xi `$ is the maximum number of emitted photons. Energy exchange in this equation is expressed in terms of the differences between the total absorption cross section (the terms with $`l>0`$ in the sum) and emission (the terms with $`l<0`$) of $`l`$ photons.

## III Computational model

Our simulations of non-relativistic electron scattering in the field of an ion and in a uniform laser field are similar to studies of other authors fraiman ; num1 . We solve Newton’s equation for the motion of a test electron

$$m\ddot{𝐫}=e𝐄\mathrm{sin}\omega tU(r),$$
(5)
where $`U(r)`$ is the electrostatic potential describing the electron-ion interaction calculated along the electron trajectory, $`𝐫(t)`$. The electron velocity $`𝐮=\dot{𝐫}`$ is divided into two components: the quiver velocity, with amplitude $`𝐯_E`$, and the drift velocity $`𝐯`$ and therefore $`𝐮=𝐯𝐯_E\mathrm{cos}\omega t`$. In the absence of an oscillatory field and for the Coulomb potential, $`U(r)=Ze^2/r`$, the electron follows the Kepler orbit that is uniquely defined by its initial drift velocity $`𝐯_0`$ and the impact parameter $`\rho `$. The scattering process for the momentum exchange is described by the Rutherford cross-section, which corresponds to the surface area of the interaction radius $`\rho _0=Ze^2/mv_0^2`$ multiplied by a logarithmic factor which accounts for small angle collisions. The limits in the logarithm are defined by the Debye screening length $`\lambda _D=\sqrt{T_e/4\pi e^2n_e}`$ at large distances (here $`T_e`$ is the electron temperature) and $`\rho _0`$. The small impact parameters correspond to large angle scattering events, including head-on electron-ion collisions, which contribute to the Rutherford cross section with finite probabilities. Similar properties describe electron-ion scattering in a weak oscillatory field $`v_E/v_01`$, the only exception being that the long distance cut-off is defined by the length $`v_0/\omega `$ if it is smaller than the Debye length.

The classical treatment of collisional absorption and particle heating gives a result, which within logarithmic accuracy, agrees with quantum mechanical calculations in the Born approximation schlessinger ; kull . The exact agreement between the quantum scattering problem and the Kepler orbit solutions is well known and is a unique feature of the Coulomb potential landau3 . In numerical studies based on classical mechanics such as molecular dynamics dufty  or test particle simulations, the quantum diffraction effect at short distances can be introduced by imposing a cut-off in the Coulomb potential deutsch  at the electron de Broglie length, $`\lambda _B=2\pi \mathrm{}/mv`$. This cut-off was first suggested by Uhlenbeck uhlenbeck . It has allowed the successful application of classical theories to weakly degenerate systems ichimaru1 . For the simplest problem of electron scattering off a bare ion, no short range cut-off is necessary in view of the agreement between quantum and classical calculations. No such general agreement exists, however, when the oscillatory field modifies the scattering orbits. For example, the quantum theory of inverse bremsstrahlung in the Born approximation kull  produces limiting expressions with logarithmic terms involving a cut-off at $`\lambda _B`$. Also, there is no quantum theory going beyond the Born approximation comment . Thus the problem of quantum corrections to the effective interaction potential is a difficult one and has no clear solution so far.

For these reasons, the potential used in our simulations accounts for cut-offs at both short and large distances:

$$U(r)=\frac{Ze^2}{r}\left[1\mathrm{exp}\left(\frac{r}{\lambda _B}\right)\right]\mathrm{exp}\left(\frac{r}{\lambda _D}\right).$$
(6)
The two lengths, $`\lambda _B`$ and $`\lambda _D`$, are considered as free parameters that allow one to make contact with known limiting cases. Figure 1 compares the effective potential (6) with the Coulomb potential. The characteristic spatial scale of the problem is the distance $`r_C=\sqrt{Ze/E}`$ fraiman , where the strength of the Coulomb field is the same of the oscillating field. Here, the Debye length was calculated for an electron temperature of 1 keV and an electron density $`10^{20}`$ cm-3. The de Broglie length, $`\lambda _B=2\pi \mathrm{}/mv_E`$, was calculated at the quiver velocity $`v_E`$ corresponding to a laser intensity of $`10^{16}`$ W/cm2 and a wavelength of 0.25 $`\mu `$m. In our study, we examine, within the classical formulation, scattering orbits that clearly display non-perturbative modifications to the Kepler trajectories. We proceed with two dynamical models involving the Coulomb potential, as in fraiman , and the effective potential (6) with $`\lambda _B`$ calculated by using the quiver electron velocity. The upper cut-off at $`\lambda _D`$ is not important, because the correlated collisions occur for relatively small impact parameters and do not demonstrate the logarithmic divergence at large distances which is the characteristic feature of scattering in the Coulomb potential.

Numerical solutions of the equations of motion, Eq. (5), require special attention regarding the choice of a suitable time step in computations of electron orbits near the scattering center. We have used the symplectic integration algorithm candy , which ensures a good numerical accuracy and an exact conservation of dynamical invariants. Even though the test particle problem described by Eq. (5) appears simple, the construction of the statistical ensemble of results for a broad range of initial conditions for the average energy gain poses considerable challenges in terms of computer memory and time requirements. The overall geometry of our simulations is shown in Fig. 2. Electrons were launched from the plane of incidence with a prescribed impact parameter $`\rho `$ and an initial drift velocity $`v_0`$ which was always normal to the plane. The interaction sphere was defined by the screening length, $`\lambda _D`$. The electron trajectories depend on the angle of $`𝐯_0`$ with respect to the high frequency field polarization and on the field phase.

The relevance of individual scattering trajectories in calculations of the average energy gain follows from the validity of the two-body collision approximation. This is also the basic assumption of all theoretical results used here for comparison. The two-body scattering model requires that the given electron can only interact with one particular ion during the scattering event. Therefore, the amplitude of electron oscillations, $`r_E=eE/m\omega ^2`$, in the laser field must be smaller than the average distance between ions $`d=(3/4\pi n_i)^{1/3}`$. This condition, $`r_E<d`$, reads

$$5.7\times 10^{14}\lambda _0^2\sqrt{I_0}n_i^{1/3}<1,$$
(7)
where the laser wavelength $`\lambda _0`$ is in $`\mu `$m, the intensity $`I_0`$ is in W/cm2 and the ion density is in cm-3. For example, a laser pulse of intensity $`I_010^{17}`$ W/cm2 and wavelength $`\lambda _0=0.25\mu `$m satisfies the condition (7) in a gaseous plasma of $`n_i=10^{18}`$ cm-3.

The scattering trajectories that give rise to an enhancement of the electron energy gain fraiman  are illustrated in Fig. 3. Figure 3a compares three cases with the same initial conditions: a hyperbolic orbit in the absence of an oscillatory field (dotted line), a trajectory in a Coulomb field and in an oscillatory field (dashed line) and a solution to the equation of motion with an effective potential (6) and an oscillatory field (continuous line). In this example, the oscillatory field changes the trajectory and causes a head-on collision with the ion. This dramatic deflection is moderated by the quantum diffraction effect accounted for in the effective potential (6). The single trajectory of Fig. 3b corresponds to the case where the initial drift velocity is parallel to the direction of the oscillatory field. It illustrates the case of quasi-capture for the motion in the effective potential. The time dependent field reduces the instantaneous electron energy to a negative value allowing for a complicate bounded motion in the vicinity of an ion. The third example, Fig. 3c, displays two trajectories of electrons launched in the transverse direction to the oscillatory field with a small difference in the initial phase. At first, both orbits are very close to each other, yet after the collision they diverge dramatically. This sensitivity to phase and a strong modification of the scattering trajectory after an almost head-on collision, are manifestations of the effect of Coulomb focusing described in Refs. corkum0 ; brabec ; ivanov .

Our calculations focus on the energy gain of a single electron $`ϵ`$ for a given impact parameter averaged over the field phase and the rate of energy gain $`dϵ/dt`$ for a single electron averaged over the impact parameters and the field phase. These quantities were calculated directly from the simulation results according to the following expressions:

$$ϵ=\frac{m}{2}v_0^2v_0^2,\frac{dϵ}{dt}=n_iv_0ϵd^2\rho ,$$
(8)
where $`𝐯_0`$ and $`𝐯_0^{}`$ are the electron drift velocities before and after the collision (cf. Fig. 2) and the angular brackets denote the average over the field phase. The electron energy gain from the oscillatory electric field (8) is compared with known analytical and numerical results.

In the following three sections we describe the results of our test particle simulations of the electron energy gain for different initial electron distribution functions. We start with the beam-like mono-energetic electron distribution function, followed by an isotropic Maxwellian distribution and finally with an anisotropic electron distribution function.

## IV Energy gain for a mono-energetic electron beam

Here, we analyse the energy exchange between the oscillatory field and electrons of a given initial velocity. The ensemble of different trajectories are constructed by varying the initial phase and the impact parameter. In particular, the parallel and perpendicular polarization of the field with respect to the initial direction of the drift velocity are considered.

### IV.1 Analytical expressions

The geometry of a parallel field polarization allows further simplification of the expressions for the electron energy gain obtained by Dawson and Oberman, Eq. (1):

$`{\displaystyle \frac{dϵ}{dt}}={\displaystyle \frac{4\pi n_iZ^2e^4}{mv_0}}{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}\{{\displaystyle \frac{2J_l^2(lv_E/v_0)}{1+l^2\mathrm{}^2\omega ^2/m^2v_0^4}}`$
$`{\displaystyle \frac{lv_E}{v_0}}J_l\left({\displaystyle \frac{lv_E}{v_0}}\right)[J_{l1}\left({\displaystyle \frac{lv_E}{v_0}}\right)J_{l+1}\left({\displaystyle \frac{lv_E}{v_0}}\right)]\times `$
$`\mathrm{ln}(1+{\displaystyle \frac{m^2v_0^4}{l^2\mathrm{}^2\omega ^2}})\},`$ (9)
and Kroll and Watson, Eq. (4):

$$\frac{dϵ}{dt}=\frac{\pi n_iZ^2e^4}{mv_0}\underset{l=l_{min}}{\overset{+\mathrm{}}{}}\underset{1\sqrt{1+\xi l}}{\overset{1+\sqrt{1+\xi l}}{}}\frac{dy\xi l}{(y+\xi l/2)^2}J_l^2\left(\frac{2v_Ey}{v_0\xi }\right).$$
(10)
Equation (IV.1) follows from Eq. (1) after integration over the transfered momentum $`𝐤`$ with the upper limit cut-off at $`k_{max}=mv_0/\mathrm{}`$, which corresponds to quantum short distance cut-off at the de Broglie length calculated at the initial electron drift velocity. The parameter $`\rho _{min}=\mathrm{}/mv_0`$. A more detailed study of the quantum mechanical expression (10) can be found in Ref. daniel , where it is shown that an enhancement of the total cross-section occurs for $`v_0v_E`$ and an oscillatory behavior is observed for $`v_0<v_E`$. Note, that in the classical limit, $`\xi l1`$, Eq. (10) gives an expression which is similar to Eq. (IV.1).

In the limit of a weak electric field, where $`v_E/v_0\xi 1,`$ one recovers from Eq. (10) the well known expression

$$\frac{dϵ}{dt}=\frac{4\pi n_iZ^2e^4v_E^2}{mv_0^3}\mathrm{ln}\left(\frac{2mv_0^2}{\mathrm{}\omega }\right)$$
(11)
which accounts for the quantum cut-off at small distances and the high frequency cut-off, $`\rho _{max}=v_0/\omega `$, at large distances.

There is no easy way of simplifying theoretical results for the perpendicular polarization of the field. We use the original expressions (1) and (4) for the comparisons with our numerical results.

### IV.2 Simulation results

The results of numerical simulations are illustrated by Figs. 4, 5, and 6. They display the average electron energy gain, $`dϵ/dt`$ (8), as a function of initial drift velocity normalized to $`v_E`$. Different $`v_0/v_E`$ values along horizontal axis were achieved by changing the initial drift velocity and keeping $`v_E`$ constant. In addition to $`v_0/v_E`$, the results also depend on the strength of the electric field. We characterize this dependence by using the parameter $`r_E/\rho _E`$, where $`\rho _E=Ze^2/mv_E^2=\rho _0v_0^2/v_E^2`$ is the distance for a large angle scattering event for an electron with the velocity $`v_E`$. It replaces $`\rho _0`$ in the definition of scattering cross-section for trajectories with $`v_0/v_E<1`$. The parameter $`r_E/\rho _E`$ can be conveniently presented as

$$r_E/\rho _E=3.4\times 10^{20}I_0^{3/2}\lambda _0^4/Z$$
(12)
where the laser intensity is in W/cm2 and the wavelength is in $`\mu `$m. It is large for typical conditions of present laser-plasma interaction experiments. However, the limit of small values of this parameter is also instructive to consider.

Figure 4 corresponds to $`r_E/\rho _E1`$ while Figs. 5, 6 describe the high intensity limit, $`r_E\rho _E`$. As can be seen from Figs. 3, that are also drawn for $`r_E/\rho _E1`$, scattering events of electrons experiencing large amplitude oscillations in the electric field are very sensitive to phase and are likely to result in large energy gains but only for a selected number of particles. Figure 5 shows numerical results for the energy gain (8) with the Coulomb potential (large dots) and the effective potential $`U(r)`$ (6) (continuous line) for the transverse launch ($`𝐯_0𝐯_E`$) case. They are compared with results of the Dawson-Oberman theory (1) (dashed line), the classical model (2) (dotted line), and the Kroll-Watson theory (4) (dash-dotted line). The energy gain (8) for the parallel launch ($`𝐯_0𝐯_E`$) is compared with results of three theoretical models, (2), (IV.1), and (10) in Fig. 6.

Two main features characterize the numerical results in Figs. 4, 5, and 6: an agreement between theoretical and numerical results for high velocities, $`v_0v_E`$, and increasing values of the energy gain for small initial velocities, $`v_0v_E`$ which is in strong disagreement with all theoretical predictions. An agreement at large drift velocities follows from the applicability of the Born approximation to these trajectories. All curves in Figs. 4, 5, and 6 are well approximated by a classical result predicting the decrease of the energy gain as $`v_0^3`$. At small drift velocities, the calculated energy gain becomes larger. The increase is more dramatic for large fields, $`r_E/\rho _E1`$, and for trajectories calculated with the Coulomb potential (large dots) as compared to the effective potential (6) (solid lines) in Figs. 5 and 6.

### IV.3 Weak field limit

We identify a weak field case by $`r_E/\rho _E1`$. According to Eq. (12), this corresponds to relatively small intensities, $`I_010^{13}Z^{2/3}\lambda _0^{8/3}`$ W/cm2 where the laser wavelength is in $`\mu `$m. The energy gain rate is shown in Fig. 4. Electrons that are launched in either the parallel or perpendicular directions with respect to the oscillatory electric field lead to similar results. In both cases we observe an enhancement of the energy gain at small velocities.

In order to interpret these numerical results, we first recall the electron scattering theory in the Coulomb field. The particle orbit is uniquely defined by the initial drift velocity, $`v_0`$, and the impact parameter, $`\rho `$. The shortest distance along the Kepler trajectory that separates the electron from the ion is $`r_{min}=\sqrt{\rho _0^2+\rho ^2}\rho _0`$. In estimates which follow, we assume that the electron can gain a large amount of energy, on the order of $`mv_E^2/2`$, if it approaches the ion at least as close as $`\rho _E`$. Then, by assuming $`r_{min}=\rho _E`$, we find the characteristic impact parameter for these events $`\rho _{int}=\sqrt{\rho _E^2+2\rho _E\rho _0}\sqrt{2\rho _E\rho _0}`$ for $`v_E>v_0`$.

During the motion inside the interaction sphere, the electron oscillates many times with a relatively small amplitude, $`r_E`$, and each time it undergoes a small angle deflection from the original hyperbolic orbit. In an analogy with the Coulomb focusing process corkum0 , each such deflection brings the electron closer to the ion and eventually the electron undergoes a large angle scattering which terminates the interaction. As this trajectory modification involves many oscillation periods, it depends weakly on the initial phase. One could estimate the energy gain rate for the case $`v_E>v_0`$ by taking the product of the number of scattering events per unit time $`n_iv_0\pi \rho _{int}^2`$ and the gained energy $`mv_E^2`$. That gives the following estimate:

$$dϵ/dt2\pi n_iZ^2e^4/mv_0.$$
(13)
The $`1/v_0`$ dependence results in large values of the energy gain at small velocities in contrast to the previous theoretical results which do not depend on $`v_0`$: $`dϵ/dt4\pi Z^2n_ie^4\mathrm{\Lambda }/mv_E`$ do . The correlation effects dominate the energy gain in the limit of small electron velocities, $`v_0<v_E/\mathrm{\Lambda }`$.

The impact of this mechanism in the enhancement of the electron energy gain is greatly reduced for small intensities where $`r_E/\rho _E1`$. Results of Fig. 4 have been obtained for the marginal case of $`r_E\rho _E`$. Still this gain enhancement is quite important as it is independent of the field polarization and the initial phase. That makes the weak field limit distinct from the regime $`r_E/\rho _E1`$ which is discussed in two next subsections.

### IV.4 Strong field, perpendicular launch

The interpretation of high intensity results ($`r_E/\rho _E1`$) is more complicated because of the strong dependence on the field phase. We first consider the energy gain of electrons launched transversely to the electric field direction, Fig. 5. Examples of relevant trajectories are presented in Fig. 3c. The energy gain is characterized by two impact parameters: $`\rho _x`$ – in the direction the oscillating field and $`\rho _y`$ – in the direction of the initial drift velocity. As one can see in Fig. 7, the energy gain is an asymmetric function: it exists over a wide range of impact parameters $`\rho _xr_E`$, while it is limited to much narrower range $`\rho _y\rho _{int}`$ in the perpendicular direction.

Modifications of electron orbits that lead to an enhancement of energy gain also depends on the field phase $`\varphi `$. This is shown in Fig. 8, which presents the dependence of the energy exchange on $`\varphi `$ for two initial drift velocities: the high velocity, $`v_0=3v_E`$ (solid line) and the slow initial motion, $`v_0=0.3v_E`$ (dashed line). Assuming that incident electrons have a uniform phase distribution, one concludes that only a small fraction of electrons $`\rho _E/(v_0/\omega )`$ interact strongly with the ion while passing it within the distance of $`\rho _E`$. A change in the electron energy can be negative or positive. This depends on the direction of electron motion, either towards or away from ion when it enters the sphere of radius $`\rho _E`$. For fast electrons, the energy loss for particular phases turns out to be exactly the same as the gain for other phases with no net change to the Born approximation results. For slow particles (dashed line in Fig. 8) the reduction of electron energy can result in quasi-trapped trajectories. This motion is irregular and produces a net energy gain when the electron finally escapes the ion. This effect is more important for slower electrons as it leads to increasing energy gain at small drift velocities as shown in Fig. 5.

An estimate of the energy gain rate $`dϵ/dt`$ involves the characteristic cross section, $`\sigma =\pi \rho _x\rho _y`$, and the characteristic energy gain $`mv_E^2`$. The impact parameter, $`\rho _x`$, along the field direction can be found by assuming that the upper bound for the closest approach and large angle scattering is $`r_E`$, that is, $`r_E=\sqrt{\rho _0^2+\rho _x^2}\rho _0`$ leading to the estimate $`\rho _x=r_E\sqrt{1+2\rho _0/r_E}`$. As discussed above, one have also to account for the fact, that only a small fraction of electrons, $`\rho _E/(v_{int}/\omega )`$, undergo collisions with this impact parameter. Here, the effective drift velocity $`v_{int}=v_0\sqrt{1+2\rho _0/r_E}`$ is calculated at the distance $`r_E`$ from the ion. By including these results in the expression for the energy gain, one obtains

$`{\displaystyle \frac{dϵ}{dt}}`$ $``$ $`\pi n_iv_0(2\rho _0\rho _E)^{1/2}r_E(1+2\rho _0/r_E)^{1/2}\times `$ (14)
$`{\displaystyle \frac{\rho _E\omega }{v_0(1+2\rho _0/r_E)^{1/2}}}{\displaystyle \frac{\pi n_iZ^2e^4}{mv_0}}.`$
This estimate is very similar to Eq. (13) for the energy gain in the weak field limit. The $`1/v_0`$ dependence explains the growth of the energy gain at small velocities. These are the large angle collisions of selected electrons that contribute to this enhancement. This is moderated by quantum diffraction effects, which are accounted for in the expression for the effective potential (solid line in Fig. 5) if the interaction radius $`\rho _{int}`$ becomes smaller than the de Broglie length.

### IV.5 Strong field, parallel launch

The perpendicular launch considered above represents a general case – electrons entering the interaction sphere at arbitrary angles with respect to the polarization vector of the oscillatory field are gaining on average the same amount of energy as those launched in the perpendicular direction. Only the electrons propagating parallel to the polarization vector make an exception. As shown in Fig. 6, these electrons are gaining much larger energies at small velocities as compared to the previous cases. The explanation of this anomalous energy gain is related to Coulomb focusing corkum0  or the ”parachute effect” fraiman ; fraimannew . Electrons launched in the parallel direction can return to the same ion several times. Each time, as the electron passes the ion, its trajectory is deflected towards the ion. The effective impact parameter decreases until a head-on collision terminates the interaction. The number of possible returns scales as $`v_E/v_0`$.

This scattering geometry has been analyzed and described in detail by Fraiman et al. fraiman ; fraimannew . A typical dependence of the energy gain averaged over the field phase $`ϵ`$ (8) on the impact parameter for the parallel launch is shown in Fig. 9 for $`v_0=0.3v_E`$. It displays three maxima in the energy exchange corresponding to three electron returns. The second curve (solid line) in Fig. 9 has been obtained with the effective potential (6). The amount of energy gained by electrons depends on the short range behavior of the interaction potential at $`r\rho _0`$ and is reduced when close distance interactions are cut-off (in the present example $`\lambda _B>\rho _0`$). Figure 10 illustrates the sensitivity of the electron energy gain to the field phase and the relative importance of head-on collisions. The latter are represented by the large maximum near $`\varphi =3\pi /2`$. The effect of quasi-capture and random trajectories do not strongly contribute to the energy gain in this case.

An estimate of the heating rate can be made by following the arguments of Fraiman et al. fraiman ; fraimannew . When passing by an ion, an electron with an impact parameter $`\rho `$ is scattered by a small angle $`\delta \theta 2\rho _E/\rho `$ in the Coulomb potential. This deflection changes the original orbit and the electron acquires a transverse velocity component $`v_{}=\delta \theta v_E`$. That electron will arrive into the head-on collision with the ion after half a field period if $`\rho =\pi v_{}/\omega `$. From the above expression one finds $`\rho =\sqrt{2\rho _Er_E}\sqrt{2}r_C`$. The electrons with larger impact parameters will also be attracted to the ion but only after several oscillations. Since the number of possible returns is $`v_E/v_0`$, the characteristic impact parameter for this strong interaction is estimated as $`r_Cv_E/v_0`$. Then the average electron energy gain for the parallel launch can be evaluated as

$$\frac{dϵ}{dt}2\pi n_iv_0\frac{\rho _E}{r_C}r_C^2\frac{v_E^2}{v_0^2}mv_E^2=2\pi n_i\frac{Z^2e^4}{mv_0}\left(\frac{r_E}{\rho _E}\right)^{1/2},$$
(15)
where we have also accounted for the fact that only fraction of electrons proportional to $`\rho _E/r_C`$ will undergo strong interactions (cf. Ref. fraiman ). Although the formula (15) demonstrates the same dependence on the electron drift velocity as (14), the energy gain is $`(r_E/\rho _E)^{1/2}1`$ times bigger due to the parachute effect discussed above. This Coulomb attraction due to multiple collisions is effective, however, only for small angles between $`𝐯_0`$ and the electric field vector, $`\theta \sqrt{\rho _E/r_E}`$, where the velocity $`v_{}`$ has a well defined sign. For larger angles, the electron-ion interaction reduces to a single collision event, as presented previously for the case of perpendicular launch.

## V Heating rate for an isotropic electron distribution

The energy gain rates (14) and (15) calculated for mono-energetic electrons can be used for evaluating heating rates for various electron distribution functions. Here we consider the common case of an isotropic Maxwellian electron distribution function which is established due to electron-electron and electron-ion collisions. Before analyzing the statistical significance of correlated collisions, we recall the basic elements of the classical theories.

By taking an average with respect to a Maxwellian distribution function with the temperature $`T_e`$, the Dawson-Oberman theory (1) gives the following formula for the heating rate per electron:

$$\frac{dT_e}{dt}=\sqrt{2\pi }\frac{8n_iZ^2e^4v_E^2}{mv_{Te}^3}\underset{l=1}{\overset{+\mathrm{}}{}}l^2_0^{q_m}\frac{dq}{q^4}I_l(q)\mathrm{exp}\left(\frac{l^2mv_E^2}{2q^2T_e}\right),$$
(16)
where $`I_l(x)=_0^xJ_l^2(y)𝑑y`$, $`v_{Te}=\sqrt{T_e/m}`$ is the electron thermal velocity, and the parameter $`q_m=2\mathrm{}\omega /T_e`$ defines the cut-off at small distances. The quantum mechanical theory (4) gives the following heating rate

$`{\displaystyle \frac{dT_e}{dt}}`$ $`=`$ $`\sqrt{2\pi }{\displaystyle \frac{16n_iZ^2e^4v_E^2}{\mathrm{}\omega v_{Te}}}{\displaystyle \underset{l=1}{\overset{+\mathrm{}}{}}}l\mathrm{sinh}{\displaystyle \frac{l\mathrm{}\omega }{2T_e}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dq}{q^4}}I_l(q)\times `$ (17)
$`\mathrm{exp}\left({\displaystyle \frac{l^2mv_E^2}{2q^2T_e}}{\displaystyle \frac{q^2\mathrm{}^2\omega ^2}{8mv_E^2T_e}}\right),`$
One can define the effective collision frequency by the relation: $`dT_e/dt=\nu _{eff}mv_E^2/2`$ and use Eqs. (16) and (17) to evaluate it. It is well-known that in the weak field limit, $`v_Ev_{Te}`$, the effective collision frequency coincides with the electron-ion transport collision frequency $`\nu _{ei}=4\sqrt{2\pi }Z^2e^4n_i\mathrm{\Lambda }/3m^2v_{Te}^3`$, where the Coulomb logarithm $`\mathrm{\Lambda }`$ is defined by Eq. (3). In the strong field limit, $`v_Ev_{Te}`$, $`\nu _{eff}`$ decreases as $`v_E^3`$ and there is a simple expression from Ref. vick ; faehl  which provides a very good interpolation for all existing analytical theories:

$$\nu _{eff}=\nu _{ei}(1+v_E^2/6v_{Te}^2)^{3/2}.$$
(18)
Figure 11a compares the variation of effective collision frequency with $`v_E/v_{Te}`$ found from our simulations and from theoretical formulas presented above. The numerical heating rates were calculated by the direct average of the energy gain rate over the isotropic Maxwellian distribution function. Numerical results have been obtained for the cases when the electron temperature is 200 eV and also for 500 eV. All models and numerical simulations generally agree for $`v_Ev_{Te}`$. Some discrepancies at small values of $`v_E`$ follow from differences between theoretical approximations, in particular, from the treatment of short range interactions.

The most dramatic difference between the known theories and our results is illustrated in Fig. 11b which shows that the calculated heating rate reaches a constant values in the large field limit, $`v_E>v_{Te}`$. Our results for the effective collision frequency can be approximated by the following formula:

$$\nu _{eff}=\nu _{ei}(1+0.3v_E^2/v_{Te}^2)^1.$$
(19)
This expression gives a $`1/v_{Te}`$ dependence in the heating rate for the limit of large values of $`v_E/v_{Te}`$. This agrees with the conclusions of Ref. fraiman . Such behavior agrees with the $`1/v_0`$ dependence of the energy gain rate for the perpendicular electron launch discussed in sections IV.3 and IV.4. The contribution of electrons with parallel velocities is not evident in Eq. (19) because small angles, $`\theta \sqrt{\rho _E/r_E}`$, make a small statistical contribution to the average quantity for the isotropic distribution function.

We recall that the enhancement of the heating rate due to the correlation effect takes place in a sufficiently rarified plasma where the electron quiver radius $`r_E`$ is smaller than the inter-particle distance $`d`$. The condition (7) ensures this limitation.

## VI Heating rates for an anisotropic electron distribution

It was shown in Sec. IV.5 that slow electrons ($`v_0<v_E`$) propagating along the oscillatory field are the main contributors to the enhanced heating rates. However this effect is confined to small angles and it disappears if electrons are distributed isotropically. Therefore for typical laser plasma interaction conditions, the heating rate might depend dramatically on the number of slow electrons and on the anisotropy of actual electron distribution function. In particular, this effect could be important for plasmas produced by the ionization of gases by high intensity laser pulses. According to Refs. corkum1 ; bychenk , the velocity distribution function of the photoelectrons displays a strong anisotropy along the direction of the linearly polarized laser beam, if tunnel ionization is the dominant mechanism. The following simple analytical expression captures the main features of electrons produced by the tunnel ionization process bychenk :

$$F_0(𝐯)=n_e(m/2\pi T_e)^{1/2}\delta (𝐯_{})\mathrm{exp}(mv_{}^2/2T_e)$$
(20)
where $`v_{}`$ is the electron velocity in the direction of laser polarization and the effective electron temperature $`T_e`$ may be a function of the ionization potential and the laser field amplitude bychenk . In our analysis we treat $`T_e`$ as a free parameter.

We compare results of various analytical models with our simulations after averaging them with respect to the distribution function (20). This averaging reduces to a simple integration of expressions for the electron energy gain in the parallel geometry (IV.1) and (10) with a one-dimensional Maxwellian velocity distribution. The main contribution to this integral comes from slow electrons. In the case where $`v_E<v_{Te}`$, all classical theories give similar results that depend on the cut-off parameters. For the quantum approach, the heating rate per electron derived from Eq. (10) has the following form:

$`{\displaystyle \frac{dT_e}{dt}}`$ $`=`$ $`n_iZ^2e^4\sqrt{{\displaystyle \frac{\pi m}{2T_e}}}{\displaystyle \underset{l=l_{min}}{\overset{+\mathrm{}}{}}}l\xi _E{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}du\mathrm{exp}({\displaystyle \frac{mv_E^2}{2T_e}}u^2)\times `$ (21)
$`{\displaystyle \underset{u\sqrt{u^2+\xi _El}}{\overset{u+\sqrt{u^2+\xi _El}}{}}}𝑑y{\displaystyle \frac{J_l^2(2y/\xi _E)}{(yu+\xi _El/2)^2}},`$
where $`\xi _E=2\mathrm{}\omega /mv_E^2`$.

Figure 12 shows the intensity dependence of the heating rates that are obtained from our numerical model and from various analytical expressions. The striking discrepancy between the analytical results and our numerical simulations is a consequence of the large number of slow electrons in the electron distribution function (20) propagating parallel to the polarization direction. Contrary to the case of a three dimensional Maxwellian distribution function, expression (20) contains a finite number of electrons at zero velocity and therefore the heating rate diverges as a logarithm of $`v_{min}`$, if one averages the energy gain (15) with the one-dimensional electron distribution function. This divergence can be resolved by accounting for the fact that the maximum impact parameter $`\rho _{max}r_Cv_E/v_0`$ cannot be greater than the inter-particle distance $`d`$. This leads to the following estimate $`v_{min}v_Er_C/d`$. Then one finds the following expression for the heating rate per electron:

$$\frac{dT_e}{dt}=\frac{\sqrt{2\pi }n_iZ^2e^4}{mv_{Te}}\left(\frac{r_E}{\rho _E}\right)^{1/2}\mathrm{ln}\left(d\sqrt{\frac{T_e}{Ze^2r_E}}\right).$$
(22)
This expression is proportional to $`v_E^{3/2}`$. It approximates well numerical results as shown with the solid line in Fig. 12.

Photo-ionized plasmas with anisotropic electron distributions provide an example of a case where enhanced heating rates are obtained due to correlated collisions as described in our paper. It is clear that the particular form (20) of the electron distribution function can be quickly altered due to rapid heating. However, inverse bremsstrahlung heating can potentially support a certain level of anisotropy in the electron distribution function. More work needs to be done in developing the self-consistent electron distribution for this case and in comparing this model with experiments.

## VII Conclusions

We have performed numerical simulations based on the test particle model which involves classical electron trajectories in the electric field of an ion and a homogeneous, high frequency laser field. We have confirmed recent results by Fraiman et al. fraiman  regarding the electron orbits that closely approach an ion several times during each scattering event and subsequently lead to a large angle deflection and to a large enhancement in the electron energy. These correlated collisions have been ignored in preceding theories which significantly underestimate the electron heating rates by neglecting large angle scattering and the important field phase dependence of electron trajectories. The inverse-bremsstrahlung heating rates have been evaluated for different electron distribution functions, including: mono-energetic beam-like electrons, Maxwellian and highly anisotropic distribution relevant to photo-ionized gases. The energy gain for electrons with small drift velocities and for large amplitude electric fields, $`v_0<v_E`$, increases as $`1/v_0`$ and exceeds the previously known results if $`v_0<v_E/\mathrm{\Lambda }`$. This behavior contradicts results of all analytical theories and it is especially important for electrons with initial velocities that are parallel to the field amplitude. When the quantum short-distance cut-off is introduced through an effective potential, the energy gains demonstrate a similar dependence on the particle drift velocity but a smaller magnitude of enhancement is obtained if the cut-off length, $`\lambda _B=2\pi \mathrm{}/mv_E`$, is larger than the characteristic radius of the Coulomb interaction $`\rho _C=Ze^2/T_e`$.

For isotropic Maxwellian distribution functions, the heating rate in the limit of $`v_E>v_{Te}`$ is independent on the laser intensity. It can be approximated by the interpolation expression (19) for the effective collision frequency. For highly anisotropic electron distribution functions such as those produced by photo-ionization, the new effect of large angle correlated scattering dominates the heating rate at high laser intensities. It leads to $`v_E^{3/2}`$ dependence of the heating rate.

Our numerical results involve calculations of the electron scattering on a single ion. They are relevant to realistic plasma conditions provided the two-body collision approximation $`r_E<d`$ (7) is satisfied. This imposes an upper limit on the laser intensity for a given ion density and typically holds for very underdense plasmas. Also, we have implicitly assumed that the plasma is weakly coupled and that the average distance between ions, $`d`$ is smaller than the screening length, $`\lambda _D`$.

The enhancement of the inverse bremsstrahlung heating rates should influence laser plasma interaction experiments. For example, the collisional heating of long preformed plasmas by high intensity short laser pulses accompanies several important applications, including X-ray lasers, laser wake-field particle accelerators, and the high-harmonics generation. Moderate ion densities in the plasma waveguide ($`n_i=10^{17}10^{19}`$ cm-3) should make it easier to satisfy condition (7). High intensity picosecond laser pulses would allow investigation of the heating rate dependence on $`v_E/v_{Te}`$ over the wide range of values. Atomic processes involved in X-ray lasers are very sensitive to plasma temperature and could provide a diagnostic method for the inverse bremsstrahlung heating rates in addition to other diagnostics techniques such as the Thomson scattering and absorption measurements.

The authors thank Dmitri Romanov for his advise concerning numerical calculations. This work was supported by the Natural Sciences and Engineering Research Council of Canada, the Alberta Ingenuity Fund and the Russian Foundation for Basic Research (grant N 00-02-16063).


format: text

text:
# Wavepacket reconstruction via local dynamics in a parabolic lattice

## Abstract

We study the dynamics of a wavepacket in a potential formed by the sum of a periodic lattice and of a parabolic potential. The dynamics of the wavepacket is essentially a superposition of “local Bloch oscillations”, whose frequency is proportional to the local slope of the parabolic potential. We show that the amplitude and the phase of the Fourier transform of a signal characterizing this dynamics contains information about the amplitude and the phase of the wavepacket at a given lattice site. Hence, complete reconstruction of the the wavepacket in the real space can be performed from the study of the dynamics of the system.

preprint: v 1.55

The reconstruction of the wavefunction of a quantum system in the real space from experimentally accessible measurements is seldom possible. Even when the spatial probability distribution is accessible, phase measurements in general require the use of a completely different, and often incompatible, technique. This is essentially due to the fact that a measurement of the position of a quantum system is usually accompanied by a wavepacket reduction that erases the phase information. In a one-atom Magneto-optical trap ref:OneAtom  (prepared in a reproductive way), repeated measurements of the fluorescence light allow reconstruction of the atom probability distribution with a precision of the order of the wavelength of the detected light, whereas a measurement of wavepacket phases needs a specific atom-interferometry setup.

In contrast, the quantum dynamics of a system is highly sensitive to both the amplitude and the phase of the wavepacket. The idea at the root of the present letter is to use information on the dynamics of the system to determine simultaneously both these quantities, with a spatial resolution at least comparable to that of the methods mentioned in the preceding paragraph. The dynamics exploited here is related to the so-called “Bloch oscillation” (BO) observed in the motion of a quantum particle in a spatially periodic lattice to which a constant slope is added, forming a “tilted lattice”. In such a potential, instead of climbing down the slope of the potential, the wavepacket performs bounded periodic oscillations in real and momentum spaces. The origin of this oscillatory motion can be retraced to quantum interference among pieces of the original wave packet that are transmitted and reflected by the successive potential barriers of the lattice, in much analogy with a chain of optical interferometers with a linearly varying effective index ref:OptBO . The coherence of the wavepacket is thus a pre-requisite to the observation of the BO. The frequency of oscillation (the “Bloch frequency”) corresponds to the energy shift from one lattice well to the next, and thus to the slope of the potential. This behavior, theoretically predicted in 1934 ref:Zener , has been observed experimentally in solid-state superlattices ref:SolidStateBO , with cold atoms ref:BlochOsc , with Bose-Einstein condensates ref:BECBO , and also with photons, in arrays of optical waveguides ref:OptBO . A theoretical study of the BO in an frame close to that adopted here can be found in Ref. ref:WSbasis .

Let us consider a quantum system, the “atom”, of mass $`M`$, placed in an one-dimensional potential formed by a spatially-periodic lattice of period $`d`$, to which a parabolic term $`ax^2/2`$ is added. With cold-atoms, this can be obtained by using a far-detuned standing wave which creates a sinusoidal lattice ref:BlochOsc ; ref:BECBO ; ref:ExpWannierStark , superposed to an independent, focalized laser beam, producing a parabolic term (a laser-atom detuning much greater than the natural width of the concerned electronic transition reduces spontaneous emission, insuring the system to be conservative). The Hamiltonian of the system is

$$H=\frac{1}{2M^{}}\frac{^2}{x^2}V_0\mathrm{cos}(2\pi x)+\frac{x^2}{2}$$
(1)
where we introduced dimensionless units in which distances are measured in units of the lattice period $`d`$, energy is measured in units of $`ad^2`$ \[or, equivalently, time is measured in units of $`\mathrm{}/(ad^2`$)\], $`M^{}Mad^4/\mathrm{}^2`$ is a reduced mass, and $`\mathrm{}=1`$. Without loss of generality, we shall use a sinusoidal lattice.

The energy spectrum of the Hamiltonian Eq. (1) is plotted in Fig. 1 and shows two distinct kinds of states. For energies below a threshold value $`\epsilon _L`$, the corresponding eigenstates are localized in symmetric lattice wells $`n`$ and $`n`$ (where the lattice site at $`n=0`$ is tied to the symmetry axis of the parabola). For each possible value of the energy, one finds a symmetric and an anti-symmetric quasi-degenerated states, which are noted resp. $`\varphi _n^S(x)`$ and $`\varphi _n^A(x)`$. We suppose that the depth of the lattice is such that there are only one symmetric and one antisymmetric state for each couple of lattice wells $`n`$ and $`n`$ ($`V_0=90`$). The corresponding eigenenergies are noted $`n^2/2+\delta _n^S`$ and $`n^2/2+\delta _n^A`$. For energies above $`\epsilon _L`$, one finds both localized and delocalized states. At high energies, the latter are similar to the eigenstates of the harmonic oscillator and their eigenenergies approximately display a linear dependence on $`n`$ (although this is not easily seen in Fig. 1).

We suppose that the atom wavepacket has been prepared in a superposition of localized states only. For high enough lattice amplitude $`V_0,`$ one can neglect the energy shifts $`\delta _n^S`$ and $`\delta _n^A`$ note:quasi-degeneracy , and consider as a new basis the left and right quasi-eigenstates $`\varphi _{\pm n}(x)\left[\varphi _n^S(x)\pm \varphi _n^A(x)\right]/\sqrt{2}`$, which are localized in the lattice site $`\pm n`$ with energy $`ϵ_n=n^2/2`$ (up to an additive constant). Moreover, the localized eigenstates $`\varphi _n(x)`$ are almost invariant under a translation by an integer number of steps of the lattice: $`\varphi _n(x)\varphi _m[x(nm)d]`$. This property is exactly verified for a tilted lattice, and is a good approximation in the present problem as long as the parabolic potential do not perturb too much the translational invariance of the lattice wells over the spatial extension corresponding to the atom wavepacket (that is, as long as $`V_01`$ in our units).

At time $`t=0`$ the wavepacket $`\mathrm{\Psi }_0(x)`$ is centered at $`n_0`$ and has a spatial extension $`\mathrm{\Delta }n`$: $`\mathrm{\Psi }_0=_nc_n\varphi _n(x)`$. At time $`t`$ one has:

$$\mathrm{\Psi }(x,t)=\underset{n}{}c_n\mathrm{exp}\left[i\frac{n^2}{2}t\right]\varphi _n(x)$$
(2)
from which one can determine the dynamics of the system ref:WSbasis . For instance, the packet average position is

$$\overline{x}=\mathrm{\Psi }|x|\mathrm{\Psi }=\underset{n}{}X_mc_nc_{nm}^{}\mathrm{exp}\left[im(n\frac{m}{2})t\right]$$
(3)
($`X_m\varphi _n\left|x\right|\varphi _{n+m}`$ is almost independent of $`n`$).

Eq. (3) shows that, given the evolution of a quantum system, its coherences $`c_nc_{nm}^{}`$ can be obtained from the complex amplitude of the oscillations at the Bohr frequency $`\omega _{nm}`$:

$$\omega _{nm}=ϵ_nϵ_m=m(nm/2)\text{ .}$$
(4)
The parabolic component of the potential produces non-degenerate Bohr frequencies. In a tilted lattice corresponding to a constant force $`F`$, the coherences between neighbor sites $`c_nc_{n1}^{}`$ correspond to an oscillation with the degenerate ($`n`$-independent) Bohr (or Bloch) frequency $`\omega _B=Fd/\mathrm{}`$ ($`\omega _B=F`$ in reduced units). The parabolic potential, on the contrary, has a local slope $`F=x`$ and the related “local” Bloch frequency depends on $`n`$. The local Bloch frequency (LBF) is by definition the Bohr frequency associated to neighbor sites at $`n`$ and $`n1`$ ($`m=1`$):

$$\omega _B(n)=ϵ_nϵ_{n1}=n\frac{1}{2}\text{.}$$
(5)
The presence of such a frequency in the dynamics is a signature of the coherence $`c_nc_{n1}^{}`$ that can be used to track the wavepacket amplitude and phase, as shown below.

Experimentally, it is more convenient to work in the momentum-representation, because the spatial amplitude of the BO is in general very small whereas, in the momentum space, the amplitude of the oscillation is of the order of $`\mathrm{}/d`$ (the extension of a Brillouin zone of the lattice), which is easy to resolve with available experimental methods. For cold atoms, for example, stimulated Raman transitions ref:QuChaos  can be used to that end. A great deal of information can be obtained by just monitoring the evolution of the probability for a given value $`p`$ of the momentum, $`P(p,t)`$. We thus work in momentum space and introduce the eigenstates $`\phi _n(p)`$ in $`p`$representation:

$`\phi _n(p)`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \mathrm{exp}(ipx)\varphi _{n_0}\left[x(nn_0)\right]𝑑x}`$ (6)
$`=`$ $`\mathrm{exp}\left[i(nn_0)p\right]\phi _{n_0}(p)`$
where we used the approximate translation invariance of the real-space eigenfunctions $`\varphi _n(x)`$. The evolution of the momentum probability-distribution is then:

$`P(p,t)`$ $`=`$ $`|\mathrm{\Psi }(p,t)|^2=|\phi _{n_0}(p)|^2\times `$ (7)
$`(1+{\displaystyle \underset{m0}{}}\mathrm{exp}(imp)`$
$`{\displaystyle \underset{nm}{}}c_nc_{nm}^{}\mathrm{exp}[im(n{\displaystyle \frac{m}{2}})t])\text{ .}`$
Fig. 2 displays the Fourier transform of the quantity $`P_0(t)P(0,t)`$ and has been obtained as follows. The time-dependent probability distribution $`P_0(t)`$ is calculated by numerical integration of the Schrödinger equation with the hamiltonian of Eq. (1) and with the initial wavepacket $`\mathrm{\Psi }_0(x)`$ shown in the inset of Fig. 2 (the phase of $`c_n`$ is arbitrarily set to $`e^{in\pi /4}`$). The spectrum of $`P_0(t)`$ presents components at the Bohr frequencies contained in Eq. (7), from which the amplitude and phase of the coherences $`c_nc_{nm}^{}`$ can be deduced. The interpretation of the spectrum may be complicated by the fact that Bohr frequencies corresponding to energy differences higher than next-neighbor wells may show up. These frequencies, corresponding to $`m>1`$ in Eq. (4), are higher than the LBF \[Eq. (5)\] typically by a factor $`m`$. In the general case, the spectrum is thus composed of “folds” described by an integer value of $`m`$. However, a judicious choice of position of the wavepacket in the parabolic potential allow to separate one fold from the others. Frequency folds $`m=1`$ and $`m=2`$ are resolved if $`n_0+\mathrm{\Delta }n/2+1/2<2\left(n_0\mathrm{\Delta }n/21\right)`$ or, roughly, $`\mathrm{\Delta }nn_0`$. In producing Fig. 2, we managed to satisfy this condition. The first frequency fold is clearly seen as the low-frequency part of the spectrum associated to the coherences $`c_nc_{n1}^{}`$. In Fig. 2, the shape of the wavepacket is reproduced by the envelopes of both the first and the second frequency folds. This is the case if the wavepacket is smooth enough at the scale of the lattice step that we can take $`|c_n||c_{n1}|`$, and then directly identify the amplitude of the $`n`$-component of the spectrum with $`|c_n|`$.

The complete reconstruction of the spatial wavefunction from the $`c_n`$ coefficients requires, in principle, the knowledge of the eigenstates of system. The detailed spatial shape of the eigenstates may depend on experimental parameters that may not be known with a sufficient precision. With our system, the shape of the localized states is almost independent of $`n`$, as it is essentially determined by the lattice, the parabolic term giving a small correction. As they are sharply localized, they can be considered as delta-function-like, probing the local value of the wavepacket. The detailed shape of the eigenstates is thus not required for the reconstruction, at the a precision corresponding to a lattice step.

We can now proceed to the complete reconstruction of a wavepacket. We consider the first fold ($`m=1`$) of the spectrum, and the zero momentum probability $`P_0(t)`$ only. The other folds have essentially the same structure; other values of the probed momentum may be used to extract complementary information or improve the accuracy of the reconstruction. We thus reduce Eq. (7) to the simpler form:

$`Q_0(t)`$ $``$ $`{\displaystyle \frac{P_0(t)}{|\phi _{n_0}(0)|^2}}1`$ (8)
$`=`$ $`{\displaystyle \underset{n}{}}c_nc_{n1}^{}\mathrm{exp}\left[i\left(n{\displaystyle \frac{1}{2}}\right)t\right]\text{ .}`$
One then sees that the Fourier component $`\stackrel{~}{Q_0}(n)=c_nc_{n1}^{}`$ of $`Q_0`$ at frequency $`n1/2`$ directly gives the coherence between neighbor sites. If the wavepacket envelope is smooth, the amplitudes can be directly obtained by taking $`|c_n|^2|\stackrel{~}{Q_0}(n)|`$. If the phase also varies smoothly, the phase shift from the site $`n1`$ to the site $`n`$ is also directly obtained from the phase of $`\stackrel{~}{Q_0}(n)`$, up to an unavoidable, but unimportant, global phase factor.

A more precise determination of the coefficients $`c_n`$ is obtained by using the second frequency fold which, as one can easily convince oneself, gives the coherences $`c_nc_{n2}^{}\stackrel{~}{Q_0^{}}(n)`$, and the identity:

$$\frac{\stackrel{~}{Q_0}(n)\stackrel{~}{Q_0}(n+1)}{\stackrel{~}{Q_0^{}}(n+1)}=|c_n|^2\text{ .}$$
(9)
The phase shifts can be straightforwardly obtained from the $`\stackrel{~}{Q_0}(n)`$. We display in Fig. 3 the coherences $`c_n`$ deduced from the spectrum of Fig. 2. The agreement with the amplitudes of the original wavepacket is excellent. The figure also shows real and imaginary parts of the complex amplitudes in very good agreement with the original wavepacket.

We numerically verified that the above method can be applied to wavepackets of arbitrary shape. The reconstructed state then reproduces the wavepacket at the scale of the lattice step, all finer features are lost. A large number of delocalized states can appear in the decomposition of the initial state, but their contribution to the Bloch-oscillation spectrum is to add a wide background whose amplitude is much smaller than the peaks associated to the coherence between localized states.

Consider now the case in which there are many atoms simultaneously interacting with the potential. If the density of atoms is low, so that there is in average one atom per site (or less), the present method reconstructs the atom distribution. If the atomic density is higher, the method measures the ensemble average of the atomic correlation $`c_nc_{nm}^{}`$. This quantity is very important because it determines the ability of the system to display coherent quantum dynamics ref:WSbasis . Finally, note that the present method can also be applied to a tridimensional lattice: adding a parabolic potential in a given direction of the lattice will probe the wavefunction along this direction.

In conclusion, we have presented a method allowing reconstruction of the spatial wavefunction of a quantum system using the Fourier spectrum of the Bloch oscillations in a periodic lattice superposed to a parabolic potential. Our method is quite general: it does not depend, e.g. on internal properties of the quantum particle being studied, or on the shape of the lattice wells. A parabolic shape is not necessary provided that (i) the Bohr (or Bloch) frequencies of the system are not degenerated and (ii) the approximate translational symmetry of the localized states is preserved. The method requires a potential shape presenting both a non-constant and non-vanishing slope in the region of interest: the local slope induces local Bloch oscillation and the spatial dependence, in turn, leads to the position-sensitivity of the method. In principle, the method is applicable to a variety of systems, including cold atoms in light potentials, Bose-Einstein condensates and electrons in an adequately constructed superlattice. Applied to a Bose-Einstein condensate, whose Bloch oscillations have been recently observed ref:BECBO , the present method may be used to probe nonlinear collective interactions. The method is also able to probe the spatial entanglement of two or more atoms, which can be of interest in the context of quantum information and quantum computing. We shall treat these interesting applications in the future.

The authors would like to acknowledge D. Delande, P. Szriftgiser and H. Lignier for fruitful discussions. This work is partially supported by an action “ACI Photonique” of the Ministère de la Recherche. Laboratoire de Physique des Lasers, Atomes et Molécules (PhLAM) is UMR 8523 du CNRS et de l’Université des Sciences et Technologies de Lille. Centre d’Etudes et Recherches Lasers et Applications (CERLA) is FR 2416 du CNRS, supported by Ministère de la Recherche, Région Nord-Pas de Calais and Fonds Européen de Développement Economique des Régions (FEDER).


text:
# References

hep-ph/0203021

RAL-TR-2002-006

29 Apr 2002

Neutrino Matter Effect Invariants and the Observables of

Neutrino Oscillations

P. F. Harrison

Physics Department, Queen Mary University of London

Mile End Rd. London E1 4NS. UK 111E-mail:p.f.harrison@qmul.ac.uk

and

W. G. Scott

Rutherford Appleton Laboratory

Chilton, Didcot, Oxon OX11 0QX. UK 222E-mail:w.g.scott@rl.ac.uk

## Abstract

We generalise our previous observation on the invariance of the Jarlskog determinant to matter effects in neutrino oscillations. Within the context of standard neutrino oscillation theory with matter effects, we present the complete set of (five) matter invariant observables for neutrino propagation in matter. We give some examples of their application.

To be published in Physics Letters B

1 Introduction Neutrino oscillations violate lepton flavour conservation. The amplitude $`A_{\alpha \beta }`$ for a neutrino of flavour $`\alpha `$ to be detected as a neutrino of flavour $`\beta `$ is given as a function of propagation distance $`L`$ by the (matrix) equation:

$$A=\mathrm{exp}(iHL)$$
(1)
where $`H`$ is the effective Hamiltonian in a flavour basis (ie. in a weak basis in which the charged lepton mass matrices are diagonal - this is defined up to an arbitrary phase redefinition of the charged lepton fields). We may take:

$$H=U\mathrm{diag}(0,\mathrm{\Delta }_{21},\mathrm{\Delta }_{31})U^{}$$
(2)
where $`U`$ is the conventional MNS matrix and the $`\mathrm{\Delta }_{ij}`$ are defined by $`\mathrm{\Delta }_{ij}=(m_i^2m_j^2)/2E`$, with $`m_i`$ the neutrino masses and $`E`$ the neutrino energy. In this paper, all calculations are assumed to be in such a flavour basis.

We will find the following theorem useful. Any function, $`f(H)`$, of a Hermitian operator $`H`$, can be expanded using Lagrange’s formula :

$$f(H)=\underset{i}{}f(\lambda _i)X^i$$
(3)
where the $`\lambda _i`$ are the eigenvalues of $`H`$ and the $`X^i`$ are hermitian (projection) operators given by:

$$X^i=\frac{_{ji}(H\lambda _j)}{_{ji}(\lambda _i\lambda _j)}.$$
(4)
Since $`H`$ satisfies its own characteristic equation, the complete product $`_i(H\lambda _i)`$ (with no factors dropped) clearly vanishes. It follows that the $`X^i`$ are both left and right eigenvectors of $`H`$ with eigenvalues $`\lambda _i`$, ie. $`(H\lambda _i)X^i=X^i(H\lambda _i)=0`$. With the normalisation defined as in Eq. 4, we have $`(X^i)^2=X^i`$, so the $`X^i`$ are indeed projection operators, representing the (pure) mass-eigenstate density-matrices.

The ancient formula Eq. 3 has been used previously for the specific case of three-family neutrino oscillations in eg. Ref. , and has recently re-appeared in the same context in Ref. . Applied to Eq. 1 it gives:

$$A=\underset{i}{}X^i\mathrm{exp}(i\lambda _iL).$$
(5)
Any $`X^i`$ is readily computed from the corresponding column of the MNS matrix :

$$X_{\alpha \beta }^i=U_{\alpha i}U_{\beta i}^{}$$
(6)
(where no summation over $`i`$ is implied). We note here that the elements of the matrices $`X^i`$ are dependent on the phase-convention chosen for the charged lepton fields and are therefore not observables (likewise, of course, neither are the MNS matrix elements themselves).

Using Eq. 3 we may compute arbitrary powers of $`H`$ as follows:

$$H^r=\underset{i}{}\lambda _i^rX^i.$$
(7)
Given the eigenvalue spectrum, Eq. 4 determines the matrices $`X^i`$ in terms of $`H`$, whilst setting $`r=1`$ in Eq. 7 yields the inverse relation, determining $`H`$ in terms of the $`X^i`$. The case $`r=0`$ of Eq. 7 expresses the unitarity of $`U`$. For an $`n\times n`$ Hamiltonian, there are only $`n`$ independent equations of the form Eq. 7, because the characteristic equation can be used to write higher powers $`rn`$ in terms of the $`n`$ lowest powers $`H^r,0r<n`$.

The squared amplitude $`|A_{\alpha \beta }|^2`$ for $`\alpha \beta `$ ($`\alpha =\beta `$) gives the appearance (survival) probabilities as a function of propagation distance (the dependence on neutrino energy is implicit in the factors $`\mathrm{\Delta }_{ij}`$):

$`P(\alpha \beta )=|A_{\alpha \beta }|^2`$ $`=`$ $`\delta _{\alpha \beta }4{\displaystyle \underset{i<j}{}}K_{\alpha \beta }^{ij}\mathrm{sin}^2(\mathrm{\Delta }_{ij}L/2)`$ (8)
$`+`$ $`8J_{\alpha \beta }\mathrm{sin}(\mathrm{\Delta }_{12}L/2)\mathrm{sin}(\mathrm{\Delta }_{23}L/2)\mathrm{sin}(\mathrm{\Delta }_{31}L/2)`$ (9)
where

$$K_{\alpha \beta }^{ij}=\mathrm{Re}(X_{\alpha \beta }^iX_{\alpha \beta }^j),$$
(10)
and

$$J_{\alpha \beta }=\mathrm{Im}(X_{\alpha \beta }^iX_{\alpha \beta }^j)=\pm J(\mathrm{for}\alpha \beta )\mathrm{or}=0(\mathrm{for}\alpha =\beta ).$$
(11)
$`J`$ is Jarlskog’s $`CP`$\- and $`T`$-violating invariant . Inspection of Eqs. 9-11 indicates that the phase-convention-independent products of the form $`X_{\alpha \beta }^iX_{\alpha \beta }^j`$ $`(ij)`$ are observables of the system , the real parts $`K_{\alpha \beta }^{ij}`$ parameterising $`CP`$-conserving oscillation amplitudes and the imaginary parts $`J_{\alpha \beta }`$ parameterising $`CP`$ violation. We note that unitarity imposes many constraints among the observables, $`K_{\alpha \beta }^{ij}`$ (and $`J`$), and that once three of these are specified (such that at least two different flavour indices are used) along with $`J`$, the remaining ones may be determined via unitarity.

When neutrinos propagate in matter, the Hamiltonian in the flavour basis is modified: $`HH^{}`$, with

$$H^{}=H+\mathrm{diag}(a,0,0)$$
(12)
where $`a=\sqrt{2}G_FN_e`$. This modifies both the eigenvalues and eigenstates of the neutrino Hamiltonian (in this paper, “primed” symbols will be used to denote the values of quantities when neutrinos propagate in matter and unprimed symbols for their vacuum counterparts).

We recently  pointed out that the determinant of the commutator of the lepton mass matrices is invariant with respect to matter effects,

$$\mathrm{Det}[D_{\mathrm{}}^2,H^{}]=\mathrm{Det}[D_{\mathrm{}}^2,H]$$
(13)
where $`D_{\mathrm{}}`$ is the charged-lepton mass matrix, which is diagonal by construction in this basis. This allowed the derivation of a matter-invariant quantity, expressible both as a function of elements of the effective Hamiltonian, and as a combination of the eigenvalue differences, $`\mathrm{\Delta }_{ij}`$, and the $`CP`$ violation parameter, $`J`$ :

$`\mathrm{Im}(H_{e\mu }^{}H_{\mu \tau }^{}H_{\tau e}^{})`$ $`=`$ $`\mathrm{Im}(H_{e\mu }H_{\mu \tau }H_{\tau e})`$ (14)
$`=\mathrm{\Delta }_{12}^{}\mathrm{\Delta }_{23}^{}\mathrm{\Delta }_{31}^{}J^{}`$ $`=`$ $`\mathrm{\Delta }_{12}\mathrm{\Delta }_{23}\mathrm{\Delta }_{31}J.`$ (15)
This quantity has found application in several instances . It is not however, the only matter-invariant combination of neutrino oscillation observables. Here we count the number of independent such combinations and derive a complete set.

2 A Complete Set of Matter Invariant Observables Clearly, in a flavour basis, the modification to the effective Hamiltonian induced by matter, Eq. 12, affects only the “$`ee`$” element:

$$H_{ee}^{}=H_{ee}+a$$
(16)
whereas the remaining eight elements of $`H`$ are matter-invariant in this basis, ie.

$$H_{\alpha \beta }^{}=H_{\alpha \beta },\alpha \beta \mathrm{`}\mathrm{`}ee\mathrm{"}.$$
(17)
The invariants of Eq. 17 are not immediately useful, however, as the Hamiltonian (vacuum or matter) is not itself observable. It transforms under (separate) phase redefinitions of the flavour eigenstates:

$$H\mathrm{\Phi }H\mathrm{\Phi }^{}$$
(18)
where

$$\mathrm{\Phi }=\mathrm{exp}\{i\mathrm{diag}(\varphi _e,\varphi _\mu ,\varphi _\tau )\}$$
(19)
so that

$$H_{\alpha \beta }e^{i(\varphi _\alpha \varphi _\beta )}H_{\alpha \beta }$$
(20)
without any observable consequence. The $`X^i`$ transform similarly. If our invariants are to be readily applicable, we must seek to relate observable quantities directly.

We can count the number of independent observables represented by $`H`$ as follows. An arbitrary $`n\times n`$ Hermitian matrix comprises $`n^2`$ independent real parameters ($`n`$ on-diagonal, $`n(n1)`$, off-diagonal). Eqs. 18-20 imply that for the effective Hamiltonian, $`n1`$ of these are unobservable relative phases, so that we would expect, at most, $`n^2n+1`$ observable quantities. In addition however, subtraction of any multiple of the identity from $`H`$ (a “trace transformation”) leaves all neutrino oscillation observables invariant, so that only eigenvalue differences are relevant, leaving finally only $`n(n1)`$ independent observables in neutrino oscillations. For the $`3\times 3`$ MNS oscillation phenomenology, we have six observables, which are conventionally considered to be three real angles, one phase and two mass-squared differences. They could equally well be taken to be three independent $`K_{\alpha \beta }^{ij}`$ (see Eq. 9), $`J`$ and two eigenvalue differences, $`\mathrm{\Delta }_{21}`$ and $`\mathrm{\Delta }_{31}`$, or alternatively six phase- and trace-invariant functions of elements of $`H`$ (see below). When neutrinos propagate in matter, one parameter is modified, Eq. 16, so that we may expect, in general, $`n(n1)1`$ matter invariances (eg. five, in the three-family case).

We seek relations analogous to Eqs. 16-17 but involving only functions of the $`H_{\alpha \beta }`$ which are invariant under both phase transformations and trace transformations. Taking Eq. 15 as a model, we seek furthermore to write them in terms of the observables of Eq. 9. Six independent such relations should exist. Relations which do not involve the “ee” flavour label will, in addition, be matter-invariant. The off-diagonal elements, $`H_{\alpha \beta }`$ $`(\alpha \beta )`$, are already invariant under trace transformations, and we can form phase-transformation invariants simply by taking their moduli-squared:

$`|H_{\alpha \beta }^{}|^2`$ $`=`$ $`|H_{\alpha \beta }|^2`$ (21)
$`=|{\displaystyle \underset{i}{}}\lambda _i^{}X_{\alpha \beta }^{}_{}{}^{}i|^2`$ $`=`$ $`|{\displaystyle \underset{i}{}}\lambda _iX_{\alpha \beta }^i|^2`$ (22)
$`={\displaystyle \underset{j<k}{}}\mathrm{\Delta }_{jk}^2K_{\alpha \beta }^{jk}`$ $`=`$ $`{\displaystyle \underset{j<k}{}}\mathrm{\Delta }_{jk}^2K_{\alpha \beta }^{jk},\alpha \beta `$ (23)
where we have used the relation, Eq. 7 with $`r=1`$ to obtain Eq. 22, and with $`r=0`$ (unitarity) to obtain Eq. 23. The quantities in Eq. 23 are all observable (see eg. Eq. 9) and the three independent combinations are each matter-invariant (only three combinations because $`|H_{\alpha \beta }|=|H_{\beta \alpha }|`$). 111The diagonal analogue of Eqs. 21-23: $`_{j<k}\mathrm{\Delta }_{jk}^2K_{\alpha \alpha }^{jk}=(H_{\alpha \alpha })^2(H^2)_{\alpha \alpha }`$ is not independent, the LHS equating through unitarity to: $`_{j<k}\mathrm{\Delta }_{jk}^2K_{\alpha \beta }^{jk}+_{j<k}\mathrm{\Delta }_{jk}^2K_{\alpha \gamma }^{jk},\alpha \beta \gamma `$, and the RHS being simply $`|H_{\alpha \beta }|^2|H_{\alpha \gamma }|^2`$

In a complementary fashion, the diagonal elements of $`H`$ are already phase-invariant and we can form trace-transformation invariants by taking differences, whereby:

$`H_{\mu \mu }^{}H_{\tau \tau }^{}`$ $`=`$ $`H_{\mu \mu }H_{\tau \tau }`$ (24)
$`=(X_{\mu \mu }^{}_{}{}^{}2X_{\tau \tau }^{}_{}{}^{}2)\mathrm{\Delta }_{21}^{}+(X_{\mu \mu }^{}_{}{}^{}3X_{\tau \tau }^{}_{}{}^{}3)\mathrm{\Delta }_{31}^{}`$ $`=`$ $`(X_{\mu \mu }^2X_{\tau \tau }^2)\mathrm{\Delta }_{21}+(X_{\mu \mu }^3X_{\tau \tau }^3)\mathrm{\Delta }_{31}`$ (25)
which is matter-invariant, and

$`H_{ee}^{}Tr(H^{})/3`$ $`=`$ $`H_{ee}Tr(H)/3+2a/3`$ (26)
$`=(X_{ee}^{}_{}{}^{}2{\displaystyle \frac{1}{3}})\mathrm{\Delta }_{21}^{}+(X_{ee}^{}_{}{}^{}3{\displaystyle \frac{1}{3}})\mathrm{\Delta }_{31}^{}`$ $`=`$ $`(X_{ee}^2{\displaystyle \frac{1}{3}})\mathrm{\Delta }_{21}+(X_{ee}^3{\displaystyle \frac{1}{3}})\mathrm{\Delta }_{31}+2a/3`$ (27)
which may be said to be matter-covariant.

The three invariances of Eq. 22 correspond to the “sum-rules” of Ref. , although they are derived in a different way here, without appealing to the commutator of the mass-matrices 222The quantities $`Z_{ij}`$ and $`\stackrel{~}{Z}_{ij}`$ defined in Ref.  are in fact simply the off-diagonal elements of the effective Hamiltonian in vacuum and in matter respectively.. Their expression in terms of our observables, Eq. 23, and the relations of Eqs. 25 and 27 are new. Our derivation demonstrates that all the above invariances follow straightforwardly from the fact that the Wolfenstein term affects only the “$`ee`$” element of the effective Hamiltonian, Eqs. 16, 17.

The five combinations Eqs. 23-27 all relate $`CP`$-conserving observables (the diagonal elements $`X_{\alpha \alpha }^i=|U_{\alpha i}|^2`$ are observable). The combination Eq. 15 parameterises $`CP`$ violation and is therefore independent. It is manifestly matter-, phase- and trace-transformation invariant, as required. We have therefore completed the derivation of six independent identities, Eqs. 15 and 23-27, relating particular combinations of matter-modified neutrino oscillation observables to their vacuum counterparts. Five of the combinations (Eqs. 15, 23 and 25) are matter-invariant, while one combination (Eq. 27) depends in a simple way on the matter-density. As regards neutrino oscillation phenomenology, these six combinations constitute a complete set, whereby any further matter-invariants (which can certainly be found, eg. $`\mathrm{Re}(H_{e\mu }H_{\mu \tau }H_{\tau e})`$) must always be expressible in terms of those given here.

3 Some Applications As a first application, we show briefly how the formalism presented in the last Section applies to the simple case of two neutrino families. In this case, there is only one independent off-diagonal element of the effective Hamiltonian, $`|H_{e\mu }|`$, which leads to the analogue of the matter-invariant in Eq. 23:

$`|H_{e\mu }^{}|^2`$ $`=`$ $`|H_{e\mu }|^2`$ (28)
$`=\mathrm{\Delta }_{12}^2K_{e\mu }^{12}`$ $`=`$ $`\mathrm{\Delta }_{12}^2K_{e\mu }^{12}`$ (29)
$`\mathrm{\Delta }_{12}^2\mathrm{sin}^22\theta ^{}`$ $`=`$ $`\mathrm{\Delta }_{12}^2\mathrm{sin}^22\theta `$ (30)
while the two-family analogue of Eq. 27 is

$`(X_{ee}^{}_{}{}^{}2{\displaystyle \frac{1}{2}})\mathrm{\Delta }_{21}^{}`$ $`=`$ $`a/2+(X_{ee}^2{\displaystyle \frac{1}{2}})\mathrm{\Delta }_{21}`$ (31)
$`(\mathrm{sin}^2\theta ^{}{\displaystyle \frac{1}{2}})\mathrm{\Delta }_{21}^{}`$ $`=`$ $`a/2+(\mathrm{sin}^2\theta {\displaystyle \frac{1}{2}})\mathrm{\Delta }_{21}`$ (32)
$`\mathrm{cos}2\theta ^{}\mathrm{\Delta }_{21}^{}`$ $`=`$ $`a\mathrm{cos}2\theta \mathrm{\Delta }_{21}`$ (33)
where we have used the definitions of the $`X^i`$ and $`K_{\alpha \beta }^{ij}`$ in terms of the $`2\times 2`$ mixing matrix elements to recover the standard results . We have found one matter-invariant combination of observables , and one which is matter-covariant, the complete set for the two-family case (compared with five and one respectively for the three family case).

In the case of three families, we consider the appearance probabilities Eq. 9 with $`\alpha \beta `$, in matter of arbitrary (uniform) density. Expanding each term in $`\mathrm{\Delta }_{ij}^{}L/2`$ for $`\mathrm{\Delta }_{ij}^{}L/2`$ small, we find:

$`P^{}(\alpha \beta )`$ $`=`$ $`{\displaystyle \underset{i<j}{}}K_{\alpha \beta }^{{}_{}{}^{}ij}(\mathrm{\Delta }_{ij}^2L^2\mathrm{\Delta }_{ij}^4L^4/24+\mathrm{})`$ (34)
$`+`$ $`J_{\alpha \beta }^{}\mathrm{\Delta }_{12}^{}\mathrm{\Delta }_{23}^{}\mathrm{\Delta }_{31}^{}L^3+\mathrm{}`$ (35)
$`=`$ $`|H_{\alpha \beta }^{}|^2L^2+\mathrm{Im}(H_{e\mu }^{}H_{\mu \tau }^{}H_{\tau e}^{})L^3+𝒪(\mathrm{\Delta }_{ij}^4L^4)+\mathrm{}`$ (36)
$`=`$ $`|H_{\alpha \beta }|^2L^2+\mathrm{Im}(H_{e\mu }H_{\mu \tau }H_{\tau e})L^3+𝒪(\mathrm{\Delta }_{ij}^4L^4)+\mathrm{}`$ (37)
We have used Eqs. 23 and 15 to identify the coefficients respectively of the $`CP`$-conserving term of order $`L^2`$ and the $`CP`$-violating term of order $`L^3`$ with manifestly matter-invariant combinations of elements of the effective Hamiltonian. Thus, they are determined directly in terms of elements of the vacuum Hamiltonian. The term of order $`L^4`$ is not matter-invariant, and is the lowest-order term through which matter effects enter. The well-known  result that matter has no observable effect for small propagation distances, is explained here by the fact that off-diagonal elements of $`H`$ are unmodified by matter. The result, Eq. 37, is valid for all mass and mixing patterns, but is, of course, useful only for small propagation distances.

A particular success of the matter invariance, Eq. 15, is that it enables the calculation of the universal $`CP`$ violating asymmetry in matter of uniform density, with no restriction to small $`L`$, knowing only the vacuum Hamiltonian and the matter-modified eigenvalue differences  (which are readily calculated in terms of vacuum quantities, ):

$$P^{}(\alpha \beta )P^{}(\beta \alpha )=\mathrm{Im}(H_{e\mu }H_{\mu \tau }H_{\tau e})\frac{\mathrm{sin}(\mathrm{\Delta }_{12}^{}L/2)\mathrm{sin}(\mathrm{\Delta }_{23}^{}L/2)\mathrm{sin}(\mathrm{\Delta }_{31}^{}L/2)}{(\mathrm{\Delta }_{12}^{}/2)(\mathrm{\Delta }_{23}^{}/2)(\mathrm{\Delta }_{31}^{}/2)}.$$
(38)
This expression does not require knowledge of the matter-modified eigenstates (ie. mixing parameters), whose dependence is contained completely within the matter-invariant factor $`\mathrm{Im}(H_{e\mu }H_{\mu \tau }H_{\tau e})`$. A comparison of the $`CP`$-conserving and $`CP`$-violating terms in Eq. 37 leads us to suspect that, perhaps more generally, the matter-invariant quantities $`|H_{\alpha \beta }|^2`$ play a similar role for $`CP`$-conserving observables to that played by $`\mathrm{Im}(H_{e\mu }H_{\mu \tau }H_{\tau e})`$ for $`CP`$-violating observables. It is interesting therefore to enquire whether a simplification similar to that of Eq. 38 exists also for the $`CP`$-conserving parts of $`P(\alpha \beta )`$, in terms of the matter-invariants $`|H_{\alpha \beta }|^2`$.

Some guidance is obtained by returning briefly to the two-family case. Consider Eq. 9 for the $`e\mu `$ appearance probability in matter of uniform density:

$`P^{}(e\mu )=|A_{e\mu }^{}|^2`$ $`=`$ $`4K_{e\mu }^{{}_{}{}^{}12}\mathrm{sin}^2(\mathrm{\Delta }_{21}^{}L/2)`$ (39)
$`=`$ $`|H_{e\mu }|^2{\displaystyle \frac{\mathrm{sin}^2(\mathrm{\Delta }_{21}^{}L/2)}{(\mathrm{\Delta }_{21}^{}/2)^2}}`$ (40)
where we have used Eq. 30 to write the probability as the product of a matter-invariant part $`(|H_{e\mu }^{}|^2=|H_{e\mu }|^2)`$ and a part with a generic dependence on eigenvalue differences of the form $`[\mathrm{sin}(\mathrm{\Delta }^{}L/2)/(\mathrm{\Delta }^{}/2)]^n`$. This exact expression, Eq. 40, is clearly analogous to Eq. 38 in the sense we require.

We do not a priori expect the three-family case to be as straightforward, but it will turn out that much of the above simplicity is retained, especially in the case of a hierarchical neutrino spectrum (as seems to be preferred experimentally). We may start again with Eq. 9, and assume a hierarchical spectrum, $`\mathrm{\Delta }_{21}^{}<<\mathrm{\Delta }_{31}^{}`$, $`\mathrm{\Delta }_{31}^{}\mathrm{\Delta }_{32}^{}\mathrm{\Delta }^{}`$ for all matter-densities considered. Then in the limit $`\mathrm{\Delta }_{21}^{}L/2<<1`$ (which is less restrictive than the limit in which Eq. 37 applies), we find for the leading oscillation of the appearance probability (disappearance probabilities are anyway readily obtained from appearance probabilities through unitarity):

$`P^{}(\alpha \beta )=|A_{\alpha \beta }^{}|^2`$ $``$ $`4(K_{\alpha \beta }^{{}_{}{}^{}31}+K_{\alpha \beta }^{{}_{}{}^{}32})\mathrm{sin}^2(\mathrm{\Delta }^{}L/2)`$ (41)
$`=`$ $`4|X_{\alpha \beta }^{}_{}{}^{}3|^2\mathrm{sin}^2(\mathrm{\Delta }^{}L/2)`$ (42)
where the last step follows from unitarity. We also require that $`|X_{\alpha \beta }^{}_{}{}^{}2|\mathrm{\Delta }_{21}^{}<<|X_{\alpha \beta }^{}_{}{}^{}3|\mathrm{\Delta }_{31}^{}`$ for all matter-densities considered, in order that the leading oscillation dominates. This will generally be true, for the hierarchical spectrum we have assumed, but may not be in some special cases eg. if the mixing matrix elements also have an extreme hierarchy, see below.

Applying the same approximations to Eq. 23 we have:

$`|X_{\alpha \beta }^{}_{}{}^{}3|^2\mathrm{\Delta }^2=(K_{\alpha \beta }^{{}_{}{}^{}31}+K_{\alpha \beta }^{{}_{}{}^{}32})\mathrm{\Delta }^2|H_{\alpha \beta }^{}|^2=|H_{\alpha \beta }|^2(K_{\alpha \beta }^{31}+K_{\alpha \beta }^{32})\mathrm{\Delta }^2=|X_{\alpha \beta }^3|^2\mathrm{\Delta }^2`$ (43)
(ie. the combination $`|X_{\alpha \beta }^{}_{}{}^{}3|\mathrm{\Delta }^{}`$ is a matter-invariant, as long as $`|X_{\alpha \beta }^{}_{}{}^{}2|\mathrm{\Delta }_{21}^{}<<|X_{\alpha \beta }^{}_{}{}^{}3|\mathrm{\Delta }_{31}^{}`$ remains valid, as already assumed). Combining Eqs. 42 and 43, we find that

$$P^{}(\alpha \beta )=|A_{\alpha \beta }^{}|^2|H_{\alpha \beta }|^2\frac{\mathrm{sin}^2(\mathrm{\Delta }^{}L/2)}{(\mathrm{\Delta }^{}/2)^2}$$
(44)
which generalises Eq. 40 to the three-family case, when the specified approximations apply. 333 The (possibly realised ) case $`|U_{e3}|^2=0`$, is a specific case which violates the assumption $`|X_{\alpha \beta }^{}_{}{}^{}2|\mathrm{\Delta }_{21}^{}<<|X_{\alpha \beta }^{}_{}{}^{}3|\mathrm{\Delta }_{31}^{}`$ for $`\alpha =e`$ or $`\beta =e`$, because then, $`X_{\alpha \beta }^{}_{}{}^{}3=0`$, and the leading oscillation vanishes. However, in this case, the long-wavelength oscillation becomes dominant, and Eqs. 42 and 44 become exactly true, with the substitutions $`|X_{\alpha \beta }^{}_{}{}^{}3||X_{\alpha \beta }^{}_{}{}^{}2|`$ and $`\mathrm{\Delta }^{}\mathrm{\Delta }_{21}^{}`$ etc.

Thus Eq. 44 comes close (as close as we believe is possible) to being the analogue of Eq. 38 for the CP-conserving case. Where valid, Eq. 44 retains all the calculational expediency previously recognised in connection with Eq. 38, for the CP-violating case. We should say that in general (for arbitrary mass-squared differences, mixing angles and propagation lengths) no such analogy exists between the $`CP`$-conserving and the $`CP`$-violating contributions. The simplification exemplified by Eq. 38 is a consequence of the universality (for three families only) of the CP-violating term, as embodied in the Jarlskog parameter, $`J`$ , but in general, no such universality exists for the CP-conserving $`K_{\alpha \beta }^{ij}`$ coefficients (the two-family case being an exception). However, for phenomenologically realistic oscillation parameters (see eg. ), Eq. 44 is certainly applicable to the leading oscillation of $`P^{}(\mu \tau )`$ (and of $`P^{}(\tau \mu )`$) for most energies of interest, and, unless it turns out that $`|X_{e\beta }^3|\mathrm{\Delta }_{31}|X_{e\beta }^2|\mathrm{\Delta }_{21}`$ for $`\beta =\mu ,\tau `$ ($`|X_{e\beta }^3|`$ depends on the presently poorly determined $`U_{e3}`$), it should also be relevant for appearance probabilities involving electron neutrinos.

4 Relation to Standard Formalism In order to relate our matter-covariant formulation to the standard one in terms of mixing angles, we first need to remove the implicit dependence of our matter-invariants on the neutrino energy (which has been left in until now in the interests of a concise discussion). We note that all the preceding equations in this paper remain unaltered, if we make the energy re-scaling $`H2EH`$, $`\mathrm{\Delta }_{ij}2E\mathrm{\Delta }_{ij}\mathrm{\Delta }m_{ij}^2`$, $`\lambda _i2E\lambda _i`$, $`a2Ea𝒜`$, $`LL/2E`$ and likewise for all primed quantities. The observables $`K_{\alpha \beta }^{ij}`$, $`J`$, and their matter-modified counterparts remain unchanged throughout, as do the oscillation probabililties.

We can now give the relationships between our energy-scaled matter-invariants (and covariant) (Eqs. 15 and 21-27) and the six parameters conventionally used to describe neutrino oscillation phenomenology, $`\theta _{12}`$, $`\theta _{23}`$, $`\theta _{13}`$, $`\delta `$, $`\mathrm{\Delta }m_{21}^2`$ and $`\mathrm{\Delta }m_{31}^2`$. Defining $`D=\mathrm{\Delta }m_{31}^2\mathrm{\Delta }m_{21}^2s_{12}^2`$ and $`F=\frac{1}{8}S_{12}S_{23}S_{13}c_{13}`$, we have:

$$_1|_{e\mu }|^2=c_{13}^2[s_{23}^2s_{13}^2D^2+\frac{1}{4}c_{23}^2S_{12}^2(\mathrm{\Delta }m_{21}^2)^2]+2FD\mathrm{\Delta }m_{21}^2c_\delta $$
(45)
$$_2|_{e\tau }|^2=c_{13}^2[c_{23}^2s_{13}^2D^2+\frac{1}{4}s_{23}^2S_{12}^2(\mathrm{\Delta }m_{21}^2)^2]2FD\mathrm{\Delta }m_{21}^2c_\delta $$
(46)
$$_3|_{\mu \tau }|^2=\frac{1}{4}\{[c_{13}^2S_{23}D(C_{12}S_{23}+S_{12}C_{23}s_{13}c_\delta )\mathrm{\Delta }m_{21}^2]^2+S_{12}^2s_{13}^2s_\delta ^2(\mathrm{\Delta }m_{21}^2)^2\}$$
(47)
$$_4\mathrm{Im}(_{e\mu }_{\mu \tau }_{\tau e})=\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{31}^2(\mathrm{\Delta }m_{31}^2\mathrm{\Delta }m_{21}^2)Fs_\delta $$
(48)
$$_5_{\mu \mu }_{\tau \tau }=C_{23}[C_{12}\mathrm{\Delta }m_{21}^2c_{13}^2D]S_{12}S_{23}s_{13}c_\delta \mathrm{\Delta }m_{21}^2$$
(49)
$$_0_{ee}\frac{1}{3}Tr=D(s_{13}^2\frac{1}{3})\frac{1}{3}C_{12}\mathrm{\Delta }m_{21}^2$$
(50)
where we have used the conventional parameterisation  and $`c_{ij}=\mathrm{cos}\theta _{ij}`$, $`s_{ij}=\mathrm{sin}\theta _{ij}`$, $`C_{ij}=\mathrm{cos}2\theta _{ij}`$, $`S_{ij}=\mathrm{sin}2\theta _{ij}`$ etc. Eqs. 45-49 remain true in matter (with primed quantities on the RHS), as well as in vacuum and show that (sums of) products of (matter-density-dependent) oscillation amplitudes and mass-squared-differences are matter-invariant, thereby accounting for the folklore “no-win” theorem , that increased mixing-angles in matter seem to be compensated by reduced masses and hence longer oscillation lengths, and vice-versa. The covariant quantity on the LHS of Eq. 50 transforms to the matter case as $`_0^{}_{ee}^{}Tr(^{})/3=_0+2𝒜/3`$, where $`𝒜=2\sqrt{2}G_FN_eE`$, while the quantities on the RHS take their density-dependent values.

The usual six parameters ($`\theta _{12},\theta _{23},\theta _{13},\delta ,\mathrm{\Delta }m_{21}^2\mathrm{and}\mathrm{\Delta }m_{31}^2`$), are popular for their intuitive relevance, but they all have complicated dependencies on the density of the matter through which the neutrinos propagate  . Different experiments operating over different baselines at different energies are all subject to different modifications of each of these conventional parameters due to matter effects. However, in each case, only one additional parameter is actually present, the (path-averaged) Wolfenstein term, $`𝒜`$. By contrast, as we have seen, the five combinations of observables, Eqs. 45-49 are all matter-invariant, while the sixth observable, Eq. 50, has a very simple dependence on the matter density. Clearly, these quantities could be used as the basis for a new parameterisation of neutrino oscillation phenomenology. This is especially true for the matter-invariants involving only off-diagonal elements: $`|_{\alpha \beta }|^2(\alpha \beta )`$ and $`\mathrm{Im}(_{e\mu }_{\mu \tau }_{\tau e})`$, which are determined directly in appearance experiments as the amplitudes of the leading $`CP`$-even oscillations (given a strong-enough mass-hierarchy) and the $`CP`$-asymmetry respectively. Not only are they matter-invariant, but they are more directly observed even than the conventional mixing angles.

To some degree then, one may say that the problem of deconvoluting matter effects on the standard mixing parameter set is a problem of man’s own making, which is solved (at least in the uniform-density case we have considered) by simply switching to a different parameter set. While it is probably unrealistic to assume that the new parameterisation proposed here will be immediately adopted universally, we believe that our parameters may have a practical value, beyond their evident conceptual importance as matter-invariant observables.

Acknowledgements It is a pleasure to thank R. Edgecock for useful comments. This work was supported by the UK Particle Physics and Astronomy Research Council (PPARC).