text:
# 1 Introduction

## 1 Introduction

For many years since the first work of Colgate & White (1966) , numerical simulations of core-collapse supernova explosions have been exciting topics. Here a “successful” explosion must explain the following observed facts:

* The explosion energy (observed in SN 1987A): 1.5 $`\pm \mathrm{\hspace{0.17em}0.5}\times 10^{51}`$ ergs.
* Asymmetry observed by polarization measurements in SN 1987A, SN 1993J, and other supernovae.
* Neutron star mass: 1.4 $`M_{}`$ .
* Explosive nucleosynthesis: 56Ni ($``$ 0.07 $`M_{}`$ for SN 1987A).

So far, no “successful” explosions have been reported in 1-D and Multi-Dimensional simulations except work by Wilson and Mayle (e.g. ; see, a review by Sato of this volume for the status of numerical simulation studies).

In the so-called “delayed explosion” scenario, an important role of convection, which is simulated with spherical neutrino radiation field, has been widely accepted. Here we study the supernova explosion energy with paying special attention on anisotropic neutrino emission and locally intense neutrino heating thereby produced. Such a new idea that anisotropic neutrino radiation may affect the explosion mechanism itself was first proposed by Shimizu et al. (1994) , one of the coauthors in the present study, and the effect of anisotropic neutrino radiation on the explosion energy was intensively studied by Shimizu et al. (2001) .

Figure 1 depicts the dimensionless entropy contour map  to explain the overall feature of this idea. First, we assume that a stronger neutrino flux is emitted in the direction of the rotation (polar) axis of a proto-neutron star than that in the equatorial plane. This anisotropic neutrino emission, for example, may be caused by the rotation of the core of a progenitor star (see  for other possibilities). Then high-entropy hot bubble is created along the polar axis behind the primary shock. Convective motion occurs, but the convection becomes very global compared with other models with spherical neutrino radiation. The final important feature of our model is that the shock wave itself is deformed as seen in Figure 1. It was emphasized in Shimizu et al.  that a reasonable degree of anisotropy, lz/lx = $``$1.05 - $``$1.2, is enough to increase the explosion energy by a large factor. In the above, lz and lx are, respectively, the local neutrino flux in the polar axis and that in the equatorial plane. (see  for details of parameterization of neutrino flux distributions). Shimizu et al. (2001) concluded that global anisotropic neutrino radiation is an alternative mechanism to revive the shock wave and lead to a successful explosion other than so far suggested “convective trigger”.

In this article we are going to show that global anisotropy of the neutrino radiation is the most effective mechanism of increasing the explosion energy than the models of spherically symmetric and fluctuated neutrino radiation when the total neutrino luminosity is given. We discuss the reason why, and demonstrate how sensitively the success of a supernova explosion depends on the neutrino luminosity.

Our simulation is performed by solving 2-dimensional hydrodynamic equations in spherical coordinates. A generalized Roe’s method is employed to solve the hydrodynamical equations with general equations of state (EOSs). Since our method is very rigorous to treat neutrinos above the neutrinosphere, we put our initial condition (t = 0) when a spherical steady shock is stalled above the neutrinosphere. The details of our numerical technique, together with the EOS and the initial condition used, are described in the previous article . It is noted here that we have improved the numerical code in the papers  and  to avoid a numerical error near the pole; although the error was not serious for the investigation of the explosion energy, this may affect the results of explosive nucleosynthesis.

## 2 Reasoning of Powerful Explosion with Anisotropic Neutrino Radiation

Figure 2 shows the evolution of the explosion energy, as well as the thermal, kinetic, and gravitational energies for models of lz/lx = 1.1 and Tν = 4.7 MeV. The quantity Tν is the neutrino temperature, and we consider here only electro-type neutrinos and the total neutrino luminosity is fixed (see  for details). In Figure 2, we observe that the difference in the evolution of the explosion energy between the globally anisotropic model and the spherically symmetric model is clearly prominent.

Two dashed lines in the explosion energy (Eexpl) in Figure 2 are for the cases when the neutrino field is fluctuated in space (see  for details). Such neutrino fluctuations may occur due to gravitational oscillation of proto-neutron stars as proposed by Burrows et al. . The model for the lower dashed line of Eexpl in Figure 2 has larger mode number of fluctuations in the neutrino flux. It is therefore seen that the explosion energy decreases as the mode number of fluctuations in the neutrino flux increases, and finally approaches that of spherical explosion.

Let us mention here the reasoning of this powerful explosion for the globally anisotropic neutrino radiation model. Figure 3 shows the accumulation of the absorbed and emitted energies due to neutrino heating and cooling for both globally anisotropic (lz/lx = 1.1) and spherical models.

The balance between heating and cooling determines the evolution of the explosion energy. The crucial point here is the following: The neutrino heating is caused by neutrino absorption due to nuclei and neutrino scattering off electrons and positrons, and hence the heating rate depends on the neutrino temperature Tν: It varies as T$`{}_{}{}^{6}{}_{\nu }{}^{}`$. On the other hand, the neutrino cooling is caused by neutrino emission due to electron captures and the emission of thermal (photo, pair, plasma) neutrinos. The cooling rate is then a very sensitive function of the matter temperature T, roughly proportional to T6.

The reasoning of the powerful explosion for the globally anisotropic neutrino radiation model can be explained as follows. The locally intense neutrino radiation in our model causes local heating along the polar axis. The thermal pressure around the heated matter becomes high, and hence it pushes part of the shock wave outward. Namely, the shock revival occurs in one direction first along the axis of rotation. Then the pressure gradient does work along the shock front, and the shock revival prevails to all directions (Note the Rankine-Hugoniot condition). This leads to the earlier shock revival than the spherical model. Thus the matter temperature behind the shock decreases rapidly due to expansion. Accordingly, the neutrino cooling rate rapidly drops. Heating rate, on the other hand, remains almost unchanged between the the anisotropic and the spherical models. Heating dominates cooling as a result (see Figure 3). Finally, the explosion energy increases because of the suppression of neutrino cooling.

Notice here that the position where the cooling rate suddenly starts to decrease in Figure 3 (i.e., t $``$ 80 ms for the anisotropic, and t $``$ 230 ms for the spherical model) corresponds to the shock revival time, which can be seen in Figure 2 as the position where the explosion energy suddenly starts to rise up.

## 3 Sensitivity of Explosions to Neutrino Luminosities

In Figure 4 we show a model case with a reduced temperature which fails to explode when the neutrino field is spherically symmetric and fluctuated, but does produce an adequate explosion when the field is reasonably anisotropic (lz/lx = 1.1); the difference between the model with global anisotropy and the other models is extremely remarkable.

We now see that two spherically symmetric models which yield very different explosion outcomes in Figures 2 and 4, but these differ by only 1 % in neutrino temperature. This is a clear demonstration of the sensitivity of the supernova problem on the neutrino luminosity and energy.

## 4 Concluding Remarks

It has been discussed that the globally anisotropic (jet-like) explosion is the most effective mechanism to increase the explosion energy than spherical and neutrino-fluctuated explosions when the total neutrino luminosity is given. While the enhancement of total neutrino luminosity is the most important ingredient for driving to an energetic explosion as concluded by Janka & Müller , the total neutrino luminosity cannot be simply increased to explain the observed explosion energy. This is because such treatment often leads to the problem of Ni overproduction, especially in the case of essentially spherical models. We have demonstrated that the enhancement of the total neutrino luminosity is not always necessary: the less luminous the neutrino radiation is, the more important the neutrino anisotropy is (compare Figures 2 and 4). This is a confirmation of the claim by Shimizu et al. , using a different neutrino flux distributions and an updated numerical code.

The degree of anisotropy assumed here, lz/lx = 1.1, can be translated into the rotation period of a proto-neutron star of $``$ 15 ms on the assumption of the Maclaurin spheroid (see ). We note that the following saturation properties of the neutrino anisotropy has been discussed in the papers ; the effect of anisotropic neutrino radiation on the explosion energy appears with a small value (e.g., lz/lx $``$ 1.05), and the effect saturates with a certain value (e.g., lz/lx $``$ 1.2; the actual value depends on the assumed form of the local neutrino flux). We emphasize here again that a large anisotropy in the neutrino radiation is NOT required for the explosion mechanism.

Finally, we also showed that the success of an explosion depends strongly on a subtle difference in the neutrino energy and luminosities (compare the results of Figures 2 and 4; the temperature difference between 4.70 and 4.65 MeV is only 1 %). What is most important is that anisotropic neutrino radiation, or local neutrino heating, has been suggested as an alternative key to the explosion mechanism of supernovae.


format: text

text:
# Heat conduction in one-dimensional Yukawa chains

## Abstract

Heat conduction in one-dimensional Yukawa chains is investigated. It is shown numerically that it has the abnormal heat conduction which is proportional to the system size. Effects of asymmetric external potential, the modified Frenkel-Kontorova one, on the heat conduction of system are also studied. It is found that the asymmetric property of external potential can induce asymmetric thermal conductivity that can be used to be a effective thermal rectifier. In certain of system parameters the heat flux are significantly different for two opposite direction. One can parametrically control the heat flux through this system by changing the potential strength and width.

Recent years many interesting works have been addressed to the problem of heat conduction lep1 ; lep2 ; hu1 ; cas1 ; pos1 ; hu2 ; sav1 ; alo1 ; cas2 . It is well known that heat conductivity of one dimensional (1D) lattice is an important problem which is related to microscopic foundation of the Fourier law. It also provides an abundant example for the microscopic origin of the macroscopic irreversibility phenomenon. Detailed review of the problem is presented in the recent paper lep1 .

Except of extensive numerical studies focusing on the validity of Fourier law $`J=\kappa dT/dx`$ for many models such as the Fermi-Pasta-Ulam (FPU) model lep2 ; hu1 , the ding-a-ling model cas1 ; pos1 , the Frenkel-Kontorova (FK) model hu2 ; sav1  and the Lorenz-gas model alo1  and so on, some other different and important problems have also attracted people attention. For example most recently Casati et al. cas2  investigated the possibility to control the energy transport inside a nonlinear 1D chain connecting with two thermostats at different temperatures. They emphasized that controlling heat conduction by nonlinearity opens a new possibility to design a thermal rectifier, i.e., a lattice that carries heat preferentially in one direction.

In this Letter we investigate heat conduction problem for a new type of model wherein there is a Yukawa-Coulomb interaction potential. It arises from dust-grain plasma which is extensively studied recently due to its new features with dust crystal or collective waves lin1 ; tho1 ; vla1 ; xie1 ; xie2 ; sam1 . Our study shows that the heat conductivity is abnormal in absence of external field as $`\kappa N`$ as in harmonic chains. Under FK field, however, it is normal which is similar to the cases of former studies but with harmonic or/and FPU interactions hu2 . Surprisingly for the dust chains with real wake potential vla1 , modified FK field, the heat conduction exhibits a mixture property of normal and abnormal behaviors. It is found that the asymmetric property of external wake potential can induce asymmetric thermal conductivity that can be used to be a thermal rectifier. In a certain of system parameters, by changing the strength and width of wake potential, the heat fluxes for two opposite direction can be different remarkably.

In a general form the Hamiltonian of the considered system is written as

$$H=\underset{i}{}h_i,h_i=\frac{p_i^2}{2}+V(x_{i1,}x_i)+U(x_i)$$
(1)
where $`V(x_{i1,}x_i)`$ stands for the interaction potential of the nearest-neighbor particles and $`U(x_i)`$ is a external on-site potential. For dust particles system we have the interaction potential (refer to xie2 )

