text:
# Differential rotation decay in the radiative envelopes of CP stars

## 1 Introduction

The present work was originally motivated by the existence of stars of spectral type A of which a subgroup is peculiar to their strong magnetic fields. These peculiar A and late B stars, which are collectively called CP stars, rotate significantly slower than their non-magnetic relatives. The dipole axes of the magnetic fields of CP stars have various orientations, with a tendency to larger tilts at faster rotation (Landstreet & Mathys 2000).

Young stars suffer from considerable rotational braking in particular during their pre-main sequence evolutionary phase. Angular momentum is lost partly through stellar winds in the T Tauri phase. Additional braking probably applies to T Tauri stars having an accretion disk (classical T Tauri stars – CTTS). These stars also form a slowly rotating subgroup among the T Tauri stars (Bouvier et al. 1993) just as the slow CP stars among ordinary A and late B stars. Magnetic fields exiting the accretion disk will couple to the stellar field and exert torques as well. The braking is efficient for solar-mass T Tauri stars as found by Cameron & Campbell (1993). Recently, Stȩpień (2000) proposed the analogy between the CP star subgroup and the CTTS among T Tauri stars and explains the slow CP stars by rotational locking with the disk in their pre-main-sequence phase.

If stars are considerably slowed down at the surface, how fast rotates their interior? The microscopic viscosity of stellar plasma is extremely small. We can estimate the time-scale of viscous decrease of the rotation in case of external braking simply by estimating the time-scale due to the gas’ viscosity of roughly $`\nu =10`$ cm2/s. For a stellar radius of $`R=10^{11}`$ cm, we get a viscous time-scale of $`R^2/\nu 10^{13}`$ yr which is 4–5 orders of magnitudes longer than the life-time of A stars and even a thousand times more than the age of the universe. Viscous decay of the differential rotation in any radiative shell of a star is thus not applicable. As long as there is no convection providing turbulent transport of angular momentum (whose efficiency may, however, not be very high), we should expect strong differential rotation on the way from the stellar surface to the deep interior.

A differential rotation in a radiative star will be prone to the magneto-rotational instability which requires only two things: an angular velocity decreasing with axis distance and a weak magnetic field (Balbus & Hawley 1991). The question what “weak” means will be addressed in Section 2 in the discussion of the initial magnetic field used for our simulations. The instability is known to evolve quickly on the time-scale of the rotation period. It has turned out to be an efficient generator of turbulence in accretion disks. See e.g. Kitchatinov & Rüdiger (1997) for a linear global analysis and Hawley (2000) and Arlt & Rüdiger (2001) for simulations. However, the instability is as well applicable to a stellar interior as long as the angular velocity decreases with axis distance in parts of the spherical domain. The mechanism is typically termed magneto-rotational instability or Balbus-Hawley instability. It is a consequence of the local and linear MHD equations; a lower limit to the magnetic field is only imposed by the magnetic diffusivity $`\eta `$ which is extremely small for stellar plasma. Field geometry is also almost irrelevant for the onset of the instability. The magneto-rotational instability must be quite ubiquitous in stellar radiative zones as soon as differential rotation emerges, likely to be caused by surface braking.

We would like to address the question of how long it takes to turn the radiative envelope of the star into a uniformly rotating shell. Is there any chance to maintain a rotation profile depending only on the axis distance? The answer on whether or not differential rotation may be expected, will have consequences for the geometry of magnetic fields or possibly for their generation in CP stars.

The questions were investigated numerically. In the second Section we will describe our setup for the simulations. The third and fourth Sections deal with the hydrodynamic and magnetohydrodynamic evolutions of a thick spherical shell. We will discuss the consequences of the computations in Section 5.

## 2 Simulations

The simulations apply the spectral, spherical MHD code of Hollerbach (2000). The computational domain covers a full spherical shell from the inner radius $`r_\mathrm{i}=0.2`$ to the outer radius $`r_\mathrm{o}=1.0`$. We start from the non-ideal MHD equations with the kinematic viscosity $`\nu `$ and the magnetic diffusivity $`\eta `$. The time-dependent, incompressible, non-dimensional equations are then

$`\left({\displaystyle \frac{}{t}}+u\right)u`$ $`=`$ $`p+\mathrm{Pm}^2u+(\times B)\times B`$ (1)
$`{\displaystyle \frac{B}{t}}`$ $`=`$ $`^2B+\times (u\times B)`$ (2)
with the usual meanings of $`u`$, $`B`$, and $`p`$ as the velocity, magnetic field, and pressure which is not explicitly calculated in this model, but eliminated by applying the curl-operator to Eq. (1). Lengths are normalized with the radius of the sphere, $`R`$, times are measured in diffusion times $`R^2/\eta `$, velocities are normalized with $`\eta /R`$ as well as magnetic fields with $`\sqrt{\mu \rho }\eta /R`$. Note that the permeability $`\mu `$ and the density $`\rho `$ are constants in our approach. This normalization leads to the magnetic Prandtl number,

$$\mathrm{Pm}=\frac{\nu }{\eta },$$
(3)
measuring the ratio of diffusive to viscous time-scales.

Instead of the physical $`u`$ and $`B`$, we integrate the potentials $`e`$, $`f`$, $`g`$, and $`h`$ which compose the physical quantities by

$`u`$ $`=`$ $`\times (e\widehat{r})+\times \times (f\widehat{r}),`$ (4)
$`B`$ $`=`$ $`\times (g\widehat{r})+\times \times (h\widehat{r}),`$ (5)
where $`\widehat{r}`$ is the unit vector in radial direction. The potentials for $`u`$ and $`B`$ are decomposed into Chebyshev polynomials for the radial dependence and into spherical harmonics for the angular dependence. The representation by potentials implies that $`u=0`$ and $`B=0`$ are always fulfilled automatically.

The initial conditions for the velocity represents a rotation profile in which the angular velocity decreases with the cylinder radius, $`s=r\mathrm{sin}\theta `$, according to

$$\mathrm{\Omega }=\frac{\mathrm{Rm}}{\sqrt{1+(2s)^{2q}}},$$
(6)
where Rm is the magnetic Reynolds number which is determined by the normalized angular velocity on the axis,

$$\mathrm{Rm}=\frac{R^2\mathrm{\Omega }_0}{\eta }.$$
(7)
This Rm will be varied in our simulations; we always put $`q=2`$. A profile depending on the axis distance appears to be a reasonable choice for the internal rotation of a star being prone to magnetic coupling with an accretion disk.

According to the Rayleigh criterion of hydrodynamic stability,

$$\frac{dj^2}{ds}>0,$$
(8)
where $`j=s^2\mathrm{\Omega }`$ is the angular momentum per unit mass, the above rotation profile (6) will provide us with a hydrodynamically stable configuration. For large axis distances $`s1`$ the profile with $`q=2`$ would be marginally stable, but within our shell of finite radius, the Rayleigh criterion is not violated.

The construction of the initial magnetic field is based on a vertical, homogeneous field, onto which we impose a non-axisymmetric perturbation of Fourier mode $`m=1`$. The total initial magnetic field can be written as

$$B=B_0[\widehat{z}+ϵ\mathrm{sin}(kz+\pi /4)\widehat{x}],$$
(9)
where $`\widehat{z}`$ is the unit vector in the direction of the rotation axis and $`\widehat{x}`$ is a unit vector in the equatorial plane. The wave number of the perturbation is $`k=4\pi `$. We added $`\pi /4`$ to the second term in (9) in order to provide mixed parity to the system. Equatorial and axial symmetry are thus broken allowing the system to develop flows and fields in all modes.

The initial magnetic field contains $`B_0=100`$ and $`ϵ=1`$ in most of our simulations. At $`\mathrm{Rm}=\mathrm{10\hspace{0.17em}000}`$, this configuration implies a magnetic energy which is two orders of magnitude smaller than the kinetic energy in the initial rotation (precisely $`E_{\mathrm{mag}}=0.017E_{\mathrm{kin}}`$). The magnetic energy is thus much smaller than the rotational energy as required for the onset of the magneto-rotational instability. In a real stellar environment, these quantities are roughly 30 orders of magnitudes apart.

We have also reduced the perturbation amplitude to $`ϵ=0.1`$ in order to estimate the sensitivity of the results to this parameter. As we will see in Section 4, the onset of the magneto-rotational instability is not affected by the smaller amplitude. The decay of differential rotation takes slightly longer with a factor of 1.6 but does not increase by an order of magnitude.

Are stars in the efficient regime of the magneto-rotational instability? For accretion disks rotating according to the Keplerian law where $`\mathrm{\Omega }s^{3/2}`$, the wavelength of the most unstable mode is

$$\lambda _{\mathrm{inst}}=2\pi \sqrt{16/15}u_\mathrm{A}/\mathrm{\Omega }$$
(10)
where $`u_\mathrm{A}=B/\sqrt{\mu \rho }`$ is the Alfvén velocity (Balbus & Hawley 1998). Let us assume that the rotation profile is partly Keplerian in a star and the rotation period be 1 day. A magnetic field of 100 kG will result – for densities between 1 and $`10^2`$ g/cm3 – in wavelengths of 10 to 100 km.

