format: text

text:
# Laser Generated Magnetic Pulses: Comparison of Hot Electron Propagation in Conducting and Dielectric Material

## Abstract

We report experimental evidence of electrostatic inhibition of fast electrons, generated in a highly resistive material upon irradiation with an intense ultra-short ($`10^{16}W/cm^2`$, $`100fmsec`$) laser pulse. The experiment involves measurement of temporal evolution of self-generated magnetic pulses using pump-probe polarimetry. A comparison is made between the temporal behaviour of magnetic pulses generated with Aluminum and Glass targets. It is found that in contrast to Aluminium, self-generated magnetic pulse decays much faster in glass. This is attributed to the absence of return shielding currents in glass, which results in build up of electrostatic field, which in turn inhibits the movement of fast electrons. Fitting of experimental measurements using a one dimensional model, yields estimate of conductivity of Aluminium and glass, and penetration depth of hot electrons in these materials.

Interaction of intense ultrashort laser pulses with a solid leads to generation of large magnetic fields. Magnetic fields up to gigagauss magnitudes have been predicted in the overdense region of solid density plasmas Sudan ; Stamper , which are created by explosive ionzation of a solid target by an incident laser. Fields comparable to stellar magnitudes ($`300`$ MG) have already been experimentally realized in the laboratory with the use of super-intense ultra-short laser pulses Tatarakis . Recent studies have also established the ultrashort duration ($`picosec.`$) of these self-generated magnetic fields under femtosecond laser excitation Sandhu ; Borghesi . The temporal behaviour of the magnetic pulse contains crucial information about the dynamics of hot electron propagation inside a solid density material / plasma. This information about the penetration and transport of hot electrons Mason  inside an overdense plasma and factors affecting it, is of utmost importance to the fast ignitor approach of fusion research Kolka ; Tabak .

The pulsed high magnetic fields generated primarily near the critical density surface, result from a combination of currents, viz. the direct hot electron current which is due to copious generation of fast electrons at the critical layer via collisionless absorption of incident laser energy , and the return shielding current which is generated in response to the hot electron current Sentoku . The propagation of hot electrons inside the solid density target depends on whether the target is a dielectric or a conductor. In case of a dielectric ( material with high classical resistivity ), the return shielding current is weak, resulting in a large charge imbalance which in turn exerts a strong retarding field on the fast electrons, resulting in electrostatic inhibition. Electrostatic inhibition of hot electrons in dense matter has been seen in some experiments by measuring K $`\alpha `$ emission from layered targets. In initial experiments by Bond et. al. Bond  inhibition of fast electrons due to resistive electric field was observed in low density gold targets. In contrast to a dielectric, in a conductor ( material with extremely low classical resistivity ), the return shielding current nearly balances the hot electron current. Here the stopping of hot electrons occurs because of turbulence induced anomalous resistivity, the chaotic magnetic turbulence effects being generated by the presence of large gradients in the plasma currents Drake ; Sentoku2 . Experimental evidence of turbulence induced anomalous resistivity is presented in a recently published experiment by Sandhu et. al. Sandhu , where they irradiate a solid aluminium target with a high intensity femtosecond laser. Since the mechanism of hot electron transport ( along with return shielding current ) and their eventual stopping is different ( it is either electrostatically induced Bond  or turbulence induced Sandhu  ) in a dielctric and a conductor, it is expected that the temporal evolution of the magnetic fields which these currents generate will also be different. Thus the complex dynamics of the hot electron propagation inside a solid density material ( or plasma ) is reflected in the time evolution of the magnetic pulse.

Here we present experimental results on transport of hot electrons and their inhibition in both high and low resistivity media via measurements on temporal evolution of megagauss magnetic pulses. Recently, there has been a lot of experimental work on the transport of hot electrons Pisani , but to our knowledge no experimental result exists which gives a clear distinction between electrostatically induced and turbulence induced inhibition of fast electrons. In our experiment, this has become possible through the use of temporal measurement of magnetic pulses which is intimately related to the propagation of hot electrons inside the material. We compare two entirely different media, one conducting (solid Aluminum) and the other insulating (BK7 Glass). We use pump-probe polarimetry to measure the ellipticity produced in a probe laser pulse due to pump-generated magnetic fields. This measured ellipticity is then analyzed numerically to deduce the magnetic field as a function of time delay between pump and probe pulses. We find that the evolution of magnetic field in aluminium and glass to be drastically different. The time required for the magnetic field to decay in glass is found to be an order of magnitude smaller than in aluminium. This is a clear indication of electrostatic inhibition of fast electrons which results from absence of return shielding currents due to very low electrical conductivity of glass. We also present an one-dimensional model ( an improvement over our zero-dimensional model Sandhu  ) for magnetic field evolution which gives an interpretation of our observations. The decay time in each case is related to the electrical conductivity ( electrostatically induced or turbulence induced ), thus providing us with estimates of conductivity of glass and aluminium under conditions of elevated temperatures and pressures.

The experimental set-up is shown in figure 1. The p-polarized pump laser (100fs, 800nm) is focused at an angle of $`45^{}`$, to an intensity of $`2\times 10^{16}Wcm^2`$ on a spot size of 20 microns. The long term contrast for the pump is $`10^5`$. The probe laser used in measurements is frequency doubled (400nm) in a thin BBO crystal, and is made $`10^310^4`$ times weaker than the pump intensity and is focused at near normal incidence. In this configuration, the probe penetrates beyond the critical density ($`n_c`$) for the pump laser (i.e. upto 4$`n_c^{(806)}`$), and thus samples high density region near and beyond the critical layer, where high magnetic fields are expected to exist. It should be noted that this high density region is not accessible using tangential polarimetric measurements of Faraday rotation. The input laser fluctuations are monitored with a photodiode (PD1). The reflected probe from the interaction region is split into two arms, one arm is used in PD2 to measure the reflectivity of the probe and the second arm is used in the quarter wave plate and polarizer combination ( kept in front of photodiode PD3 ) to measure the polarization parameters of the probe. The reflectivity measurement serves to establish the zero delay and is used to normalize other polarimetric measurements. All the measurements are made as a function of time delay between pump and probe. Diffuse scatter is also monitored with another PD and is observed to be insignificant and constant in our temporal range of interest. In our experiment, we have used two type of targets, first one an aluminium target i.e. a conductor and second one a glass target i.e. an insulator. The BK7 glass target was coated with about half a micron thick aluminium on the front surface so that laser hitting the target generates identical hot electron temperatures in both cases. To establish that this is indeed so, we have measured the hard X-rays arising from hot electrons in the energy range from 50 to 500 keV for both the cases, under similar laser conditions with a NaI(Tl) detector. The hot electron spectra are shown in figure 2(a). The hot electron temperature obtained from the fit (shown as solid line) is observed to be approximately same ($`30keV`$) in both the cases (Al layer on Glass and Al metal alone). This indicates that we have identical hot electron source in both cases arising from interaction of the laser pulse with the front surface of the target. The magnetic field induced ellipticity in the probe pulse is measured using different combinations of quarter wave plate and polarizer in analyzer setup i.e. polarizer at $`0,45,90^{}`$ and $`\lambda /4`$ plate at $`45^{}`$ and polarizer at $`90^{}`$ Segre . The measurement of stokes parameters in this method rules out the presence of any random depolarization. Experimental curve showing ellipticity vs time delay is also shown in figure 2(b). The significant difference is observed in time scale when we compare the two cases, that of glass and aluminum. We now numerically analyze the measured ellipticity in order to deduce the magnetic field present in the interaction region.

An electromagnetic wave travelling through a magnetized plasma with wave vector perpendicular to the quasistatic magnetic field can be analyzed in terms of two characteristic modes, the ordinary mode (O-wave) with electric field polarization vector parallel to the magnetic field and the extraordinary mode (X-wave) with electric field polarization vector perpendicular to the magnetic field Chen . These two modes have different refractive indices, which are

$$\mu _O=\left[1\frac{\omega _{pe}^2}{\omega ^2}\right]^{1/2}$$
(1)
$$\mu _X=\left[1\frac{\omega _{pe}^2(\omega ^2\omega _{pe}^2)}{\omega ^2(\omega ^2\omega _{pe}^2\omega _{ce}^2)}\right]^{1/2}$$
(2)
where $`\mu _O`$ represents refractive index for $`O`$-mode and $`\mu _X`$ represents refractive index for $`X`$-mode. Here $`\omega `$, $`\omega _{pe}`$ and $`\omega _{ce}`$ are respectively the laser frequency, the electron plasma frequency, and the electron cyclotron frequency i.e. $`eB/mc`$. The difference in refractive indices results in accumulation of phase difference between O- and X- waves as they propagate through the plasma. This difference in phase is observed as induced ellipticity. If we represent the laser polarization in terms of stokes vector $`s=(s_1,s_2,s_3)`$ , $`s_1=cos(2\chi )cos(2\psi )`$, $`s_2=cos(2\chi )sin(2\psi )`$ , $`s_3=sin(2\chi )`$, where $`\psi `$ is orientation of polarization ellipse, then the ellipticity $`ϵ`$ is given by $`ϵ=tan(\chi )`$. The evolution of stokes vector inside the magnetized plasma is given by the difference in refractive indices of the characteristic waves i.e. $`ds/dz=\mathrm{\Omega }(z)\times s(z)`$, where $`\mathrm{\Omega }=(\omega /c)(\mu _O\mu _X)`$. The solution of this polarization evolution equation is written as $`s(z)=M.s_{input}`$, where $`M_{ij}`$ is the transition matrix which obeys the equation $`dM_i/dz=\mathrm{\Omega }\times M_i`$, where $`M_i=(M_{1i},M_{2i},M_{3i})`$ Segre . We integrate this equation numerically inside the plasma by dividing it into small slabs, where within each slab, the plasma parameters are assumed to be constant.

For the purpose of integration, the required plasma density profile is modelled by assuming that the plasma expands into vacuum at ion sound speed in a self similar fashion Mora ; Kruer , which yields an exponential density profile of the form $`n_e=n_sexp(z/c_st)`$, where $`n_s`$ is the solid density and $`c_s=\sqrt{(ZkT_e/m_i)}`$ is the ion sound speed. We have chosen the background electron temperature $`T_e=100eV`$ Multi  and the average ionization state $`Z=5`$ Milchberg . At the vacuum end, which is chosen at a distance $`L_s`$ from the critical layer, we put a cut off on the density ‘$`n_{end}`$’ which is chosen sufficiently low so that it makes practically no difference to the results. The plasma slab to be integrated over is given by $`L_s=c_st[ln(n_{critical}/n_{end})]`$. It is found that the results are not sensitive to the chosen density profile; even the use of linear density profile yields results within acceptable error bars. In addition to the density profile, we have assumed uniform magnetic field inside the plasma column. Although, this assumption is not strictly correct, it makes little difference to the results, and is a reasonable one as shown later. The solution $`M_{ij}(z)`$, and $`s(z)=M.s_{input}`$ can now be determined for any small plasma slab ( where the density is assumed to be constant ) for a given value of magnetic field B in definition of $`\mathrm{\Omega }`$. Thus recursively solving in small steps and maintaining unitarity of M one can find the full transition matrix for whole plasma column, for any chosen density profile. From the transition matrix, the final stokes vector can be computed which in turn gives the ellipticity of the reflected probe.