$$V(r_{ij})=\frac{q^2}{r_{ij}}\mathrm{exp}(\frac{r_{ij}}{\lambda _d})\alpha \frac{q^2}{r_{ij}}$$
(2)
where $`\alpha `$ is the parameter that ranges from $`0`$ to $`1`$, $`\lambda _d`$ is the plasma Debye length and $`q`$ is the dust grain charges, while $`r_{ij}`$ is the distance between particles $`i`$-th and $`j`$-th. For simplicity we have assumed that the masses $`m_d`$ and charges of all the grains $`q`$ are the same and constant, while for convenience we shall choose parameter $`\alpha `$ that makes the neighbor dust particle is located at the bottom of $`V`$ (see more clear in the following). Note that the plasma electrons and ions do not appear explicitly in the simulations but their effects, which result in an attractive force to dust, are taken into account in the effective interaction potential as the second term of Eq.(2). On the other hand in dusty plasma we know that there exist wake potential acting the dusts along downstream vla1

$$U(r_{ti})=\frac{qq_t}{r_{ti}}\frac{2}{1M^2}\mathrm{cos}(\frac{r_{ti}}{\lambda _d\sqrt{M^21}})$$
(3)
where $`q_t`$ is the upstream test dust particle charge, $`M>1`$ is the Mach number of ion flow and $`r_{ti}`$ is the distance of $`i`$-th particle from the test particle. The potential (3) is named as the FK-modified one in this paper because it is different from FK potential with varying amplitude.

Let us turn to the side of the motion of system. A general approach is employed here that are used by many authors like Lepri et al. and Hu et al., namely, two Nose-Hoover thermostats are put on the first and last particle, keeping the temperature at $`T_+`$ and $`T_{}`$, respectively. For simplicity and convenience, firstly we need to normalize the basic physical quantities: the mass by the particle mass $`m_d`$, the space distance by the dust inter-particle distance $`a_d`$ and the time by the $`\sqrt{ma_d^3/q^2}`$ respectively. Then the equations of motion for particles are

$$\stackrel{}{x}_1=\xi _+\stackrel{}{x}_1+f_1f_2,\xi _+=\frac{\stackrel{}{x}_1^2}{T_+}1$$

$$\stackrel{}{x}_i=f_if_{i+1},i=2,\mathrm{},N1$$
(4)
$$\stackrel{}{x}_N=\xi _{}\stackrel{}{x}_N+f_Nf_{N1},\xi _{}=\frac{\stackrel{}{x}_N^2}{T_{}}1$$

where $`f_i=V^{^{}}(x_ix_{i1})U^{^{}}(x_i)`$ is the force. Usually the fixed boundary conditions $`x_0\mathrm{const}`$ and $`x_{N+1}\mathrm{const}`$ are assumed. Now the interaction potential in normalized form can be written

$$V(x_{i1,}x_i)=\left(\mathrm{exp}[\chi (x_ix_{i1})]\alpha \right)/(x_ix_{i1})$$
(5)
where parameter $`\chi =a_d/\lambda _d`$ is the ratio of inter-particle distance to the plasma Debye length. Many studies show that when $`\chi `$ is large the nearest neighbor interaction is appropriate instead of particle-pair interaction sam1 . In our following numerical studies $`\chi =5`$ is generally used. The particles are arranged initially at equilibrium positions of $`V`$ and the normalized inter-distance equals to 1, therefore, it means that parameter $`\alpha `$, adjusting the attraction of plasma to dust particle, must satisfy $`\alpha =(1+\chi )\mathrm{exp}(\chi )`$. For the external FK-modified potential Eq.(3) we have its normalized form

$$U(x_i)=\gamma _2\mathrm{cos}(\gamma _1x_i)/x_i,$$
(6)
where $`\gamma _1=\chi /\sqrt{M^21}`$ and $`\gamma _2=2\beta /(1M^2)`$ with $`\beta =q_t/q`$ the ratio of upstream test particle charge to downstream particles charge. Now we denote $`x_i`$ as the position of $`i`$-th dust particle in the chain from the origin of test dust particle $`x_{0t}0`$. In correspondence for the normalized external FK potential one has $`U(x_i)=\gamma _2\mathrm{cos}(\gamma _1x_i)`$. The different temperature profiles for different external potentials of field-free, FK and FK-modified are shown in Fig.1. The imposed temperature are $`T_+=0.12`$ and $`T_{}=0.08`$ and $`\gamma _1=1`$. The potential strengths are $`\gamma _2=0.3`$ and $`2`$ for FK and FK-modified, respectively. Obviously the temperature gradients are constructed globally only in the FK case.

Usually the heat flux at the $`i`$-th position is given by $`J_i(t)=\frac{1}{2}(\stackrel{}{x}_i+\stackrel{}{x}_{i+1})f_{i+1}`$. Numerically the time average $`J=J_i(t)`$ is independent of the index for long time enough. We plot in Fig.2 the heat flux $`JN`$ vs the system size $`N`$ for three different cases of field-free (a), FK (b) and FK-modified (c). Obviously the abnormal heat conductivity $`\kappa N`$ is obtained for purely Yukawa interaction chains, which is similar to the case of harmonic situation. In the Yukawa chains for standard FK external potential, although it seems unclear about whether there exists phase transition of heat conduction between normal and abnormal hu2 ; sav1 , our results indicates that the normal heat conductivity exists, which support the conclusion of normal heat conduction in Ref. hu2 . From Fig.2 (c), however, we can see that it is abnormal globally on heat conduction in the Yukawa chains under FK-modified on-site external potential while the local temperature gradient appeared in the left region of system, see Fig.1(c). The dependence of $`J_L/J_R`$, the heat flux ratio of opposite direction for the Yukawa-Coulomb model, on parameter $`\gamma _2`$, the strength of the FK-modified on-site external potential is plotted in Fig.3, where the power law $`J_L/J_R10^{0.24\gamma _2}`$ is fitted in our numerical parameter regime. Our numerical experiments shows the validity of this power law until to $`\gamma _2=4`$ which indicates that the difference of two directional heat current can possibly over one order. However in practice of dusty plasmas it can only be realized by a relatively small $`\gamma _2`$ due to the test dust charge equals or less than the charge of particles in chains. For example, for given $`\chi =5`$, $`\gamma _1=1`$ and $`\beta =q_t/q1`$ it corresponds to the practice ion flow Mach number $`M=\sqrt{26}`$ (usually in the sheath $`M>1`$), and $`\gamma _22`$ which indicates $`J_L/J_R3`$.

In order to see how the single kicked particle transfer energy to the others, in Fig.4 we plotted the momentum excitations with different external potentials of three cases of field-free, FK and FK-modified. We can see that in the case of field-free there exists a solitary-like wave which can carry energy independently and not transform it to others. In case of FK there exist phonons which can transfer heat to the others. And in the case of FK-modified there seems a combination of two kind waves mentioned above: solitary-like and scatter phonons. Moreover it seems that there are always more solitary-like waves which propagates in the left-direction than that in the right direction. For example, there are three solitons in Fig.3(f) but only one in Fig.3(e). Maybe this is one reason why in this case the left-directional heat flux, $`J_L`$, is larger than that of right directional, $`J_R`$.

On the other hand according to the expression of heat flux, the ratio of them of two-direction is approximately

$$\frac{J_L}{J_R}\sqrt{\frac{T_+}{T_{}}}\frac{r_>}{r_<}\mathrm{exp}\chi (r_>r_<)$$
(7)
where $`r_>`$ and $`r_<`$ are the average inter-particle distance in the region where the temperature is almost constant for $`[T_+,T_{}]`$, right-direction heat-flow $`J_R`$, and $`[T_{},T_+]`$, left-direction heat-flow, respectively. In general it is shown that $`r_>`$ is a little larger than $`r_<`$ by numerical experiments. It is not surprising that the more phonons are scattered in the case of $`[T_+,T_{}]`$ than that of vise versa. Therefore even if $`r_>`$ is just a very little greater than $`r_<`$, $`J_L/J_R`$ would be greatly larger than one. For example, $`r_>1.08`$ and $`r_<0.92`$, then we have $`J_L/J_R3`$. However in the case of Harmonic interaction potential because of its symmetric quadratic property about equilibrium position the ratio of opposite directional heat flux is in order of $`\sqrt{T_+/T_{}}`$, for example by choosing $`T_+=0.12`$, $`T_{}=0.08`$, the ratio of $`J_L/J_R1.23`$. Our numerical results confirm this simple approximation analysis. Indeed if we instead of harmonic interaction chains the heat flux difference for two direction can not be over $`23\%`$ under asymmetric potential.

Finally we calculate heat flux in different potential width $`\gamma _1^1`$ (see Fig.5). It shows that we can also control the energy flow by changing the potential width $`\gamma _1^1=\sqrt{M^21}/\chi `$ through adjusting the ion flow velocity, which is associated with $`M`$ or and the plasma particles density which is associated with $`\chi `$. Obviously the great difference of heat flux at low temperature caused by the width of potential smear out as the temperature increases to over $`0.5`$. It is noted that the ratio of heat flux $`J_R(\gamma _1=1)/J_R(\gamma _1=\pi )`$ is about over two orders difference.

In summary we have studied the heat conduction problem in the system of 1D Yukawa interaction chains with asymmetric on-site modified FK wake potential. We can parametrically control the heat flux through this system by changing the potential strength and width. Moreover it is also possible to adjust the heat flux by changing the potential from FK to modified FK in the Yukawa chains, thus we provide a simple method to change the properties of the system, from a normal conductor obeying Fourier law, down to an almost perfect insulator. Some theoretical consideration has also been attempted to understand for our results and it is found that in the case of Yukawa interaction potential, the asymmetric property about equilibrium position has played a crucial role to improve the ratio of opposite directional heat flux. It seems that there exist some modes which are combinations of solitary-like and scattering phonons. Certainly more formidable theoretical works is still challenging our research on this problem in the future.

Authors would like to thank Drs. L. Wang, L. Yang and H. Zhang for their useful discussions. This work was supported by the grants of the Hong Kong Research Grants Council (RGC) and the Hong Kong Baptist University Faculty Research Grant (FRG). B. S. X also acknowledge the support of funding in part by the National Natural Science Foundation of China (NSFC) under Grant Nos. 10275007 and 10135010.


text:
# Intermittency is a consequence of turbulent transport in nonlinear systems

## Abstract

Intermittent high-amplitude structures emerge in a damped and driven discrete nonlinear Schrödinger equation whose solutions transport both energy and particles from sources to sinks. These coherent structures are necessary for any solution that has statistically stationary transport properties.