At this point, we can see that the magneto-rotational instability will set in much less promptly if $`\lambda _{\mathrm{inst}}>R`$, thus imposing an upper limit to the magnetic field. In principle, Eq. (10) reflects the aforementioned fact that the magnetic energy must be significantly smaller than the kinetic energy in order to be called a “weak field”. For stars, this is no problem as $`\lambda _{\mathrm{inst}}`$ and $`B`$ are proportional, and the limiting magnetic fields for which $`\lambda _{\mathrm{inst}}>R`$ are many orders of magnitudes larger than the fields assumed in the stellar interior.

There is no lower limit for $`B`$ in ideal MHD, but the non-vanishing diffusivity $`\eta `$ leads to a minimum magnetic field necessary for the magneto-rotational instability. Essentially, the growth rate of a perturbation with $`\lambda _{\mathrm{inst}}`$ must be larger than the decay rate of a structure of the same wavelength. The former is independent of $`\lambda `$, whereas the latter changes with $`\lambda ^2`$. Balbus & Hawley (1998) give estimates of the minimum field for the limits of $`\mathrm{Pm}1`$ and $`\mathrm{Pm}1`$. Both cases lead to $`B1`$ G in stellar interiors, assuming $`\eta =1000`$ cm2/s.

Even though magnetic fields in Ap stars are well in the suitable range for instability, we have to check the applicability of our initial $`B`$ for the numerical model which requires a much larger magnetic diffusion. This strong diffusion makes the instability window much narrower than it is in reality. Also the stellar $`\lambda _{\mathrm{inst}}`$ is invisible in the simulations due to limited resolution.

Using the angular velocity of about 6000 at $`s=0.6`$, the lower, diffusive limit is $`B=75`$, while the upper constraint from the size $`R`$ of the domain is $`B=920`$. In the computational setup with $`\mathrm{Rm}=\mathrm{10\hspace{0.17em}000}`$ and $`B_0=100`$, we obtain a wavelength of the most unstable mode of $`\lambda _{\mathrm{inst}}=0.1`$ near $`s=0.6`$. Under real conditions, the range of suitable fields spans many orders of magnitudes, since $`R`$ is huge and $`\eta `$ is extremely small.

The values of $`\mathrm{Pm}`$, $`\mathrm{Rm}`$, the initial magnetic field strength $`B_0`$ and the amplitude of the perturbation, $`ϵ`$, are the free parameters in the equations. Stellar gases possess magnetic Prandtl numbers of $`\mathrm{Pm}0.01`$. Values different from unity are typically difficult to achieve by numerical schemes. Values vastly different from unity mean that the time-scales for the diffusive processes in velocity and magnetic fields differ very much and are thus hard to cover appropriately by one simulation. We will vary Pm in our MHD simulations to evaluate the reasonableness of using Pm near unity.

The velocity and magnetic field in our computational domain are decomposed into 50 Chebyshev polynomials, 100 Legendre polynomials, and 30 Fourier modes. This resolution is sufficient to resolve the above wavelength of $`\lambda _{\mathrm{inst}}=0.1`$ in our numerical setup. The nonlinear terms $`(u)u+(\times B)\times B`$ and $`\times (u\times B)`$ are computed on a suitable number of collocation points in real space, and the spectral decomposition of these “right sides” are fed into the implicit time-stepping scheme of the linear part of Eqs. (1) and (2). First simulations with lower resolutions of $`(K,L,M)=(25,30,30)`$ and $`(K,L,M)=(50,50,30)`$ led to the same results within 3%.

The boundary condition for the flow is stress-free at the innermost and outermost radius. Vacuum conditions are imposed to the magnetic field at the inner and outer boundaries.

## 3 Hydrodynamic evolution

Before studying the magnetohydrodynamic case, we have to assess the evolution of the rotation flow without magnetic fields. This is of particular interest since the stress-free boundary conditions are not compatible with the initial azimuthal velocity profile $`u_\varphi (r\mathrm{sin}\theta )`$. The rotation profile will lead to meridional circulations which equalize the differential rotation on the viscous time-scale. We thus determine the purely hydrodynamic decay time before we can turn to magnetic configurations and their instabilities.

We measure a decay time with $`\mathrm{\Omega }(r,\varphi )`$ in the equatorial plane. The quantity is averaged over $`\varphi `$, providing a one-dimensional function of $`r`$. The equatorial plane was chosen just for simplicity. Function (6) is fitted to that profile varying Rm and $`q`$. We determine the time when $`q(t)=1`$ and call it the decay time $`\tau _{\mathrm{decay}}`$.

A model with $`\mathrm{Rm}=\mathrm{10\hspace{0.17em}000}`$ gradually decays and reaches $`q=1`$ after 65 rotation periods $`\tau _{\mathrm{rot}}`$. We will later see that this time is much longer than the decay time caused by the magneto-rotational instability. The kinetic energy in the meridional circulation is less than 1/2000 of the energy of the azimuthal velocities. At $`\mathrm{Rm}=\mathrm{50\hspace{0.17em}000}`$ the viscous decay lasts much longer than 200 rotation periods (Fig. 1). The ratio of azimuthal to meridional kinetic energies now exceeds $`10^4`$. Viscosities as low as stellar microscopic values are of course not achievable by numerical simulations. We can turn to magnetohydrodynamic computations knowing that the desired mechanism of differential-rotation decay should work on a time-scale of order 10 rotational periods or less to be essentially unaffected by viscosity.

## 4 MHD evolution

### 4.1 Simulation results

The following simulations including magnetic fields regard the enhanced decay rate of differential rotation as expected from the magneto-rotational instability. An initial magnetic field as described in Eq. (9) in Section 2 causes poloidal flows with vortices of the same size as the perturbation of the magnetic field. They form in all places, except where the gradient of the angular velocity is negligible, i.e. close to the rotation axis. A series of vertical cuts through the velocity field at four equidistant times is shown in Fig. 2. The projected velocity vectors indicate that the problem is numerically resolved. Note that only every second vector of our collocation point grid is plotted for the sake of visibility; at radii smaller than 0.45 only every third vector. The same series is shown for the magnetic field in Fig. 3. The vector lengths in both graphs are not comparable among the four slices; they are scaled for best visibility.

The first velocity snap-shot after 1.6 orbital revolutions shows the development of cells as the direct consequence of the Lorentz forces resulting from the magnetic-field perturbation. Two counter-rotating vortices represent one wave length of the perturbation from the second term in (9). Roughly four of these waves fit into the sphere according to the wavenumber $`k=4\pi `$. The second slice shows the emergence of relatively thin sheets of strong radial flows in up and down directions. These features become very prominent in the third figure after $`4.8\tau _{\mathrm{rot}}`$. They are actually quite extended over several tens of degrees in azimuthal direction. A detailed plot with the full resolution of our collocation point grid is shown in Fig. 4 magnifying a localized upstream. The fourth velocity slice of Fig. 2 shows an almost equalized rotation profile and a decay of small-scale features in the flow.

The latitudinal resolution of the model shown in Fig. 3 is plotted in Fig. 5 as a series of three Legendre spectra at $`t=0.001`$, 0.003, and 0.010, corresponding to $`1.6\tau _{\mathrm{rot}}`$, $`4.8\tau _{\mathrm{rot}}`$, and $`16\tau _{\mathrm{rot}}`$ resp. Maximum and minimum power span 2.5 orders of magnitudes in the most turbulent case (solid line). The initial power contrast is $`10^7`$ (dotted line); after the redistribution of angular momentum, the contrast quickly reaches the same order of magnitude again (dashed line in Fig. 5). The spectra of the velocity fields are very similar. The Fourier spectra as given in Fig. 6 show very satisfying power contrast all through the simulation.

The change of specific angular momentum as a function of axis distance is plotted in Fig. 7. The profiles of $`s^2\mathrm{\Omega }`$ are – for simplicity – again taken from the equatorial plane and averaged over the $`\varphi `$-direction. The initial differential rotation profile of (6) with $`\mathrm{Rm}=\mathrm{10\hspace{0.17em}000}`$ and $`q=2`$ is shown as a flat curve, while the steepest, nearly parabolic lines are the final distribution of angular momentum, corresponding to a nearly uniform rotation. We plotted dashed lines for $`t4.8\tau _{\mathrm{rot}}`$; a gap between these and the solid lines marks the transition period when strongest transport of angular momentum is found.

The redistribution of angular momentum is a combined result of stresses from velocity and magnetic field fluctuations. The averages $`u_r^{}u_\varphi ^{}`$ and $`B_r^{}B_\varphi ^{}`$, which are again taken in the equatorial plane only, show a clear domination of magnetic stresses over kinetic stresses, occasionally by a factor 10.