In addition to the phase difference introduced by differences in refractive indices, there is an extra phase difference which aries due to slightly different turning points for O and X-wave. The X-mode has a turning point i.e. right cut-off in the dispersion diagram, earlier than the O-mode, which reflects at usual critical density. This extra path travelled by O-mode as compared to X-mode leads to an additional phase difference, which we have included in our calculation of net ellipticity. Thus the effect of magnetic field manifests itself in the form of ellipticity of the reflected probe light by both of the above mentioned processes i.e. difference in refractive indices and difference in turning points.

Using the above described procedure, for each time delay, we iteratively calculate the magnetic field B required to yield the experimentally observed ellipticity. The final plots of magnetic field for solid aluminium target and BK7 glass coated with aluminium layer are shown in fig. 3(a) and fig. 3(b) respectively. The squares represent experimental points and the solid curves, which show a reasonable fit to the experimental data are obtained using a one-dimensional model developed below. It is observed that the peak magnetic field reached in both cases are of the same order ( $`32MG`$ in case of BK7 glass and $`22MG`$ in case of solid aluminium ). The most striking difference between the two cases is seen in the time durations of the self-generated magnetic pulses (note different time axis of fig 4(a) and (b)). In case of BK7 glass coated with alumunium, the time duration of the magnetic pulse ($`1picosec`$) is an order of magnitude smaller than the time duration ($`10picosec`$) of the magnetic pulse for pure aluminium. This is a clear indication of high background resistivity of glass, which curtails the return currents and hence results in inhibition of hot electrons because of electrostatic fields. This leads to much faster resistivite decay of magnetic fields in case of glass sample as compared to solid aluminium sample, where the decay is mainly due to turbulence induced anomalous resistivity Sandhu .

We now model the temporal evolution of magnetic field using the following equation, which describes the mechanism of quasistatic magnetic field generation under EMHD approximation.

$$\frac{\stackrel{}{B}}{t}=\frac{c^2}{4\pi \sigma }^2\stackrel{}{B}+\frac{c}{\sigma }\left(\stackrel{}{}\times \stackrel{}{j}_{hot}\right)$$
(3)
The above equation is derived by taking the curl of the equation of motion of background plasma electrons that carry the return shielding currents, $`\stackrel{}{j}_p=(c/4\pi )(\times \stackrel{}{B})\stackrel{}{j}_{hot}`$, where $`\stackrel{}{j}_p=\sigma \stackrel{}{E}`$, $`\sigma `$ being the conductivity of the background plasma. The first term describes the magnetic field decay due to resistive damping of the plasma shielding currents and the second term describes the magnetic field generation due to hot electron currents. Assuming the magnetic field to be in the azimuthal direction, equation (3) can be approximated as

$$\frac{B}{t}\frac{B_\varphi }{\tau }+S(z,t)$$
(4)
where the diffusion term is approximated as $`B/\tau `$ with $`\tau (4\pi \sigma /c^2)(\mathrm{\Delta }r)^2`$, and the source term is approximated as $`S(z,t)(c/\sigma (\mathrm{\Delta }r))j_{hot}(z,t)`$. Here $`\mathrm{\Delta }r`$ is the laser spot radius $`10\mu `$. Our main interest lies in understanding the evolution of the magnetic field after the laser is switched off. This is because the experimental results show that, although the laser is on only for $`100fmsec.`$, the magnetic field peaks and remains on for picoseconds time duration. Taking $`B=B_{las}`$ at $`t=\tau _{laser}`$, the solution of equation (4) is given by

$`B`$ $`=`$ $`B_{las}e^{(t\tau _{laser})/\tau }+e^{t/\tau }{\displaystyle _{\tau _{laser}}^t}S(t,z,r)e^{t/\tau }𝑑t`$ (5)
$``$ $`B_{las}e^{(t\tau _{laser})/\tau }e^{t/\tau }{\displaystyle _{\tau _{laser}}^t}{\displaystyle \frac{c}{\sigma \mathrm{\Delta }r}}j_{hot}(z,t)e^{t/\tau }𝑑t`$
where $`j_{hot}=en_hv_h`$, $`n_h`$ is the hot electron density and $`v_h`$ is the velocity of the hot electron fluid. To make an estimate of $`n_h`$ and $`v_h`$, we use the formalism given by Bell et. al. Bell . According to Bell , the evolution of hot electron density ($`n_h`$) is governed by the following nonlinear diffusion equation

$$\frac{n_h}{t}=\frac{}{z}\left(\frac{\sigma T_h}{e^2n_h}\frac{n_h}{z}\right)$$
(6)
Here $`T_h`$ is the hot electron temperature. It can be shown that the above equation is valid even without the restricted assumption $`j_{total}=j_{hot}+j_p0`$, used in ref. Bell . Since our interest lies in $`t>\tau _{laser}`$, we use the solution of equation (6) in this temporal regime, as given in ref. Bell

$$n_h=\frac{2n_0z_0}{\pi }\frac{L}{z^2+L^2}$$
(7)
with

$`L(t)`$ $`=`$ $`z_0\left[{\displaystyle \frac{5\pi \sigma T_h}{3e^2n_0z_0^2}}\left(t\tau _{laser}\right)+1\right]^{3/5}`$ (8)
$`n_0`$ $`=`$ $`{\displaystyle \frac{2}{9}}{\displaystyle \frac{I_{abs}^2\tau _{laser}e^2}{\sigma T_h^3}}`$ (9)
$`z_0`$ $`=`$ $`{\displaystyle \frac{3\sigma T_h^2}{e^2I_{abs}}}`$ (10)
where the absorbed intensity $`I_{abs}=fI_{incident}`$, $`f`$ being the fraction absorbed. Here $`n_0`$ is the density of hot electrons at $`z=0`$, at time $`t=\tau _{laser}`$ and $`z_0`$ is the characteristic stopping length such that $`n_0z_0`$ is the total number of hot electrons produced at time $`t=\tau _{laser}`$. The constants $`n_0`$ and $`z_0`$ have been derived by equating the absorbed laser energy “$`I_{abs}\tau _{laser}`$” to the hot electron kinetic energy (see ref. Bell ). The above solution (7) is a self-similar solution of equation (6) in which the spatial shape remains the same but it expands in time with a scale length $`L(t)`$. Using the above expression for $`n_h`$ and $`L(t)`$, we estimate $`j_{hot}`$ as

$`j_{hot}`$ $`=`$ $`en_hv_h`$ (11)
$`=`$ $`e{\displaystyle \frac{2n_0z_0}{\pi }}{\displaystyle \frac{L}{z^2+L^2}}\left[\alpha {\displaystyle \frac{dL}{dt}}\right]`$
where $`v_h`$, the hot electron velocity is taken to be proportional to $`\frac{dL}{dt}`$, $`\alpha `$ being the proportionality constant. Substituting the expression for $`\frac{dL}{dt}`$ in $`j_{hot}`$ and using it in equation (5) along with $`\sigma \frac{c^2}{4\pi (\mathrm{\Delta }r)^2}\tau `$, we get

$$B=B_{las}e^y+Ae^y\tau _0^y\frac{(py\tau +1)^{1/5}e^y}{z^2+z_0^2(py\tau +1)^{6/5}}𝑑y$$
(12)
where $`y=\frac{(t\tau _{laser})}{\tau }`$, $`A=\frac{2cz_0\alpha T_h}{e\mathrm{\Delta }r}`$ and $`p=\frac{5\pi \sigma T_h}{3e^2n_0z_0^2}`$. We now use the above expression for $`B(t,z)`$ at $`z=0`$ to model the magnetic field evolution as a function of time, for both BK7 glass coated with aluminium and solid aluminium using $`\tau `$ ( which is related to conductivity $`\sigma `$ ) and $`f`$ ( fraction of light absorbed ) as free parameters. $`z=0`$ is the point where hot electrons are generated; so it is actually the critical density point. The proportionality contstant $`\alpha `$ is taken to be unity. $`B_{las}`$, the magnetic field at time $`t=\tau _{laser}`$ is $`15MG`$ and $`18MG`$ respectively for aluminium and glass. The best fit curves ( solid lines ) are shown in fig. 3(a) and fig. 3(b) respectively. From the best fit, the relevant parameters ( $`f`$ and $`\tau `$ ) for aluminium and glass are:

Aluminum: f = 0.5; $`\tau `$ = 3.5 ps

Glass: f = 0.2; $`\tau `$ = 0.2 ps

Before presenting a discussion of our results obtained from one-dimensional modelling and numerical fit, we first justify the assumption we have made in our numerical analysis of ellipticity viz. the assumption of uniform magnetic field inside the plasma slab. Fig. 4 shows the typical magnetic field as a function of $`z`$ and $`t`$ obtained from our one-dimensional model. This plot is for aluminium with $`f`$ and $`\tau `$ as calculated above. As expected, the magnetic field diffuses in $`z`$ and decreases in amplitude as time increases. Although the magnetic field $`B`$ varies with $`z`$, we shall see below that its variation is not significant over the range of scale lengths which are relevant for our experiments. In Fig. 5, using equation (12) we have plotted the magnetic field as a function of position at a time corresponding to the peak of the magnetic field, for both aluminium and glass. In this plot, we have used the same value of $`f`$ and $`\tau `$ as obtained above. This plot clearly shows that the variation of magnetic field over the plasma slab is negligible. This justifies that the assumption of uniform magnetic field in our analysis of ellipticity is a reasonable one.

We now discuss our results. From the results of the fit, it is clear that the time for magnetic field decay in aluminium is an order of magnitude larger than in glass. Calculation of conductivity using $`\sigma =(c^2/4\pi (\mathrm{\Delta }r)^2)\tau `$, gives $`\sigma _{Al}=2.5\times 10^{14}s^1`$ and $`\sigma _{glass}=1.43\times 10^{13}s^1`$. As of now there are no measurements of conductivity at elevated temperatures for glass and very few in aluminium. We see that resistivity of aluminium as deduced from these measurements is an order of magnitude higher than the value of resistivity deduced from reflectivity measurements by Milchberg et. al. Milchberg . This, we believe, is due to turbulence induced anomalous resistivity which affects the return currents. This conclusion is also supported by simulation where similar effect of anamalous resistivity and stopping of hot electrons have been seen Sentoku . In the case of glass, the situation is more complex as at low temperatures glass is almost non-conducting $`(\sigma 0)`$. However, at elevated temperatures, as a result of target heating by collisional effects and large electric fields exceeding breakdown threshold, conductivity becomes finite. There exists observation of $`8\mu `$ wide collimated ionization tracks extending upto $`150300\mu `$ and lasting for $`60ps`$ in glass with a $`10^{17}Wcm^2`$ intensity laser Teng . Although we cannot estimate here the amount of target heating or ionization, we do get an estimate of resitivity, which we assume to be uniform in space and time. Our measurements show that neutralization of hot electron current is clearly not as effective in glass as it is in aluminium, due to high background resistivity. This results in inhibition of hot electron propagation in glass through generation large of electrostatic fields. We now make an estimate of the penetration depth ($`L_f`$) of hot electrons using equation (8). Using the parameters $`z_0`$ and $`\sigma `$ from the fit and $`T_{hot}`$ from observations, we get, for glass, $`L_f=4.7\times 10^5cm`$ in $`1ps`$ and $`2.1\times 10^4cm`$ in $`10ps`$, and for aluminium $`L_f=1.6\times 10^3cm`$ in $`10ps`$. These numbers indicate that hot electrons are penetrating an order of magnitude more distance in a conducting background. Similar results on electrostatic inhibition have also been measured in plastic targets using $`K_\alpha `$ emission from layered targets. This technique has its limitations due to its reliance on x-ray measurements which is a secondary diagnostic (with limitations of sensitivity, calibration, and reabsorption) and on numerical codes for interpreting X-ray data.