Turbulent flows transferring energy from a stirring range at large scales to the dissipation range at small scales consist of dissimilar components, a broad spectrum of eddies and randomly occurring, intermittent coherent structures. The first cascade leads to Kolmogorov-like finite-flux power spectra. The second component is particularly visible in the anomalous short-scale behavior of the higher order moments of velocity differences. Compared to fully developed three dimensional turbulence, it is easier to quantify these two components in turbulent systems of weakly coupled dispersive waves. In this case, the Kolmogorov-Zakharov spectrum is a stationary finite flux solution of kinetic equations that follow from three and four wave resonances of the weak interactions. The second component emerges since the wave turbulence approximation nnb  is almost never valid at very low and very high wavenumbers where the ’weak’ coupling approximation breaks down, leading to the emergence of fully nonlinear structures dya . In short, despite the fact that the amplitudes are, on the average, small, the weakly nonlinear dynamics can lead to intermittent and localized high amplitude events and anomalies in high order moments. It may also lead to a contamination of low order moments and to power spectra which, at least in some wavenumber ranges, are dominated by strongly nonlinear events. Nowhere is this more evident than in the illuminating studies of Cai, Majda, McLaughlin, and Tabak (CMMT) mmt  (later confirmed by Zakharov, et al. zak ) on damped, driven and on freely decaying weakly nonlinear dispersive one dimensional wave systems. Indeed, CMMT found that, in damped and driven systems, there were some situations in which the nonlinear solutions dominated at almost all scales. In the freely decaying state, they found the Kolmogorov-Zakharov spectra were much more likely to appear. This of course was disappointing because the strength of the Kolmogorov-Zakharov solution of the undamped, undriven kinetic equation is that it describes exactly what would be expected if an energy source at $`k=0`$ feeds at a constant flux rate through an inertial range to a viscous sink at $`k=\mathrm{}`$.
The purpose of this Rapid Communication is to demonstrate in a simple but representative model that in driven, damped systems in which there are fluxes of two conserved densities (energy and particle number), the realization of a statistically steady  state demands the existence of high amplitude coherent structures because they play a crucial role in balancing the energy flux budget. We show that a single cascade of weakly interacting waves transporting particles from a source to a sink leads to a steady net loss of energy of the system. Any steady solution requires contributions from the nonlinear terms in order to offset this energy loss. This is achieved by the intermittent formation and destruction of high-amplitude structures. Fig.1 shows the intermittent emergence of peaks for a

discrete nonlinear Schrödinger equation

$$i\dot{\varphi }_n+\varphi _{n+1}+\varphi _{n1}2\varphi _n+\varphi _n^2\varphi _n^{}=(\varphi ,t)𝒟(\varphi ,t)$$
(1)
for a chain of complex oscillators. $`𝒟`$ is a short-wavelength damping term while $``$ drives the system on long space scales. Fig.1(a) displays the typical long-time behavior of $`|\varphi _n|^2`$ in a sector of 30 lattice sites from a chain of $`N=512`$ oscillators with periodic boundary conditions over 30 time units. $`|\varphi _n|^2`$ is small on average, but high-amplitude structures emerge locally. The peak shows a ’breathing’ behavior, decreasing and increasing irregularly. Fig.1(b) shows the average density $`\rho `$ of oscillators with an amplitude $`|\varphi |`$ for weak and for strong damping forces. The density for small $`|\varphi |`$ is Gaussian, but for larger $`|\varphi |`$ the intermittent high-amplitude structures lead to a slower non-Gaussian decay. In the weakly damped system, the density increases slightly for high amplitudes so that there is a small hump near $`|\varphi |=2`$. The number of these intermittent events increases with the driving force, but the hump remains in the vicinity of $`|\varphi |=2`$.
In the simulations we use time-periodic $`\delta `$-kicks as damping and driving forces $`𝒟`$ and $``$, which allows a simple control of the energy input and output: The driving increases the homogeneous mode $`c_0(1+\lambda _{})c_0`$, and the damping decreases modes on a short space scale as $`c_k(1\lambda _𝒟)c_k`$ where $`c_k=(1/N)_n\varphi _nexp(ikn)`$ and the wavenumber $`k`$ is in the Brillouin-zone $`[0,2\pi ]`$. In numerical studies we apply the damping to the short-wavelength modes with $`7\pi /8k9\pi /8`$, so that the impact of the damping does not decrease with the system’s size. Analytical studies are simplified by restricting the damping to the mode $`k=\pi `$. $`𝒟`$ and $``$ are zero in the intervals of one time unit between the time-periodic kicks. No important changes are found for shorter intervals between the kicks, or for continuous driving and damping.
The dynamics between the kicks is governed by the nonlinear Schrödinger equation that derives as $`i\dot{\varphi }_n=H/\varphi _n^{}`$ from the Hamiltonian

$$\begin{array}{ccc}\hfill H& =& H_2+H_4\hfill \\ & =& _n2\varphi _n\varphi _n^{}\varphi _n\varphi _{n+1}^{}\varphi _n^{}\varphi _{n+1}\hfill \\ & & \frac{1}{2}_n\varphi _n^2\varphi _n^2\hfill \end{array}$$
(2)
The modulus square norm or particle number

$$A=\underset{n}{}\varphi _n\varphi _n^{}=N\underset{k}{}c_kc_k^{}$$
(3)
is a second conserved quantity. In recent studies rune  we have shown that this isolated Hamiltonian system with no damping and driving force forms localized high-amplitude structures as a statistical consequence from the thermalization under the constraint of its two conserved quantities. Such peaks are generated in a self-focusing process of low-amplitude waves with long wavelengths and a low energy. A typical initial condition from which high-amplitude peaks emerge is the Rayleigh-Jeans distribution of the power $`|c_k|^2=(\beta (\omega (k)+\gamma ))^1`$ with a positive temperature $`\beta ^1`$. $`\gamma `$ is the chemical potential, $`\omega (k)=22\mathrm{cos}k`$ is the frequency. As the system approaches its state of maximum entropy, the spectrum of low-amplitude fluctuations becomes flat so that the power is equipartitioned on all modes ($`\beta 0,\gamma \mathrm{}`$). During this transformation of the low-energy spectrum to an equipartitioned spectrum, stable high-amplitude peaks with a negative energy emerge as a by-product of the production of entropy in the low-ampitude waves. No such peaks emerge from low-amplitude short waves with a high energy.

In contrast, the long time behavior of solutions of the damped and driven system is governed, not by the values of $`H`$ and $`A`$, but by the fluxes  of both quantities. Each driving step $``$ changes the amplitudes $`\varphi _n`$ by $`\mathrm{\Delta }\varphi _n=\lambda _{}c_0`$, so that the total number of particles increases by $`\mathrm{\Delta }_{}A=_n\varphi _n\mathrm{\Delta }\varphi _n^{}+\mathrm{\Delta }\varphi _n\varphi _n^{}=2\lambda _{}N|c_0|^2`$ while the damping changes the particle number by $`\mathrm{\Delta }_𝒟A=2\lambda _𝒟N|c_\pi |^2`$ for small values of $`\lambda _{,𝒟}`$. Consequently, there is a flow of particles from the source at $`k=0`$ to the sink at $`k=\pi `$, and the power spectrum (Fig.2(a)) decays with the wavenumber. With most power gathered at low wavenumbers, the spectrum is similar to a thermodynamic spectrum (Fig.2(a)) of the corresponding undamped and undriven Hamiltonian system with a low energy. The high-amplitude structures are generated just by the same mechanism in such an isolated Hamiltonian system and in the damped and driven system where the permanent particle flow maintains the bias of the spectrum. The spectrum decays even exponentially for strong damping forces Fig.2(b). Particle loss and gain are balanced when the flow in a sufficiently large system is constant so that

$$\mathrm{\Delta }_{}A+\mathrm{\Delta }_𝒟A=0$$
(4)
or $`\lambda _{}|c_0|^2=\lambda _𝒟|c_\pi |^2`$. The driving parameter $`\lambda _{}`$ and the damping parameter $`\lambda _𝒟`$ regulate the particle number and the particle flux of the system. By choosing the driving parameter as $`\lambda _{}=\frac{\mathrm{\Delta }A}{2N|c_0|^2}`$ we obtain a constant flow of particles $`\mathrm{\Delta }_{}A=\mathrm{\Delta }A`$ into the system in the numerical simulations. The damping parameter $`\lambda _𝒟`$ is fixed. For a strong damping to forcing ratio, $`\lambda _𝒟=ϵ^2\lambda _{}`$ with $`ϵ^21`$, it follows from (4) that $`|c_\pi |=ϵ|c_0|`$ and

$$\lambda _{}|c_0|=ϵ\lambda _𝒟|c_\pi |$$
(5)
The input and output of particles also changes the quadratic coupling energy and the quartic energy. Gains and losses again have to match in a stationary nonequilibrium state so that

$$\mathrm{\Delta }_{}H_2+\mathrm{\Delta }_{}H_4+\mathrm{\Delta }_𝒟H_2+\mathrm{\Delta }_𝒟H_4=0$$
(6)
In order to understand the role of the strongly nonlinear terms in the energy flow, we analyze how it divides up among these four terms. The change of the quadratic coupling energy $`H_2=N\omega (k)c_kc_k^{}`$ is given by the flux of particles times the frequency $`\omega (k)`$. The influx of particles through the driving force leads to a zero energy influx $`\mathrm{\Delta }_{}H_2=\omega (0)\mathrm{\Delta }A=0`$ since $`\omega (0)=0`$. The damping leads to a loss of coupling energy $`\mathrm{\Delta }_𝒟H_2=\omega (\pi )\mathrm{\Delta }A`$ with $`\omega (\pi )=4`$. Consequently, the particle flux leads to a net loss of coupling energy $`\mathrm{\Delta }H_2=\mathrm{\Delta }_{}H_2+\mathrm{\Delta }_𝒟H_2=4\mathrm{\Delta }A`$.

This energy loss must be balanced by an energy gain in the nonlinear term $`H_4`$ for maintaining a stationary nonequilibrium state. The viscosity changes the nonlinear energy $`H_4=(1/2)\varphi _n^2\varphi _{n}^{}{}_{}{}^{2}`$ as

$$\mathrm{\Delta }_𝒟H_4=2\lambda _𝒟NRe(d_\pi c_\pi ^{})$$

where $`d_k=(1/N)_n|\varphi _n|^2\varphi _ne^{ikn}`$. The driving changes this energy contribution by

$$\mathrm{\Delta }_{}H_4=2\lambda _{}NRe(d_0c_0^{})$$

Fig.3(a) shows $`\mathrm{\Delta }_𝒟H_2`$, $`\mathrm{\Delta }_{}H_4`$ and $`\mathrm{\Delta }_𝒟H_4`$ as functions of time for a chain of 512 oscillators for a total particle flux per time unit $`\mathrm{\Delta }A=0.512`$ and a strong damping to forcing ratio where $`\lambda _𝒟=0.2`$. This corresponds to the simulation of Fig.1(a), the strong damping case of Fig.1(b), Fig.2(b) and to Figs.4(c),(d). The fluctuations indicate single collapse events which have a smaller relative strength for larger system sizes. Fig.3(b) shows the time average of $`\mathrm{\Delta }_𝒟H_2`$, $`\mathrm{\Delta }_{}H_4`$, $`\mathrm{\Delta }_𝒟H_4`$ as functions of the damping coefficient $`\lambda _𝒟`$ and the particle flux $`\mathrm{\Delta }A`$, where the weakly driven case ($`\mathrm{\Delta }A=0.512`$) corresponds to all other simulations. The impact of the driving force on the energy $`\mathrm{\Delta }_{}H_40`$ is negligible compared to the two effects of the damping $`\mathrm{\Delta }_𝒟H_2`$ and $`\mathrm{\Delta }_𝒟H_4`$. This predominance of $`\mathrm{\Delta }_𝒟H_4`$ over $`\mathrm{\Delta }_{}H_4`$ (i.e. $`\lambda _{}Re(d_0c_0^{})\lambda _𝒟Re(d_\pi c_\pi ^{})`$ ) follows from $`\lambda _{}|c_0|\lambda _𝒟|c_\pi |`$ by equation (5) and from $`|d_\pi |𝒪(|d_0|)`$. The second relation is due to the spiky shape of $`|\varphi _n|^2\varphi _n`$ that leads to a flat power spectrum. In addition, the phases of $`d_\pi `$ and $`c_\pi `$ are correlated by the phase velocity of the peaks while $`d_0c_0^{}`$ oscillates randomly as it is governed by low-amplitude fluctuations. $`\mathrm{\Delta }_{}H_4`$ yields a significant negative contribution only for high particle fluxes $`\mathrm{\Delta }A`$ and weak damping forces $`\lambda _𝒟`$ (at $`\lambda _𝒟=0.02`$, $`\mathrm{\Delta }A=2.5`$ in Fig.3b).
For all other cases, both $`\mathrm{\Delta }_{}H_2`$ and $`\mathrm{\Delta }_{}H_4`$ are zero or close to zero, so that (6) reduces to $`\mathrm{\Delta }_𝒟H_2\mathrm{\Delta }_𝒟H_4`$. This yields $`Re(c_\pi ^{}d_\pi )4c_\pi ^{}c_\pi `$ or the inequality