In the same way as in Section 3 we measure the decay time of the differential rotation by the time it takes the system to cross a $`q=1`$ profile. Despite the enormous flows emerging, the azimuthally averaged angular velocity provides us with profiles $`\mathrm{\Omega }(r)`$ in the equatorial plane which are fairly compatible with the two-parameter function (6). Fig. 8 shows examples of such fits for two velocity snap-shots of the simulation illustrated by snap-shots in Figs. 2 and 3. It is thus still reasonable to use the steepness $`q`$ for defining the decay time even for the simulations where magneto-rotational instability emerges.

The decay of $`q`$ with time is given by the solid line in Fig. 9. A short transition period between $`4.0\tau _{\mathrm{rot}}`$ and $`5.5\tau _{\mathrm{rot}}`$ can be seen. The third velocity slice of Fig. 2 with strong radial up and down flow sheets falls right in the middle of this period. Towards the end of the computation, the system oscillates around an equilibrium state with $`q=0`$, and magnetic and kinetic energies decay exponentially.

The transition period is less marked when we go to lower diffusivities, i.e. to higher magnetic Reynolds numbers Rm. The steepness $`q`$ diminishes more gradually, but still on a time-scale which is an order of magnitude shorter than viscous decay in purely hydrodynamic simulations. (Note, however, that the difference to viscous decay is expected to be much larger in a real star.) We have added the non-magnetic model with $`\mathrm{Rm}=\mathrm{10\hspace{0.17em}000}`$ as a dashed line in Fig. 9 illustrating the marginal influence of viscous decay on the MHD simulations.

### 4.2 Application to stellar parameters

The rotation profile seems to decay on the rotational time-scale. This is apparently way too fast for any trace of $`\mathrm{\Omega }(s)`$-rotation in stars with radiative envelopes. We will later see that this is not quite true. There are the real physical quantities which are hard to match in a computer simulation. The diffusive time-scale is orders of magnitudes longer than the rotational time-scale. We achieved to make them more than four orders of magnitudes different and may obtain an extrapolation towards real stellar parameters. The magnetic Reynolds number $`\mathrm{Rm}=R^2\mathrm{\Omega }/\eta `$ in stellar radiative zones is about $`10^{13}`$$`10^{14}`$. The highest Rm achieved numerically is 50 000 in this presentation. The dependence of the decay time $`\tau _{\mathrm{decay}}`$ on the magnetic Reynolds number and magnetic Prandtl number is shown in Fig. 11. The decay times are given in rotation periods which is a few days for CP stars. Fortunately, we found no significant dependence of the decay times on Pm. As the true magnetic Prandtl number will be of the order of $`10^2`$ for stars, we may assume that our Pm near unity will not imply severe differences from the real physics.

Also the amplitude $`ϵ`$ of the perturbation of the initial vertical magnetic field was changed to $`ϵ=0.1`$. The resulting decay of differential rotation for the model with $`\mathrm{Rm}=\mathrm{10\hspace{0.17em}000}`$ and $`\mathrm{Pm}=1`$ is shown in Fig. 10 in terms of the steepness $`q`$ of the rotation profile. The decay time increased from 5 to 8 rotation periods. Also the ratio of magnetic to kinetic stresses as the constituents of the transport of angular momentum is not changed and reaches 10:1.

With a series of computations with various Rm, we can make an attempt in extrapolating the decay time to stellar conditions. It is found that the decay time scales as the magnetic Reynolds number, in particular for $`\mathrm{Rm}\mathrm{10\hspace{0.17em}000}`$ which is the interesting interval. Measured in diffusion times, the decay time depends only slightly on Rm. We derive the relation

$$\frac{\tau _{\mathrm{decay}}}{\tau _{\mathrm{rot}}}=\frac{\mathrm{Rm}}{2000}$$
(11)
where $`\tau _{\mathrm{rot}}`$ is the rotation period of the star. If the stellar parameters $`R=2R_{}`$, $`\tau _{\mathrm{rot}}=1`$ day, and $`\eta =1000`$ cm2/s (Spitzer 1956) are applied, the magnetic Reynolds number is $`310^{14}`$. The diffusivity is the variable ingredient here; our diagram in Fig. 11 can thus be annotated with $`\eta `$ in stead of Rm (see upper abscissa). If Rm enters (11) with an exponent unity, the decay time is actually locked to the diffusion time and not to the rotation period. If we adopt a stellar $`\mathrm{Pm}=10^2`$ and rewrite Eq. 11 as $`\tau _{\mathrm{decay}}/\tau _{\mathrm{rot}}=0.0005\mathrm{Pm}\mathrm{Re}`$, we see that the differential-rotation decay as a consequence of the magneto-rotational instability is of $`O(10^5)`$ times faster than the viscous decay scaling with Re.

Relation (11) delivers a decay time of $`310^8`$ yr for the stellar parameters given above. This is of the order of the life-time of an A star. This extrapolation bears of course a wide uncertainty. If we look at the graph for $`\mathrm{Pm}=10`$ the power of Rm is slightly lower, $`\tau _{\mathrm{decay}}/\tau _{\mathrm{rot}}=0.002\mathrm{Rm}^{0.85}`$. The resulting decay time for $`\mathrm{Rm}=310^{14}`$ is $`\tau _{\mathrm{decay}}=10^7`$ yr.

## 5 Summary

Since differential rotation in radiative stellar zones cannot be damped by viscosity in a life-time of a star, we investigate the magnetic evolution which is likely to imply the onset of the magneto-rotational instability providing efficient angular-momentum transport. MHD simulations of spherical shells were performed showing that the instability emerges indeed, and quick equalization of differential rotation is found. An extrapolation to stellar parameters gives decay times of differential rotation of the order of 10–100 million years. This is the time-scale on which redistribution of angular momentum is taking place. The magneto-rotational instability grows on a scale of rotation periods, but the nett efficiency of angular momentum transport in the fully nonlinear regime is obviously a different one.

While solar-type stars have long enough life-times to rotate uniformly in their radiative cores, the life-time of stars of spectral type A is of the same order of magnitude as the decay time of about 100 million years. CP stars have magnetic fields providing the conditions for the magneto-rotational instability – along with an initial differential rotation caused by interactions with the accretion disk of the pre-main sequence life. Because of the comparable stellar life-time and decay time, differential rotation should thus be present in CP stars during a considerable period of their life. As long as the total angular momentum is conserved, we should expect a slight increase of surface angular momentum with age. Since an A star roughly doubles its radius during the presence on the main sequence, the resulting decrease of surface rotation may balance with the increase due to the mechanism described in this paper. Age data and precise rotation periods of CP stars will be needed to test this result.


format: text

text:
# Self-avoiding knots

## Abstract

Scaling arguments are used to analyze the size of topologically constrained closed ring polymer with excluded volume. It is found that there exists a finite range of polymer thickness (excluded volume) in which self-avoidance is unimportant and polymer swelling compared to the Gaussian size is entirely due to the topological constraints.

The ability to knot is the most obvious and the least understood property of long strings, such as ropes, threads, spaghetti, flux lines, or polymers. Knots in DNA have been directly observed Cozarelli ; Wang . The compelling, albeit circumstantial, proof of their importance can be seen in the fact that Nature had undertaken to design special cell machinery (called topological enzymes) spending energy to simplify DNA topology Rybenkov . Knots and other topological constraints are also of great importance in equilibrium and dynamics properties of polymer materials, such as melts, gels, networks, etc. Despite this ample motivation, current understanding of polymer topology is very limited.

Apart from mathematical theory of knot invariants Kaufman , the study of knots in polymers was dominated for a long time by computer simulations Maxim ; Muthukumar ; Deguchi ; Deutch . The incorporation of sophisticated topological invariants to compute entropy of the topologically constrained polymers remains a difficult challenge Nechaev . More recently, some insights into the statistical properties of knots were developed based on various scaling arguments tube\_inflation ; Stasiak\_book ; 3\_5\_uzla ; Kardar1 ; Kardar2 ; Kardar5 . As in other cases in soft condensed matter physics deGennes\_book ; deGennes\_Views , scaling is the theory of choice where more sophisticated approaches are too difficult.

The simplest question of polymer physics is that about the size (say, root-mean-squared gyration radius) of a single macromolecule in a dilute solution. In general, this question is not understood for the topologically constrained polymer ring. As a first step in this direction, it was recently shown 3\_5\_uzla  (confirming the earlier conjecture desCloizeaux ) that topological constraints in the knot effectively lead to self-avoidance on the large enough length scales. As a result, the size of a polymer ring with the topology of a trivial knot (an unknot) scales as $`R\mathrm{}N^\nu `$ even for the polymer with negligible excluded volume. Here, $`N`$ is the number of effective segments, $`\mathrm{}`$ is the segment length, and $`\nu 0.5883/5`$ is the well known critical exponent of self-avoiding walks deGennes\_book . This result has been also extended for the non-trivial knots using the aspect ratio of maximally inflated tube, $`p`$, as a topological invariant tube\_inflation . Overall, the result of the work 3\_5\_uzla  reads