In conclusion, the penetration of fast electrons into a target can result in very high electric field, which is not quenched by thermal return currents. This can lead to severe impediment in the flow of fast electron currents. We provide evidence of inhibition of fast electrons in an insulating media. Our results on magnetic pulse measurement and modelling in terms of hot electron currents as source, yield measurements of conductivity of conducting and dielectric media under extreme conditions. These results are very important for laser fusion schemes, which rely on ignition of hot spot by energy deposition by fast electrons.

We would like to thank D. Mathur, M. Krishnamurthy,and P. P. Rajeev for discussions and useful suggestions.


text:
# Optimizing the ensemble for equilibration in broad-histogram Monte Carlo simulations

## I Introduction

At first-order phase transitions and in systems with many local minima of the free energy such as frustrated magnets or spin glasses, conventional Monte Carlo methods simulating canonical ensembles have very long equilibration times. Several simulation methods have been developed to speed up such systems, including the multicanonical method Multicanonical , broad histograms BroadHistograms , parallel tempering Tempering , the Wang-Landau algorithm WangLandau  and variations thereof WLVariations . Most of these methods simulate a flat-histogram ensemble: Instead of sampling a configuration of energy $`E`$ with Boltzmann weight $`w(E)\mathrm{exp}(\beta E)`$, they use weights $`w(E)1/g(E)`$, where $`g(E)`$ is the density of states. The probability distribution of the energy, $`n(E)=w(E)g(E)`$, then becomes constant, producing a flat energy histogram. Naively, one might assume that sampling all energies equally often produces an unbiased random walk in energy. However, it was recently shown Limitations  that the growth with the number of spins $`N`$ of the “tunneling times” between low and high energy in any local-update flat-histogram method is stronger than the naive $`N^2`$ of an unbiased random walk in energy for various 2D Ising models: as $`N^{2.4}`$ for the ferromagnetic and $`N^{2.9}`$ for the fully-frustrated models. For the 2D $`\pm J`$ spin glass, exponential growth was observed Limitations .

In view of these results for flat-histogram simulations Limitations  we have asked how this type of simulation can be improved, both in terms of computer time and statistical errors. The general type of application we have in mind is to the equilibrium behavior of a system that is very slow to equilibrate in a conventional simulation, such as domain walls in ordered phases, low-energy configurations of frustrated systems or a spin-glass ordered phase. Our new algorithm instead simulates a broad-histogram ensemble, where the system can, at “equilibrium” in this ensemble, wander to part of its phase diagram where equilibration is rapid. We look specifically at histograms that are broad in energy, but in general another variable other than the energy could be used. To minimize the statistical errors of measurements in the energy range of interest, one maximizes the number of statistically-independent visits. For a glassy phase, the system will relax very little as long as it remains in that phase, so to get a new statistically-independent visit to the phase the system has to leave it and equilibrate elsewhere (usually at high energies). Thus the quantity we want our simulation to maximize is the number of round trips – between low and high energy – per unit computer time. This should minimize both the simulation’s equilibration time and the statistical errors in the low-energy regime of interest.

In this paper, we present an algorithm that systematically optimizes the ensemble simulated to maximize the rate of round trips in energy. We use a feedback loop that reweights the ensemble based on preceding measurements of the local diffusivity of the total energy. This detects the “bottlenecks” in the simulation as minima in the diffusivity (at critical points in the cases we study), and reallocates resources to those energies in order to minimize the slowdown. We find that the resulting statistical errors in the density of states as estimated by this new algorithm are nearly uniform in energy, in strong contrast to flat-histogram simulations where the errors are much larger at low energy than at high energy. While our algorithm is rather general and should be widely applicable to study complex systems, we have developed and tested it on ferromagnetic (FMI) and fully-frustrated (FFI) square-lattice Ising models.

## II Feedback of local diffusivity

In our simulations, the system’s energy does a random walk in the energy range between two extremal energies, $`E_{}EE_+`$ where we take the lowest energy $`E_{}`$ to be the system’s ground state, although this is not necessary for our approach. Consider a general ensemble, with weights $`w(E)`$, which define the acceptance probabilities for moves based on the Metropolis scheme

$$p(EE^{})=\mathrm{min}(1,\frac{w(E^{})}{w(E)}).$$
(1)
Our algorithm iteratively collects data from batch runs which simulate with a fixed ensemble. During a simulation detailed balance is strictly satisfied at all times. For a reasonably large number of sweeps we can thus measure the equilibrium distribution of the energy in this ensemble which is $`n_w(E)w(E)g(E)`$. The simulated system does a biased and Markovian random walk in configuration space. Since we bias this walk based only on total energy, the projection of this random walk onto that variable is what we will discuss. This projection, which ignores all properties of the state other than its total energy, results in a random walker that is non-Markovian, with its memory stored in the system’s configuration. Thus the simulation may be viewed as a biased non-Markovian random walker moving along the allowed energy range between the two extremal energies.

To measure the round trips we add a label to the walker that says which of the two extremal energies it has visited most recently. The two extrema act as “reflecting” and absorbing boundaries for the labeled walker: e.g., if the label is plus, a visit to $`E_+`$ does not change the label, so this is a “reflecting” boundary. However, a visit to $`E_{}`$ does change the label, so the plus walker is absorbed at that boundary. The steady-state distributions of the labeled walkers satisfy $`n_{}(E)+n_+(E)=n_w(E)`$. It is important to note that the behavior of the labeled walker is not affected by its label except when it visits one of the extrema and the label changes. When the unlabeled walker is at equilibrium, the labeled walker is in a nonequilibrium steady state. Let $`f(E)=n_+(E)/n_w(E)`$ be the fraction of the walkers at $`E`$ that have label plus, so they have most recently visited $`E_+`$. The above-discussed boundary conditions dictate $`f(E_{})=0`$ and $`f(E_+)=1`$.

To calculate the rate of round trips, we note that in steady state the current $`j`$ of the labeled walkers is independent of $`E`$. The plus and minus walkers drift in opposite directions and the equilibrium unlabeled walker has no net current. We first examine the case of a continuous energy $`E`$. The steady-state current from $`E_+`$ to $`E_{}`$ to first order in the derivative is

$$j=D(E)n_w(E)\frac{df}{dE},$$
(2)
where $`D(E)`$ is the walker’s diffusivity at energy $`E`$. There is no current if $`f`$ is constant, since this is equilibrium; this is why the current is to leading order proportional to $`df/dE`$. If one rearranges the above equation and integrates on both sides, noting that $`j`$ is a constant and $`f`$ runs from 0 to 1, one obtains

$$\frac{1}{j}=_E_{}^{E_+}\frac{dE}{D(E)n_w(E)}.$$
(3)
In the following we separately discuss how we can maximize the rate of round-trips for Metropolis and $`N`$-fold way dynamics based on this estimate of the current.

### II.1 Metropolis dynamics

For Metropolis dynamics the rate of round-trips is simply proportional to the current. To maximize the round-trip rate, the above integral, Eq. (3), must be minimized. However, there is a constraint: $`n_w(E)`$ is a probability distribution and must remain normalized. We do this by adding a Lagrange multiplier:

$$_E_{}^{E_+}𝑑E\left(\frac{1}{D(E)n_w(E)}+\lambda n_w(E)\right).$$
(4)
To minimize this integrand, the ensemble, that is the weights $`w(E)`$ and thus $`n_w(E)`$ are varied. At this point we assume that the dependence of $`D(E)`$ on the weights can be neglected. This is justified by noting that the rates of transitions between configurations depend only on the ratios of weights, so the diffusivity $`D(E)`$ is unchanged when the weights are multiplied by an energy-independent constant. By ignoring the variation of $`D(E)`$ with the weights, we are assuming that the adjustments to the weights are slowly-varying in energy, which is true for most systems, particularly for large systems where the energy range being studied is large. With this assumption, the optimal weighting which minimizes the above integrand is

$$n_w^{(opt)}=\frac{1}{\sqrt{D(E)\lambda }}=\frac{df^{(opt)}}{dE}.$$
(5)
Thus for the optimal ensemble with Metropolis dynamics, the probability distribution is simply inversely proportional to the square root of the local diffusivity.

### II.2 $`N`$-fold way dynamics

Since Metropolis dynamics can be slowed down by high rejection rates of singular moves, e.g. in the vicinity of the fully polarized ground state of the FMI, or the occurrence of multiple, generally accepted zero-energy moves it can be advantageous to introduce rejection-free single-spin flip updates such as the $`N`$-fold way NFoldWay . $`N`$-fold way dynamics involve two time scales, the walker’s time and the computer time. At a given energy level the two time scales differ by the (energy-dependent) lifetime of a given spin configuration. The random walk with $`N`$-fold way dynamics is an equilibrium process when measured in walker’s time, that is the equilibrium distribution, $`n_w(E)`$, is proportional to the amount of walker’s time the walker spends at $`E`$. However, for the ensemble optimization with $`N`$-fold way dynamics we want to speedup the random walk measured in computer time. This setup with two clocks requires a slightly different reweighting procedure than is presented above for Metropolis dynamics.

As for the Metropolis dynamics the amount of walker’s time it takes to make a round trip is proportional to $`1/j`$ given in Eq. (3). However, we are interested in minimizing the amount of computer time spent, so we need to multiply this by the ratio of computer time to walker’s time at $`E`$ which we denote as $`t(E)`$. Let us assume the distribution $`n_w(E)`$ is normalized to integrate to one. Then for one unit of walker’s time, the fraction spent at $`E`$ in $`dE`$ is $`n_w(E)dE`$. The amount of computer time used per unit walker’s time is thus

$$T=_E_{}^{E_+}n_w(E)t(E)𝑑E.$$
(6)
To find the weights that minimize the round-trip time as measured in computer time, the full quantity we want to minimize is thus

$$_E_{}^{E_+}\frac{dE}{D(E)n_w(E)}_E_{}^{E_+}t(E^{})n_w(E^{})𝑑E^{}.$$
(7)
Since the probability distribution $`n_w(E)`$ occurs in both the numerator and denominator of the integrand there is no need to enforce its normalization by a Lagrange multiplier. To extremize the integrand, we will again vary the weights $`w(E)`$, which gives the following condition for the optimum:

$$\frac{T}{D(E)n_w^2(E)}=\frac{t(E)}{j},$$
(8)
so the weights should be chosen to give the optimal distribution

$$n_w^{(opt)}(E)=\sqrt{\frac{jT}{D(E)t(E)}}.$$
(9)
For the optimal ensemble with $`N`$-fold way dynamics, the probability distribution is thus larger at the points with smaller $`t(E)`$ (since they do not cost a lot of computer time) and smaller diffusivity $`D(E)`$.

### II.3 Feedback iteration

To feed back we simulate with some trial weights $`w(E)`$, get steady-state data for $`n_w(E)`$ and $`f(E)`$ and thus obtain estimates for the diffusivity via

$$D(E)=\frac{j}{n_w(E)\frac{df}{dE}}.$$
(10)
For Metropolis dynamics chose new weights $`w^{}(E)`$ so that