$$|d_\pi |4|c_\pi |$$
(7)
The importance of the nonlinearity $`\mathrm{\Delta }_𝒟H_4`$ in the energy balance is reflected in the system’s first three moments: Defining $`\psi _n=\varphi _n(1)^n`$, the inequality (7) reads $`||\psi _n|^2\psi _n|4|\psi _n|`$. The first moment $`|\psi _n|\sqrt{\mathrm{\Delta }A}`$ increases with the particle flux. The second moment $`|\psi _n|^2=A`$ is the particle number, which is small if the system is strongly damped and weakly driven. Equation (7) shows that the third moment exceeds the first moment, so that the distribution is strongly non-Gaussian. High-amplitude excitations (Fig.1(a)) leading to the non-Gaussian tails of Fig.1(b) are necessary to maintain a stationary nonequilibrium state under the flow of energy and particles.
To understand how the high-amplitude structures balance the energy flow, we compute the input and output of energy and particles as functions of the amplitude $`|\varphi |`$ of the oscillators. The driving kick feeds a number $`\mathrm{\Delta }_{}a(|\varphi |)d|\varphi |=_{n:|\varphi ||\varphi _n||\varphi |+d|\varphi |}2Re(\varphi _n\mathrm{\Delta }_{}\varphi _n^{})`$ of particles to those lattice sites where where the amplitude is between $`|\varphi |`$ and $`|\varphi |+d|\varphi |`$. The total particle gain as used in (4) is $`\mathrm{\Delta }_{}A=_0^{\mathrm{}}\mathrm{\Delta }_{}a(|\varphi |)d|\varphi |`$. Similarly, the forcing changes the nonlinear energy at these lattice sites by $`\mathrm{\Delta }_{}h_4(|\varphi |)d|\varphi |=_{n:|\varphi ||\varphi _n||\varphi |+d|\varphi |}2|\varphi _n|^2Re(\varphi _n^{}\mathrm{\Delta }_{}\varphi _n)`$. The increments of the damping $`\mathrm{\Delta }_𝒟a(|\varphi |)`$ and $`\mathrm{\Delta }_𝒟h_4(|\varphi |)`$ are defined analogously. $`\mathrm{\Delta }_{}h_2(|\varphi |)`$ is again zero, and $`\mathrm{\Delta }_𝒟h_2(|\varphi |)=\omega (\pi )\mathrm{\Delta }_𝒟a(|\varphi |)`$. Fig.4 shows the time average of these influxes and outfluxes of particles and energy as functions of $`|\varphi |`$. In the weakly damped case (Fig.4(a),(b) with $`\lambda _𝒟=0.02`$, $`\mathrm{\Delta }A=0.512`$), the particles are fed into the system mainly at lattice sites with moderate amplitudes $`|\varphi |0.5`$ (Fig.4(a)), while particles are removed mainly at sites with high amplitudes $`|\varphi |2`$. While the input of particles has a negligible effect on the balance of nonlinear energy, the output of particles leads to a gain of energy peaked at $`|\varphi |2`$. Remarkably, a very small number of sites with high amplitudes is responsible for practically the total outflux of particles and the corresponding feedback of nonlinear energy.
This energy gain can be seen for the most simple case of an isolated real peak $`\varphi _n=\chi \delta _{nl}`$ whose height decreases by $`\mathrm{\Delta }_𝒟\varphi =\lambda _𝒟c_\pi \chi `$ during the damping kick, so that the coupling energy changes by $`\mathrm{\Delta }_𝒟H_2=8\chi \mathrm{\Delta }_𝒟\varphi `$. This increases the nonlinear energy by $`\mathrm{\Delta }_𝒟H_4=2\chi ^3\mathrm{\Delta }_𝒟\varphi `$. These fluxes are balanced if the peak amplitude is $`\chi =2`$ which agrees with the results of Figs.4(a),(b).

Particles removed from the tip of the peak have to be replaced constantly through the driving mode at $`k=0`$. During the Hamiltonian evolution, the focusing process increases the peak again while short-wave fluctuations are radiated.
The input of particles in the strongly damped equation (Fig.4(a),(b) with $`\lambda _𝒟=0.2`$, $`\mathrm{\Delta }A=0.512`$) is maximal for low amplitudes $`|\varphi |0.4`$ similarly to the weakly damped case of Fig.4(a). The damping leads to a significant loss of particles at $`|\varphi |1.5`$, and to a tiny loss of particles in the domain of low-amplitude fluctuations at $`|\varphi |0.2`$, but, paradoxically, also to a particle gain at $`|\varphi |0.7`$. While there is still a net loss of particles, this partial recycling of particles increases the ratio of energy gain per particle loss. This allows the system to satisfy the balance condition (7) with a peak amplitude less than $`|\varphi |=2`$. As a simple model for this mechanism we assume two peaks at the sites $`l`$ and $`l+3`$ where $`\varphi _l=\chi >\xi =\varphi _{l+3}`$ are real and positive. A damping kick decreases the higher peak to $`\varphi _l=\chi \mathrm{\Delta }_𝒟\varphi `$ and increases the lower peak to $`\varphi _{l+3}=\xi +\mathrm{\Delta }_𝒟\varphi `$. The decrease of the higher peak $`\varphi _l`$ models the particle loss at $`|\varphi |1.5`$ in Fig.4(c) and increase of the lower peak $`\varphi _{l+3}`$ represents the gain of particles at $`|\varphi |0.7`$. The coupling energy loss is $`\mathrm{\Delta }_𝒟H_2=8(\chi \xi )\mathrm{\Delta }_𝒟\varphi `$ and the nonlinear energy gain is $`\mathrm{\Delta }_𝒟H_4=2(\chi ^3\xi ^3)\mathrm{\Delta }_𝒟\varphi `$. The gain of nonlinear energy can outweigh the loss of coupling energy for $`\chi >2/\sqrt{3}`$.
In summary, coherent structures are an essential component of the transport of particles from the large scales to small dissipation scales. The particle flow leads to a steady loss of coupling energy, so that the fluctuations have a low ratio of energy per particle. Such long-wave fluctuations form high-amplitude structures rune , and it is the role of the coherent structures to balance the energy input and output by exploiting the nonlinear component of the energy. To offset this energy loss, high amplitude structures containing a significant amount of negative nonlinear energy must be formed and destroyed. Energy is transfered from the nonlinearity $`H_4`$ to the coupling $`H_2`$ during the peak growth, and pruning or destroying the peaks increases the systems total energy. The peaks are essentially the pipes through which energy and particles flow under the constraints of the balance conditions.
Obviously, this effect follows purely from the dispersion and the nonlinearity, and it is not restricted to the one dimensional discrete system of our simulations. A similar behavior was found dya  for the two-dimensional continuous focusing nonlinear Schrödinger equation when the forcing is applied at a midscale wavenumber and the damping at large wavenumbers. This leads to an energy flux to high wavenumbers, and to an inverse particle flux to low wavenumbers. The particle flux builds condensates from which collapsing solitons emerge that carry particles to the dissipation scale and feed both particles and energy back into the wave field at these high wavenumbers. With an almost Gaussian statistics of the waves and a Poissonian distribution of the intermittent coherent events, the system resembles a two species gas. The presence of collapses does not appear to affect low moment statistics of two- or three-dimensional systems if the damping is sufficiently strong, while the collapses contaminate the power spectrum generated by wave-wave interactions in nonintegrable one dimensional systems. The important point to stress is that no particles would be dissipated and no steady state would be achieved without strongly nonlinear collapses.


text:
# Transition radiation by matter-wave solitons in optical lattices

## Abstract

We demonstrate that matter-wave solitary pulses formed from Bose condensed atoms moving inside optical lattices continuously radiate dispersive matter waves with prescribed momentum. Our analytical results for the radiation parameters and the soliton decay rate are found to be in excellent agreement with numerical modelling performed for experimentally relevant parameters.

Recent observations of matter-wave solitons  have clearly been amongst the most breaking achievements in the burgeoning field of Bose-Einstein condensation of dilute atomic gases. Balance between the spatial dispersion of matter waves and repulsive or attractive interatomic interactions ensures existence of dark  or bright  solitons, respectively. Dispersion of the atomic condensates can, however, be reversed by embedding the condensate into a periodic potential created by standing light waves, i.e. optical lattice . The idea of changing the dispersion sign and of the possible observation of the bright matter-wave solitons in the condensates with repulsive interatomic interaction has been around for a while, see  and references therein, and more theoretical results have been produced recently, see, e.g., , in the view of the rapid maturing of the experimental techniques . The concept of the dispersion control by periodic potentials is also well known in solid-state physics  and a very active topic of research in nonlinear optics, see, e.g., .

To understand the initial motivation which leads to the results described below, it is instructive to recall the effect of the transition radiation known from classical electrodynamics . Transition radiation is a continuous emission of electromagnetic waves by a charged particle moving with a constant velocity in spatially inhomogeneous medium. This radiation is emitted because the field created by the particle has different characteristics in different parts of the medium. When the particle moves, the field reorganizes itself continuously and shakes off some of its parts in form of the radiation. We expect that similar phenomenon should take place with matter-wave solitons in optical lattices.

In order to change the dispersion sign of matter waves forming the initially resting one-dimensional packet embedded into the optical lattice, one needs to position it in momentum space somewhere between the inflexion point of the energy-momentum, i.e., dispersion, characteristic and the edge of the Brillouin zone, see Fig. 1(a) . Only exactly at the edge of the zone does the group velocity go to zero. Solitons with a spread of quasi-momenta centered at the edge will therefore be the only resting bright solitons in the condensates with repulsive inter-atomic interaction. Our primary interest below is, however, moving solitons. Once a solitonic wave packet moves through the periodic potential one can expect that its structure will not be able to instantly readjust itself to perfectly fit the conditions that the local density maxima are positioned at the center of the local potential minima. Therefore the moving soliton has to continuously overcome the energy barrier, which leads to the continuous emission of matter waves. This energy barrier is analogous to the Peierls-Nabarro potential known for the solitons in discrete systems . Using analogy with the electromagnetic case we term this effect as transition radiation of matter waves.