$$R\{\begin{array}{ccc}\mathrm{}N^{1/2}p^{1/6}\hfill & \mathrm{when}& \hfill N<pN_{00}\\ \mathrm{}N^\nu p^{\nu +1/3}N_{00}^{\nu +1/2}\hfill & \mathrm{when}& \hfill N>pN_{00}\end{array}.$$
(1)
Here, $`N_{00}`$ is the characteristic length of knotting; it comes from the observation Muthukumar ; Deguchi  and mathematical proof theorem  that the probability to have an unknot in a phantom chain (which freely crosses itself) goes like $`e^{N/N_{00}}`$. Numerically, $`N_{00}300`$ Muthukumar ; Deguchi . Thus, the theory predicts that $`RN^{3/5}`$, if $`N`$ is large enough, independently of the knot type. This prediction was scrutinized using computer simulations ShimamuraDeguchiDA ; ShimamuraDeguchiNET ; StasiakNeScience ; the most recent work StasiakNeScience  claims complete agreement with the predicted scaling.

In this paper, our goal is to study the interplay of topological constraints and excluded volume interactions - two factors both of which are inevitably present in any real polymer. We imagine a polymer ring of $`N`$ effective segments of the length $`\mathrm{}`$ and diameter $`d`$ each. We would like to emphasize that end closing itself would not change the polymer size scaling, affecting only a numerical prefactor. This means, if we imagine a ring with excluded volume, but phantom (for which all topological classes of conformations are available, like in the presence of topological enzyme Sikorav ), then the usual prediction deGennes\_book  holds:

$$R_g\{\begin{array}{ccc}\mathrm{}N^{1/2}\hfill & \mathrm{when}\hfill & N<1/r^2\hfill \\ \mathrm{}N^\nu r^{2\nu 1}\hfill & \mathrm{when}\hfill & N>1/r^2\hfill \end{array}.$$
(2)
Thus, our goal if to generalize together the Eqs. (1) and (2). Our results are summarized by the Fig. 1, in which the diagram of scaling regimes is presented in terms of variables $`N`$ and $`r=d/\mathrm{}`$. The paper is organized around the discussion of this diagram.

As an important input of our analysis, we shall rely on the following insight gained by computer experiments Klenin ; Deguchi\_excluded  regarding the probability of knotting for the excluded volume polymer. As one might have guessed, excluded volume enhances the probability of a trivial knot, somewhat suppressing all non-trivial knots. What is much more difficult to guess is that the probability of a trivial knot remains exponential in $`N`$, with characteristic length very sharply depending on the excluded volume parameter. According to Klenin ; Deguchi\_excluded ,

$`p_{\mathrm{unknot}}(N,r)`$ $``$ $`\mathrm{exp}\left[N/N_0(r)\right],`$
$`N_0(r)`$ $``$ $`N_{00}\mathrm{exp}\left[\beta r^\upsilon \right],`$
$`N_{00}`$ $``$ $`300,\beta 30,\upsilon 1.`$ (3)
In fact, the authors of the work Deguchi\_excluded  found that according to their numerical data, the value of $`\upsilon `$ is slightly below unity, $`\upsilon 0.85`$. They also claimed that the coefficient $`A`$ me depend on the model. Neither of these details is important for our purposes here, we will only use the fact that the characteristic length of knotting, $`N_0(r)`$ is sharply increasing with $`r`$.

We are now in a position to start discussing the scaling diagram, Fig. 1. Let us begin with a trivial knot, or an unknot, in which case $`p`$ (in the scaling sense) is of order unity ($`p=2\pi `$).

First of all, there is obviously a regime when neither excluded volume nor topological constraints are important. This range is labeled as Gaussian in Fig. 1, because chain size in this regime scales as that of a Gaussian polymer, $`RN^{1/2}`$. This happens if the conditions $`N<1/r^2`$ (2) and $`N<N_0(r)`$ (3) are met simultaneously. In particular, the latter condition, according to the formula (3) ensures that even for the phantom chain the part of conformation space corresponding to all non-trivial knots is very small; therefore, set of conformations of phantom chain is practically the same as that of a trivial knot.

When we cross over the $`GS`$ boundary $`N1/r^2`$, excluded volume effect takes over, but topological constraints may still remain unimportant if $`N<N_0(r)`$. In this regime, labeled as dominated by self-avoidance in Fig. 1, chain size is given by $`R\mathrm{}N^\nu r^{2\nu 1}`$.

Another regime is obtained from Gaussian if we cross over the $`N=N_0(r)`$ line $`GT`$, remaining at $`N<1/r^2`$. In this case, we should imagine the chain as consisting of $`N/N_0(r)`$ blobs, of $`N_0(r)`$ monomers each. Inside each blob, as it length is smaller than the characteristic knotting length, the topological constraints are insignificant, as well as the excluded volume; therefore, the blob size is about $`R_{\mathrm{blob}}\mathrm{}N_0(r)^{1/2}`$. On the other hand, the chain of blobs is topologically constrained, and, therefore, according to 3\_5\_uzla , we get for the chain size $`RR_{\mathrm{blob}}\left(N/N_0(r)\right)^\nu \mathrm{}N^\nu N_0(r)^{\nu +1/2}`$.

At this stage, it is useful to check the self-consistency of our result. As we know tube\_inflation ; 3\_5\_uzla , the trivial knot size results from the balancing of two entropic factors. On the one hand, chain has entropic elasticity $`R^2/N\mathrm{}^2`$ which resists swelling. On the other hand, topologically constrained blobs resist interpenetration. At the equilibrium $`R`$, both entropic contributions are of the same order. We must now check that the second virial term corresponding to the excluded volume interactions is smaller than either of the above mentioned terms. Given that the second virial coefficient of elongated (rod-like) segments is about $`\mathrm{}^2d`$, the self-consistency condition reads

$$\mathrm{}^2d\frac{N^2}{R^3}<\frac{R^2}{\mathrm{}^2N}.$$
(4)
Note, that this condition is equivalent to $`R>R_{\mathrm{Flory}}`$, where $`R_{\mathrm{Flory}}\mathrm{}N^{3/5}r^{1/5}`$ is the Flory estimate of polymer size without topological constraints deGennes\_book . Thus, this condition means simply that excluded volume is unimportant if the chain swells more for topological reasons than it would due to the excluded volume. Simple algebra indicates that this condition is met if $`N_0(r)<1/r^2`$, which is equivalent to $`r<r^{}`$. Numerically, according to eq. (3), $`r^{}0.03`$.

Finally, let us see what happens if $`N_0(r)>1/r^2`$ and $`N>N_0(r)`$. Diagram Fig. 1 indicates that in this case we cross over the $`TS`$ boundary and enter again the same self-avoidance regime which we discussed previously n terms of crossing over the $`GS`$ line. Let us show why this is indeed the same regime. We should imgine once again that our chain consists of $`N/N_0(r)`$ blobs, with $`N_0(r)`$ monomers in each blob. In the present case, since $`N_0(r)>1/r^2`$, excluded volume is important for each blob, so that the blob size is about $`R_{\mathrm{blob}}\mathrm{}N_0(r)^\nu r^{2\nu 1}`$. Since the chain of blobs is topologically constrained, overall chain size scales as $`RR_{\mathrm{blob}}\left(N/N_0(r)\right)^\nu \mathrm{}N^\nu r^{2\nu 1}`$, which does not involve the topological quantity of $`N_0(r)`$ and coincides with our result for the self-avoidance dominated regime.

This completes the analysis for the case of an unknot (trivial knot).

Let us now consider some non-trivial knot, characterized by the length-to-diameter ratio of maximally inflated tube, also called ideal knot representation. The specific values of $`p`$ for many knots were computed by Pieranskii et al Pieranskii . Following the ideology of tube\_inflation ; 3\_5\_uzla , we imagine the polymer as confined in the maximally inflated tube such that it is unknotted inside the tube. This corresponds to viewing the polymer as consisting of $`p`$ roughly spherical blobs, with $`N/p`$ monomers in each blob. Furthermore, since the tube is maximally inflated, the blobs fill completely (in the scaling sense) the entire volume of the coil, which means $`R^3pR_{\mathrm{blob}}^3`$, or

$$Rp^{1/3}R_{\mathrm{blob}}.$$
(5)
Therefore, all we have to do is to address all possible regimes of the chain inside the blob.

The simplest case is when neither excluded volume nor topology are of importance inside the blob. This happens when both $`N/p<1/r^2`$, which is below the $`GS`$ line in Fig. 1, and $`N/p<N_0(r)`$, which is above or to the left of the $`GT`$ line. In this regime every blob is Gaussian, yielding $`R_{\mathrm{blob}}\mathrm{}(N/p)^{1/2}`$ and, according to formula (5),

$`R\mathrm{}N^{1/2}p^{1/6},\mathrm{when}N`$ $`<`$ $`p/r^2`$
$`\mathrm{and}N`$ $`<`$ $`pN_0(r).`$ (6)
This is Gaussian regime of Fig. 1.

When $`N/p`$ exceeds $`1/r^2`$, we cross over the $`GS`$ line into the regime dominated by self-avoidance. In this case, the blob size is about $`R_{\mathrm{blob}}\mathrm{}(N/p)^\nu r^{2\nu 1}`$. Therefore, formula (5) yields

$`R\mathrm{}N^\nu p^{\nu +1/3}r^{2\nu 1},\mathrm{when}N`$ $`>`$ $`p/r^2`$
$`\mathrm{and}r`$ $`>`$ $`r^{}.`$ (7)
Similar to the case of the trivial knot, this regime, which we call self-avoidance dominated, is valid all the way from the $`GS`$ boundary $`Np/r^2`$ to the $`TS`$ boundary $`r=r^{}`$.

Finally, if $`r`$ is small, but the blob length $`N/p`$ exceeds $`N_0(r)`$, then blobs become topologically constrained. This means, it becomes important that the chain inside the tube is unknotted, and so must be every blob. In this regime, blob size scales as $`R_{\mathrm{blob}}\mathrm{}(N/p)^\nu N_0(r)^{\nu +1/2}`$. Plugging this into formula (5) leads to

$`R\mathrm{}N^\nu p^{\nu +1/3}N_0(r)^{\nu +1/2},\mathrm{when}N`$ $`>`$ $`pN_0(r)`$
$`\mathrm{and}r`$ $`<`$ $`r^{}.`$ (8)
This is topology dominated regime.

It is instructive to summarize the results by looking at the ratio of the given knot size $`R`$ to the size of a phantom ring polymer (2), with the same length $`N`$, segment length $`\mathrm{}`$, and excluded volume $`r`$:

$$\xi =\frac{R_{\mathrm{knot}}(p)}{R_{\mathrm{phantom}}}.$$
(9)
This quantity is useful, because $`R_{\mathrm{phantom}}`$ can be also understood as the average size over all knot types, and such quantity can be easily extracted from computer simulations data ShimamuraDeguchiDA ; ShimamuraDeguchiNET . Various regimes for this quantity are outlined in the Fig. 2, where the cross-over lines of the scaling diagram Fig. 1 are superimposed with the line $`r1/\sqrt{N}`$ separating Gaussian and swollen regimes for the phantom polymer.

Somewhat counterintuitive result is that there is a parameter range (shaded in Fig. 2) in which topologically constrained knot is more expanded than phantom average over all knots. This happens because for the phantom chain (with annealed topological disorder) in this range more complex knots (with $`p`$ larger than the given value) are prevailing, they are more compact, and they dominate the average $`R_{\mathrm{phantom}}`$. This result is at least qualitatively in agreement with the simulation data ShimamuraDeguchiDA .

Metzler et al Kardar1 ; Kardar2  recently reexamined the possibility of ”knot segregation,” which was briefly discussed at the end of tube\_inflation  and later analyzed numerically KatrichStasiak\_segregation . The idea is that the entropy of a knotted polymer may be maximaized by the conformations in which knot is tightened within a short part of the chain, while the rest of the polymer fluctuates as a free unknotted loop. While the tendency to tighten the knot is quite strong for charged polymers Kardar5 , for neutral system it is governed by sub-linear in $`N`$ terms of free energy and, therefore, is very sensitive to the details of structure and interactions. In this paper, we assumed that topological constraints are delocalized along the polymer chain backbone. In the context of our work, this is justified for the following reasons. First, this issue is altogether irrelevant for the important case of the trivial knot, which is the first case considered in this paper. Second, for the non-trivial knots, it is necessary to understand the delocalized states in order to be able to push the question of knot segregation beyond the slip-link model of Kardar2 .

It should also be noted that in this paper we always assumed knots simple enough such that polymer chain is far from fully stretched state in the tube. The regime of such strong knotting requires special analysis Rabin .

To conclude, we presented the scaling theory determining the size of a polymer ring with quenched topology of a certain knot.

###### Acknowledgements.

The author acknowledges useful discussion and correspondence with M. Shimamura and T. Deguchi.


text:
# A NEW ANALYSIS OF THE 𝜋⁺𝜋⁻ AND 𝜋⁰𝜋⁰ DATAPresented by L. Leśniak at the Meson 2002, Seventh International Workshop on Production, Properties and Interaction of Mesons, Cracow, Poland, May 24-28, 2002

## Abstract

A joint analysis of the $`\pi ^+\pi ^{}`$ and $`\pi ^0\pi ^0`$ S-wave amplitudes including results of three experiments has been performed. Using our theoretical model of the $`\pi \pi `$ production in which not only the pion exchange but also the $`a_1`$ meson exchange were taken into account we have obtained a unique solution (called ”down-flat”) for the scalar-isoscalar phase shifts. Thus a long-standing problem of ambiguities in the S-wave isoscalar amplitude has been finally resolved and the second (”up-flat”) solution was eliminated.

PACS 13.60.Le, 13.75Lb, 14.40Cs

Analyses of experiments on the $`\pi \pi `$ production processes are important sources of information about the properties of the scalar mesons. Our knowledge of the basic parameters for the scalar-isoscalar mesons like $`f_0(500)`$, $`f_0(980)`$, $`f_0(1370)`$ or $`f_0(1500)`$ is still very limited. This makes tests of the QCD predictions about a possible existence of the scalar glueballs in the 1 to 2 GeV range to be very difficult.

The existence of ambiguities in the $`\pi \pi `$ S-wave isoscalar amplitude has been a serious problem since more than thirty years. Two experiments on the $`\pi ^{}p\pi ^+\pi ^{}n`$ reaction, both using the 17.2 GeV/c $`\pi ^{}`$ beam, have been performed at CERN already more than 20 years ago . The second experiment was done by the CERN-Cracow-Munich Collaboration on a polarized target. In 1997 we have shown  that the data obtained in those two high statistics experiments could be equally well described by at least two different solutions (called ”up–flat” and ”down–flat”) for the $`\pi \pi `$ S-wave phase shifts corresponding to the isospin 0 amplitude. Fortunately in 2001 new results from the E852 experiment performed at Brookhaven (USA) on the reaction $`\pi ^{}p\pi ^0\pi ^0n`$ at 18.3 GeV/c have been published . These $`\pi ^0\pi ^0`$ data combined with the previous $`\pi ^+\pi ^{}`$ results were used to eliminate one of the ambiguous solutions. Below we shall briefly describe a method to achieve this aim using our model of the $`\pi \pi `$ production . A more complete description of this analysis can be found in .

First of all we have merged two low momentum transfer squared $`|t|`$ bins measured by the E852 Group into a single one which was quite close to the $`(0.005÷0.20)`$GeV2/c2 bin used in the experiments performed at CERN. The $`m_{\pi \pi }`$ effective mass range chosen in this joint analysis was between 600 and 1460 MeV. Since the E852 collaboration has provided only the numbers of events and not the absolute cross sections we have determined the normalization comparing the intensities of the strongest $`D_0`$ wave saturated by the $`f_2(1270)`$ resonance. The data on the reaction $`\pi ^+p\pi ^+\pi ^+n`$ have been used to calculate the $`\pi \pi `$ scalar isospin amplitude.

At high energies of incoming pions and at small momentum transfers the above three reactions are dominated by the one pion exchange mechanism. However the pseudovector $`a_1(1260)`$ exchange, although smaller than the one pion exchange, is also present and cannot be neglected. We consider two independent $`S`$–wave transversity amplitudes $`g`$ and $`h`$ for the reaction $`\pi ^{}p_{}\pi ^+\pi ^{}n`$ which differ by a proton or neutron spin projection on the axis perpendicular to the production plane. By $`g_0`$ and $`h_0`$ we denote similar amplitudes of the reaction $`\pi ^{}p\pi ^0\pi ^0n`$. All these amplitudes have two components proportional to the off-shell $`T(\pi \pi \pi \pi )`$ and $`t(a_1\pi \pi \pi )`$ amplitudes. Both $`T`$ and $`t`$ amplitudes have two terms corresponding to isospin 0 or 2 of the pion pair. We use the isospin symmetry to relate different production amplitudes of the final $`\pi ^+\pi ^{}`$ or $`\pi ^0\pi ^0`$ pairs. New E852 results for the S-wave intensity $`I_0=|g_0|^2+|h_0|^2`$ provide us with only one additional constraint at each $`m_{\pi \pi }`$. This information combined with the knowledge of the amplitudes $`g`$ and $`h`$ from  is not sufficient to calculate four independent isospin amplitudes $`T_0`$, $`T_2`$ and $`t_0`$, $`t_2`$. To proceed further we have neglected the isotensor $`a_1`$ exchange amplitude $`t_2`$ which is presumably the smallest amplitude at low momentum transfers. Then one can directly relate the amplitudes $`g_0`$ to $`g`$ and $`h_0`$ to $`h`$, namely $`g_0=0.5C_\pi UT_2g`$ and $`h_0=0.5C_\pi U^{}T_2h`$, where $`C_\pi `$ and U are known functions . Thus using the amplitudes $`g`$ and $`h`$ found in the analysis of the $`\pi ^+\pi ^{}`$ data and knowing the amplitude $`T_2`$ we can calculate the $`I_0`$ intensity. The results in the critical $`m_{\pi \pi }`$ range below the $`K\overline{K}`$ threshold are shown in Fig. 1. One notices that it is the ”down-flat” solution which reasonably agrees with the E852 data. On the other hand the ”up-flat” points are too low, especially in the region of $`m_{\pi \pi }=(870÷970)`$ MeV.

In  the phases of the transversity amplitudes $`g`$ and $`h`$ have been obtained essentially from the assumed resonant $`\rho (770)`$, $`f_2(1270)`$ and $`\rho _3(1690)`$ phases. Now we can directly calculate one phase independently for each $`m_{\pi \pi }`$ using the new $`\pi ^0\pi ^0`$ results for $`I_0`$ in addition to the $`\pi ^+\pi ^{}`$ data. In this way one can adjust the phases of $`g`$ and $`h`$ to fit simultaneously the $`\pi ^+\pi ^{}`$ and the $`\pi ^0\pi ^0`$ cross sections. Then the $`\pi \pi `$ scalar-isoscalar phase shifts and inelasticities can be recalculated.

Again for the ”down-flat” solution to the $`\pi ^+\pi ^{}`$ data we have obtained a reasonable behaviour of the inelasticity parameter which, as one should expect, is close to unity below the $`K\overline{K}`$ threshold. However, the inelasticity corresponding to the ”up-flat” solution strongly deviates from unity. Also the ”up-flat” phase shifts have an unphysical minimum at 900 MeV. This is a clear signal that the ”up-flat” solution should be rejected. Since the ”down-flat” phase shifts and inelasticities behave regularly we accept the ”down-flat” solution, resulting from the joint analysis of the $`\pi ^+\pi ^{}`$ and $`\pi ^0\pi ^0`$data as a unique physical one.


text:
# Intrinsic Josephson effect and nonequilibrium soliton structures in two-gap superconductors

## Abstract

We predict a new dynamic state in current-carrying superconductors with multicomponent order parameter. If the current density J exceeds a critical value $`J_t`$, an interband breakdown caused by charge imbalance of nonequilibrium quasiparticles occurs. For $`J>J_t`$, the electric field penetrating from current leads gives rise to various static and dynamic soliton phase textures, and voltage oscillations similar to the nonstationary Josephson effect. We propose experiments to observe these effects which would probe the multicomponent nature of the superconducting order parameter.

\]

There are experimental evidences that several superconductors, including $`MgB_2`$ maz ; 2gap , $`NbSe_2`$ nbse2 , the heavy-fermion $`UPt_3`$ upt3 , the organic $`(TMTSF)_2X`$ gerome  and $`\kappa BEDT`$ org  compounds may have a multicomponent order parameter $`\psi `$ with internal degrees of freedom legget ; su . For instance, for two weakly coupled s-wave order parameters $`\psi _1=\mathrm{\Delta }_1e^{i\theta _1}`$ and $`\psi _2=\mathrm{\Delta }_2e^{i\theta _2}`$ corresponding to two different parts of the Fermi surface (like in $`MgB_2`$), the internal degree of freedom is the interband phase difference $`\theta (𝐫,t)=\theta _1\theta _2`$. In this case, besides the usual phase-locked states ($`\theta =0`$ or $`\pi `$), intrinsic phase textures $`\theta (𝐫,t)`$ can occur. The soft modes associated with the fluctuations of $`\theta (𝐫,t)`$ are nearly decoupled from the gap fluctuations and behave like the Anderson plasmons in Josephson junctions legget . These modes may manifest themselves as resonances in the ac Josephson effect agt , or static $`2\pi `$ phase kinks in $`\theta (x)`$ tanaka .

In this Letter we show that, although static metastable textures $`\theta (𝐫)`$ do not manifest themselves in the equilibrium magnetic response, the situation radically changes for nonequilibrium current states where the charge imbalance near normal leads gives rise to dynamic phase slip structures $`\theta (𝐫,t)`$ propagating into a superconductor (Fig. 1). We predict a new dynamic state above the critical current density $`J>J_t`$ which marks the onset of the current-induced breakdown of the superconducting state. For the multicomponent $`\psi `$, such breakdown is a two-stage process. First, at $`J=J_t`$, an interband phase breakdown occurs, resulting in spontaneous phase solitons in $`\theta (x,t)`$ and ac voltage oscillations. At higher currents, $`J=J_d`$, both gaps $`\mathrm{\Delta }_{1,2}`$ are suppressed by the pairbreaking effects. For weak interband coupling, $`J_t`$ is much smaller than the depairing current density $`J_d`$.

We obtain the equations of motion for $`\theta `$ and the electric field $`𝐄`$ near the critical temperature $`TT_c`$, using the time-dependent Ginzburg-Landau (TDGL) equations kopnin  generalized to a two-gap superconductor:

$$\mathrm{\Gamma }_\mu (_t\psi _\mu 2\pi ci\phi /\varphi _0)=\delta F/\delta \psi _\mu ,$$
(1)
where $`\mu `$ runs from 1 to 2, $`\phi `$ is the electric potential, $`\varphi _0`$ is the flux quantum, $`c`$ is the speed of light, $`\mathrm{\Gamma }_\mu `$ are damping constants, and the free energy $`F=d^3𝐫(f_1+f_2+f_m+f_{int})`$ contains the magnetic part $`f_m=|\times 𝐀|^\mathrm{𝟐}/\mathrm{𝟖}\pi `$, the GL intraband part $`f_\mu `$, and the interband interaction $`f_{int}`$

$`f_\mu =\alpha _\mu |\psi _\mu |^2+{\displaystyle \frac{\beta _\mu }{2}}|\psi _\mu |^4+g_\mu \left|\left(+{\displaystyle \frac{2\pi i}{\varphi _0}}𝐀\right)\psi _\mu \right|^2,`$ (2)
$`f_{int}=\gamma (\psi _1\psi _2^{}+\psi _1^{}\psi _2)=2\gamma \mathrm{\Delta }_1\mathrm{\Delta }_2\mathrm{cos}\theta .`$ (3)
Here $`𝐀`$ is the vector potential, $`\gamma `$ is a coupling constant, the asterisk means complex conjugation. The qualitative results of this work remain valid for any other periodic dependencies of $`f_{int}(\theta )`$ on $`\theta `$, for example, $`f_{int}=2\mathrm{\Delta }_1^2\mathrm{\Delta }_2^2(\gamma _1\mathrm{cos}2\theta +\gamma _2)`$, suggested for heavy-fermion superconductors upt3 . We focus here on two-gap superconductors for which two cases are possible: (i) The Cooper instability occurs at the same temperature in both bands, thus $`\alpha _\mu (T)TT_c`$, $`\gamma TT_c`$, while $`g_\mu `$, and $`\beta _\mu `$ are constants. (ii) The critical temperature $`T_{c1}`$ in the band 1 is higher than $`T_{c2}`$, in the band 2, so $`\psi _1`$ appears spontaneously below $`T_c`$, while $`\psi _2`$ at $`T_{c1}<T<T_c`$ is induced by the interband interaction, because $`\alpha _1TT_{c1}`$, $`\alpha _2TT_{c2}`$, and $`\gamma `$, $`g_\mu `$ and $`\beta _\mu `$ are constants suhl . In this paper we consider the first scenario, assuming that $`\gamma \alpha _{1,2}`$, which is likely to occur in $`MgB_2`$ maz ; 2gap .

Eqs. (1) yield the following gap equations

$$\mathrm{\Gamma }_\mu \dot{\mathrm{\Delta }}_\mu =g_\mu ^2\mathrm{\Delta }_\mu \stackrel{~}{\alpha }_\mu \mathrm{\Delta }_\mu \beta _\mu \mathrm{\Delta }_\mu ^3\gamma \mathrm{\Delta }_{\overline{\mu }}\mathrm{cos}\theta ,$$
(4)
where $`\overline{\mu }=2`$ if $`\mu =1`$, and $`\overline{\mu }=1`$ if $`\mu =2`$, $`\theta =\theta _1\theta _2`$, $`\stackrel{~}{\alpha }_\mu =\alpha _\mu +(2\pi Q_\mu /\varphi _0)^2g_\mu `$, and $`𝐐_\mu =𝐀\varphi _0\theta _\mu /2\pi `$. The imaginary part of Eq. (1) gives

$`\mathrm{\Gamma }_\mu \left(\dot{\theta }_\mu {\displaystyle \frac{2\pi c}{\varphi _0}}\phi \right)\mathrm{\Delta }_\mu ^2={\displaystyle \frac{2\pi g_\mu }{\varphi _0}}(\mathrm{\Delta }_\mu ^2𝐐_\mu )\pm \gamma \mathrm{\Delta }_1\mathrm{\Delta }_2\mathrm{sin}\theta `$ (5)
$`\times \times 𝐀=(4\pi /c)(\sigma 𝐄+𝐉_s)`$ (6)
with the plus sign corresponding to $`\mu =1`$. The current density $`𝐉=\sigma 𝐄+𝐉_s`$ has the ohmic component $`𝐉_n=\sigma 𝐄`$ proportional to the electric field $`𝐄=\phi \dot{𝐀}/c`$, and the normal conductivity $`\sigma `$. The supercurrent

$$𝐉_s=8\pi ^2c(g_1\mathrm{\Delta }_1^2𝐐_1+g_2\mathrm{\Delta }_2^2𝐐_2)/\varphi _0^2$$
(7)
is a sum of independent intraband contributions. The boundary conditions, $`(i_n+2\pi A_n/\varphi _0)\psi _\mu =i\psi _\mu /l_\mu `$, between a superconductor and a normal metal imply zero perpendicular supercurrents for both $`\psi _1`$ and $`\psi _2`$. Excluding $`𝐐_{1,2}=(\varphi _0𝐉_s\pm 4\pi cg_{2,1}\mathrm{\Delta }_{2,1}^2\theta )\varphi _0/[8\pi ^2c(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2)]`$ from Eqs. (2) and (7), we obtain the energy density $`f=f_\theta +f_e+_\mu [\alpha _\mu \mathrm{\Delta }_\mu ^2+g_\mu (\mathrm{\Delta }_\mu )^2+\beta _\mu \mathrm{\Delta }_\mu ^4/2]`$, where the magnetic and phase dependent parts, $`f_m`$ and $`f_\theta `$ turn out to be decoupled from each other:

$`f_m={\displaystyle \frac{H^2}{8\pi }}+{\displaystyle \frac{\varphi _0^2J_s^2}{(4\pi c)^2(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2)}}`$ (8)
$`f_\theta ={\displaystyle \frac{g_1g_2\mathrm{\Delta }_1^2\mathrm{\Delta }_2^2(\theta )^2}{(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2)}}+2\gamma \mathrm{\Delta }_1\mathrm{\Delta }_2\mathrm{cos}\theta `$ (9)
For constant gaps $`\mathrm{\Delta }_\mu `$, variation of Eq. (8) with respect to $`𝐇`$ yields the London equation $`\lambda ^2^2𝐇=𝐇`$ with the magnetic penetration depth $`\lambda ^2=\varphi _0^2/32\pi ^3(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2)`$. If $`\gamma \alpha _{1,2}`$, the gaps $`\mathrm{\Delta }_\mu `$ are decoupled from $`\theta `$, so variation of Eq. (9) with respect to $`\theta `$ yields the sine-Gordon equation $`L_\theta ^2^2\theta =\pm \mathrm{sin}\theta `$ which has a single-soliton solution tanaka , or periodic solutions analogous to the Josephson vortices in a magnetic field barone . However, the interband phase difference $`\theta `$ is decoupled from static magnetic fields, while the Gibbs free energy of a Josephson contact contains the interaction term $`(𝐇\theta _𝐉)`$. Thus, equilibrium nonuniform solutions $`\theta (x)`$ are energetically unfavorable as compared to the uniform phase-locked state $`\theta =0`$ for $`\gamma <0`$ or $`\theta =\pi `$ for $`\gamma >0`$. By contrast, various dynamic or quenched phase textures can be generated during current-induced interband breakdown.

We consider the weak coupling limit $`\gamma \alpha _\mu `$, for which the uniform gaps $`\mathrm{\Delta }_{1,2}`$ are unaffected by $`\theta (𝐫,t)`$, so Eqs. (5)-(6) give a time-dependent London equation with the account of dissipation and charge imbalance effects imb . Expressing $`𝐐_{1,2}`$ in Eqs. (5) in terms of $`𝐉`$ and $`\theta `$, and then subtracting the equations for $`\theta _1`$ and $`\theta _2`$ from each other, we obtain the following equation for $`\theta `$

$$\tau _\theta \dot{\theta }=L_\theta ^2^2\theta \pm \mathrm{sin}\theta +\alpha _\theta \text{div}𝐉_s,$$
(10)
where the relaxation time $`\tau _\theta `$, the decay length $`L_\theta `$, and the charge coupling parameter $`\alpha _\theta `$ are given by

$`\tau _\theta =\mathrm{\Gamma }_1\mathrm{\Gamma }_2\mathrm{\Delta }_1\mathrm{\Delta }_2/|\gamma |(\mathrm{\Gamma }_1\mathrm{\Delta }_1^2+\mathrm{\Gamma }_2\mathrm{\Delta }_2^2),`$ (11)
$`L_\theta ^2=g_1g_2\mathrm{\Delta }_1\mathrm{\Delta }_2/|\gamma |(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2),`$ (12)
$`\alpha _\theta ={\displaystyle \frac{\varphi _0\mathrm{\Delta }_1\mathrm{\Delta }_2(\mathrm{\Gamma }_1g_2\mathrm{\Gamma }_2g_1)}{4\pi c|\gamma |(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2)(\mathrm{\Gamma }_1\mathrm{\Delta }_1^2+\mathrm{\Gamma }_2\mathrm{\Delta }_2^2)}}`$ (13)
The signs in Eq. (10) correspond to the sign of $`\gamma `$. As seen from Eq. (10), the phase mode is unaffected by any distribution of bulk supercurrents, unless there is a nonequilibrium charge imbalance $`\text{div}𝐉_s=\sigma \text{div}𝐄`$ caused by the electric field penetrating from normal current leads. The excess charge density provides a driving term in the sine-Gordon equation for $`\theta `$, where the coupling constant $`\alpha _\theta `$ is nonzero if $`g_1\mathrm{\Gamma }_2g_2\mathrm{\Gamma }_1`$, which implies different electron diffusivities in the bands 1 and 2.

To obtain the equation for $`𝐄`$, we add Eqs. (5) for $`\theta _1`$ and $`\theta _2`$, then take the gradient of the sum and express $`𝐐_{1,2}`$ in terms of $`𝐉`$ and $`\theta `$. This yields

$$\tau _e\dot{𝐄}+𝐄L_e^2\text{grad}\text{div}𝐄+\alpha _e\dot{\theta }=\tau _e\dot{𝐉}/\sigma ,$$
(14)
where $`𝐉(t)`$ is the driving current density, $`L_e`$ is the electric field penetration depth, $`\tau _e`$ is the charging time constant, and the coupling term $`\alpha _e\dot{\theta }`$ describes an electric field caused by moving phase textures:

$`\tau _e=\sigma \varphi _0^2/8\pi ^2c^2(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2),`$ (15)
$`L_e^2=\sigma \varphi _0^2/8\pi ^2c^2(\mathrm{\Gamma }_1\mathrm{\Delta }_1^2+\mathrm{\Gamma }_2\mathrm{\Delta }_2^2),`$ (16)
$`\alpha _e=2|\gamma |\mathrm{\Delta }_1\mathrm{\Delta }_2\alpha _\theta `$ (17)
We first use Eqs. (10) and (14) to calculate the phase textures in a microbridge of length $`2a`$ (Fig. 1a) for which Eq. (10) at $`\gamma <0`$ takes the form

$$L_\theta ^2\theta ^{\prime \prime }\mathrm{sin}\theta +\alpha _\theta J_s^{}=0$$
(18)
Here $`\theta ^{}(x)`$ and $`J_s(x)`$ vanishes at the bridge edges, $`x=\pm a`$, where $`E(x,t)=E_0\mathrm{cosh}(x/L_e)/\mathrm{cosh}(a/L_e)`$ equals the value $`E_0=J/\sigma `$ provided by an external power source, $`\theta ^{}(0)=E^{}(0)=0`$, and the prime denotes differentiation over x. Now we obtain the condition under which Eq. (18) has a stable solution $`\theta (x,\theta _a)`$ localized at the edge of a long ($`aL_\theta `$) bridge, where $`\theta _a=\theta (a)`$. We first consider the limit $`L_\theta L_e`$, for which $`J_s^{}(x)=\sigma E^{}(x)`$ is essential only in a narrow region $`aL_e<x<a`$. Multiplying Eq. (18) by $`\theta ^{}`$ and integrating from 0 to $`aL_e`$, we obtain $`\theta ^{}(aL_e)=2\mathrm{sin}\theta _a/2`$. Here $`\theta ^{}(aL_e)=\alpha _\theta J/L_\theta ^2`$ is obtained by integrating Eq. (18) over $`aL_e<x<a`$, where $`\theta (xL_e)\theta _a`$, and $`\mathrm{sin}\theta `$ can be neglected. Excluding $`\theta ^{}(aL_e)`$, we arrive at the equation $`\alpha _\theta J=2L_\theta |\mathrm{sin}\theta _a/2|`$, which has solutions for $`\theta _a`$ only below the critical current density $`J_t=2L_\theta /\alpha _\theta `$:

$$J_t=\frac{8\pi c(\mathrm{\Gamma }_1\mathrm{\Delta }_1^2+\mathrm{\Gamma }_2\mathrm{\Delta }_2^2)\sqrt{|\gamma |g_1g_2(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2)}}{\varphi _0\sqrt{\mathrm{\Delta }_1\mathrm{\Delta }_2}|g_1\mathrm{\Gamma }_2g_2\mathrm{\Gamma }_1|}$$
(19)
In the opposite limit $`L_\theta L_e`$, the Laplacian in Eq. (10) can be neglected. Then the steady-state solutions $`\mathrm{sin}\theta (x)=\sigma \alpha _eE^{}(x)`$, and $`E(x)=E_0\mathrm{cosh}(x/L_e)/\mathrm{cosh}(a/L_e)`$, exist only if $`\sigma \alpha _eE^{}(a)<1`$, or $`J<J_t=L_e/\alpha _e\mathrm{tanh}(a/L_e)`$, where

$$J_t=\frac{|\gamma |[\sigma (\mathrm{\Gamma }_1\mathrm{\Delta }_1^2+\mathrm{\Gamma }_2\mathrm{\Delta }_2^2)]^{1/2}(g_1\mathrm{\Delta }_1^2+g_2\mathrm{\Delta }_2^2)}{\sqrt{2}\mathrm{\Delta }_1\mathrm{\Delta }_2|g_1\mathrm{\Gamma }_2g_2\mathrm{\Gamma }_1|\mathrm{tanh}(a/L_e)}$$
(20)
If $`\gamma (T)`$ and $`\alpha _\mu (T)`$ linearly pass through zero at $`T_c`$, Eqs. (19) and (20) yield $`J_t(T_cT)^{3/2}`$.

For $`J>J_t`$, the charge-induced interband breakdown gives rise to a striking dynamic phase texture in which solitons periodically appear near the current leads and then propagate to the bulk. The actual soliton dynamics depends on the particular boundary conditions, so we performed numerical simulations of Eqs. (10)-(17) for the different geometries shown in Fig. 1. We consider the limit $`\alpha _e\alpha _\theta \text{min}(\tau _\theta L_e^2,\tau _eL_\theta ^2)`$ of weak coupling between $`\theta `$ and $`E`$, for which the driving term divE in Eq. (10) is mostly determined by the static $`E(x)`$ near the current lead, while the last term in Eq. (17) gives only a small ac correction to $`E`$ of the second order in $`\alpha _e`$stat .

Fig. 2 shows the soliton formation near the current lead in the thin microbridge (Fig 1a) as J was instantaneously turned on from zero to a value below $`J_t`$. Such current step produces a phase soliton localized near the edge, while the bulk of the bridge remains in the phase-locked state with either $`\theta =0`$ or $`\pi `$, depending on the sign of $`\gamma `$. For $`J<J_t`$, this behavior is characteristic of all geometries shown in Fig. 1, whereas the dynamic phase texture for $`J>J_t`$ can be very different.

We start with the bridge geometry (Fig. 1a), for which $`E(x,t)`$ and $`\theta (x,t)`$ are even and odd functions of x, respectively, so the boundary conditions are: $`E(\pm a,t)=E_0`$, $`E^{}(0,t)=0`$, $`\theta (0)=0`$, and $`\theta ^{}(\pm a,t)=0`$. The latter condition ensures that supercurrents in both bands vanish at the normal electrodes, where $`J=\sigma E_0`$. In this case a soliton first appears at the bridge edge, but for $`J>J_t`$, it is pushed to the bulk by the strong gradient of $`E(x)`$. Then the next soliton forms near the edge and the process repeats periodically, resulting in the propagation of the soliton chain into the bulk, as shown in Fig. 3. After the first soliton reaches the center of the bridge, it stops (because $`\theta (0)=0`$), while new solitons keep entering the bridge from the current lead. During this soliton pileup, the mean slope $`\overline{\theta }^{}(t)`$ increases, reaching a critical value $`\overline{\theta }^{}`$ at which the soliton generation at the edge stops and a static phase texture forms. For $`JJ_c`$, the maximum soliton density $`n=\overline{\theta }_c^{}/2\pi `$ can be estimated from the static Eq. (10) in which the rapidly oscillating $`\mathrm{sin}\theta `$ can be neglected, and $`E(x)=J\mathrm{sinh}(x/L_e)/\sigma \mathrm{sinh}(a/L_e)`$. This yields $`\theta ^{}(x)=\alpha _\theta JL_e[\mathrm{cosh}(a/L_e)\mathrm{cosh}(x/L_e)]/L_\theta ^2\mathrm{sinh}(a/L_e)`$, giving for the long bridge $`(aL_e)`$:

$$\overline{\theta }_c^{}=JL_e\varphi _0(\mathrm{\Gamma }_1g_2\mathrm{\Gamma }_2g_1)/4\pi cg_1g_2(\mathrm{\Gamma }_1\mathrm{\Delta }_1^2+\mathrm{\Gamma }_2\mathrm{\Delta }_2^2)$$
(21)
During the formation of the soliton chain, $`t<t_c\tau _\theta a\overline{\theta }_c^{}/2\pi `$, voltage oscillations are generated on the bridge. A similar transient behavior occurs at the point contact (Fig. 1b), for which div$`𝐄=2I\mathrm{exp}(r/L_e)r_0/\pi L_e(r_0+L_e)r`$, where $`I`$ is the current through the semi-spherical contact of radius $`r_0`$. In this case the increase of $`I`$ above $`I_c`$ causes propagation of concentric soliton shells until the static structure forms as the critical phase gradient $`\theta _c^{}`$ is reached.

A very different soliton dynamics occurs in the 4-terminal geometry (Fig. 1c), for which the supercurrents make $`90^{}`$ turns around the central stagnation point $`(x=0)`$ where $`\theta =0`$. As a result, $`E(x)`$ and $`\theta (x,t)`$ become odd and even functions of x, respectively, which radically changes the dynamics of the phase textures as compared to the bridge, since now the solitons do not stop as they reach the center, but keep moving until they reach the opposite current lead. This is due to the fact that the driving charge density div$`𝐄`$ does not change sign all the way along the horizontal leg of the cross in Fig. 1c. The corresponding total extra charge (which plays the same role here as the transport current in a Josephson junction) is exactly compensated by the opposite charge distributed along the vertical leg. Thus, the mean phase gradient $`\overline{\theta }^{}`$ remains constant, while the phase slippage at the normal leads results in a continuous soliton shuttle and periodic voltage oscillations between the current lead and the center as shown in Fig. 4. For $`L_eL_\theta `$, the ac voltage $`V_\omega (t)=_{m=1}^{\mathrm{}}V_m\mathrm{cos}(m\omega t+\psi _m)`$ can be estimated by integrating Eq. (14) from $`0`$ to $`a`$, which yields $`\tau _e\dot{V}_\omega +V_\omega =\alpha _e[\dot{\theta }_a(t)\dot{\theta }(0,t)]`$. Here the values $`\theta _a(t)\theta (0,t)`$ obey the RSJ equation $`\tau _\theta \dot{\theta }_a+\mathrm{sin}\theta _a=\beta `$, for which $`\tau _\theta \dot{\theta }_a=\sqrt{\beta ^21}[1+2_{m=1}^{\mathrm{}}(\beta \sqrt{\beta ^21})^m\mathrm{cos}m\omega t]`$, $`\beta =J/J_t`$, and $`\omega =\tau _\theta ^1\sqrt{\beta ^21}`$ al . As a result, we obtain $`\mathrm{tan}\psi _m=m\omega \tau _\theta `$, and

$$V_m\frac{2\alpha _e}{\tau _\theta }\sqrt{\beta ^21}\frac{(\beta \sqrt{\beta ^21})^m}{\sqrt{1+(m\omega \tau _e)^2}}.$$
(22)
$`V_m(\beta )`$ is maximum at $`\beta =\beta _m`$, where $`\beta _m1+0.5(\tau _\theta /m\tau _e)^{2/3}`$ for $`\tau _\theta \tau _e`$. Thus, a two-gap superconducting bridge between normal banks exhibits the behavior of an overdamped Josephson junction. If the widths of the superconducting legs in Fig. 1c are different, either a unidirectional and mutually orthogonal soliton motion may occur, if $`J>J_t`$ in one or both legs, respectively.

The soft interband phase mode could also manifest itself in rf absorption at frequencies below the small gap $`\mathrm{\Delta }_2`$. Interaction of the $`\theta `$-mode with the rf electric field depends on the polarization of $`𝐄`$: if $`𝐄(t)`$ is parallel to the sample surface, then div$`𝐄=0`$, so the phase mode is not excited by the rf field. If, however the rf field has a component perpendicular to the sample surface, then the phase mode contributes to the rf impedance.

In conclusion, we predict an interband breakdown induced by nonequilibrium quasiparticles in two-gap superconductors. It results in spontaneous generation of static and/or dynamic phase textures, causing the Josephson voltage oscillations. The observation of these oscillations would unambiguously indicate the multicomponent nature of the superconducting order parameter.

This work was supported by the NSF MRSEC (DMR 9214707), AFOSR MURI (F49620-01-1-0464) (AG), and by the US DOE Office of Science under contract No. W31-109-ENG-38 (AG and VV).