$$n_w^{}(E)=A\sqrt{n_w(E)\frac{df}{dE}},$$
(11)
where $`A`$ is a normalization constant whose value is not needed to run the next “batch” of the simulation with the new weights $`w^{}(E)`$. For $`N`$-fold way dynamics the new weights $`w^{}(E)`$ are chosen to be

$$n_w^{}(E)=\sqrt{n_w(E)\frac{df}{dE}\frac{T}{t(E)}}.$$
(12)
In practice we work with the logarithms of the weights, so the reweighting becomes

$`\mathrm{ln}w^{}(E)`$ $`=`$ $`\mathrm{ln}w(E)`$
$`+{\displaystyle \frac{1}{2}}\left(\mathrm{ln}\left\{{\displaystyle \frac{df}{dE}}\right\}\mathrm{ln}n_w(E)\mathrm{ln}t(E)\right),`$
where all energy independent terms have been dropped as they introduce a constant shift only. For Metropolis updates the last term $`\mathrm{ln}t(E)`$ can also be dropped. Each subsequent batch should be run significantly longer than the previous one – in our implementation we double the number of sweeps – in order to get better statistics, and fed back to improve the estimates of the optimal weights.

## III Implementation and Applications

We implemented this algorithm for 2D Ising models with single-spin-flip Metropolis and $`N`$-fold way dynamics, found the optimal ensembles for the FMI and FFI models, obtained the scaling of round-trip times and calculated the density of states and its statistical errors for both models. In both cases we used the ground state, $`E_{}=E_0`$, and zero energy, $`E_+=0`$, as the energy limits.

In the initial batch mode step we simulated a flat-histogram ensemble for small system sizes using either the exact density of states footnote  or a rough estimate thereof calculated with the Wang-Landau algorithm WangLandau . For larger systems ($`N>64\times 64`$) where the round-trip times for the flat-histogram ensemble are more than a magnitude larger than for the optimized ensemble, we produced an initial estimate of the optimal weights by extrapolating the optimized weights of smaller systems. For all batch mode steps the fraction of labeled walkers, $`f(E)`$, was determined by recording two histograms, one for the equilibrium (unlabeled) walker and one for the labeled walker which most recently visited $`E_{}`$. The derivative $`df/dE`$ was then estimated by a linear regression of several neighboring points at each energy level. The number of points used for the regression can be reduced for subsequent batch mode steps as the estimate of $`f(E)`$ becomes increasingly accurate due to better statistics. In the final batch mode steps of our calculations the regression was performed using a minimum of three points. In general, there is a trade-off between the accuracy in the measurement of the local diffusivity and the number of feedback steps. For the Ising models we study we found rapid convergence to the optimal ensemble. For small systems, $`N32\times 32`$, an initial batch mode step with some $`10^5`$ to $`10^6`$ sweeps was sufficient to find the optimal weights after the first feedback step. Since the possible energy levels are discrete for the Ising model, special care is taken when applying the reweighting derived for the continuum limit, particularly at the boundaries of the energy interval $`[E_{},E_+]`$. However, we did not encounter any subtlety for either model.

### III.1 Fully frustated Ising model

We first present our results for the fully-frustated model (FFI), which has a critical point at its ground state, and shows rather simple scaling of our algorithm’s behavior with energy and system size. For the optimized ensemble of the FFI the histogram of the equilibrium random walker is no longer flat, but exhibits a power-law divergence at its ground state, as shown in Fig. 1. This divergence reflects the behavior, $`D(E)[(EE_0)/N]^2`$ of the diffusivity, as is seen in the inset of Fig. 1. These power-law behaviors extend from the first few points, $`EE_0=O(1)`$, up nearly to zero energy, $`EE_0=O(N)`$. If we accept that the critical exponent for the diffusivity is indeed 2, then the optimal distribution scales as $`n_w1/[(EE_0)\mathrm{ln}N]`$, and the round-trip time as $`\tau (N\mathrm{ln}N)^2`$, consistent with our results shown in Figs. 1 and 2.

Noting that for our optimized ensemble the system spends a large fraction of its time near the ground state where many Metropolis moves are rejected, we applied a version of our algorithm that instead uses single-spin-flip rejection-free $`N`$-fold way updates. We find the $`N`$-fold way updates do give a significant speedup compared to Metropolis dynamics, but do not change the $`\tau (N\mathrm{ln}N)^2`$ scaling of the round-trip time. In comparison to the performance of flat-histogram sampling we find a substantial speedup up to a factor of around $`50`$ for the largest simulated system with $`N=128\times 128`$ spins, see Fig. 2.

### III.2 Ferromagnetic Ising model

We now turn to the results for the ferromagnetic Ising model (FMI) which exhibits a finite-temperature second order phase transition. After applying the feedback, we find a peak in the histogram near the critical energy, as shown in Fig. 3. For Metropolis updates a second divergence close to the fully polarized ground state appears which is eliminated by changing the dynamics to rejection-free $`N`$-fold way moves. However, the minimum in the diffusivity at the critical point remains with $`N`$-fold way dynamics and the resulting peak in the histogram is not suppressed. With increasing system size this power-law divergence moves towards the critical energy of the infinite system, $`E_c/2N0.71`$, as illustrated in the inset of Fig. 4. For both types of single-spin-flip moves we find that the rate of round trips between the magnetically ordered and disordered phases of the ferromagnet appear to scale as $`\tau (N\mathrm{ln}N)^2`$ as for the FFI model, see Fig. 4.

## IV Statistical errors

Finally we address the statistical errors of measurements performed during the simulation. Standard tools can be used for the error analysis as the simulated random walk in configuration space is a conventional Markov chain Monte Carlo simulation. Only the projection of this random walk onto energy space becomes non-Markovian which is irrelevant for the measurements.

For each batch mode step simulating a fixed statistical ensemble $`w(E)`$ we can measure the density of states, $`\mathrm{ln}g(E)=\mathrm{ln}n_w(E)\mathrm{ln}w(E)`$, from the recorded equilibrium distribution $`n_w(E)`$. Comparing our results with the exact density of states we find perfect agreement within the statistical errors as illustrated for the FMI in the inset of Fig. 5. The observed distribution of statistical errors is nearly flat in energy, which is a further improvement compared to flat-histogram simulations where the errors can be orders of magnitude larger at low energy than at high energy WangLandau . The statistical error is found to scale as $`\mathrm{\Delta }\mathrm{ln}g(E)1/\sqrt{\text{round trips}}`$ with the number of round trips in energy which is shown in the main panel of Fig. 5. For different system sizes we find the statistical errors to collapse onto a single $`1/\sqrt{\text{round trips}}`$ dependence which a posteriori validates our goal of maximizing the rate of round trips.

## V Conclusions

The presented algorithm should be widely applicable to study the equilibrium behavior of complex systems, such as glasses, dense fluids or polymers. To speed up the system’s equilibration the rate of round trips in energy is maximized by systematically optimizing the statistical ensemble based on measurements of the local diffusivity. We find that the relative statistical error in the density of states as calculated with the new method scales as $`O(1/\sqrt{\text{round trips}})`$. For the 2D ferromagnetic and fully frustrated Ising models the round-trip time from the ground state to the maximum entropy state scales like $`O([N\mathrm{log}N]^2)`$ which is a significant speedup compared to the power law behavior $`O(N^{2+z})`$ of flat-histogram algorithms.

The idea of performing round-trips in energy is similar to the parallel tempering algorithm Tempering  which simulates replicas of the system at various temperatures. The swapping of replicas at neighboring temperatures can be viewed as a random walk of the replicas along the temperature. In order to maximize the round-trips in temperature one can use our algorithm to systematically optimize the simulated temperature set which we will discuss in a forthcoming publication FutureTempering .

## VI Acknowledgments

We thank S. Sabhapandit for providing the exact density of states for the FFI and R. H. Swendsen for helpful discussions. ST acknowledges support by the Swiss National Science Foundation. DAH is supported by the NSF through MRSEC grant DMR-0213706.


text:
# ’Superluminal’ tunnelling, quantum measurements and the speed of information transfer

## Abstract

We exploit the analogy between the transfer of a pulse across a scattering medium and Aharonov’s weak measurements to resolve the long standing paradox between the impossibility to exceed the speed of light and the seemingly ’superluminal’ behaviour of a tunnelling particle in the barrier or a photon in a ’fast-light’ medium. We demonstrate that ’superluminality’ occurs when the value of the duration $`\tau `$ spent in the barrier is uncertain, whereas when $`\tau `$ is known accurately, no ’superluminal’ behavior is observed. In all cases only subluminal durations contribute to the transmission which precludes faster-than-light information transfer, as observed in a recent experiment.

Recent experiments  on transmitting information-containing features of an optical pulse across the ’fast-light’ medium, in which the group velocity exceeds the vacuum speed of light $`c`$, have renewed the interest in the so-called ’superluminal’ propagation phenomenon. It is well known that a wavepacket transmitted across the potential barrier, undersized wave guide or fast-light medium may arrive in a detector ahead of the one that propagates freely, as if it has crossed the scatter infinitely fast. This phenomenon, often referred to as ’apparent superluminality of quantum tunnelling’ was first noticed more than seventy years ago  and often discussed since (for Reviews see , a recent selection of different views on the subject can be found in ). It seems to raise the question about the possibility of faster-than-light travel, e.g., in classically forbidden regions, inaccessible to a classical particle. It is, however, broadly understood that the paradox results from an incorrect identification of the transmitted peak with the incident one, since the incident pulse undergoes severe reshaping in the barrier region. Although it has been suggested , that a barrier may, in some manner, carve the transmitted pulse from the forward tail of the incident one, the question of how exactly this is achieved remains open. The precise mechanism of this reshaping, its relation to Aharonov’s ’weak’ measurements  of the time delay and the implications for the information transfer, studied in , are the subjects of this Letter. Earlier work on the weak nature of the tunnelling times can be found in -, and a more recent approach relating ’superluminality’ to superoscillations is given in .
We start by revisiting  the analysis of a quantum system prepared at $`t=0`$ in the initial state $`|I>`$ and then post-selected (observed) at $`t=\mathrm{T}`$ in the final state $`|F>,`$

$$|I>=\underset{\nu }{}a_\nu |\nu >|F>=\underset{\nu }{}b_\nu |\nu >,\widehat{A}|\nu >=A_\nu |\nu >.$$
(1)
If at some $`t`$, $`0<t<\mathrm{T}`$, the system is subjected to a von Neumann-type measurement  of an operator $`\widehat{A}`$, the state of the pointer with position $`\tau `$ after post-selection is given by

$$<\tau |M>=\underset{\nu }{}G(\tau A_\nu )\eta _\nu \eta _\nu b_\nu ^{}a_\nu $$
(2)
where $`G(\tau )`$ is the initial (e.g., Gaussian) state of the meter at $`t=0`$ and $`_\nu `$ must be replaced by an integral $`𝑑\nu `$ if $`\widehat{A}`$ has a continuous spectrum. The meter is then read, i.e., the pointer position is accurately determined. For an accurate ’strong’ measurement, the width $`\mathrm{\Delta }\tau `$ is small (compared to the separation between the eigenvalues or, if the spectrum is continuous, to the scale on which $`\eta _\nu `$ varies considerably) and the meter’s readings occur close to the eigenvalues $`A_\nu `$. For an inaccurate, or ’weak’, measurement, $`\mathrm{\Delta }\tau `$ is large, and the interference between overlapping Gaussians in Eq.(2) may, for a special choice of $`|I>`$ and $`|F>`$, produce anomalous readings in the regions where $`\widehat{A}`$ has no eigenvalues . Next we show that when using the coordinate of the tunnelled particle to estimate the time delay it has experienced in the barrier, $`\tau `$, one, in fact, performs a weak measurement of $`\tau `$. The anomalously small value of $`\tau `$ that is obtained, results, just as an Aharonov’s weak value, from the quantum uncertainty inherent to the procedure.
Consider a one-dimensional wave packet transmitted across a short-ranged potential barrier $`V(x)`$, contained inside the region $`0<x<b`$. At $`t=0`$ the particle is prepared as an incident wave packet with a mean momentum $`k_0`$, centred at $`x=x_I<0`$