We start our analysis from the Gross-Pitaevskii (GP) equation describing evolution of the macroscopic wave function of the zero-temperature Bose-Einstein condensate (BEC) interacting with off-resonant standing light wave. We assume that the condensate is tightly confined by the external harmonic potential along the $`Y`$ and $`Z`$ directions having the trap frequencies $`\omega _{Y,Z}=2\pi \times 400s^1`$ and that any deviations of the potential along $`X`$-axis from the $`\mathrm{sin}^2k_lX`$ produced by the intensity of the standing laser field with wavenumber $`k_l`$ can be disregarded. We take for our estimates that $`k_l=2\pi /800nm^1`$ and consider BEC made of $`{}_{}{}^{87}Rb`$ atoms with a two-body scattering length $`a5.4\times 10^9m`$. The characteristic transverse width of this BEC is then given by $`w=\sqrt{\mathrm{}/(m\omega _Y)}0.5\mu m`$. Assuming that the condensate profile along the $`Y,Z`$ directions is given by the lowest mode of the harmonic potential, we can derive the one-dimensional GP equation, which describes dynamics of the $`X`$ dependent part of the full wave function. The dimensionless normalized form of this equation is :

$$i_t\psi =_x^2\psi \beta \psi \mathrm{cos}2x+|\psi |^2\psi .$$
(1)
Here the dimensionless time $`t`$ and spatial coordinate $`x`$ are measured respectively in the units of $`T_0=2m/(\mathrm{}k_l^2)5\times 10^5s`$ and $`1/k_l0.13\mu m`$. Meaning of the dimensionless parameter $`\beta `$ is easily inferred from the expression for the lattice potential in the physical units, which is taken as $`4\beta E_r\mathrm{sin}^2k_lX`$, where $`E_r=\mathrm{}^2k_l^2/(2m)\mathrm{}\times 20kHz`$ is the recoil energy. To estimate number of atoms $`N`$ in the condensate we introduce the effective area in the $`(Y,Z)`$-plane, $`A_{eff}=2\pi w^21.5\mu m^2`$, and the atom density $`n10^{14}cm^3`$. Then $`N=|\psi |^2𝑑x\times A_{eff}^2k_ln^{2/3}/(8\pi |a|)4\times 10^3\times |\psi |^2𝑑x`$.

We proceed by expanding $`\psi `$ over the Bloch functions, $`b(x,k)`$ : $`\psi (x,t)=𝑑k\stackrel{~}{\psi }(k,t)b(x,k)`$. Here $`b(x,k)`$ are eigenfunctions of the operator $`\widehat{L}=_x^2+\beta \mathrm{cos}2x`$, such that $`\widehat{L}b=ϵ(k)b`$, where $`k`$ is the quasi-momentum introduced as $`b(x,k)=g(x,k)e^{ixk}`$ , $`ϵ(k)`$ is the energy of the non-interacting, i.e. linear, matter-waves, and $`g(x,k)`$ is the function with the spatial period $`\pi `$. The first allowed, first forbidden (gap) and a small part of the 2nd allowed energy bands are shown in the energy-momentum plot in Fig. 1(a).

Choosing $`k=k_s[1,1]`$, we expand $`g(x,k)`$ in a Taylor series around $`k=k_s`$ and demonstrate that

$$\psi (x,t)=e^{ik_sx}\widehat{G}_sA(x,t).$$
(2)
Here, $`A(x,t)=𝑑k\stackrel{~}{\psi }(k,t)e^{i(kk_s)x}`$ and $`\widehat{G}_s`$ is the linear differential operator: $`\widehat{G}_s=_{n=0}^{\mathrm{}}(1/n!)_k^ng(x,k_s)(i_x)^n`$. Similarly $`\widehat{L}\psi `$ can be represented as

$$\widehat{L}\psi =e^{ik_sx}\widehat{G}_s\widehat{}_s[i_x]A(x,t),$$
(3)
where $`\widehat{}_s[i_x]=_{m=0}^{\mathrm{}}(1/m!)_k^mϵ(k_s)(i_x)^m`$ is the energy operator. Subscript ’$`s`$’ refers to quantities and functions calculated for $`k=k_s`$.

After substitution of Eqs. (2), (3) into Eq. (1), we then replace $`\widehat{G}_s`$ with its 1st order approximation $`g(x,k_s)g_s`$ and $`\widehat{}_s`$ with $`\widehat{}_{2s}ϵ_siϵ_s^{}_x1/2ϵ_s^{\prime \prime }_x^2`$. Here $`ϵ_s^{}`$ and $`ϵ_s^{\prime \prime }`$ are, respectively, group velocity and group velocity dispersion of matter-waves. Assuming smallness of the nonlinearity and averaging the resulting equation over $`g_s`$ we derive the renowned nonlinear Schrödinger (NLS) equation $`i_tA=\widehat{}_{2s}A+\alpha _s|A|^2A`$, where $`\alpha _s=𝑑x|g_s|^4`$. Below we are interested in the bright solitons with $`ϵ_s`$ belonging to the first allowed energy band. These are given by $`A_s=R_s(\xi )\mathrm{exp}\{i\kappa tiϵ_st\}`$, where

$$R_s(\xi )=\sqrt{\frac{2\kappa }{\alpha _s}}sech\left\{\xi \sqrt{\frac{2\kappa }{ϵ_s^{\prime \prime }}}\right\},\xi =xϵ_s^{}t,$$
(4)
$`sgn\kappa =sgn\alpha _s=sgnϵ_s^{\prime \prime }`$ and $`\kappa 0`$ is the nonlinearity induced energy shift. Details of the derivation of NLS equation from Eq. (1) has been previously published using the method of multiple scales . Approximations made above and in  imply that $`\psi `$ has sufficiently narrow spread of quasi-momenta around $`k=k_s`$. However, the approach introduced here is readily adaptable to give an access to the small amplitude corrections having quasi-momenta detuned far from $`k_s`$, see Eqs. (8), (9) below. Note, here that mobile envelope solitons (4) are very different in their properties from practically immobile solitons occupying primarily one or few lattice sites and considered, e.g., in .

To introduce the effect studied and explained in this work, we first present the results of numerical modelling of Eq. (1) with initial conditions in the form $`g_s(x)R_s(x)e^{ik_sx}`$ for values of $`k_s(k_0,1]`$, $`ϵ^{\prime \prime }(k_0)=0`$, ensuring that the effective mass, $`1/ϵ_s^{\prime \prime }`$, is negative, see Fig. 1(a). Taking $`k_s=1`$, i.e. fixing soliton parameters at the point corresponding to the zero of the group velocity, we observe formation of the ideal solitary pulse . The Fourier spectrum of this solution contains series of the equidistant peaks, the location of which is determined by the spectrum of the corresponding Bloch function $`b(x,k_s=1)`$. Note, that the wave-numbers $`q`$, parameterizing the Fourier transform of $`\psi `$, $`\psi =𝑑x\stackrel{~}{\mathrm{\Psi }}(q)e^{iqx}`$, are linked with quasi-momenta $`k`$ as $`q=k\pm 2n`$. Here and below $`n=0,1,2,\mathrm{}`$ numbers the Brillouin zones. For values of $`k_s1`$ we have observed the quasi-solitonic pulses, which, while travelling, leave behind the trail of small amplitude radiation, Fig. 2(a). The radiation effect becomes noticeably stronger for $`k_sk_0`$. Spectra of the radiating solitons have a distinct peak, see Fig. 2(b), which is absent in the spectra of the ideal resting solitons. The overall results of the extensive series of numerical experiments unambiguously indicate that solitary pulses moving through the optical lattice continuously emit radiation with certain spectrally localized quasi-momenta. By analogy with electromagnetic case we term this radiation as transition radiation of matter waves.

The initial conditions used above in the form of the $`sech`$ envelope superimposed on the Bloch function $`g_s(x)e^{ik_sx}`$ are difficult, though probably not impossible to prepare in the real experiments. However, current experimental techniques allow straightforward preparation of the gaussian matter-wave packets with spectrum centered around the zero quasi-momentum and setting lattices in motion. The lattice moving with velocity $`2k_s`$ will then effectively shift the central momentum of the wave packet to $`k_s`$. $`2k_s`$ equals to the group velocity of the free, i.e. without the lattice, matter waves with quasi-momentum $`k_s`$. In turn, the backward scattering of matter waves from the moving lattice is expected to create the second strong and other peaks in the spectrum of the wave-packet. The results of modelling of Eq. (1) with moving potential $`\mathrm{cos}2(x2k_st)`$ and initial conditions in the form of the gaussian packet are shown in Figs. 2(c,d). Taking into account shift of the axes, one can see that these results are in remarkable agreement with those obtained using the solitonic initial conditions, see Figs. 2(a,b). Good resolution of the radiation peaks in Figs. 2(b,d) suggests even that they can possibly be recorded experimentally.

To understand and give analytical interpretation of the observed radiation we develop a perturbative approach, allowing us to predict quasi-momenta and the amplitude of the emitted wave and, thereby allowing an estimate of the decay rate of a solitary pulse. To proceed we form the ansatz

$$\psi =\psi _s(x,\xi ,t)+\phi (t,x),$$
(5)
where the first term approximates the solitonic part of the wave function, $`\psi _s=g_s(x)A_s(\xi ,t)e^{ik_sx}`$, and the second one is an arbitrary perturbation. Substituting (5) into (1) and assuming that $`|\phi |`$ is small we find that evolution of $`\phi `$ is governed by

$$i_t\phi +\widehat{L}\phi 2|\psi _s|^2\phi \psi _s^2\phi ^{}=S(x,t),$$
(6)
where $`S(x,t)=i_t\psi _s\widehat{L}\psi _s+|\psi _s|^2\psi _s`$ is the source term, which is different from zero because $`\psi _s`$ is not an exact solution of Eq. (1).

The structure of the radiation tail observed in the numerical modelling corresponds to the spatially extended lattice eigenmode. The natural mechanism for excitation of a selected eigenmode from the continuum is the energy and wave-number resonance with one Fourier component of the source term $`S`$. The resonance condition ensures that there exists a lattice mode which is always in-phase and therefore interferes constructively with one Fourier component of the source term. The latter can be represented in the form $`S=e^{ik_sxiϵ_sti\kappa t}_jh_j(x)f_j(\xi )`$, where the sum is taken over all the terms appearing in the right-hand side after substitution of the explicit expression for the soliton and $`h_j(x)`$ are some functions with a period of $`\pi `$. Replacing $`f_j(\xi )`$ through their Fourier integrals $`f_j(\xi )=𝑑Q\stackrel{~}{f}_j(Q)e^{iQ\xi }`$ and recalling that $`\xi =xϵ_s^{}t`$, one can easily find that the lattice modes satisfying the resonance condition are those, which have quasi-momenta $`k=k_r`$ obeying the condition $`ϵ(k_r)=\kappa +ϵ_s+ϵ_s^{}Q_r`$. $`Q_r=Q(k_r)`$ is the detuning of the wave-number of the resonant wave $`q_r=k_r\pm 2n`$ from the soliton quasi-momentum $`k_s`$, i.e. $`Q_r=k_r\pm 2nk_s`$. Using that $`ϵ(k_r)=ϵ(q_r)`$ the resonance condition can be rewritten in the form