$$\mathrm{\Psi }_0(x)<x|I>=_{\mathrm{}}^{\mathrm{}}C(kk_0)exp(ikx)𝑑k.$$
(3)
where the factor $`C(kk_0)`$ insures that only positive momenta contribute to the integral. At some large $`t=\mathrm{T}`$, it is post-selected in the transmitted state ($`\mathrm{}=1`$)

$$\mathrm{\Psi }_F(x)<x|F>=_{\mathrm{}}^{\mathrm{}}T(k)C(kk_0)exp[ikxiE(k)\mathrm{T}]\mathrm{dk},$$
(4)
where $`T(k)`$ is the barrier transmission amplitude. Using the convolution property of the Fourier integral, it is convenient to rewrite Eq.(4) as

$$<x|F>=\mathrm{\Psi }_\mathrm{T}(xx^{})\xi (x^{})𝑑x^{}$$
(5)
where

$$\mathrm{\Psi }_\mathrm{T}(x)=_{\mathrm{}}^{\mathrm{}}C(kk_0)exp[ikxiE(k)\mathrm{T}]\mathrm{dk},$$
(6)
is the state that would evolve from the initial one under free propagation, and $`\xi (x^{})`$ is the Fourier transform of $`T(k)`$. We can continue the discussion in terms of the ’tunnelling times’ by identifying ($`v_0`$ is the velocity corresponding to the wavevector $`k_0`$)

$$\tau (x^{})x^{}/v_0$$
(7)
with the delay experienced by the particle in the barrier. Changing the variables in Eq.(5) and separating the inessential phase associated with the free motion, we obtain

$$\mathrm{\Psi }_F=\mathrm{exp}[ik_0xiE(k_0)\mathrm{T}]\mathrm{G}(\tau (\mathrm{x})\tau )\eta (\tau )d\tau $$
(8)
where $`G`$ is the envelope of $`\mathrm{\Psi }_T`$,

$$G(\tau )\mathrm{exp}[iE(k_0)\mathrm{T}+\mathrm{ik}_0\mathrm{v}_0\tau ]\mathrm{\Psi }_\mathrm{T}(\mathrm{v}_0\tau )$$
(9)
and

$$\eta (\tau )(2\pi )^1v_0T(k)\mathrm{exp}[i(k_0k)v_0\tau ]𝑑k.$$
(10)
Comparing Eq.(8) with Eq.(2) shows that the relation between the time delay $`\tau `$ and the particle’s position $`x`$ is that between the measured quantity, whose amplitude distribution is $`\eta (\tau )`$ and the position of the pointer, whose initial state is determined by the envelope of the initial pulse. As in the original Aharonov’s approach, the particle is post-selected in its transmitted state. However, unlike in Ref.-, no external pointer variable is employed and its role is played by the particle’s own position. Equivalently, registering the transmitted particle at a location $`x`$ amounts to measuring the time delay $`\tau `$ of a particle with the momentum $`k_0`$. Finding the particle roughly the the width of the barrier $`b`$ ahead of the free one, as it happens in tunnelling , corresponds to a negative time delay of $`b/v_0`$. Importantly, the accuracy to which the time delayed is evaluated is limited (if the spreading of the wavepacket is neglected ) by the uncertainty $`\mathrm{\Delta }x`$ of the particle’s position in its initial state , $`\mathrm{\Delta }\tau \mathrm{\Delta }x/v_0`$.
In the case of tunnelling, the momentum spread of the initial wavepacket must be at least  small enough for all its components to tunnel, rather than to pass over the barrier. Such wavepackets, broad in the coordinate space, correspond to inaccurate weak measurements which may produce anomalous ’superluminal’ readings, even when the amplitude distribution $`\eta (\tau )`$ contains only non-negative time delays. Indeed, in the complex $`k`$-plane, $`T(k)`$ may only have poles on the positive imaginary axis and in the lower half-plane . The poles of the first kind (I) correspond to the bound states supported by $`V(x)`$, while those of the second kind (II) correspond to scattering resonances. Closing the contour of integration in Eq.(10) as appropriate, we obtain

$`\eta (\tau )=2\pi i{\displaystyle \underset{I}{}}Res_nT\mathrm{exp}(ik_nv_0\tau ),\tau <0`$ (11)
$`=2\pi i{\displaystyle \underset{II}{}}Res_nT\mathrm{exp}(ik_nv_0\tau ),\tau >0.`$ (12)
where $`Res_nT`$ denotes the residue of $`T(k)`$ at the $`n`$-th pole. Thus, the poles (I), if present, produce negative time delays and are responsible, in the classical limit, for the speed up of a particle passing above a potential well. It is interesting to note that a potential well too shallow to support a bound state would not speed up a passing wavepacket. If, on the other hand, no bound states are present, then

$$\eta (\tau )0\text{ for }\tau <0$$
(13)
and the ’spectrum’ of the time delays in Eq.(8) is confined in the $`0\tau <\mathrm{}`$ semi-axis. Note that Eq.(13) demonstrates the causal nature of the scattering process, since the condition $`Imk_n<0`$, used in its derivation, also ensures that $`ImE(k_n)<0`$ for $`Rek_n>0`$ and the resonance states containing outgoing waves, $`Rek_n>0`$, are emptied, rather than filled up, as the time increases .
It is clear now that with a careful choice of $`T(k)`$, and, therefore, $`\eta (\tau )`$ the ’superluminal’ pulse can be produced, in an explicitly causal manner, from the front tails of $`\mathrm{\Psi }_T(xx^{})`$, all delayed relative to free propagation. One such system is a particle tunnelling across a potential barrier. However, our approach only relies on the analyticity of the transmission amplitude $`T(k)`$ and can also be applied to evanescent propagation in waveguides  and optical propagation through ’fast-light’ media (see ,), which have additional advantage of dealing with wavepackets not subject to spreading in vacuum.
Information transfer is often associated with propagation of non-analytic features, such as cut-offs  and next we will show that Eq.(13) ensures that it cannot be transferred faster than light. For a simple example, consider the propagation of an electromagnetic pulse across two narrow semi-transparent mirrors, broadly similar to the setup studied in . If the mirrors are modelled by $`\delta `$-functions of magnitude $`\mathrm{\Omega }`$ located at $`x=0`$ and $`x=b`$, respectively, $`T(k)`$ is given by the multiple scattering expansion

$$T(k)=\underset{m=0}{\overset{\mathrm{}}{}}T^{(m)}(k)(1+R(k))\underset{m=0}{\overset{\mathrm{}}{}}R(k)^{2m}\mathrm{exp}(2imkb)$$
(14)
where

$$R(k)=i\mathrm{\Omega }/(2k+i\mathrm{\Omega })$$
(15)
is the reflection amplitude for a single $`\delta `$-function placed at the origin $`x=0`$. Accordingly, the distribution $`\eta (\tau )`$ is decomposed into subamplitudes $`\eta _m(\tau )`$, $`\eta _m(\tau )0`$ for $`\tau <2mb/c`$, each peaked near $`\tau _m=2mb/c`$ (Fig.1a). For a large $`\mathrm{\Omega }`$,$`\mathrm{\Omega }b1`$, the widths can be neglected and the incident pulse is split into a number of discrete path modes corresponding to $`2m`$, $`m=0,1,2\mathrm{}`$ additional reflections experienced by the ray between $`x=0`$ and $`x=b`$ , and Eq.(4) reduces to Eq.(2) for a variable with a discrete spectrum $`\{\tau _m\}`$,

$$\mathrm{\Psi }_F(x)\underset{m}{}G(\tau (x)\tau _m)T^{(n)}(k_0).$$
(16)
A numerical evaluation of $`\mathrm{\Psi }_F(x)`$, for $`\mathrm{\Omega }b=100`$, corresponding to a weak Gaussian measurement,

$$G(x)=\mathrm{exp}(x^2/\mathrm{\Delta }x^2),\mathrm{\Delta }x>b,$$

shows (see Fig1.b) how a set of nearly Gaussian shapes, each delayed by $`\tau _m>0`$, interfere to produce a negative time delay $`b/c`$. Note that this is the best speed up which can be achieved with the model Eq.(14) .
It is now straightforward to show that this speed up effect cannot be used to send information faster than light. Aharonov and co-workers have already demonstrated  that a weak von Neumann pointer cannot be used for this purpose. Rather, they argued, the meter acts as a filter, extracting, in a non-trivial manner, the signal, otherwise hidden by a noise. The same argument applies to ’superluminal’ propagation. If the incident wave packet is chosen to be the rear half of a Gaussian with a sharp front, $`G_{}(x)=[1\theta (x)]\mathrm{exp}(x^2/\mathrm{\Delta }x^2)`$, Eqs.(8) and (13) show that the transmitted field will vanish outside the causal boundary $`x_B=c\mathrm{T}+\mathrm{x}_\mathrm{I}`$, or for $`\tau <0`$, as shown in Fig.1c. Note that in Fig.1c the narrow spikes near $`\tau =\tau _m`$ result from large oscillations of $`\eta _m(\tau )`$ clearly visible in Fig1.a. Equally unsuccessful would be an attempt at superluminal transfer of information encoded in a sharp cut-off at the rear of the incident pulse, $`G_+(x)=\theta (x)\mathrm{exp}(x^2/\mathrm{\Delta }x^2)`$. Figure 1d shows that by inspecting the field at $`x>c\mathrm{T}`$, an observer cannot decide whether the whole Gaussian, or only its front half was incident on the barrier, and must await the arrival of the information-carrying part of the signal.
Recent experiments, in which a detector was to distinguish between a cut ($`G_+`$) and an uncut ($`G`$) signals,  have shown that the information detection time for pulses propagating through the fast-light medium is somewhat longer than that in vacuum, even though the group velocity in the medium is in the highly superluminal regime. As in Fig.1d, the absence of superluminal information transfer has a simple explanation. Since the medium (in this case, the potassium vapor) does not bind photons, the transmission amplitude cannot have poles in the upper half of the $`k`$-plane and $`\xi (x^{})`$ in Eq.(5) must vanish for $`x^{}>0`$. The most advanced term in Eq.(5), $`\mathrm{\Psi }_T(x)`$ contains a cutoff at $`x=x_BcT+x_I`$, and the ’superluminal’ ($`x>x_B`$) part of the transmitted field in Fig.1d builds up from the front tails of $`\mathrm{\Psi }_T(xx^{})`$, $`x^{}0`$, unaffected by the cut-off, and for this reason is the same as the advanced uncut field (cf. Fig.2a of Ref.). For $`x<x_B`$ the field has a complicated shape and the actual delay in detecting the back face of the pulse is likely to be caused by its deformation while propagating through the medium. To advance the cut-off beyond $`x_B`$, and achieve a truly superluminal information transfer, one requires a one-dimensional system capable of supporting both bound and scattering states of a photon, in the way a finite-depth potential well supports bound states of an electron. At present, we are not aware of the existence of such systems,
In summary, the notion of ’superluminality’ in wavepacket propagation is based on relating the final position $`x`$ of the transmitted particle to the time $`\tau `$ it is supposed to have spent in the scatterer. Quantally, $`\tau `$ and $`x`$ are related as the measured quantity and the pointer position in a measurement, whose accuracy is determined by the coordinate width of the pulse. In cases where apparent ’superluminality’ is observed, e.g., in tunnelling or optical propagation through ’fast-light’ media, the measurement is inevitably weak. Even if no negative time delays contribute to the transmission ($`\eta (\tau )0`$ for $`\tau <0`$), it may, therefore, produce an anomalous reading by constructing a ’superluminal’ pulse from the front tails of the components of the transmitted pulse, all delayed relative to free propagation. Just as an improvement in the accuracy destroys anomalous weak values , a choice of a narrow incident pulse destroys superluminal propagation by making higher incident momenta pass over the barrier, or spill outside the anomalous dispersion region. As a result, a ’strong’ measurement registers only the ’subluminal’ time delays, as illustrated in Fig.1e for the simple model described above. The contradiction between the impossibility of faster-than-light travel and observing an apparently ’superluminal’ pulse is, therefore, resolved in a typically quantum mechanical fashion: when ’superluminality’ is present, one does not know the delay, and cannot claim that the duration spent in the scatterer is shorter than $`b/c`$. Conversely, when the delay is known, no ’superluminal’ transmission is observed. The absence of negative virtual delays in optical propagation does, however, limit the speed of information transfer to $`c`$ and below, as a non-analytic feature (e.g., a cut-off) of the initial pulse may only travel as far as the most advanced component in Eq.(5), i.e., at most by $`c\mathrm{T}`$. Beyond this point the field will either vanish if the front part of the pulse was discarded, or remain identical to the uncut field if its the rear tail has been removed as the recent experiments by Stenner et al  show.