$$ϵ(q_r)=\kappa +ϵ_s+(q_rk_s)ϵ_s^{}.$$
(7)
The geometrical meaning of Eq. (7) is clear. The right-hand side of (7) equals to the energy of the Fourier component of the soliton (4), while its left hand side is simply the energy of the linear dispersive wave. Eq. (7) can be solved for $`q_r`$ by plotting the tangent line to the periodically extended dispersion characteristic at the point $`ϵ=ϵ_s`$. Then one should make the parallel upshift of this tangent by $`\kappa `$ and find points of the intersection with the dispersion characteristic itself. Topologically it is self-evident that only resting solitons do not produce any real roots of Eq. (7) and do not radiate into the lattice modes. For $`ϵ_s^{}0`$ Eq. (7) has infinitely many real roots corresponding to different values of $`n`$. Root $`q_r`$ for $`n=0`$, i.e. when $`q_r=k_r`$, as a function of $`k_s`$ is shown in Fig. 1(b). Dots on this graph show $`(k_r,k_s)`$ pairs measured from the direct numerical modelling of Eq. (1). Modelling of Eq. (1) has not revealed any resonances in the Brillouin zones with $`n0`$. It indicates that coupling into the higher-order resonances is negligible, primarily because spectral strength of the source term is very weak for the values of $`Q_r`$ with $`n0`$. Note, that the secondary less intense spectral peaks, seen in Figs. 2(b,d), are described by the peaks in the Fourier spectra of $`b(x,k_{r,s})`$, and are given by $`q_{r,s}=k_{r,s}\pm 2n`$. At the same time, multiple roots of Eq. (7) are not linked by any simple algebraic expression.

In order to calculate the amplitude of the radiated wave we form the anzats

$$\psi \psi _s(x,\xi ,t)+g_r(x)e^{ik_rxiϵ_rt}W(t,\xi ),$$
(8)
where subscript ’$`r`$’ refers to the quantities calculated at $`k=k_r`$. Substituting (8) into Eq. (1) we take into account that $`\widehat{L}g_r(x)e^{ik_rx}W(t,\xi )g_r(x)e^{ik_rx}\widehat{}_rW(\xi ,t)`$. The latter expression can be easily inferred by comparison with Eq. (3). Assuming that $`\delta k_s|k_sk_r|2`$, where $`\delta k_s`$ is the width of the soliton in the quasi-momentum space and $`2`$ is the width of the Brillouin zone, and disregarding terms nonlinear in W, one gets the averaged equation for the amplitude of the radiation field: $`\stackrel{}{s}=i_t\stackrel{}{W}+\widehat{}_r\stackrel{}{W}`$, where $`\stackrel{}{s}=(P,P^{})^T`$, $`\stackrel{}{W}=(W,W^{})^T`$, $`P=e^{i\xi (k_sk_r)}[\alpha _4\{\widehat{}_s\widehat{}_{2s}\}R_s+\alpha _3|R_s|^2R_s]`$, $`\widehat{D}_r[i_\xi ]=ϵ_r+iϵ_s^{}_\xi +\widehat{}_r[i_\xi ]`$,

$`\widehat{}_r=\left[\begin{array}{cc}\widehat{D}_r+\alpha _2|R_s|^2& \alpha _1R_s^2e^{2i\xi (k_sk_r)}\\ \alpha _1^{}R_s^2e^{2i\xi (k_rk_s)}& \widehat{D}_r^{}\alpha _2^{}|R_s|^2\end{array}\right],`$ (11)
$`\alpha _1=𝑑xg_s^2g_{r}^{}{}_{}{}^{2}`$, $`\alpha _2=2𝑑xg_s^2|g_r|^2`$, $`\alpha _3=\alpha _s𝑑xg_sg_r^{}𝑑x|g_s|^2g_sg_r^{}`$, and $`\alpha _4=𝑑xg_sg_r^{}`$.

The approximate solution we obtain for $`W`$ is

$$W(t,\xi )iC\left\{\mathrm{\Theta }(\zeta \xi )\mathrm{\Theta }(\zeta \xi +\zeta t(ϵ_s^{}ϵ_r^{}))\right\}.$$
(12)
Here $`\zeta =sgn(ϵ_r^{}ϵ_s^{})`$, $`\mathrm{\Theta }`$ is the Heaviside function and $`C`$ is the amplitude, which can not be generally expressed in a closed analytical form and, therefore, was calculated numerically. $`C`$ characterizes the spectral intensity of the source term for $`k=k_r`$. Heaviside functions in Eq. (12) describe the tail of the radiation field having the length $`t|ϵ_r^{}ϵ_s^{}|`$. The tail starts at the soliton ($`\xi =0`$) and extends beyond ($`ϵ_r^{}<ϵ_s^{}`$) or in front ($`ϵ_r^{}>ϵ_s^{}`$) of it.

Substituting the ansatz (8) into the conservation law $`_t𝑑x|\psi |^2=0`$, we can estimate the rate, $`\gamma `$, of the transfer of particles from the soliton to the radiation: $`\gamma =C^2|ϵ_s^{}ϵ_r^{}|`$. Plot of $`\mathrm{log}_{10}\gamma `$ vs $`k_s`$ for two values of $`\kappa `$ is shown in Fig. 3(a). Naturally, the rate of transfer increases, when detuning $`|k_rk_s|`$ decreases for $`k_sk_0`$, see Fig. 1(b). This is because for $`k_sk_0`$ the lattice mode is resonant with the most intense central part of the soliton spectrum. Contrary, $`|k_rk_s|`$ increases for $`k_s1`$, and the radiation amplitude decays almost exponentially. For example, semi-analytical Eq. (12) gives that the initial condition with $`k_s=0.87`$ used to generate Fig. 2 has the initial decay rate $`3\times 10^4`$ particles per second. Providing that this rate is constant in time, the soliton life-time would be $`0.033`$s. However, radiation carries away both density and momentum from the solitonic part of the field. Therefore solitonic parameters $`\kappa `$ and $`k_s`$ are in fact functions of time. In particular, we have found that the radiation results in convergence of the soliton quasi-momentum to some limit value, which is always closer to $`1`$ than the initial $`k_s`$. Thus, radiation emission slows the soliton. For example, taking initial conditions with $`\kappa =0.01`$ and $`k_s=0.87`$ we have observed that over $`0.05`$s ($`1000`$ dimensionless time units) $`k_s`$ shifts to $`0.887`$, see Fig. 3(b), and $`\kappa `$ to $`0.006`$. The relatively small increase of $`k_s`$ is accompanied by the decrease of the soliton decay rate, which drops, accordingly with plots shown in Fig. 3(a), from $`3\times 10^4`$ to $`10^2`$ particles per second. Let us stress again that the soliton life-time tends to infinity when $`k_s1`$. Correspondingly, the estimated soliton life-time increases dramatically to $`10`$s. Fig. 3(b) also shows the numerically computed slowdown in the growth of the number of particles in the radiation component of the condensate. Soliton decay rate inferred from these data and corresponding theoretical points calculated for instant values of the soliton parameters $`\kappa `$ and $`k_s`$ are in good agreement, see inset in Fig. 3(b), which confirms validity of our theoretical method.

In summary, we have reported transition radiation by matter-wave solitons moving through the optical lattice. This effect extends family of the already known and related quantum radiative effects such as, e.g., sound emission by precessing quantized vortices  and by dark matter-wave solitons oscillating in a harmonic trap . Note, also that our main conclusions and techniques can be used to predict and analyze radiation by spatial optical solitons moving in nonlinear photonic crystals . This work was partially supported by the INTAS project 211-855.


text:
# Suppression of Superconductivity in YBCO/LCMO Superlattices

## Abstract

The competition of superconductivity and magnetism in superlattices composed of alternating YBa2Cu3O7-d and La0.67Ca0.33MnO3 thin films is investigated using low-energy optical spectroscopy. The thickness of the superconducting YBCO layers is varied from 30 nm to 20 nm while the thickness of the magnetic LCMO layers is kept constant at 20 nm. We clearly observe that the superconducting condensate density in the superconducting state of superlattice is drastically reduced by the magnetic subsystem which may be connected with proximity effects that distort the gap symmetry and thus suppress superconductivity.

preprint: Superlattices/PRB/Brief Report

Since the early 1960s the interplay of superconductivity (SC) and ferromagnetism (FM) has continuously drawn attention Fisher73 . At first glance SC and FM exclude each other because the magnetic exchange field breaks the Cooper pairs. Nevertheless the possible coexistence of both phenomena due to a spatially modulated order parameter was suggested by Fulde, Ferrell, Larkin and Ovchinnikov Fulde64 ; the correspondent state was never observed in conventional SC. A good way to study the effect is to bring the nano-sized layers of the two materials into contact and to view the SC-FM interaction on the level of the proximity effects. The exploration of multilayer systems composed of alternating conventional SC and metallic FM layers was subject of recent reviews Chien99 ; Izyumov02 .

In the last years artificially grown superlattices (SLs) consisting of high-$`T_c`$ superconductors and manganese oxides which exhibit the phenomenon of colossal magnetoresistance have attracted increasing interests for various reasons. On one hand, SLs were developed as an important tool to explore the interplay between the two antagonistic SC and FM ground states, and on the other hand the injection of spin-polarized carriers can lead to new SC switching devices Goldman01 . Due to the structural similarity of the two classes of perovskite compounds, it is possible to construct a unique combination of SC cuprate and strongly FM manganese layers. It has been demonstrated by different groups Jakob95 ; Przyslupski99 ; Habermeier01 ; Prieto01  that although SC and FM are preserved in each subsystem, the superconducting transition temperature $`T_c`$ and the Curie temperature of the magnetic ordering $`T_{\mathrm{mag}}`$ are considerably suppressed. Transport and magnetization measurements clearly indicate that spin-polarized quasiparticle injection in YBa2Cu3O7-d/La0.67Ca0.33MnO3 (YBCO/LCMO) SL can strongly decrease $`T_c`$ and the critical current $`J_c`$ Vasko97 ; Yeh99 . Particularly interesting is the almost complete spin polarization in LCMO and the low carrier density in YBCO which seems to be suppressed even further in the case of SL, as was shown by optical spectroscopy Holden03 . The advantage of optical techniques is that one can directly monitor the temperature dependence of the “strength” of the conducting and superconducting condensates expressed by the spectral weight, i.e. the area under the dynamical conductivity spectrum. Up to now, however, optical measurements failed to even detect a sign of (weakened) superconductivity in SL. In this paper, we report on the first measurements of optical properties of YBCO/LCMO superlattices at the radiation quantum energies low enough to allow for observation of SC condensate response. We find strongly reduced SC superfluid density and an enhanced London penetration depth in the superlattices and ascribe it to proximity effects.

The samples are fabricated by pulsed laser deposition techniques on LaSrGaO4 (LSGO) (001) substrates which has low dielectric losses ($`\mathrm{tan}\delta =5\times 10^5`$) in the used frequency range and a favorable in-plane lattice match between LSGO ($`a=b=0.3840`$ nm) and YBCO ($`a=0.3836`$ nm, $`b=0.3883`$ nm). The growth conditions of the $`c`$-axis oriented YBCO and LCMO layers were optimized in the way to achieve optimal superconductivity and ferromagnetism. Detailed information about the fabrication of SLs can be found in Ref. Habermeier01 . Alternating LCMO and YBCO layers (five of each) are deposited as sketched in the inset of Fig. 1. For all SLs used, the thickness $`d`$ of the LCMO layers is kept at 20 nm, while the thickness of the YBCO layers changes from 20 nm to 30 nm. We have studied two superlattices, which will be referred to as 20/20 and 30/20 - the ratio of thicknesses of the YBCO to the LCMO layers. Parameters of the samples used in our study are listed in Table 1. Only the corresponding ($`00h`$) peaks are observed in X-ray diffraction $`\theta /2`$ measurements for YBCO, LCMO, and the LSGO substrate. The SC transition temperature $`T_c`$ and the magnetic transition temperature $`T_{\mathrm{mag}}`$ are determined by dc-resistance $`R(T)`$ and SQUID magnetization measurements, respectively (Fig. 1); the values are close to what we have obtained for SLs grown on SrTiO3 substrate Habermeier01 ; Holden03 . In addition, single LCMO film ($`d=20`$ nm) and single YBCO film ($`d=30`$ nm) were also prepared on LSGO substrates so that we could study optical properties of thin YBCO and LSCO films not incorporated in the SL. For a single thin YBCO film of 30 nm we find a reduced SC transition temperature of approximately 70 to 80 K.