Acknowledgements: Two of us (DS and VS) gratefully acknowledge the support of the CTSPS, Clark Atlanta University. AZM was supported by the U.S. DoE, Division of Chemical Sciences, Office of Basic Research.

Figure captions

Fig.1:
a) Real (solid) and imaginary (dashed) part of the time delay amplitude distribution, $`\eta (\tau )`$, for $`\mathrm{\Omega }b=100`$;
b) Transmitted field $`\mathrm{\Psi }_F(x)`$ (thick solid) and its components $`\mathrm{\Psi }_m`$, corresponding to different terms in Eq.(14), for $`\mathrm{\Omega }b=100`$, and a Gaussian incident pulse $`k_0b=1.4\pi `$ and $`\mathrm{\Delta }x/b=2.85`$;
c) Same as b)but for an incident Gaussian pulse truncated at front.
d) Same as b)but for an incident Gaussian pulse truncated at rear.
e) Same as b) but for $`k_0b=7\pi `$ and $`\mathrm{\Delta }x=0.57`$.


text:
# Effects of color superconductivity on the structure and formation of compact stars

## Abstract

We show that if color superconducting quark matter forms in hybrid or quark stars it is possible to satisfy most of recent observational boundaries on masses and radii of compact stellar objects. An energy of the order of $`10^{53}`$ erg is released in the conversion from a (metastable) hadronic star into a (stable) hybrid or quark star in presence of a color superconducting phase. If the conversion occurs immediately after the deleptonization of the proto-neutron star, the released energy can help Supernovae to explode. If the conversion is delayed the energy released can power a Gamma Ray Burst. A delay between the Supernova and the subsequent Gamma Ray Burst is possible, in agreement with the delay proposed in recent analysis of astrophysical data.

The new accumulating data from X-ray satellites provide important information on the structure and formation of compact stellar objects. Concerning the structure, these data, are at first sight difficult to interpret in a unique and self-consistent theoretical scenario, since some of the observations are indicating rather small radii and other observations are indicating large values for the mass of the star. Concerning the formation scenario, crucial information are provided by the very recent observations of Gamma-Ray Bursts (GRB), indicating the possibility that some of the GRBs are associated with a Supernova (SN) explosion. It has not yet been clarified if the two explosions are always simultaneous or if, at least in a few cases, a time delay can exist, with the SN preceding the GRB Amati00 ; Reeves02 ; hjorth ; rutledge ; reeves2003 .

The effect of the transition to deconfined Quark Matter (QM) on explosive processes like SNs and GRBs has been discussed by many authors. In particular, the possibility that deconfinement takes place during the core-collapse of massive stars at the moment of the bounce, has been discussed e.g. in Refs.noiJP ; sannino  and this mechanism could help the SN to explode by increasing the mechanical energy associated with the bounce. However, it seems more plausible that deconfinement takes place only when the proto-neutron star (PNS) has deleptonized and cooled down to a temperature of a few MeV Benvenuto ; Pons . The energy released in the conversion to QM produces a refreshed neutrino flux which can help the supernova to explode in a neutrino-driven scheme. Finally, another scenario is possible in which neutron stars having a small enough mass can exist as metastable Hadronic Star (HS) if a non-vanishing surface tension is present at the interface between Hadronic Matter (HM) and QM. The process of quark deconfinement can then be a powerful source for GRBs and it can also explain the possible delay between a SN explosion and the subsequent GRB noiapj .

In recent years, many theoretical works have investigated the possible formation of a diquark condensate in QM, at densities reachable in the core of a compact star alf1 ; alf5 ; schafer . The formation of this condensate can deeply modify the structure of the star alf4 ; baldo1 ; blas1 ; lugo1 . We present here an extension of the previous works, showing that it is possible to satisfy the existing boundaries on mass and radius of a compact stellar object if a diquark condensate forms in a Hybrid Star (HyS) or a Quark Star (QS). Moreover, we show that the formation of diquark condensate can significantly increase the energy released in the conversion from a purely HS into a more stable star containing deconfined QM.

To describe the high density Equation of State (EOS) of matter we adopt standard models in the various density ranges. Concerning the hadronic phase we use relativistic non-linear models glen2 ; ms96 . At very low density we use the standard EOSs of Refs.negele ; baym . For the QM phase we adopt a MIT-bag like model in which the formation of a diquark condensate is taken into account. To connect the two phases of our EOS, we impose Gibbs equilibrium conditions.

It is widely accepted that the Color-Flavor Locking phase (CFL) is the real ground state of QCD at asymptoticly large densities. We are interested in the bulk properties of a compact star and we adopt the simple scheme proposed in Refs.alf4 ; lugo1  where the thermodynamic potential is given by the sum of two contributions. The first term corresponds to a “fictional” state of unpaired QM in which all quarks have a common Fermi momentum chosen to minimize the thermodynamic potential. The other term is the binding energy $`\mathrm{\Delta }`$ of the diquark condensate expanded up to order $`(\mathrm{\Delta }/\mu )^2`$. In Ref.alf4  the gap is assumed to be independent on the chemical potential $`\mu `$. In the present calculation we consider a $`\mu `$ dependent gap resulting from the solution of the gap equation alf1 . The resulting QM EOS reads:

$`\mathrm{\Omega }_{CFL}`$ $`=`$ $`{\displaystyle \frac{6}{\pi ^2}}{\displaystyle _0^\nu }k^2(k\mu )dk`$ (1)
$`+`$ $`{\displaystyle \frac{3}{\pi ^2}}{\displaystyle _0^\nu }k^2(\sqrt{k^2+m_s^2}\mu )dk{\displaystyle \frac{3\mathrm{\Delta }^2\mu ^2}{\pi ^2}}`$
with $`\nu =2\mu \sqrt{\mu ^2+\frac{m_s^2}{3}},`$ and the quark density $`\rho `$ is calculated numerically by deriving the thermodynamic potential respect to $`\mu `$. Pressure and energy density read:

$`P`$ $`=`$ $`\mathrm{\Omega }_{CFL}(\mu )B\mathrm{\Omega }^e(\mu _e)`$ (2)
$`E/V`$ $`=`$ $`\mathrm{\Omega }_{CFL}(\mu )+\mu \rho +B+\mathrm{\Omega }^e(\mu _e)+\mu _e\rho _e.`$ (3)
In Fig.1 we have collected most of the analysis of data from X-ray satellites, concerning masses and radii of compact stellar objects sanwal ; cottam ; quaint ; heinke ; dey ; li ; burwitz . Observing Fig.1, we notice that the constraints coming from a few data sets (labeled “e”, “f”111A very recent reanalysis of the data of the pulsar SAX J1808.4-3658, discussed in Ref. li , seems to indicate slightly larger radii, of the order of 9-10 km for a star having a mass of 1.4-1.5 $`M_{}`$ Poutanen . “g” and maybe also constraint “h”222In Ref.zane  an indication for an even more compact stellar object can be found. Anyway, the so-called thermal radius obtained in these analysis could be significantly smaller than the total radius of the star) indicate rather unambiguously the existence of very compact stellar objects, having a radius smaller than $`10`$ km. At the contrary, at least in one case (“a” in the figure), the analysis of the data suggests the existence of stellar objects having radii of the order of 12 km or larger, if their mass is of the order of 1.4 $`M_{}`$. We recall that it is difficult from an astrophysical viewpoint to generate compact stellar objects having a mass smaller than 1 $`M_{}`$. Therefore the most likely interpretation of constraint “a” is that the corresponding stellar object does not belong to the same class of objects which have a radius smaller than $`10`$ km. Concerning constraint “b”, it can be satisfied both with a very compact star or with a star having a larger radius. The apparent contradiction between the constraints “e”, “f”, “g” and the constraint “a” can be easily accommodated in our scheme, since it can be the signal of the existence of metastable purely HS which can collapse into a stable configuration when deconfined QM forms inside the star.

Finally, constraints (“c”333The result of Ref.quaint  is $`M/M_{}=1.88\pm 0.13`$. In Fig.1 only the lower limit is displayed. and “d”444If the observed X-ray emission is due to continuing accretion, a smaller mass is allowed, $`M/M_{}=1.4`$. ) do not provide stringent limits on the radius of the star, but they put strong constraints on the lower value of its mass. It is in general not easy to obtain stellar configurations having both large masses and very small radii. As we will see, the existence of an energy gap associated with the diquark condensate helps in circumventing this difficulty, since the effect of the gap is to increase the maximum mass of stars having a huge content of pure QM.

In Fig.1 we show a few theoretical M-R relations which correspond to the scenario we are proposing. More precisely, we show a thick-dashed line corresponding to HSs (GM1), a thick dot-dashed line corresponding to HySs (GM1, $`B^{1/4}=170`$ MeV, $`\mathrm{\Delta }_2`$) and a thick solid line corresponding to QSs ($`B^{1/4}=170`$ MeV, $`\mathrm{\Delta }_4`$). Similar shapes can be obtained using the EOS of Ref.ms96 . Both the HyS and the QS lines can satisfy essentially all the constraints derived from observations. The shapes of the gaps $`\mathrm{\Delta }_i`$ are shown in Fig. 2. In conclusion, in our scheme most of the compact stars are either HySs or QSs having a mass in the range $`1.21.8M_{}`$ and a radius $`8.510`$ km. Stars having a significantly larger radius (as the one suggested by constraint “a”) correspond in our scheme to metastable HSs which can exist if their mass is not too large, as we show in the following.

Let us now discuss $`\mathrm{\Delta }E`$, the energy released in the conversion from HS to HyS or QS. $`\mathrm{\Delta }E`$ is the difference between the gravitational mass of the HS and that of the final HyS or QS having the same baryonic mass. As mentioned in the introduction, a possibility is that deconfinement takes place a few seconds after the bounce, when the PNS has deleptonized and its temperature has dropped down Benvenuto ; Pons . In particular, for stars having a small mass the formation of QM takes place only at $`T`$ few MeV. Notice that for a star having a mass of order 1.4$`M_{}`$ and using the relativistic EOSs discussed in this paper, hyperons are present in the initial configuration, since the typical mass at which hyperons start forming is $`1M_{}`$. The energy released during quark deconfinement powers a new neutrino flux which can be useful to make the supernova explode. $`\mathrm{\Delta }E`$ is shown in Table I, for a PNS having a mass of 1.4 $`M_{}`$. As it can be seen, $`\mathrm{\Delta }E`$ can be as large as 1053 erg, if the final configuration corresponds to an HyS and three times as large if a QS is obtained. The effect of the gap is to increase the energy released and to allow QS configurations in cases where an HyS would be obtained in the absence of quark pairing. Let us now remark that the deconfinement transition can be delayed if a non-vanishing surface tension at the interface between HM and QM exists and if the mass of the HS is not too large. This possibility was not discussed in Ref.Pons  and it is the main ingredient of our model. To compute the time needed to form QM we use the technique of quantum tunneling nucleation. We can assume that the temperature has no effect in our scheme because, as discussed above, when QM forms the temperature is so low that only quantum tunneling is a practicable mechanism.

In Ref.noiapj  it was proposed that the central density of a pure HS (containing hyperons) can increase, due to spin down or mass accretion, until its value approaches the deconfinement critical density. At this point a spherical virtual drop of QM can form. The potential energy for fluctuations of the drop radius $`R`$ has the form lif :

$$U(R)=\frac{4}{3}\pi R^3n_q(\mu _q\mu _h)+4\pi \sigma R^2+8\pi \gamma R$$
(4)
where $`n_q`$ is the quark baryon density, $`\mu _h`$ and $`\mu _q`$ are the hadronic and quark chemical potentials, all computed at a fixed pressure $`P`$, and $`\sigma `$ is the surface tension for the interface separating quarks from hadrons. Finally, the term containing $`\gamma `$ is the so called curvature energy. For $`\sigma `$ we use standard values from 10 to 40 MeV/fm2 and we assume that it takes into account, in a effective way, also the curvature energy. The value of $`\sigma `$ was estimated in Ref.jaffe  to be $`10`$ MeV/fm2. Values for $`\sigma `$ larger than $`30`$ MeV/fm2 are probably not useful at the light of the results of Refs.voskresensky ; alf2 .

The calculation proceeds in the usual way: after the computation (in WKB approximation) of the ground state energy $`E_0`$ and of the oscillation frequency $`\nu _0`$ of the virtual QM drop in the potential well $`U(R)`$, it is possible to calculate in a relativistic frame the probability of tunneling as iida

$$p_0=\mathrm{exp}[2\frac{A(E_0)}{\mathrm{}}]$$
(5)
where

$$A(E)=_R_{}^{R_+}𝑑R\sqrt{[2M(R)+EU(R)][U(R)E]}$$
(6)
Here $`M(R)=4\pi \rho _h\left(1\frac{n_q}{n_h}\right)^2R^3`$, $`\rho _h`$ is the hadronic energy density and $`n_h`$, $`n_q`$ are the baryonic densities at a same and given pressure in the hadronic and quark phase, respectively. Finally, $`R_\pm `$ are the classical turning points. The nucleation time is then equal to $`\tau =(\nu _0p_0N_c)^1`$, where $`N_c`$ is the number of centers of droplet formation in the star, and it is of the order of $`10^{48}`$ iida . $`\tau `$ can be extremely long if the mass of the metastable star is small enough but, via mass accretion, it can be reduced from values of the order of the age of the universe down to a value of the order of days or years. We can therefore determine the critical mass $`M_{cr}`$ of the metastable HS for which the nucleation time corresponds to a fixed small value (1 year in Tab. 1).

In Table II we show the value of $`M_{cr}`$ for various sets of model parameters. In the conversion process from a metastable HS into an HyS or a QS a huge amount of energy $`\mathrm{\Delta }E`$ is released. We see in Table II that the formation of a CFL phase allows to obtain values for $`\mathrm{\Delta }E`$ which are one order of magnitude larger than the corresponding $`\mathrm{\Delta }E`$ of the unpaired QM case ($`\mathrm{\Delta }=0`$). Moreover, we can observe that $`\mathrm{\Delta }E`$ depends both on magnitude and position of the gap.

In the model we are presenting, the GRB is due to the cooling of the justly formed HyS or QS via neutrino - antineutrino emission. The subsequent neutrino-antineutrino annihilation generates the GRB. In our scenario the duration of the prompt emission of the GRB is therefore regulated by two mechanisms: 1) the time needed for the conversion of the HS into a HyS or QS, once a critical-size droplet is formed and 2) the cooling time of the justly formed HyS or QS. Concerning the time needed for the conversion into QM of at least a fraction of the star, the seminal work by Oli87  has been reconsidered by HB88 , where it has been shown that the stellar conversion is a very fast process, having a duration much shorter than 1s. On the other hand, the neutrino trapping time, which provides the cooling time of a compact object, is of the order of a few ten seconds ignazio , and it gives the typical duration of the GRB in our model.

In conclusion, comparing the theoretical mass-radius curves with recent observational data, we find that color superconductivity is a crucial ingredient in order to satisfy all the constraints coming from observations. The difficult problem posed by astrophysical data indicating the existence of stars which are both very compact and rather massive can be solved either with hybrid or quark stars. Concerning hybrid stars, the gap increases the maximum mass of the stable configuration, while keeping the corresponding radius $``$ 10 km.

The superconducting gap affects also deeply the energy released in the conversion from hadronic star into hybrid or quark star. We assume that the deconfinement transition only takes place when the star has deleptonized and cooled down, in agreement with the results of Ref.Benvenuto ; Pons . If deconfinement occurs immediately after deleptonization, the energy released can help the SN to explode. If, at variance, the transition is delayed, a metastable hadronic star can form. Its subsequent transition to a stable configuration, containing deconfined quark matter, can power a GRB via the annihilation of neutrinos and antineutrinos emitted during the cooling of the newly formed compact star. The energy released is significantly increased by the effect of the chemical-potential dependent superconducting gap and it can reach a value of the order of $`10^{53}`$ erg. The proposed mechanism could explain recent observations indicating a possible delay between a Supernova and the subsequent Gamma Ray Burst Amati00 ; Reeves02 ; hjorth .


text:
# References

Stratification of Moving Multicomponent Bose-Einstein Condensates

V.I. Yukalov1 and E.P. Yukalova2

1Bogolubov Laboratory of Theoretical Physics
Joint Institute for Nuclear Research, Dubna 141980, Russia

2Department of Computational Physics
Laboratory of Information Technologies
Joint Institute for Nuclear Research, Dubna 141980, Russia

## Abstract

A mixture of the multicomponent Bose-Einstein condensate is considered, where each component moves with its own velocity. As a result of the relative motion, the mixture stratifies when the relative velocity reaches a critical value. Stability conditions for a binary moving mixture are derived and the critical velocity is found.

Ultracold atomic gases, forming Bose-Einstein condensates, are nowadays an active field of research (see reviews \[1–5\]). Dilute Bose gas can be trapped and cooled down to temperatures when practically all atoms are in Bose-Einstein condensate. One can create mixtures of several Bose-condensed species. The stability of such mixtures depends on the mutual atomic interactions. The condition of stability for binary mixtures was, first, found by Nepomnyashchy , considering the collective excitation spectrum in the Bogolubov and long-wave approximations. When this stability condition is not valid, the mixture is not stable and stratifies onto separate components. The analysis of Nepomnyashchy  was based on the Green functions approach for a superfluid mixture with broken gauge symmetry. As was shown later , the same stability condition holds true for a fluid mixture with conserved gauge symmetry.

Recently, it has become feasible to create moving wave packets of Bose-Einstein condensates, which could collide and move through each other \[8–12\]. Therefore, it is interesting to study the interplay between the possibility of different Bose components to mix and the influence on their stability of their relative motion through each other. Such a relative motion of two fluids could change the stability condition of mixing . Bose-Einstein condensates, formed in modern traps, are made of dilute gases, for which the Gross-Pitaevskii equation provides a very good description \[1–5\]. Actually, this equation is an exact equation for the coherent system at zero temperature . The aim of the present communication is to analyze a mixture of Bose-Einstein condensates, whose components move through each other, under conditions of low temperature, when the system is described by the Gross-Pitaevskii equation. This situation corresponds to current experiments with Bose-Einstein condensates of dilute gases.

An atomic gas, of $`N`$ atoms in volume $`V`$, is dilute when its density $`\rho N/V`$ and the $`s`$-wave scattering length $`a`$ are connected by the inequality $`\rho a^31`$. This can be generalized to a mixture of species enumerated by the index $`i=1,2,\mathrm{}`$ as the inequality $`\sqrt{\rho _i\rho _j}a_{ij}^31`$, in which $`a_{ii}a_i`$ and $`a_{ij}`$ are the corresponding scattering lengths, and $`\rho _iN_i/V`$. The Gross-Pitaevskii equation for a multicomponent mixture is

$$i\mathrm{}\frac{}{t}\eta _j(𝐫,t)=\widehat{H}_j[\eta ]\eta _j(𝐫,t),$$
(1)
where $`\eta _i(𝐫,t)`$ is the coherent field of the $`i`$-component; in the nonlinear Hamiltonian

$$\widehat{H}_i[\eta ]=\frac{\widehat{𝐩}^2}{2m_i}+U_i+\underset{j}{}A_{ij}|\eta _j|^2$$
(2)
the first term is the kinetic-energy operator, with $`\widehat{𝐩}i\mathrm{}\stackrel{}{}`$, and $`m_i`$ being atomic mass, $`U_i`$ is an external (trapping) potential, and the third term describes contact atomic interactions with

$$A_{ij}2\pi \mathrm{}^2a_{ij}\frac{m_i+m_j}{m_im_j}.$$

The normalization condition for each component reads

$$\eta _i^2|\eta _i(𝐫,t)|^2𝑑𝐫=N_i.$$

Elementary collective excitations of the system can be obtained by linearizing the Gross-Pitaevskii equation (1) around a given stationary solution. Following here this way, we should keep in mind that our aim here is to study a mixture of moving components. Assume that each component moves with a fixed velocity $`𝐬_i`$. Then taking account of this motion can be done by means of the Galilean transformation replacing $`\eta _j`$ by $`\eta _j\mathrm{exp}\{\frac{i}{\mathrm{}}m_j𝐬_j𝐫\}`$. To find collective excitations in the moving mixture, we write the coherent field $`\eta _j`$ in the form

$$\eta _j(𝐫,t)=\left[\phi _j(𝐫)+u_j(𝐫)e^{i\omega t}+v_j^{}(𝐫)e^{i\omega t}\right]\mathrm{exp}\left\{\frac{i}{\mathrm{}}\left(m_j𝐬_j𝐫E_jt\right)\right\},$$
(3)
where $`\phi _j(𝐫)`$ satisfies the stationary equation