Optical measurements were done in transmission mode at the lower end of the far-infrared spectral range (6 \- 30 cm-1) at temperatures $`5\mathrm{K}<T<300`$ K, with the help of the millimeter-submillimeter coherent-source spectrometer utilizing backward-wave oscillators as monochromatic frequency-tunable radiation sources, as described elsewhere in detail Kozlov98 . In a Mach-Zehnder arrangement we measure the transmission coefficient Tr and the phase shift $`\phi _t`$ of the radiation passing through the plane-parallel substrate with the SL on it. Spectra of both, conductivity $`\sigma (\omega )`$ and dielectric permittivity $`ϵ(\omega )`$, of single films or of composite SL are obtained directly (no Kramers-Kronig relations used) from Tr and $`\phi _t`$ on the basis of the general Fresnel expressions for the layered system DresselGruner  (the dielectric parameters of the substrate are measured beforehand). In these measurements the probing radiation passes through the superlattice which is considered as a single layer and whose effective optical properties are thus averaged over those of SC and FM layers composing it, in conditions when the radiation wavelength is much larger compared to the layer thickness; we note that this transmission technique is thus more sensitive to the optical properties of the films when compared to the optical reflectivity techniques.

Fig. 2 shows the frequency dependence of the conductivity $`\sigma (\omega )`$ and dielectric permittivity $`ϵ(\omega )`$ of single LCMO and YBCO films recorded at different temperatures. For the LCMO film, basically no frequency dependence is observed for $`\sigma `$ whose value increases at lower temperatures, the spectra of $`ϵ`$ show a slight increase to high frequencies and do not significantly depend on the temperature - both observations indicate good metallic (Drude-like) properties. In the case of YBCO, displayed in Fig. 2c and d, both $`\sigma `$ and $`ϵ`$ are frequency independent at high temperatures (metallic behavior) but reveal a pronounced dispersion at 5 K. A strong increase of $`\sigma `$ is seen towards the lowest frequencies, similar to that observed in all high-Tc’s and connected to uncondensed (normal) carriers Tanner92 ; Bonn96 ; such dispersion leads to a peak in the temperature dependence of the conductivity, in case of our film - around 50 K for the frequency of 25 cm-1. The dielectric permittivity $`ϵ`$ shows in the SC state a divergent dispersion $`ϵ(1/\omega )^2`$ which represents the dielectric response of the zero-frequency delta-function in the conductivity spectrum, responsible for the infinite DC conductivity DresselGruner . Both behaviors of $`\sigma `$ and $`ϵ`$ of the YBCO film clearly display the appearance of a SC state. We have analyzed the 5 K spectra of YBCO film modeling the response of the condensed SC carriers as $`ϵ_s(\omega )=(\omega _p^s/\omega )^2`$ and a Drude expressions DresselGruner  for $`\sigma (\omega )`$ and $`ϵ(\omega )`$ for uncondensed carriers contribution; the results are listed in Table 1. For the plasma frequency of the superconducting carriers at $`T=5`$ K we obtain $`\omega _p^s/2\pi =1780`$ cm-1 from which we get the value of the London penetration depth $`\lambda _L=c/\omega _p^s=0.9\mu `$m ($`c`$ is the velocity of light). Both values are significantly different from those obtained for single crystals and thicker films Hardy93 ; Djordjevic98 : the value of $`\lambda _L`$ is strongly enhanced and $`\omega _p^s`$ decreased, probably due to the effect of the film thickness, which will be discussed in a subsequent publication.

We now turn to the dielectric responses of the superlattices for which the frequency and temperature dependences are shown in Figs. 3 and 4. Similar to the single YBCO film, the dielectric permittivity spectra of SLs also reveal a signature of a superconducting response (Fig. 3) – a diverging behavior at the lowest temperature, $`ϵ(\omega )(1/\omega )^2`$. At the same time the absolute values of $`ϵ`$, which represent a measure of the SC condensate density, are noticeably, by 6 to 10 times, reduced compared to a single YBCO film. This means that while the superconductivity is still maintained in SLs, its “strength” is strongly suppressed. By fitting the spectra of $`\sigma (\omega )`$ and $`ϵ(\omega )`$ for a SLs as described above we obtain the plasma frequencies of the SC condensate and the London penetration depths, the values given in Table 1 remark2 . We find a strong decrease of the SC carriers density ($`\omega _p^s`$) and a correspondent increase of $`\lambda _L`$ when going from the single YBCO film to the 30/20 and further to the 20/20 superlattice.

The conductivity spectra of SLs also show behaviors different from that in the single YBCO film. The absolute values of the conductivities for the 20/20 and 30/20 SLs are reduced and no significant dispersion is observed, both indicating an increase of the scattering rate for the normal (uncondensed) carriers.

The observed suppression of SC in the superlattices cannot be caused by such effects as incomplete oxygenation or surface disorder, as was shown in Habermeier01  and Sefrioui2002 . We also note that we were not able to fit the measured spectra of transmission coefficient and phase shift by considering the SL as consisting of 5 YBCO and 5 LCMO layers each having properties of single YBCO and LCMO films. This means that these films, when incorporated in a SL, cannot be considered as independent layers. We thus conclude that it is the interaction of FM and SC condensates that is responsible for the strong reduction of the observed SC carrier density. Our findings are in accordance with estimations of the spin diffusion length of approximetely 30 nm Soltan03 .

In the study of impurity doped YBCO, the normalized critical temperature $`T_c/T_{c0}`$ ($`T_{c0}`$ corresponds to the transition temperature of a high quality YBCO) is found to show the dependence on the normalized zero temperature superfluid density $`n_s/n_0`$ Ulm95 ; Franz97  which is quite different from that predicted by a standard Abrikosov-Gor’kov theory Abrikosov . Using $`T_{c0}`$=92 K (YBCO single crystal), we found that the ratio $`T_c/T_{c0}`$ for our superlattices shows behaviors similar to that in Zn- and Ni-doped YBCO thin films ; those results were explained by a gapless superconductivity. In fact for low-$`T_c`$ SC/FM superlattices Sun et al. Sun02  recently proposed that gapless superconductivity appears in both FM and SC regions near the interface. In the superlattices, the influence of the FM layers upon the SC layers is expressed via a strong plane perturbation, magnetic exchange field and proximity effect. As first pointed out by Abrikosov and Gor’kov Abrikosov , gapless superconductivity is a common feature in SC in the presence of strong perturbation. In the gapless region, the electronic properties are drastically modified and lead, for example, to a linearly temperature dependent specific heat, reduced $`T_c`$ and enhanced London penetration depth Maki69 , in accordance with our optical observations.

In summary, from our low-energy optical experiments we find a strong weakening of the superconducting response of superconducting/ferromagnetic YBCO/LCMO superlattices when the thickness of the YBCO SC layer in the superlattice decreases: considering the superlattice as an effective medium we find a strong decrease of the effective density of the superconducting condensate and correspondent enhancement of the London penetration depth. These results indicate a strong suppression of superconductivity in YBCO/LCMO superlattices which can be attributed to the pair-breaking effect of the ferromagnetic layers.


text:
# Advertising, consensus, and ageing in multilayer Sznajd model

e-mail: ab127@uni-koeln.de

Abstract: In the Sznajd consensus model on the square lattice, two people who agree in their opinions convince their neighbours of this opinion. We generalize it to many layers representing many age levels, and check if still a consensus among all layers is possible. Advertising sometimes but not always produces a consensus on the advertised opinion.

Keywords: Sociophysics, phase transition, ageing, Monte Carlo simulation

Many recent models tried to simulate how in a society a consensus may emerge . Particularly well studied is the Sznajd model  where two neighbouring people (on a one-dimensional chain or higher-dimensional lattice) convince their neighbours if and only if these two people agree in their opinions. We generalize here the square lattice to a multilayer model, where each layer correponds to a different age group of the simulated people. We check if the phase transition known for the simple square lattice is also valid in the multilayer model and ask how many different initial opinions are allowed if a consensus is still desired. Finally, we check the effects of advertising  favouring the first of the several possible opinions.

Each site on an $`L\times L`$ square lattice has one of $`m`$ possible opinions, initially random. A pair of nearest neighbours, randomly selected in sequential updating, may share the same opinion $`j`$. In this case, all six neighbours of this pair are convinced of that opinion, that means they are forced to adopt the same value $`j`$. This model has a phase transition  as a function of the initial concentration $`p`$: If for $`m=2`$ the initial distribution of opinions is random and if the fraction $`p`$ of sites has the first and the remaining fraction $`1p`$ the other opinion, then for large lattices everybody ends up with the first opinion if $`p<1/2`$ and everybody with the other opinion if $`p>1/2`$. If $`m>2`$ we assume bounded confidence , which means that no neighbours can be convinced which differ by more than one opinion unit from the central pair. In that case  $`m=3`$ usually still leads to a consensus while $`m=4`$ seldomly allows a consensus: It is difficult to form a stable government coalition with four and more parties.

Now we generalize this planar Sznajd model to a multilayer model of $`N`$ square lattices on top of each other, representing $`N`$ age intervals in the life of an individual. With a certain probability $`q`$, called mortality, an individual dies, the younger ones on the same two-dimensional site all move up one time unit in their age, and a baby is born on the newly freed site, having the opinion of the parent. Within each layer only, and not between layers, the usual Sznajd convincing process takes place. First, we found that for two layers we still have a phase transition. Below we will check the chances of finding a consensus as a function of the number $`N`$ of layers in four cases: with an age-independent mortality (Fig.1), with a mortality increasing exponentially with age (Fig.2), with a probability to be convinced inversely proportional to mortality (Fig.3), and finally again with age-independent mortality and convincing probability but a constant advertising effort in favour of the first of $`m`$ opinions (Fig.4).

In this multilayer model we continue the simulations until the basic layer (the babies) reaches a complete consensus within that layer only. Then we find the majority opinion in each of the higher layers, and we declare a consensus if and only if the majority opinion in each of the higher layers agrees with the one surviving opinion in the basic layer.

With $`m3`$ possible opinions, a consensus is always formed for large lattices and long enough times, while for $`m4`$ a consensus is rare, just as for the single layer. We used $`L=13,17,23,31,61`$ and sometimes 101. Fig.1 shows for a $`31\times 31`$ lattice that for small mortalities $`q=0.001`$ the probability to reach a consensus, already below 7 percent for a single layer, decays exponentially towards zero with increasing number $`N`$ of layers. In contrast, for large mortalities, $`q=0.5`$, it stays at its single-layer value for all $`2N10`$. The higher mortality leads to a stronger coupling between layers and thus to easier overall consensus; note that no convincing takes place between different layers.

This trend of Fig.1 is well reproduced in Fig.2 if we include ageing, that means a mortality $`q`$ increasing exponentially with age $`n=1,2,\mathrm{},N`$ according to the Gompertz law of the 19th century: $`q=\gamma \mathrm{exp}[10(nN)/(N1)]`$. Not much changes if we include that young people are more easily convinced than old people and thus take a convincing probability $`\mathrm{exp}[10(Nn)/(N1)]=\gamma /q`$ in Fig.3.