$$\left(\widehat{H}_i[\phi ]+\widehat{𝐩}𝐬_i+\frac{m_is_i^2}{2}\right)\phi _i(𝐫)=E_i\phi _i(𝐫),$$
(4)
with the normalization $`\phi _i^2=N_i`$. Substituting Eq. (3) into the Gross-Pitaevskii equation (1) and linearizing it with respect to the functions $`u_j(𝐫)`$ and $`v_j(𝐫)`$, we obtain

$$\left(\mathrm{}\omega \widehat{𝐩}𝐬_i\widehat{H}[\phi ]\frac{m_is_i^2}{2}+E_i\right)u_i\phi _i\underset{j}{}A_{ij}\left(\phi _j^{}u_j+\phi _jv_j\right)=0,$$

$$\left(\mathrm{}\omega \widehat{𝐩}𝐬_i+\widehat{H}[\phi ]+\frac{m_is_i^2}{2}E_i\right)v_i+\phi _i^{}\underset{j}{}A_{ij}\left(\phi _j^{}u_j+\phi _jv_j\right)=0.$$
(5)
This is a generalization of the Bogolubov - de-Gennes equations to the case of a mixture of moving components.

Being interested in collective excitations above the ground state of trapped atoms, we keep in mind that the wave functions, defined by Eq. (4), are real-valued. Then we can resort to the local-density approximation, setting $`\phi _i(𝐫)\sqrt{\rho _i(𝐫)}`$ and assuming that $`u_i`$ and $`v_i`$ are proportional to $`e^{i𝐤𝐫}`$. Supposing that atomic interactions are sufficiently strong, so that mean potential energies are essentially larger that mean kinetic energies, we may omit the terms containing the spatial differentiation of $`\rho _i(𝐫)`$. In this approximation,

$$\left(\widehat{H}_i[\phi ]+\frac{m_is_i^2}{2}\right)u_i=\left(E_i+\frac{\mathrm{}^2k^2}{2m_i}\right)u_i.$$
(6)
Then Eqs. (5) reduce to

$$\mathrm{}\left(\omega 𝐤𝐬_i\frac{\mathrm{}k^2}{2m_i}\right)u_i\underset{j}{}A_{ij}\sqrt{\rho _i\rho _j}(u_j+v_j)=0,$$

$$\mathrm{}\left(\omega 𝐤𝐬_i+\frac{\mathrm{}k^2}{2m_i}\right)v_i+\underset{j}{}A_{ij}\sqrt{\rho _i\rho _j}(u_j+v_j)=0.$$
(7)
The latter equations possess a nontrivial solution if the related determinant is zero. To write down the resulting equation in a compact form, we introduce the notation

$$\omega _i^2c_i^2k^2+\left(\frac{\mathrm{}k^2}{2m_i}\right)^2,c_i^2\frac{\rho _i}{m_i}A_i,$$
(8)
where $`A_iA_{ii}`$. If there would be no interactions between different species, then each component would have the Bogolubov spectrum $`\omega _i(k)`$ given by Eq. (8). Note that the existence of the Bogolubov spectrum for a single-component Bose-Einstein condensate of trapped atoms is well confirmed by experiment . Define also

$$\omega _{ij}c_{ij}k,c_{ij}^2\left(\frac{\rho _i\rho _j}{m_im_j}\right)^{1/2}A_{ij}.$$
(9)
And let us consider the case of a binary mixture. The spectrum of collective excitations, following from Eq. (7), is given by the equation

$$\left[(\omega \epsilon _1)^2\omega _1^2\right]\left[(\omega \epsilon _2)^2\omega _2^2\right]=\omega _{12}^4,$$
(10)
in which $`\epsilon _i𝐤𝐬_i`$.

It is convenient to connect the system of coordinates with one of the components, say with the first one, so that $`𝐬_1=0`$. Then $`𝐬_2𝐬`$ is the velocity of the second component with respect to the first component. In that case, the spectral equation (10) yields

$$\omega ^42\epsilon \omega ^3\left(\omega _1^2+\omega _2^2\epsilon ^2\right)\omega ^2+2\epsilon \omega _1^2\omega +\omega _1^2\left(\omega _2^2\epsilon ^2\right)\omega _{12}^4=0,$$
(11)
where

$$\epsilon 𝐤𝐬=ks\mathrm{cos}\vartheta .$$
(12)
Equation (11) has four solutions for $`\omega =\omega (𝐤)`$. But only two of them can be stable, being positively defined for all $`\vartheta [0,2\pi ]`$. The necessary condition for Eq. (11) to possess two positive solutions for arbitrary $`\epsilon [ks,ks]`$, by Descartes theorem, is

$$\omega _1^2\omega _2^2\omega _{12}^4>\omega _1^2k^2s^2.$$
(13)
In the long-wave limit, when $`k0`$, we have the stability condition

$$c_1^2c_2^2c_{12}^4>c_1^2s^2.$$
(14)
This shows that, when the immovable mixture is stable, satisfying the inequality $`c_1^2c_2^2>c_{12}^4`$, it will retain its stability only up to the critical velocity

$$s_c=\frac{1}{c_1}\sqrt{c_1^2c_2^2c_{12}^4}.$$
(15)
Reaching this relative velocity of one component through another, the binary mixture becomes unstable and stratifies onto two spatially separated components. The stratification starts inside the cone of the angle

$$\vartheta _c=\mathrm{arccos}\frac{s_c}{s}.$$
(16)
The stability condition (14), with Eqs. (8) and (9), transforms to

$$A_1A_2A_{12}^2>\frac{m_2}{\rho _2}A_1s^2.$$
(17)
Respectively, the critical velocity (15) can be written as

$$s_c=\sqrt{\frac{\rho _2(A_1A_2A_{12}^2)}{m_2A_1}}.$$
(18)
Taking into account the relation between $`A_{ij}`$ and the scattering length $`a_{ij}`$, the stability condition (17) can be presented in the form

$$\frac{a_1a_2}{a_{12}^2}>\frac{(m_1+m_2)^2}{4m_1m_2}+\frac{m_2^2a_1s^2}{4\pi \mathrm{}^2\rho _2a_{12}^2}.$$
(19)
Finally, for colliding wave packets with nonuniform densities $`\rho _i(𝐫)`$, the stratification would occur not everywhere in space, where the moving packets overlap, but only there, where the inequality

$$\rho _2(𝐫)>\frac{m_2A_1s^2}{A_1A_2A_{12}^2}$$
(20)
becomes invalid, provided that $`A_1A_2>A_{12}^2`$.

The stability conditions for the mixture of components with relative motion have been obtained here by considering the stability of the spectrum of collective excitations. This is what is called the dynamic stability . There exists as well the notion of thermodynamic stability, when one compares the corresponding free energies, or, at zero temperature, the related internal energies. For immovable mixtures, these two types of stability result in equivalent stability conditions. However, for the mixtures of moving components, these two kinds of stability are not equivalent. To show the latter, we may, for simplicity, employ the uniform approximation, writing the internal energy of the mixture as

$$E_{mix}=\frac{1}{2V}\left(N_1^2A_1+N_2^2A_2+2N_1N_2A_{12}\right)+K,$$
(21)
where $`V=V_1+V_2`$ and the last term is the kinetic energy

$$K=\frac{1}{2}\left(N_1m_1s_1^2+N_2m_2s_2^2\right).$$
(22)
If the components are spatially separated, that is, stratified, their total internal energy is

$$E_{str}=\frac{N_1^2A_1}{2V_1}+\frac{N_2^2A_2}{2V_2}+K.$$
(23)
Because of the additivity of kinetic energy, the comparison of $`E_{mix}`$ and $`E_{str}`$ does not involve the velocities of the moving components, so that the thermodynamic stability condition $`E_{mix}<E_{str}`$ does not depend on the relative motion of these components. Moreover, employing the thermodynamic stability condition, one commonly invokes the condition of mechanical equilibrium, according to which the pressures of the stratified components,

$$p_i\frac{E_{str}}{V_i}=\frac{N_i^2A_i}{2V_i},$$

are equal, so that $`p_1=p_2`$. The latter yields the equality

$$\frac{N_1^2}{V_1^2}A_1=\frac{N_2^2}{V_2^2}A_2,$$

from where, recalling that $`V=V_1+V_2`$, we have the relations

$$\frac{V}{V_1}=1+\sqrt{\frac{N_2^2A_2}{N_1^2A_1}},\frac{V}{V_2}=1+\sqrt{\frac{N_1^2A_1}{N_2^2A_2}}.$$

This makes it possible to present Eq. (23) as

$$E_{str}=\frac{1}{2V}\left(N_1^2A_1+N_2^2A_2+2N_1N_2\sqrt{A_1A_2}\right)+K.$$
(24)
Now, comparing Eqs. (21) and (24), one notices that $`E_{mix}<E_{str}`$ if and only if $`\sqrt{A_1A_2}>A_{12}`$. However, the usage of the condition of mechanical equilibrium for moving separated components does not look appropriate. The thermodynamic stability condition can be applied only to immovable mixtures.

The method of studying the stability of mixtures with a relative motion of their components by means of the dynamic stability condition, requiring the positive definiteness for the spectrum of collective excitations, is quite general and can be applied to arbitrary mixtures. In particular, this can also be an atomic-molecular mixture. Such mixtures are formed by means of the Feshbach resonance technique . Two Bose atoms can form a quasi-molecule. Then, it is possible to create a two-component mixture consisting of atoms themselves and of the molecules formed by these atoms. When both atoms and molecules are in the Bose-condensed state, they again can be described by the system of the Gross-Pitaevskii equations for the corresponding coherent fields. Thus, if the interaction amplitude for the reaction ”atom+atom $``$ molecule” is denoted by $`B`$, then the related Gross-Pitaevskii equations are

$$i\mathrm{}\frac{}{t}\eta _1=\widehat{H}_1[\eta ]\eta _1+B\eta _1^{}\eta _2,i\mathrm{}\frac{}{t}\eta _2=\widehat{H}_2[\eta ]\eta _2+\frac{1}{2}B\eta _1^2,$$
(25)
where $`\widehat{H}_i[\eta ]`$ are the same Hamiltonians (2). When atoms as well as molecules are immovable, their mixture is stable under the condition

$$\sqrt{A_1A_2}>A_{12}+\frac{B}{\sqrt{\rho _2}},$$

where $`\rho _2N_2/V`$. But if one is able to organize the relative motion of the atomic and molecular components with respect to each other, then, to find the stability condition for such a moving mixture, one has to follow the same procedure as described above. That is, the form (3) of the coherent fields has to be substituted into Eqs. (25), after which the spectrum of collective excitations can be obtained. Requiring the positiveness of this spectrum yields the stability condition under the relative motion of components.

Concluding, if an immovable mixture of Bose-condensed species is stable, it can stratify if the components move with respect to each other. The stratification occurs when the relative velocity of motion reaches a critical value. The conditions of stability and, respectively, stratification can be derived by linearizing the Gross-Pitaevskii equations for the components, with the inclusion of the Galilean transformation for the coherent fields, as is done in Eq. (3).

Corresponding author

Prof. V.I. Yukalov

Bogolubov Laboratory of Theoretical Physics

Joint Institute for Nuclear Research

Dubna 141980, Russia

E-mail: yukalov@thsun1.jinr.ru

Tel: 7 (096) 216-3947 (office); 7 (096) 213-3824 (home)

Fax: 7 (096) 216-5084