Advertising has the effect that independent of the whole convincing process, a site assumes the first opinion with probability 0.1, at every iteration. The mortality is taken again as 0.01 independent of age, and convincing always takes place. For $`m=2`$ opinions, advertising always wins for large enough lattices, just as in a single layer ; also for $`m=3`$ advertising reaches a complete consensus for the advertised opinion. Fig.4 shows a more complicated behaviour for $`m=4`$: Now with increasing lattice size the fraction of advertising successes shrinks, possibly towards zero. It is quite independent of the number of layers which varies from 2 to 10. Fig.4 for small lattices like $`L=13`$ also shows that sometimes a consensus is reached but not for the advertised opinion.

In summary, the multilayer model is similar to the single-layer Sznajd model in allowing a consensus for two and three opinions but not for four (and more). This effect remains if advertising, a mortality increasing with age, and a convincing probability decreasing with age are introduced. The low probability of consensus, in the case of four opinions without advertising, decreases exponentially with the number of age layers.

Thanks are due to D. Stauffer for help.


text:
# 1 Introduction

## 1 Introduction

The Thomas precession is purely kinematical in origin . If a component of acceleration $`(\stackrel{}{a})`$ exists perpendicular to the velocity $`\stackrel{}{v}`$, for whatever reason, then there is a Thomas Precession, independent of other effects . When the acceleration is caused by a gravitational force field, the corresponding Thomas Precession is reasonably referred to as the Gravitational Thomas Precession (GTP). The question “Is there a gravitational Thomas Precession ?” was raised in  without an answer. However, once the physics involved in the Thomas Precession is accepted, the possibility of the existence of the GTP in planetary motion can not be ruled out in principle. The GTP was recently discussed in  wherein by introducing the GTP together with the notion of gravitational Lienard-Wiechart potentials (i.e., gravitational analogues of the electromagnetic Lienard-Wiechart potentials), the so called non-Newtonian perihelion shift of a planet (say Mercury) is seemingly explained, from a new angle, as a relativistic effect in flat-space time. Surprisingly the estimated value for the perihelion advance in this new approach happens to coincide with Einstein’s expression for the same effect, when a spin-orbit factor contributing insignificantly towards the observed effect, is reasonably dropped from the new expression. In the new approach the contribution of the GTP towards the perihelion advance came out as half of the Einstein’s predicted value and the other half seemingly originate from what we call the GLWP( gravitational Lienard-Wiechart potentials). The aim of this note is to clarify and confirm our earlier results by illustrating the GTP effect as nothing but a gravitomagnetic effect. To this end a hypothetical situation involving the motion of a comparatively low mass point particle (call it a planet) around a massive central object (call it the Sun) might be considered as under.

## 2  The GTP as a Gravitomagnetic effect

Let the massive Sun with its rest mass $`M_{}`$ be presumed fixed at the origin of some frame S, and let a relatively light planet with rest mass $`m_0`$ move under the force of gravitational interaction between them. In the frame S, the motion of the planet (treated as a point particle) is described by its position as a function of time $`\stackrel{}{r}(t)`$. At each instant $`t`$ the planet will have a velocity $`\stackrel{}{v}(t)=d\stackrel{}{r}/dt`$, an acceleration $`\stackrel{}{a}(t)=d\stackrel{}{v}/dt`$ as well as higher derivatives of position, $`\stackrel{}{b}(t)=d\stackrel{}{a}/dt`$, etc. Over a very small time interval $`t_1<t<t_2`$, the planet is thought to be moving under uniform acceleration. Then following Rohrlich  and Hill , the condition of uniform acceleration can be stated as :

$$\stackrel{}{b}+\mathrm{\hspace{0.17em}3}\gamma ^2(\stackrel{}{a}\stackrel{}{v})\stackrel{}{a}=\mathrm{\hspace{0.17em}0}$$
(1)
where

$$\gamma =(\mathrm{\hspace{0.17em}1}v^2/c^2)^{1/2}$$
(2)
The most general motion $`\stackrel{}{r}(t)`$ in S with uniform acceleration is given by the general solution of Eq. (1). To fond this solution we note first that Eq.(1) can be written as

$$d(\gamma ^3\stackrel{}{a})/dt=\mathrm{\hspace{0.17em}0}$$
(3)
and can therefore be integrated immediately, $`\gamma \mathrm{\hspace{0.17em}1}`$, to yield

$$\gamma ^3\stackrel{}{a}=\stackrel{}{g}$$
(4)
where $`\stackrel{}{g}`$ is a constant vector ( independent of $`t`$ )over a very small time interval. A simple calculation shows that Eq.(4) can be put into the following form

$$d(\gamma \stackrel{}{v})/dt\gamma (\gamma ^21)\stackrel{}{v}\times (\stackrel{}{v}\times \stackrel{}{a})/v^2=\stackrel{}{g}.$$
(5)
This equation can be written as an equation of motion of the planet, using

$$\stackrel{}{P}=m_0\gamma \stackrel{}{v},\stackrel{}{F}_g=m_0\stackrel{}{g}$$
(6)
where

$$\stackrel{}{g}=GM_{}\stackrel{}{r}/r^3$$
(7)
provided one introduces “a new force”

$$\stackrel{}{F}_{Th}=(\gamma +\mathrm{\hspace{0.17em}1})\stackrel{}{P}\times \omega _{Th}\text{where}$$
(8)
$$\stackrel{}{\omega }_{Th}=(\gamma \mathrm{\hspace{0.17em}1})\frac{\stackrel{}{v}\times \stackrel{}{a}}{v^2}$$
(9)
is the angular velocity associated with the Thomas precession. The force $`\stackrel{}{F}_{Th}`$ may be designated as the Thomas force. Thus the equation of the planet under uniform acceleration is

$$\stackrel{}{F}=\frac{d\stackrel{}{P}}{dt}=\stackrel{}{F}_g+\stackrel{}{F}_{Th}.$$
(10)
From Eqs.(2), (4), (7),(8) and (9) it can be shown that

$$\stackrel{}{F}_{Th}=m_0\stackrel{}{v}\times (\stackrel{}{v}\times \stackrel{}{g})/c^2.$$
(11)
so that the relativistic equation of motion of the planet becomes

$$\stackrel{}{F}=\frac{d\stackrel{}{P}}{dt}=m_0\stackrel{}{g}+m_0\stackrel{}{v}\times (\stackrel{}{v}\times \stackrel{}{g})/c^2=m_0\stackrel{}{g}+m_0\stackrel{}{v}\times \stackrel{}{B}_g$$
(12)
where

$$\stackrel{}{B}_g=(\stackrel{}{v}\times \stackrel{}{g})/c^2=\frac{GM_{}}{c^2r^3}(\stackrel{}{r}\times \stackrel{}{v})$$
(13)
is what we may call the gravitomagnetic induction field felt by the planet moving around the Sun in analogy with electromagnetism. It is to be noted that the origin of this gravitomagnetic field lies in the Thomas Precession in a gravitational field as per our derivation presented here. Now introducing the angular momentum of the planet as $`\stackrel{}{L}=m(\stackrel{}{r}\times \stackrel{}{v})`$, Eq.(13) can be rewritten as

$$\stackrel{}{B}_g=\frac{GM_{}}{mc^2r^3}\stackrel{}{L}$$
(14)
If the planet has an orbital gravitomagnetic moment defined by

$$\stackrel{}{\mu }_g=\left(\frac{m_0}{2m}\right)\stackrel{}{L},$$
(15)
then this would interact with the gravitomagnetic field (14) and the corresponding interaction energy according to Maxwellian Gravity  would be

$$U_{gm}=\stackrel{}{\mu }_g\stackrel{}{B}_g=\frac{GM_{}m_0}{2c^2r^3}\frac{L^2}{m^2}$$
(16)
But it can be shown that $`\frac{L^2}{m^2}=\frac{L_{0}^{}{}_{}{}^{2}}{m_{0}^{}{}_{}{}^{2}}`$, where $`\stackrel{}{L}_0=m_0(\stackrel{}{r}\times \stackrel{}{v})`$,( see the Appendix-I ), so that Eq.(16) can be reduced to

$$U_{gm}=\stackrel{}{\mu }_g\stackrel{}{B}_g=\frac{GM_{}}{2m_0c^2r^3}L_{0}^{}{}_{}{}^{2}$$
(17)
This is what represents the $`\stackrel{}{L}\stackrel{}{\omega }_{gT}`$ term in Eq.(29) of , where $`\stackrel{}{\omega }_{gT}`$ represents the gravitational Thomas Precession frequency corresponding to Eq.(24) of  with the non-relativistic assumption that $`\stackrel{}{L}\stackrel{}{L}_0`$. Then the rest of the arguments regarding the Thomas Precession’s contribution to the perihelion advance follows as given in . The result is that this contribution amounts to $`\mathrm{\hspace{0.17em}21.49}`$ arc-seconds per century for Mercury. This seems to be a clarified version of our earlier results in .The above discussion seems to look at the Gravitational Thomas Precession as a Gravitomagnetic effect in flat space-time.

Acknowledgments
The author is very much grateful to Prof. N. Barik and Prof. L. P. Singh, both of Utkal University, Vanivihar, Bhubaneswar and Dr. P. C. Naik, D. D. College, Keonjhar for their valuable suggestions for the improvement of this letter. The author also acknowledges the help received from Institute of Physics, Bhubaneswar for using its Library and Computer Centre for this work.
Appendix-I

Let us start with the force Eq.(12),viz.,

$$\stackrel{}{F}=\frac{d\stackrel{}{P}}{dt}=m_0\stackrel{}{g}+m_0\stackrel{}{v}\times (\stackrel{}{v}\times \stackrel{}{g})/c^2$$
(18)
This equation can be reduced to

$$\stackrel{}{F}=\frac{d\stackrel{}{P}}{dt}=m_0\stackrel{}{g}(\mathrm{\hspace{0.17em}1}v^2/c^2)+m_0\stackrel{}{v}(\stackrel{}{v}\stackrel{}{g})/c^2$$
(19)
Thus in relativistic central force Kepler problem we find the force is not strictly central , since it is evident from Eq.(19) that a component of the force exists in the direction of the velocity vector $`\stackrel{}{v}`$. Therefore the torque acting on the planet in the central force Kepler problem would be given by

$$\stackrel{}{N}=\frac{d\stackrel{}{L}}{dt}=\stackrel{}{r}\times \stackrel{}{F}=m_0(\stackrel{}{v}\stackrel{}{g})(\stackrel{}{r}\times \stackrel{}{v})/c^2\mathrm{\hspace{0.17em}0}$$
(20)
for non-radial motion. With $`\stackrel{}{L}=m(\stackrel{}{r}\times \stackrel{}{v})`$ and knowing that

$$c^2\frac{dm}{dt}=m_0(\stackrel{}{v}\stackrel{}{g})$$
(21)
Eq.(20) can be written as

$$\frac{d\stackrel{}{L}}{dt}=\frac{dm}{dt}\frac{\stackrel{}{L}}{m}$$
(22)
from which we can get

$$_{L_0}^L\frac{dL}{L}=_{m_0}^m\frac{dm}{m}.$$
(23)
This equation on integration yields

$$\mathrm{ln}\frac{L}{L_0}=\mathrm{ln}\frac{m}{m_0}$$
(24)
implying that $`\frac{L}{L_0}=\frac{m}{m_0}`$, or

$$\frac{L}{m}=\frac{L_0}{m_0}=constant$$
(25)
which is our required equation.