text:
# Option Pricing Formulas based on a non-Gaussian Stock Price Model

## Abstract

Options are financial instruments that depend on the underlying stock. We explain their non-Gaussian fluctuations using the nonextensive thermodynamics parameter $`q`$. A generalized form of the Black-Scholes (B-S) partial differential equation, and some closed-form solutions are obtained. The standard B-S equation ($`q=1`$) which is used by economists to calculate option prices requires multiple values of the stock volatility (known as the volatility smile). Using $`q=1.5`$ which well models the empirical distribution of returns, we get a good description of option prices using a single volatility.

Although empirical stock price returns clearly do not follow the log-normal distribution, many of the most famous results of mathematical finance are based on that distribution. For example, Black and Scholes (B-S)  were able to derive the prices of options and other derivatives of the underlying stock based on such a model. An option is the right to buy or sell the underlying stock at some set price (called the strike) at some time in the future. While of great importance and widely used, such theoretical option prices do not quite match the observed ones. In particular, the B-S model underestimates the prices of options in situations when the stock price at the time of exercise is different from the strike. In order to match the observed market values, the B-S model would need to use a different value of the volatility for each value of the strike. Such “implied volatilities” of options of various strike prices form a convex function known as the “volatility smile”.

Indeed, attempts have been made to modify the B-S model in ways that can correct for the smile effect (cf  or more recently ). However, those approaches are often very complicated or rather ad-hoc, and do not result in managable closed form solutions, which is the forte of the B-S approach. In this paper we do however succeed in developing a theory of non-Gaussian option pricing which allows for closed form solutions for European options, which are such that can be exercised exclusively on a fixed day of expiration and not before (as is the case for American options).

Our approach uses stochastic processes with statistical feedback  as a model for stock prices. Such processes were recently developed within the Tsallis generalized thermostatistics . The driving noise can be interpreted as a generalized Wiener process governed by a Tsallis distribution of entropic index $`q`$. In the limit $`q1`$ the standard model is recovered. For $`q1.5`$, this model closely fits the empirically observed distribution for many financial time series, such as stock prices  (Figure 1), SP500 index,), FX rates, etc. This is consistent with a cumulutive distribution having power tails of index 3 . We derive closed form option pricing formulas, reproducing prices which, relative to the standard B-S model, exhibit volatility smiles very close to those observed empirically (Figure 4). Note that $`q=1.5`$ well models hydrodynamic turbulence on small scales , reinforcing notions of a possible analogy between these two systems.

The standard model for stock prices is $`S=S_0e^{Y(t)}`$ where $`Y(t)=\mathrm{ln}(S(t+t_0)/\mathrm{ln}S(t_0)`$ follows

$$dY=\mu dt+\sigma d\omega $$
(1)
The drift $`\mu `$ is the mean rate of return and $`\sigma ^2`$ is the variance of the stock logarithmic return. The noise $`\omega `$ is a Brownian motion defined with respect to a probability measure $`F`$. It is a Wiener process and satisfies $`E^F[d\omega (t)d\omega (t^{})]=dtdt^{}\delta (tt^{})`$ where the notation $`E^F[]`$ means the expectation value with respect to the measure $`F`$. This model yields a Gaussian distribution for $`Y`$ resulting in a log-normal distribution for $`S`$. Within this framework, Black and Scholes were able to establish a pricing model to obtain the fair value of options on the underlying stock $`S`$.

In this paper we assume that the log returns $`Y(t)=\mathrm{ln}S(t+t_0)/\mathrm{ln}S(t_0)`$ follow

$$dY=\mu dt+\sigma d\mathrm{\Omega }$$
(2)
with respect to the timescale $`t`$. Here $`\mathrm{\Omega }`$ evolves according to the statistical feedback process

$$d\mathrm{\Omega }=P(\mathrm{\Omega })^{\frac{1q}{2}}d\omega $$
(3)
The probability distribution $`P`$ satisfies the nonlinear Fokker-Planck equation

$$\frac{}{t}P(\mathrm{\Omega },t\mathrm{\Omega }^{},t^{})=\frac{1}{2}\frac{}{\mathrm{\Omega }^2}P^{2q}(\mathrm{\Omega },t\mathrm{\Omega }^{},t^{})$$
(4)
Explicit solutions for $`P`$ are given by Tsallis distributions

$$P_q(\mathrm{\Omega },t\mathrm{\Omega }(0),0)=\frac{1}{Z(t)}\{1\beta (t)(1q)[\mathrm{\Omega }(t)\mathrm{\Omega }(0)]^2\}^{\frac{1}{1q}}$$
(5)
Choosing $`\beta (t)=c^{\frac{1q}{3q}}((2q)(3q)t)^{2/(3q)}`$ and $`Z(t)=((2q)(3q)ct)^{\frac{1}{3q}}`$ ensures that the initial condition $`P_q=\delta (\mathrm{\Omega }(t)\mathrm{\Omega }(0))`$ is satisfied. The $`q`$-dependent constant $`c`$ is given by $`c=\beta Z^2`$ with $`Z=_{\mathrm{}}^{\mathrm{}}(1(1q)\beta \mathrm{\Omega }^2)^{\frac{1}{1q}}𝑑\mathrm{\Omega }`$ for any $`\beta `$. With $`\mathrm{\Omega }(0)=0`$, we obtain a generalized Wiener process, distributed according to a zero-mean Tsallis distribution In the limit $`q1`$ the standard theory Eq(1) is recovered, and $`P_q`$ becomes a Gaussian. We are concerned with the range $`1q<5/3`$ in which positive tails and finite variances are found . The distribution for $`\mathrm{ln}S`$ becomes

$$P_q(\mathrm{ln}S(t+t_0),t+t_0\mathrm{ln}S(t_0),t_0)=\frac{1}{Z(t)}\{1\stackrel{~}{\beta }(t)(1q)[\mathrm{ln}\frac{S(t+t_0)}{S(t_0)}\mu t]^2\}^{\frac{1}{1q}}$$
(6)
with $`\stackrel{~}{\beta }=\beta (t)/\sigma ^2`$. This implies that log-returns $`\mathrm{ln}[S(t+t_0)/S(t_0)]`$ over the timescale $`t`$ follow a Tsallis distribution, consistent with empirical evidence for several markets, e.g. S&P500 (Figure 1 ) , with $`q1.5`$.

Our model exhibits a feedback from the macroscopic level characterised by $`P`$, to the microscopic level characterised by $`\mathrm{\Omega }`$. We can imagine that this is really due to the interactions of many individual traders whose actions all contribute to shocks to the stock price which keep it in equilibrium. Their collective behaviour yields a nonhomogenous reaction to returns: rare events (i.e. extreme returns) will be accompanied by large reactions, and will tend to be followed by large returns in either direction.

Using Ito calculus , the equation for $`S`$ follows from Eq(2) as

$$dS=\stackrel{~}{\mu }Sdt+\sigma Sd\mathrm{\Omega }$$
(7)
where $`\stackrel{~}{\mu }=\mu +\frac{\sigma ^2}{2}P_q^{1q}.`$ The term $`\frac{\sigma ^2}{2}P_q^{1q}`$ is the noise induced drift. Remember that $`P_q`$ is a function of $`\mathrm{\Omega }`$ with

$$\mathrm{\Omega }(t)=\frac{\mathrm{ln}S(t)/\mathrm{ln}S(0)\mu t}{\sigma }$$
(8)
(with $`t_0=0`$ for simplicity.) As in the standard case (cf ), the noise term driving $`S`$ is the same as that driving the price $`f(S)`$ of a derivative of the underlying stock. It should be possible to invest one’s wealth in a portfolio of shares and derivatives such that the noise terms cancel each other, yielding a risk-free portfolio, the return on which must be the risk-free rate $`r`$. This results in a generalized B-S PDE

$$\frac{f}{t}+rS\frac{f}{S}+\frac{1}{2}\frac{^2f}{S^2}\sigma ^2S^2P_q^{1q}=rf$$
(9)
where $`P_q(\mathrm{\Omega }(t))`$ evolves according to Eq(4). For $`q1`$ the standard case is recovered. This PDE depends explicitly only on the risk-free rate and the variance, not on $`\mu `$, but it does depend implicitly on $`\mu `$ through $`P_q(\mathrm{\Omega })`$, with $`\mathrm{\Omega }`$ given by Eq(8). Therefore, to be consistent with risk-free pricing theory, we should first transform our original stochastic equation for $`S`$ into a martingale before we apply the above analysis. This will not affect our results other than that $`\stackrel{~}{\mu }`$ will be replaced by the risk-free rate $`r`$, ultimately eliminating the dependency on $`\mu `$. We now show how this is done.

The discounted stock price $`G=e^{rt}S`$ follows $`dG=(\stackrel{~}{\mu }r)Gdt+\sigma Gd\mathrm{\Omega }`$ where $`d\mathrm{\Omega }`$ follows Eq(3). For there to be no arbitrage opportunities, risk-free asset pricing theory requires that this process be a martingale, which it is not due to the drift term $`(\stackrel{~}{\mu }r)Gdt`$. One can however define an alternative driving noise $`z`$ associated with an equivalent probability measure $`Q`$ so that, with respect to the new noise measure, the discounted stock price has zero drift and is thereby a martingale. Explicitly,

$$dG=(\stackrel{~}{\mu }r)Gdt+\sigma GP^{\frac{1q}{2}}d\omega $$
(10)
Here, $`P`$ is a non-vanishing bounded function of $`\mathrm{\Omega }`$. With respect to the initial noise $`\omega `$, $`\mathrm{\Omega }`$ relates to $`S`$ via Eq(8). That is why for all means and purposes, $`P`$ in Eq(10) is simply a function of $`S`$ (or $`G`$), and the stochastic process can be seen as a standard state-dependent Brownian one. As a consequence, both the Girsanov theorem (which specifies the conditions under which we can transform from the measure $`F`$ to $`Q`$) and the Radon-Nikodym theorem (which relates the measure $`F`$ to $`Q`$) are valid, and we can formulate equivalent martingale measures much as in the standard case . We rewrite Eq(10) as

$`dG`$ $`=`$ $`\sigma GP^{\frac{1q}{2}}dz`$ (11)
where the new driving noise term $`z`$ is related to $`\omega `$ through

$$dz=\frac{(\stackrel{~}{\mu }r)}{\sigma P_q^{\frac{1q}{2}}}dt+d\omega $$
(12)
With respect to $`z`$, we thus obtain $`dG=\sigma Gd\mathrm{\Omega }`$ with $`d\mathrm{\Omega }=P_q^{\frac{1q}{2}}dz`$ which is non other than a zero-mean Tsallis distributed generalized Wiener process, completely analogous to the one defined in Eq(3). Transforming back to $`S`$ yields $`dS=rdt+\sigma Sd\mathrm{\Omega }.`$ Compared with Eq(7), the rate of return $`\stackrel{~}{\mu }`$ has been replaced with the risk-free rate $`r`$. This recovers the same result as in the standard asset pricing theory. Consequently, in the risk-free representation, Eq(8) becomes

$$\mathrm{\Omega }(t)=\frac{1}{\sigma }(\frac{\mathrm{ln}S(t)}{\mathrm{ln}S(0)}rt+\frac{\sigma ^2}{2}_0^tP_q^{1q}(\mathrm{\Omega }(s))𝑑s)$$
(13)
This eliminates the dependency on $`\mu `$ which we alluded to in the discussion of Eq(9). As discussed later on, by standardizing the distributions $`P_q(\mathrm{\Omega }(s))`$ we can explicitly solve for $`\mathrm{\Omega }(t)`$ as a function of $`S(t)`$ and $`r`$.

Suppose that we have a European option $`C`$ which depends on $`S(t)`$, whose price $`f`$ is given by its expectation value in a risk-free (martingale) world as $`f(C)=E^Q[e^{rT}C]`$. We assume the payoff on this option depends on the stock price at time $`T`$ so that $`C=h(S(T))`$. After stochastic integration of Eq(11) to obtain $`S(T)`$ we get

$$f=e^{rT}E^Q\left[h\left(S(0)\mathrm{exp}\left(_0^T\sigma P_q^{\frac{1q}{2}}𝑑z_s+_0^T(r\frac{\sigma ^2}{2}P_q^{1q})𝑑s\right)\right)\right]$$
(14)
The key point is that the random variable $`_0^TP_q^{\frac{1q}{2}}𝑑z_s=_0^T𝑑\mathrm{\Omega }(s)=\mathrm{\Omega }(T)`$ follows the Tsallis distribution Eq(5). This gives

$`f`$ $`=`$ $`{\displaystyle \frac{e^{rT}}{Z(T)}}{\displaystyle _R}h\left[S(0)\mathrm{exp}(\sigma \mathrm{\Omega }(T)+rT{\displaystyle \frac{\sigma ^2}{2}}\alpha T^{\frac{2}{3q}}+(1q)\alpha T^{\frac{2}{3q}}{\displaystyle \frac{\beta (T)}{2}}\sigma ^2\mathrm{\Omega }^2(T))\right]`$ (15)
$`(1\beta (T)(1q)\mathrm{\Omega }(T)^2)^{\frac{1}{1q}}d\mathrm{\Omega }_T`$
with $`\alpha =\frac{1}{2}(3q)((2q)(3q))c)^{\frac{q1}{3q}}.`$ We have utilized the fact that each of the distributions $`P(\mathrm{\Omega }(s))`$ occuring in the latter term of Eq(14) can be mapped onto the distribution of $`\mathrm{\Omega }(T)`$ at time $`T`$ via the appropriate variable transformations $`\mathrm{\Omega }(s)=\sqrt{\beta (T)/\beta (s)}\mathrm{\Omega }(T)`$. A major difference to the standard case is the $`\mathrm{\Omega }^2(T)`$-term which is a result of the noise induced drift. With $`q=1`$, the standard option price is recovered .

Eq(15) is valid for an arbitrary payoff $`h`$. We shall evaluate it explicitly for a European call option, which gives the holder the right to buy the stock $`S`$ at the strike price $`K`$, on the day of expiration $`T`$. The payoff is $`C=\mathrm{max}[S(T)K,0]`$. Only if $`S(T)>K`$ will the option have value at expiration $`T`$ (it will be in-the-money). The price $`c`$ of such an option becomes

$`c`$ $`=`$ $`E^Q[e^{rT}C]=E^Q[e^{rT}S(T)]_DE^Q[e^{rT}K]_D=J_1J_2`$ (16)
where the subscript $`D`$ stands for the set $`\{S(T)>K\}`$. This condition is met if $`\frac{\sigma ^2}{2}\alpha T^{\frac{2}{3q}}+(1q)\alpha T^{\frac{2}{3q}}\frac{\beta (T)}{2}\sigma ^2\mathrm{\Omega }^2+\sigma \mathrm{\Omega }+rT>\mathrm{ln}K/S(0),`$ which is satisfied for $`\mathrm{\Omega }`$ between the two roots $`s_1`$ and $`s_2`$ of the corresponding quadratic equation. This is a very different situation from the standard case, where the inequality is linear and the condition $`S(T)>K`$ is satisfied for all values of the random variable greater than a threshold. In our case, due to the noise induced drift, values of $`S(T)`$ in the risk-neutral world are not monotonically increasing as a function of the noise. As $`q1`$, the larger root $`s_2`$ goes toward $`\mathrm{}`$, recovering the standard case. But as $`q`$ increases, the tails of the noise distribution get larger, as does the noise induced drift which tends to pull the system back. As a result we obtain

$`J_1`$ $`=`$ $`S(0){\displaystyle \frac{1}{Z(T)}}{\displaystyle _{s_1}^{s_2}}\mathrm{exp}(\sigma \mathrm{\Omega }{\displaystyle \frac{\sigma ^2}{2}}\alpha T^{\frac{2}{3q}}(1q)\alpha T^{\frac{2}{3q}}{\displaystyle \frac{\beta (T)}{2}}\sigma ^2\mathrm{\Omega }^2)`$ (17)
$`(1(1q)\beta (T)\mathrm{\Omega }^2)^{\frac{1}{1q}}d\mathrm{\Omega }`$
$`J_2`$ $`=`$ $`e^{rT}K{\displaystyle \frac{1}{Z(T)}}{\displaystyle _{s_1}^{s_2}}(1(1q)\beta (T)\mathrm{\Omega }^2)^{\frac{1}{1q}}𝑑\mathrm{\Omega }`$ (18)
The equation Eq(16) with Eq(17) and Eq(18) constititutes a closed form expression for the price of a European call. We calculated option prices for different values of the index $`q`$, and studied their properties as a function of the relevant variables such as the current stock price $`S(0)`$, the strike price $`K`$, time to expiration $`T`$, the risk free rate $`r`$ and $`\sigma `$. The results obtained by our closed form pricing formula were confirmed both by implicitly solving the generalized B-S PDE Eq(9) and via Monte Carlo simulations of Eq(11). Note that American option prices can be solved numerically via Eq(9).

We compare results of the standard model ($`q=1`$) with those obtained for $`q=1.5`$, which fits well to real stock returns. Figure 2 shows the difference in call price. In Figure 3, the B-S implied volatilities (which make the $`q=1`$ model match the $`q=1.5`$ results) are plotted as a function of $`K`$. The assymetric smile shape, which is more pronounced for shorter times, reproduces well-known systematic features of the “volatility smile” that appears when using the standard $`q=1`$ model to price real options. In Figure 4, the volatility smiles for actual traded options on BP and S&P 500 futures is shown together with those resulting from our model using $`q=1.5`$. These results are encouraging, and we are currently studying a larger sample of options data. Empirical work is required to see if arbitrage opportunities can be uncovered that do not appear when the standard model is used. Another potential application will be with respect to option replication and hedging.

Acknowledgements: Fruitful discussions with Roberto Osorio and Jeremy Evnine are gratefully acknowledged.


format: text

text:
# 1 Introduction

## 1 Introduction

A great deal of effort has been expended on studying compactifications of string theory, that is, geometrical backgrounds in which four spacetime dimensions are assumed to be large while the others are unobservably small. In order to develop a complete and viable physical model, though, we must not only describe the nature of the compactification, but also understand the dynamical processes which produced it. While a compactification can be, and usually is, imposed by hand, doing so goes against the grain of the superstring program.

We are already familiar with a number of fundamental symmetries which are broken as the universe evolves. A common and reasonable assumption that the division between large and small spatial dimensions is one of the asymmetries arising naturally during the course of cosmological evolution. For example, in the very early universe all of the dimensions may have been compact with small radius and some dynamical process may have subsequently driven significant expansion in only three of the spatial dimensions. Some time ago, Brandenberger and Vafa  proposed just such a mechanism for a universe in which all spatial dimensions were toroidal. Here, wound strings impede cosmological expansion if they fall out of thermal equilibrium. Since two-dimensional string worldsheets generically fail to intersect in more than four spacetime dimensions, the interactions required for thermal equilibrium are not robust in more than four spacetime dimensions, giving a qualitative explanation of why spatial regions with more than three dimensions do not expand. This is an attractive proposal, since it uses the simple and intrinsically stringy winding modes. Moreover, aspects of the model were verified  using a numerical algorithm originally devised by Smith and Vilenkin  and this line of argument was recently extended to include a spectrum of branes .

Unfortunately, toroidal compactifications do not lead to realistic particle physics models whereas toroidal orbifolds or Calabi-Yau compactifications of string theory (or perhaps $`G_2`$ compactifications of M-theory) have a better chance of yielding the particle phenomenology we observe. Obviously, it is worth rethinking the Brandenberger-Vafa mechanism in these contexts. In this paper we take a preliminary step in this direction.

A crucial ingredient of the Brandenberger-Vafa mechanism is that tori have nontrivial cycles in every dimension, and each of these cycles is a member of a family that moves throughout the torus. This is not true of more general compactifications. In particular, many toroidal orbifolds and Calabi-Yau spaces have trivial fundamental groups and thus no notion of wound strings. Without incontractible 1-cycles it would seem that one cannot expect an analog of the Brandenberger-Vafa mechanism to hold in these spacetimes. However, these spaces still support “pseudo-wound” strings – contractible strings that wrap the full spatial extent in a given direction. If these strings contract in time scales that are less than the cosmological Hubble time, they will be cosmologically irrelevant and cannot support a modified version of the Brandenberger-Vafa scenario. However, if their contraction time scale is significantly longer than the Hubble time this raises the possibility that pseudo-wound modes can play the same role as their stable, truly wound cousins.

As an illustrative example, consider a compactification involving an eliptically fibered $`K_3`$ manifold. $`K_3`$ has a trivial fundamental group and yet in the neighborhood of a non-degenerate fiber the space looks locally like $`U\times T^2`$ with $`U`$ an open set in $`CP^1`$. A string that winds around a cycle of the $`T^2`$ is not topologically stable, but in order to decay it must move over the base to one of 24 singular fibers on which it can shrink. If the $`CP^1`$ base is large or asymmetric in shape (e.g. the cigar shape discussed in ), it is conceivable that the growth of the toroidal fibers might be inhibited by the pseudo-wound string modes over cosmologically interesting timescales.

Similarly, even though toroidal orbifolds generally have trivial fundamental groups, string configurations that would be topologically wound on the covering space must find their way to a fixed point of the group action in order to unwind. This process introduces a time scale and here we study the lifetime of such pseudo-wound strings on orbifolds. We employ the same numerical techniques used in toroidal spacetimes  but adapt them for use on orbifolds. The $`_4`$ symmetry of the basic 2D lattice used in  allows us to construct several orbifolds which respect the usual supersymmetry conditions , and we consider:

$$T^4/_2,T^4/_4,T^6/_2,T^6/_4,T^6/_4^{},T^6/(_2\times _2)$$

where $`_4^{}`$ acts differently than $`_4`$, as explained below. For purposes of comparison, we also consider the “factored,” non-supersymmetric orbifolds

$$T^2/_2,(T^2/_2)^2,(T^2/_2)^3,T^2/_4,(T^2/_4)^2,(T^2/_4)^3.$$

In this way we study the dimensional and topological dependence of various scenarios.

Our numerical simulations assume that the background spacetime is not evolving, so while our results shed light on the stability and persistence of string networks on toroidal orbifolds, they do not directly describe the dynamical evolution of these spaces. Moreover, the Brandenberger-Vafa scenario requires that the dilaton is evolving freely , so the spacetime dynamics cannot be governed by general relativity. If the dilaton is fixed, the low energy limit of string theory reduces to general relativity and string networks actually enhance the expansion rate since they have negative pressure. However, as shown in Appendix A of Ref. , string networks do not induce rapid expansion in dilaton gravity, and the subspaces where the strings do not annihilate will expand more slowly than those in which they do. The mechanism that fixes the dilaton is not well understood, but it is reasonable to assume that it is related to the physics which clamps the values of the moduli that set the size of the internal space, which must happen if the compactification is to be stable. This leads to a consistent picture: in addition to a running dilaton, the Brandenberger-Vafa scenario requires that the moduli be free to evolve if it is to induce the anisotropy between the three large directions we observe and the remaining compact directions, and both the moduli and dilaton are fixed after the Brandenberger-Vafa mechanism has had time to work.

In Section 2 we review the numerical techniques used in the simulation, including the representation of the orbifolds we study. In Section 3 we present the results of the simulations and discuss their significance. In Section 4 we go on to discuss some general features of string and brane gas models in the context of Calabi-Yau compactifications, and we identify possible extensions to this work.

## 2 Numerical simulation

### 2.1 The Smith-Vilenkin Lattice

Our simulations use the lattice model constructed by Smith and Vilenkin . The worldsheet gauge is chosen, so we must satisfy the conditions

$`𝐗^{}\dot{𝐗}`$ $`=`$ $`0,`$ (2.1)
$`𝐗{}_{}{}^{}{}_{}{}^{2}+\dot{𝐗}^2`$ $`=`$ $`1,`$ (2.2)
as well as the equation of motion

$$\ddot{𝐗}𝐗^{\prime \prime }=0,$$
(2.3)
where primes and overdots denote derivatives with respect to the worldsheet coordinates $`\sigma `$ and $`\tau `$.

By discretizing the worldsheet coordinates as in Figure 1, Eq. (2.1) and Eq. (2.2) become

$`𝐮^{}\dot{𝐯}`$ $`=`$ $`0,`$ (2.4)
$`𝐮{}_{}{}^{}{}_{}{}^{2}+\dot{𝐯}^2`$ $`=`$ $`1`$ (2.5)
where we have defined

$`𝐮(\sigma ,\tau )`$ $``$ $`[𝐗(\sigma +\delta ,\tau )𝐗(\sigma \delta ,\tau )]/2\delta ,`$ (2.6)
$`𝐯(\sigma ,\tau )`$ $``$ $`\{[𝐗(\sigma ,\tau +\delta ){\displaystyle \frac{1}{2}}[𝐗(\sigma +\delta ,\tau )+𝐗(\sigma \delta ,\tau )]\}/\delta .`$ (2.7)
The solutions to these equations are

$`𝐗(\sigma ,\tau +\delta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}[𝐗(\sigma +\delta ,\tau )+𝐗(\sigma \delta ,\tau )]+𝐯(\sigma ,\tau )\delta ,`$ (2.8)
$`𝐯(\sigma ,\tau +\delta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}[𝐯(\sigma +\delta ,\tau )+𝐯(\sigma \delta ,\tau )]`$
$`+`$ $`[𝐗(\sigma +2\delta ,\tau )2𝐗(\sigma ,\tau )+𝐗(\sigma 2\delta ,\tau )]/4\delta .`$
These are gauge preserving and become exact solutions to the differential equations Eq. (2.1) and Eq. (2.2) as $`\delta 0`$.

It is essential that the initial conditions respect the gauge. This is ensured by also discretizing the target space and allowing only three types of links, as shown in Figure 2 (our target space grid has been rotated from that of ). The first type has $`|𝐮|=1`$ and $`𝐯=0`$. Its end points are separated by $`|\mathrm{\Delta }𝐗|=2\delta `$, which implies that two components of $`\mathrm{\Delta }𝐗`$ are $`\pm \sqrt{2}\delta `$ and the others are zero. Such links are “fully stretched” and are at rest. The second type has $`|𝐮|=|𝐯|=1/\sqrt{2}`$, implying that one component of $`\mathrm{\Delta }𝐗`$ is $`\pm \sqrt{2}\delta `$ and the others are zero, and one component of $`𝐯`$ is $`\pm \frac{1}{\sqrt{2}}`$ and the others are zero. Such partially contracted links move at 0.707 light speed perpendicular to themselves. The third type has $`𝐮=0`$ and $`|𝐯|=1`$. Its end points are degenerate, with $`\mathrm{\Delta }𝐗=0`$. Two components of $`𝐯`$ are $`\pm \frac{1}{\sqrt{2}}`$ and the others are zero. Such links are fully contracted and move at light speed diagonal to the coordinate axes. A 2D lattice is shown here but the generalization to higher dimensions is obvious.

String interactions are achieved via the intercommutation of two strings which pass through the same lattice site in the target space, as shown in Figure 3. This preserves the gauge constraint exactly, as well as energy, momentum and angular momentum. We specify the coupling constant by assigning a probability $`P`$ that strings in such a configuration will in fact intercommute.

### 2.2 Simulation of strings on orbifolds

In Figure 4 we show how $`G=_2`$ acts on the manifold $`M=T^2`$ to produce the pillow-shaped orbifold $`\mathrm{\Gamma }=T^2/_2`$ with fixed points $`O`$, $`P`$, $`Q`$, and $`R`$. Construction of $`T^2/_4`$ is similar, producing a space with three fixed points. Higher-dimensional orbifolds are of course more difficult to visualize, but are mathematically straightforward. We construct $`\mathrm{\Gamma }=T^4/_2`$, $`T^4/_4`$, $`T^6/_2`$, and $`T^6/_4`$ by letting $`z_i,`$ $`i=1,\mathrm{},D/2`$, be the variables on a complex torus defined by $`z_iz_i+1z_i+i`$. The group acts on these as follows:

$$(z_1,z_2)(e^{i\theta }z_1,e^{i\theta }z_2)$$

with $`\theta =\pi `$ for $`G=_2`$ and $`\theta =\pi /2`$ for $`G=_4`$. For $`D=6`$, $`z_3`$ is left invariant under the group action. We construct $`T^6/_4^{}`$ by acting as

$$(z_1,z_2,z_3)(e^{i\pi /2}z_1,e^{i\pi /2}z_2,e^{i\pi }z_3)$$

and finally $`T^6/(_2\times _2)`$ by

$`(z_1,z_2,z_3)`$ $``$ $`(e^{i\pi }z_1,e^{i\pi }z_2,z_3)`$
$`(z_1,z_2,z_3)`$ $``$ $`(z_1,e^{i\pi }z_2,e^{i\pi }z_3).`$
We refer to these as the “supersymmetric orbifolds”. For the “factored” cases $`(T^2/_2)^{D/2}`$ and $`(T^2/_4)^{D/2}`$ each $`z_i`$ is acted on by a separate copy of $`G`$.

We simulate string evolution on these spaces as follows. When identifying interactions, the group acts on the position of each string link, producing an orbit. Identical orbits (rather than identical positions) then give rise to intercommutation, with the group action information included in the string’s link for future string evolution.

## 3 Simulations and Results

We begin by checking our simulation algorithm against previous codes, and then perform two sets of computations. In the first we keep the intercommutation probability at $`P=1`$, and vary the topology. In the second, we focus on topologies which had potential for thermalization in more than 3+1 dimensions, and study the impact of varying $`P`$.

### 3.1 Critical density computation

Previous studies have calculated the critical density in string thermodynamics, or the maximum allowed short string density . We first determine this density as a computational check. We define a “long” string in the same way as previous analyses: a string is “long” if its extent in at least one dimension is approximately the length of that dimension. Using a string network consisting of 100,000 string links and units such that energy is measured in links and length in $`\delta `$, we obtain the following critical densities:

| $`D=3:`$ | $`0.1715\pm 0.0083`$ |
| --- | --- |
| $`D=4:`$ | $`0.0560\pm 0.0009`$ |
| $`D=5:`$ | $`0.0269\pm 0.0003`$ |

These agree to within a few percent of the results found in the literature . While the difference exceeds the uncertainty in the fit, the results also depend weakly on the overall lattice size, or the specific details of the interaction code when more than two strings intersect at a point, and these systematic effects are most likely the origin of the discrepency.

### 3.2 Unwinding rate at constant coupling

Our first investigation compares the rates of string unwinding in different topologies in the “strong coupling” regime of $`P=1`$. This allows the most efficient interactions, and hence the greatest possible change to some initially wound state. We chose the dimensions so that each space has volume of $`(200\sqrt{2}\delta )^D`$, and a total of 40,000 wound string links. For example, on $`T^2`$ we choose $`L=200\sqrt{2}\delta `$ and 100 wound strings, whereas on $`T^6/_4`$ we choose $`L=252\sqrt{2}\delta `$ and 80 wound strings. As shown in Figure 5, the strings are initially wrapped at random around the dimensions, with half of the string links being fully extended in a zig-zag pattern and half the links being fully kinetic with velocity directions chosen at random. This allows us to construct initial configurations with random velocities, so that the strings have a fair chance of moving around and interacting. This kind of winding configuration was one of several proposed in , and suited our purposes best. The string density was kept low so that complications arising from the Hagedorn transition (see ) were avoided and there was no obstacle to all strings becoming small.

To determine if the result was density-dependent, we then performed simulations where the lattice length was reduced by factors of 3 and 7, keeping the number of wound strings the same (note that the number of links in a wound string is proportional to the length of the dimension). This significantly increases the density: when the 6D lengths are reduced by a factor of 3, the density changes by a factor of $`3^5=243`$.

The simulations were run for 20 Hubble times444We use the term ‘Hubble time’ to denote the time required for light to propagate across the size of the universe. This usage is slightly loose since we are not working with an expanding background., which in the lowest density simulation, was 4000 time steps. In Figure 6 we plot the ratio of small string energy density $`\rho _{small}`$ to the total energy density $`\rho _{total}`$ for each configuration. Since all strings began large, this ratio measures the efficiency with which the strings interact to form small strings and the speed with which this happens determines whether the “pseudo-wound” strings can have an impact on the overall cosmological evolution. For these simulations, we could not use the “extent” criteria of Refs , as orbifolds do not allow simple definitions of length in each dimension. Instead, our criterion for a string to be “long” is that the total string length is at least $`L/2`$. This is not as descriptive as measuring the string’s extent in each dimension, but since thermodynamics creates a large gap between “small” and “large” strings we believe this serves our purpose .

Our initial configuration, with half of the links purely kinetic, allows half of the string to be “slack”. Thus it can lose up to half of its links and still be wound around the torus or wrap around the largest cycle of the orbifold. Since we expect that strings can shed this slack quickly, the real indication of unwinding is how quickly the short string density rises after reaching $`\rho _{small}/\rho _{total}0.5`$.

The toroidal cases behaved as expected, in agreement with Ref. . With $`D=2`$ the long strings immediately fragmented to small strings, whereas with $`D=4`$ and $`D=6`$ the long strings persisted, as expected from naive phase-space arguments based on the dimensionality of strings.

The results for the orbifold simulations are similar. In all supersymmetric cases we studied, the unwinding rate was nearly identical to the toroidal case. This apparent indifference to topology is presumably due to the high codimension of the fixed point sets, which makes it difficult for the pseudo-wound strings to find a path that results in their unwinding. In both of the $`4`$ and $`6`$ dimensional orbifolds we study, the singularities are codimension 4, providing a significant phase space to support long pseudo-wound strings.

To test this idea, we simulated examples where the topology can be factored into products of 2-dimensional orbifolds: $`(T^2/_2)^{D/2}`$ and $`(T^2/_4)^{D/2}`$, yielding lower codimension fixed point sets. As explained previously, these were constructed analagously to the supersymmetric orbifolds, except that each complex torus now has its own group action. In these examples the string can unwind within the individual 2D orbifolds around which it was initially wound. Intuitively, we expect these models to produce different results from the analogous toroidal cases, except for $`D=2`$ where the strings will overlap regardless of topology. We emphasize that these spaces are not realistic candidates for compactification.

The $`T^2/_2`$ case is very similar to the torus, immediately achieving $`\rho _{small}/\rho _{total}0.5`$ and thus unwinding quickly. For other $`(T^2/_2)^{D/2}`$ spaces, the unwinding is more efficient than for the torus, but the overall dynamics do not differ significantly. The $`(T^2/_4)^{D/2}`$ cases do differ, however. For $`D=2`$, the story is identical to $`T^2`$. In $`D=4`$, however, the behavior is nearly identical to that of $`T^4`$ for only about 1 Hubble time. After this, the strings continue to contract, taking only $``$ 3 Hubble times to achieve $`\rho _{small}/\rho _{total}0.5`$ and eventually producing completely short strings. In a cosmological setting, this rapid unwinding would not be compatible with the Brandenberger-Vafa scheme. The larger configuration space associated with the $`(T^2/_4)^3`$ prevents rapid string unwinding, although string contraction is more efficient than that on $`T^6`$ (see Figure 7).

The intuitive picture of the decay of long strings as a function of dimensionality is built around two extremes. A low dimensionality has high string density and thus allows quick but not necessarily complete fragmentation, due to the Hagedorn transition favoring the creation of long (and thus likely wound) strings. A high dimensionality has slow (if any) fragmentation, though once fragmentation takes place there is no obstacle to all strings remaining small. How the system behaves in between these extremes depends on the topology. The 2-dimensional simulations were nearly identical (immediate, but not complete, unwinding) as were the 6-dimensional cases (slow or nonexistent unwinding). All 4-dimensional cases begin with moderate unwinding. After approximately one Hubble time the $`T^4`$ and $`(T^2/_2)^2`$ fragment very slowly, whereas $`(T^2/_4)^2`$ continues to fragment moderately until all strings are small. This is because during the first Hubble time fragmentation is only possible due to mutual interaction. After this, once a string has had enough time to contract around a fixed point due to its tension, it can self-interact. This difference is explained by the topology: $`T^2/_4`$ demands that strings wound in the covering space self-intersect, whereas the others do not.

### 3.3 Varied probability of intercommutation

Of the topologies examined, $`(T^2/_4)^{D/2}`$ has the greatest potential for string unwinding, possibly allowing thermalization in more than $`3+1`$ dimensions. Even though this space is not a realistic candidate for compactification, it allows us to perform one other useful test. The previous calculations were performed at an intercommutation probability of 1, but perhaps for sufficiently low interaction rates the strings would not interact enough to unwind. To explore this, we repeated the simulation on $`(T^2/_4)^2`$ for up to 5 Hubble times with $`0P1`$, effectively varying the value of the string dilaton, which determines the intercommutation probability. As shown in Figure 8, for small coupling the unwinding rate sharply rises with $`P`$ until $`P0.4`$, beyond which it stabilizes.

## 4 Branes and General Compactifications

Our simulations show that pseudo-wound strings on supersymmetry preserving orbifold compactifications have long lifetimes. Consequently, they may play an important cosmological role by inducing anisotropic cosmological expansion similar to that seen with toroidal compactifications. So far we have focused solely on strings, but brane wrappings of higher dimensional cycles can be important cosmologically, as emphasized in Ref. . However, Ref.  restricts its analysis to toroidal spaces and it is worth briefly considering such ideas in more general compactifications.

In Refs , the authors employ the moduli space approximation, studying the cosmological evolution of the radial moduli of the assumed spatial nine-torus topology. Although this approximation has limited validity in cosmological settings, we also use it to get a qualitative understanding of the dynamics when branes are wrapped on more complex compact spaces. Imagine that we compactify nine spatial dimensions and allow the moduli of the compactification to be time dependent. For example, assume that the compactification involves the simplest non-toroidal Calabi-Yau space, namely $`K_3\times T^5`$ in type IIA string theory.555Here we are specifying the asymmetric topology by hand and asking if there is a natural mechanism that yields an asymmetric geometry. Just as in , the higher dimensional branes with $`p=8,6,4`$ generically intersect in nine space dimensions, and quickly thermalize, leaving 2-branes and 1-branes. We expect wrapped 2-branes to dominate over wrapped 1-branes in the energy-momentum tensor, so 2-branes will initially dominate the cosmology.

The holomorphic 2-cycles in this product space are the 2-cycles on the $`K_3`$ and the 2-cycles on the $`T^5`$, since $`K_3`$ has no 1-cycles to hook up with the 1-cycles on the $`T^5`$. Naive dimension counting then suggests that thermal fluctuations in the early universe cause some 5-dimensional submanifold to grow large and that efficient 2-brane/anti-2-brane annihilations in this submanifold allow the geometrical asymmetry to grow. A natural scenario, independently suggested in , is that the $`T^5`$ fluctuates to a large size, leaving the $`K_3`$ small, and continues to grow ever larger while the 2-branes wrapping the 2-cycles in the $`K_3`$ keep it small. However, this scenario is not unique. Return to the example in the Introduction, and imagine that the $`K_3`$ is eliptically fibered over a $`CP^1`$ base and that the initial thermal fluctuation involves $`CP^1\times T^3`$. Wrapped 2-branes efficiently annihilate in this subspace allowing it to continue to grow while 2-branes and 1-branes wrapping the smaller dimensions do not efficiently annihilate and hence impede their growth. Within the 5-manifold $`CP^1\times T^3`$ we then expect a thermal fluctuation to drive a larger 3-dimensional subspace within which wrapped strings can efficiently annihilate. Since $`CP^1`$ has no 1-cycles, the preferred choice would be a large $`CP^1\times S^1`$ so that wrapped strings keep the remaining $`T^2`$ small. But if the $`T^3`$ should fluctuate large first, it is hard to see what would subsequently stabilize the $`CP^1`$ at small size.

If, instead, we chose the initial topology to be a Calabi-Yau three-fold $`M`$ times, say, a $`T^3`$, how would things change? At first sight, one might imagine that the story is essentially the same. Some 5-dimensional manifold is driven large by a thermal fluctuation, and within this space a 3-dimensional subspace grows larger still, while wound 2-branes and 1-branes keep all other dimensions small. However, the nontrivial holomorphic cycles in a generic Calabi-Yau space $`M`$ are isolated. This means that the annihilation of the branes wrapping these cycles is enormously more efficient than in the toroidal cases. For example, if a thermal fluctuation should drive $`M`$ to grow momentarily larger than the $`T^3`$, dimension counting shows that 2-brane annihilation in 6-dimensional space is not efficient and thus leads us to believe that further expansion would be impeded. However, if the cycles are isolated, the branes and anti-branes which we assume to be evenly distributed with no net winding number, immediately annihilate, allowing the space to continue growing.

A potential way out of this conclusion is to note that although the holomorphic homology generators are generically isolated, sufficiently high multiples can be part of continuous families that sweep through the Calabi-Yau manifold. If branes can wrap these cycles as independent states, their annihilation would once again be inefficient, constricting expansion. Even so, a variety of geometrical evolutions might follow. For example, imagine that $`M`$ is $`K_3`$ fibered and that the $`K_3`$ fibers are themselves elliptically fibered (over a twisted product of two copies of $`CP^1`$). Various combinations of, say, the $`CP^1`$s in the base and the $`T^3`$ or its submanifolds could receive the thermal fluctuation driving a large 5-dimensional submanifold while wrapped branes keep all other dimensions small. For instance, if $`CP^1\times T^3`$ gets large we are in the same situation as in the last paragraph.

In some of the above toy scenarios, the geometrical route we illustrated involved including part of the Calabi-Yau space or the $`K_3`$ into the large spacetime manifold, partially erasing the motivation for the initial choice of vacuum topology. If the Calabi-Yau space is multiply connected or if, using the results of the last section, it is simply connected but has a sufficiently robust spectrum of pseudo-wound strings, then we can imagine a scenario in which the Calabi-Yau stays small. Namely, in $`M\times T^3`$ if a two-dimensional subspace of $`M`$ together with the $`T^3`$ should initially become larger than the other dimensions and if the $`T^3`$ then receives an additional thermal fluctuation to large size while wound or pseudo-wound strings keep the somewhat asymmetric Calabi-Yau space small, we would evolve into the more familiar situation of a small Calabi-Yau space for the extra dimensions. Moreover, in the moduli space approximation, some Calabi-Yau spaces simply do not have the necessary flexibility for all independent subspaces to get large (e.g. examples with one Kahler moduli) and hence we would then expect a thermal fluctuation to make either the $`T^3`$ large from the outset, or the whole of $`M`$ large from the outset. If wound or pseudo-wound strings can impede the growth of $`M`$, this would favor $`T^3`$ getting large, returning us to the usual scenario of a small Calabi-Yau space.

## 5 Conclusion

We have studied the unwinding time scale for pseudo-wound strings on toroidal orbifolds to decay to short string states and found that in cases relevant for supersymmetric string compactifications the time scales can be many times the Hubble time. This suggests that such strings may play an important cosmological role, a prospect that is worthy of further study. Ideally, one would tackle this problem in a dynamical spacetime, fully incorporating the general relativistic backreaction of the strings on the geometry. This would then permit us to take up the question of geometrical homogeneity within the compact spatial dimensions. Additionally, while some work has been done on extending the Brandenberger-Vafa mechanism from string gases to brane gases, there is much more to be done. As we briefly mentioned, Calabi-Yau manifolds have numerous non-trivial cycles around which branes can wrap but since these cycles are often isolated, the subsequent dynamics can be quite different from what one expects on the torus and it would be interesting to determine the resulting cosmological implications. Finally, while this scenario can produce an initial anisotropy between the large and small directions, it does not guarantee that the compactification is stable at late times, and this problem needs to be addressed separately.

## Acknowledgments

We are grateful to R. Brandenberger, D. Easson and A. Vilenkin for useful discussions. The work of R.E. is supported by the Columbia University Academic Quality Fund and the Department of Energy. The work of B.R.G. is supported in part by the DOE grant DE-FG02-92ER40699B. M.G.J. is supported by a GAANN Fellowship from the United States Department of Education and by a Pfister Fellowship. ISCAP gratefully acknowledges the generous support of Ohrstrom Foundation.


text:
# Influence of the passive region on Zero Field Steps for window Josephson junctions

## I  Introduction

The present design of low temperature superconducting devices based on semi-conductor technology enables to integrate Josephson junctions into superconducting strip-lines to realize complex devices. A simple example is the window design shown in Figure 1 where a single rectangular junction is surrounded by a uniform passive region. In other systems like integrated receivers as0 (August 2001), the local flux flow oscillator is coupled directly to a strip-line containing a few junctions which will realize the mixing of the local signal with the external signal coming from an antenna.

The electrodynamic behavior of such window junctions is given by the phase difference between the top and bottom superconducting layers. It obeys an inhomogeneous sine-Gordon equation which is equivalent to a homogeneous sine-Gordon equation for the junction domain coupled to a wave equation for the passive region.
In the static case, a small passive region causes a rescaling of the Josephson length $`\lambda _J`$ into $`\lambda _{\mathrm{eff}}>\lambda _J`$ Caputo et al. (1996a) and this renormalization was seen in the experiments by the group of Ustinov Franz et al. (2001). On the other hand a large passive region causes the destruction of the kink due to the finite length of the junction Caputo et al. (1994).

In the dynamical case two limiting systems can be considered, first the passive region can be present only along the longitudinal direction of the junction. Then one can assume a transverse profile that propagates in the long direction. Using these ideas, Lee et al derived the dispersion relation for linear superconducting strip-lines Lee (1991); Lee and Barfknecht (1992). The motion of kinks in such a window junction was studied by one of the authors using periodic boundary conditions in $`x`$ Flytzanis et al. (2000). It was shown that the speed of the kink depends on the extension of the passive region and that Cerenkov resonances occur between the cavity modes and the kink. These show up as steps in the zero field steps (ZFS) in the IV characteristics. The experiments conducted in this geometry, for a fixed set of electric parameters Barbara et al. (1996) confirm these findings.

In the second limiting case, which we will call the 1D model, the passive region exists only at each end of the junction as shown in the bottom left panel of Figure 1. Then the passive region acts in a different way, and fixes the speed of the kink Benabdallah et al. (2002). In Benabdallah et al. (2002) we also studied the 2D problem i.e. a homogeneous passive region surrounding the junction, using an adapted numerical method, the finite volume approximation together with soliton perturbation theory. The 1D model showed that interfaces can act as fluxon traps. Numerical results in the two dimensional case indicated that the sheet inductance miss-match $`L_I/L_J`$ between the junction and the passive region is the main cause of instability of the kink motion.

Here we confirm these results, give the details of the characteristic IV curves and ZFS for the 1D and 2D models and compare these two situations for different geometrical and electrical parameters. We also give a complete derivation of the continuum model used in Benabdallah et al. (2002) from the Resistive Shunted Junction (RSJ) approximation and explain discreteness effects observed in the solution of the 1D model. For the 2D model we justify the fact that the kink adapts its velocity to the one in the passive region. Finally we provide insight into the kink stability for large sheet inductance in the passive region $`L_I`$ or capacitance per unit surface $`C_I`$. In the 2D case we give specifically the different IV curves which can be seen for the different values of the sheet inductance in the passive region $`L_I`$.

The paper is organized as follows. After deriving the continuum equations for the RSJ model in section 1, we introduce the 1D model and study its zero field steps in section 2. In Section 3 we study the ZFS for the 2D model and compare them to the ones for the 1D model. We give our concluding remarks in section 4.

## II  The model

### II.1 RSJ model for window Josephson junction

A simple mathematical model of the window Josephson junction is to describe each superconductor by an array of inductances $`L`$ (see Fig. 2). The coupling elements between two adjacent nodes in each array are, a capacitor $`C`$, a resistance $`R`$ and a Josephson current $`I_c`$ Likharev (1986); Devoret . The Kirchoff laws at each couple of nodes $`(i_b,i_t)`$ in the bottom and top superconducting layers can be combined to give the relation expressing the conservation of currents at node $`i`$ in the Josephson junction

$$C_J\ddot{\mathrm{\Phi }}_i+\underset{j}{}\frac{\mathrm{\Phi }_i\mathrm{\Phi }_j}{L_{ij}}+I_i^c\mathrm{sin}\frac{\mathrm{\Phi }_i}{\mathrm{\Phi }_0}+\frac{\dot{\mathrm{\Phi }}_i}{R_i}=0,$$
(1)
and in the passive region

$$C_I\ddot{\mathrm{\Psi }}_i+\underset{j}{}\frac{\mathrm{\Psi }_i\mathrm{\Psi }_j}{L_{ij}}=0,$$
(2)
where $`\mathrm{\Phi }=\mathrm{\Phi }_t\mathrm{\Phi }_b`$ (resp. $`\mathrm{\Psi }=\mathrm{\Psi }_t\mathrm{\Psi }_b`$), is the phase difference between the two superconductors in the junction (resp. passive) part, the summation $`_j`$ is applied to the nearest neighbors. Also note that equations (1) and (2) are discretisations of Maxwell’s equations (wave equation part) and Josephson constitutive equations (sinus term), assuming an electric field normal to the plates, a magnetic field in the junction plane and perfect symmetry between the top and bottom superconducting layers. We can now obtain the model in the continuum limit, more suitable for analysis.

### II.2 Continuum limit

The continuum version of the system (1)-(2) can be derived by introducing the following quantities per unit area ($`a^2`$) of elementary cells of length $`a`$.

$$\overline{C_{I,J}}=\frac{C_{I,J}}{a^2},j_c=\frac{I_i^c}{a^2},\overline{r}=Ra^2.$$
(3)
We normalize the phases by the flux quantum $`\mathrm{\Phi }_0`$,

$$\varphi _i=\frac{\mathrm{\Phi }_i}{\mathrm{\Phi }_0},\psi _i=\frac{\mathrm{\Psi }_i}{\mathrm{\Phi }_0},$$
(4)
and introduce the Josephson characteristic length $`\lambda _J=a\stackrel{~}{\lambda _J}`$

$$\lambda _J^2=\frac{\varphi _0}{2\pi j_cL_J}=\frac{\mathrm{\Phi }_0}{j_cL_J}=\stackrel{~}{\lambda _J}^2a^2.$$
(5)
Notice that in this 2D problem the inductance associated to each cell is equal to the branch inductance $`L_J`$ (resp. $`L_I`$) in the junction (resp. passive region). This is not the case for the 1D model for which these inductances are in series, giving a total inductance proportional to the mesh size $`dx`$.

We substitute relations (3) - (5) in (1) and (2) and obtain in the junction

$$a^2\mathrm{\Phi }_0\overline{C_J}\ddot{\varphi }_i+\mathrm{\Phi }_0\underset{j}{}\frac{\varphi _i\varphi _j}{L_{ij}}+a^2j_c\mathrm{sin}\varphi _i+\frac{a^2}{\overline{r}}\mathrm{\Phi }_0\dot{\varphi }_i=0,$$
(6)
and in the passive part

$$a^2\mathrm{\Phi }_0\overline{C_I}\ddot{\psi }_i+\mathrm{\Phi }_0\underset{j}{}\frac{\psi _i\psi _j}{L_{ij}}=0.$$
(7)
After simplification of Eqs. (6)-(7) and using Eq. (5) we obtain

$$\overline{C_J}L_J\ddot{\varphi }_i+\underset{j}{}\frac{\varphi _i\varphi _j}{a^2}+\frac{1}{\lambda _J^2}\mathrm{sin}\varphi _i+\frac{L_J}{\overline{r}}\dot{\varphi }_i=0,$$
(8)
$$\overline{C_I}L_I\ddot{\psi }_i+\underset{j}{}\frac{\psi _i\psi _j}{a^2}=0.$$
(9)
Now we take the limit $`a0`$, so that
$`{\displaystyle \underset{j}{}}{\displaystyle \frac{\varphi _i\varphi _j}{a^2}}\mathrm{\Delta }\varphi `$ and $`{\displaystyle \underset{j}{}}{\displaystyle \frac{\psi _i\psi _j}{a^2}}\mathrm{\Delta }\psi `$.
and obtain the following system of partial differential equations

$$\overline{C_J}L_J\frac{^2\varphi }{t^2}\mathrm{\Delta }\varphi +\frac{1}{\lambda _J^2}\mathrm{sin}\varphi +\frac{L_J}{\overline{r}}\frac{\varphi }{t}=0,$$
(10)
$$\overline{C_I}L_I\frac{^2\psi }{t^2}\mathrm{\Delta }\psi =0.$$
(11)
We introduce the plasma frequency in the junction $`\omega _J^2=\lambda _J^2v_I^2,`$ where $`v_I^2=\overline{C_J}L_J`$ is the Swihart velocity Swihart (1961). Similarly we define the electromagnetic wave frequency in the passive part $`\omega _I^2=\lambda _J^2v^2,`$ where $`v^2=\overline{C_I}L_I`$. In the following we fix the couple $`(L_J,C_J)`$ and $`\omega _J`$ and vary the inductance $`L_I`$ and capacity $`C_I`$ in the passive region.
We normalize equations (10) - (11) in space by $`\lambda _J`$ and in time by $`\omega _J^1`$ and obtain in the junction ($`\mathrm{\Omega }_J`$)

$$\frac{^2\varphi }{t^2}\mathrm{\Delta }\varphi +\mathrm{sin}\varphi +\alpha \frac{\varphi }{t}=0,$$
(12)
and in the passive region $`(\mathrm{\Omega }\backslash \mathrm{\Omega }_J)`$

$$\left(\frac{\omega _J}{\omega _I}\right)^2\frac{^2\psi }{t^2}\mathrm{\Delta }\psi =0,$$
(13)
where the dimensionless parameter $`\alpha =\frac{\stackrel{~}{\lambda _J}}{R}\sqrt{\frac{L_J}{\overline{C_J}}}`$ is a damping coefficient which depends on the spatial discretization via $`\stackrel{~}{\lambda _J}`$.

The boundary conditions of equations (12) and (13) are of inhomogeneous Neumann type which physically indicates a lateral injection of current or/and an external magnetic field

$$\frac{1}{L_I}\psi n=J_{ext}\mathrm{on}\mathrm{\Omega },$$
(14)
where $`n`$ is the exterior normal. To this we add the interface conditions for the phase and its normal gradient, the surface current on the junction boundary $`\mathrm{\Omega }_J`$

$$\psi =\varphi \mathrm{and}\frac{1}{L_I}\frac{\psi }{n}=\frac{\varphi }{n}$$
(15)
Note that the jump condition (2nd relation in (15)) can be obtained by integrating (12) - (13) on a small surface overlapping the junction domain $`\mathrm{\Omega }_J`$.
In the rest of the paper we assume a rectangular window of length $`l=10`$ and width $`w=1`$ embedded in a rectangular passive region of extension $`w^{}`$ as shown in Figure 1. We will not consider the influence of an external magnetic field and will assume the external current feed to be of overlap type so that the boundary conditions (14) become

$`{\displaystyle \frac{1}{L_I}}{\displaystyle \frac{\psi }{y}}(x,0)`$ $`=`$ $`{\displaystyle \frac{I}{2}},{\displaystyle \frac{1}{L_I}}{\displaystyle \frac{\psi }{y}}(x,𝒲)={\displaystyle \frac{I}{2}},`$
$`{\displaystyle \frac{\psi }{x}}(0,y)`$ $`=`$ $`{\displaystyle \frac{\psi }{x}}(,y)=0,`$ (16)
where $`=2w^{}+l\mathrm{and}𝒲=2w^{}+w`$.

We assumed throughout the study a small damping $`\alpha =0.01`$ which is typical of under damped Josephson junctions.

## III A 1D window Josephson junction

### III.1 The model

We first consider the simplified situation where the passive region is present only at the two lateral ends of the device as shown in the bottom left panel of Fig. 1. For this system it is possible to derive the effective 1D equation Benabdallah et al. (2002)

$`C(x)\phi _{tt}_x\left({\displaystyle \frac{1}{L(x)}}\phi _x\right)+ϵ(x)(`$ $`\mathrm{sin}\phi +\alpha \phi _t)\gamma (x)=0`$ (17)
$`\mathrm{for}x]0,[,`$
where $`\phi `$ is the phase averaged in y across the junction and $`\gamma (x)=\frac{I}{L(x)s}`$ where $`s`$ is the junction area.
In equation (17) $`C,L`$ and $`ϵ`$ are discontinuous functions of $`x`$ defined by

$$(C(x),L(x),ϵ(x))=\{\begin{array}{cc}(1,1,1)& \mathrm{if}x]w^{},w^{}+l[\\ & \\ (C_I,L_I,0)& \mathrm{if}x]0,w^{}[]w^{}+l,[,\end{array}$$
(18)
The boundary conditions (16) reduce then to

$$\phi _x(0)=\phi _x()=0.$$
(19)
To compare different values of $`w^{}`$ we normalize the current $`I`$ by the maximum current that the junction can carry $`I_{\mathrm{max}}=s`$. To simplify the calculations we have assumed a constant current density $`\gamma `$.

We choose the values of the inductance $`L_I`$ and capacity $`C_I`$ according to two strategies. First we assume $`L_IC_I=v_I^2`$ is constant which gives a uniform velocity in the passive part. In a second set of calculations we vary $`v_I`$ by moving orthogonally to the set of hyperbolas in the $`(L_I,C_I)`$ plane.

The numerical method to obtain the (IV) curve is to start with a static kink like solution obtained by solving the static problem derived from (17) and to increase progressively the current to obtain a moving fluxon. After obtaining a stable solution in time, we compute the voltage in the center of the junction by

$`V{\displaystyle \frac{1}{T_2T_1}}{\displaystyle _{T_1}^{T_2}}\phi _t(w^{}+l/2,t)dt=`$
$`{\displaystyle \frac{\phi (w^{}+l/2,T_2)\phi (w^{}+l/2,T_1)}{T_2T_1}},`$ (20)
where the times $`T_1`$ et $`T_2`$ are given such that the voltage is stable.

### III.2 Discreteness effects

In this section we consider a window junction with homogeneous electric properties ($`L_I=C_I=1)`$ and fix the width of the passive part $`w^{}=2`$. We have chosen the number of discretization points in the $`x`$ direction as $`N=253,503,703`$ and $`1003`$. In this case we compute the IV characteristic as shown in the top panel of Fig. 3.

For small values of $`N`$, we observe resonances in the IV curve which decrease and disappear as the number of mesh points is increased from $`N=253`$ to 703. These disappear altogether for $`N=1003`$. To see if these resonances are due to the presence of the passive region, we have calculated a IV curve for a homogeneous junction ($`w^{}=0`$) and indeed they are absent. To understand the mechanism of this fine structure in the IV curve, we have plotted on the bottom panel of Fig.3 the instantaneous voltage $`\varphi _t`$ as a function of $`x`$ for a fixed time and $`N=703,1003`$ and $`1403`$. From the comparison of the two panels of Fig. 3, it is clear that the radiation is responsible for the fine structure. In fact the wave length of the excited radiation corresponds to a standing mode $`V=2n\pi /4`$ with $`n10`$.

Let us now explain how such radiation is excited by the kink as it crosses an interface from a linear to a nonlinear medium. Fig. 4 shows the phase velocity $`v_\varphi =\omega /k`$ for the continuum sine-Gordon system (where $`\mathrm{sin}\varphi \varphi `$) and the discrete system used numerically for which the second derivative is approximated by a three points difference $`\varphi _{xx}=(\varphi _{n+1}+\varphi _{n1}2\varphi _n)`$. For the former $`v_\varphi =\sqrt{1+1/k^2}>1`$ so that nonlinear modes cannot be excited via the long wave short wave Benney resonance mechanism Dodd et al. (1982). This is not true for the discrete system for which

$$v_\varphi =\frac{1}{k}\sqrt{1+\frac{4}{h^2}\mathrm{sin}^2(kh/2)}$$

can be smaller than 1 for $`h0`$ and small wave lengths. As expected the threshold $`k`$ decreases as is increased. Notice that the scale in $`k`$ corresponds to $`9\times \pi /2k12\times \pi /2`$.

We then explain quantitatively the resonances in the IV curve shown in Fig. 3 by the fact that as the kink velocity increases, it can lock with a given cavity mode. For example $`N=703`$ gives voltage (resp. velocity) steps $`V=0.4457,0.44645,0.447,0.4475`$ (resp. $`v=0.9935,0.9952,0.996,0.997)`$ which correspond to the cavity mode index $`n=14,13,12`$ and $`11.`$ For $`N=503`$ we observe steps at voltages $`V=0.4451,0.446,0.4467`$ (resp. $`v=0.9922,0.99426,0.9958)`$ for which $`n=11,10`$ and $`9`$ respectively. The lower values of $`n`$ explain the larger amplitudes observed in the radiation.

An intrinsic fine structure in the IV characteristics has been seen in experiments with window junctions with a lateral passive region. For example Thyssen et al. in Parmentier and Pedersen (1995) show a slight shift in the position (and therefore limiting fluxon velocity) of the IV curves for different values of $`w^{}`$. Ustinov et al Barbara et al. (1996) also display IV curves which depend on $`w^{}`$ and large substructures. Recent numerical work Flytzanis et al. (2000) on a Josephson window junction with a lateral passive region and periodic boundary conditions along the propagation direction confirm these resonances due to Cerenkov radiation between a soliton travelling faster than $`v=1`$ and the radiation of phase speed $`v_\varphi =\sqrt{1+1/k^2}>1`$. Here we show that when the passive region exists only in the propagation direction, resonances disappear at large resolutions. This can be expected from the calculation of the emitted power of radiation for a kink as it crosses an interface Kivshar and Chubykalo (1991); Kivshar and Malomed (1989). This quantity drops fast to zero as the speed approaches 1.

In the rest of the study we have chosen $`N=1403`$ so that numerical resonances are eliminated from the IV characteristic.

### III.3 Zero field steps (ZFS)

We now proceed to give an estimate of the time average of the voltage for the window junction. For this, we compute the velocity in the junction using the McLaughlin-Scott soliton perturbation equations McLaughlin and Scott (1978) and get for the kink velocity

$$v=\frac{1}{\sqrt{1+\left(\frac{4\alpha }{\pi \gamma }\right)^2}}.$$

When the current $`\gamma `$ increases and becomes larger than $`\alpha `$, the velocity $`v`$ gets close to 1.
In the passive part the velocity $`v_I`$ is given by (17)

$$v_I=\frac{1}{\sqrt{L_IC_I}},$$

so that the average of the voltage in time in the limit of large current can be approximated by

$$V_{1\mathrm{D}}<\phi _t>=\frac{\mathrm{\Delta }\phi }{\mathrm{\Delta }t}=\frac{2\pi }{2w^{}/v_I+l/v},$$
(21)
where the denominator is just the time take by the kink (phase jump of $`2\pi `$) to travel across the device. We now confirm this estimation first when the electric properties of the device are homogeneous and then when they are different in the junction and the passive region so that $`v_I1`$.

### III.4 Influence of geometry

First we consider a situation where the electric properties are homogeneous in the whole device so that $`(L_I=C_I=1)`$ and change the extension of the passive region $`w^{}`$.
In Fig. 5 we plot the (IV) characteristic curves for four values of $`w^{}=0,1,2`$ and $`3`$. This figure. show that the position of the ZFS moves towards the left when $`w^{}`$ increases. This agrees with formula (21). Table (1) reports the limiting values of the zero field steps and the estimates from (21), which are in excellent agreement. For large value of $`w^{}`$ (typically $`w^{}=12`$), the kink becomes unstable. This seems to be due to the fixed points that exist on each junction/passive region interface Benabdallah et al. (2002). The kink can be trapped easier by one of them as $`w^{}`$ increases because the driving force due to the current gets weaker.

### III.5 Influence of the electrical parameters

The second case tested is when the velocity in the passive region $`v_I1`$. To simplify things we fix $`w^{}`$ to $`2`$ in all the presented runs. We show in Fig. 6 the IV characteristics for four values of the velocity in the passive region $`v_I=0.33,0.5,1`$ and $`2`$. Notice that the positions of the ZFS of the 1D effective model follow (21) as shown in table (2). Fig. 6 also shows that the position of the ZFS moves from the right to the left as $`v_I`$ decreases from $`2`$ to $`0.33`$. When $`v_I`$ goes to zero so that $`L_I`$ or $`C_I`$ goes to infinity, we expect the ZFS to disappear. We observe regions of instability for $`L_I>>1`$ or $`C_I>>1`$ and the zero field steps exist only in small intervals in the $`(L_I,C_I)`$ plane. In Fig. 8 we present a detailed numerical exploration of the $`(L_I,C_I)`$ plane. The hyperbolas corresponding to $`L_IC_I=v_{I}^{}{}_{}{}^{2}`$ constant are shown. We indicate by the symbols $`(),`$ $`(\times )`$ and $`(+)`$ the points where the ZFS exist. We note that for these points the solution of the window 1D problem is a kink. The region of stability is concentrated along the diagonal. Note that the results of Fig. 8 are based only on the numerics.
When $`C_I`$ or $`L_I`$ become large we obtain particular types of solutions which we now discuss.

### III.6 Solutions for large $`L_I`$ or $`C_I`$

Instabilities of the kink motion occur when the extension of the passive region $`w^{}`$ is large, typically $`w^{}12`$. In that case the current density becomes small and cannot drive the kink out of the potential well created by the interface Benabdallah et al. (2002). Other instability factors are the electrical parameters $`(L_I,C_I)`$. We have shown in Benabdallah et al. (2002); Benabdallah (1999) that the electrical parameters act differently from the geometry. For $`C_I>>1`$ the solution is static while it is dynamic for $`L_I>>1`$. In this latter case, kink motion is possible only in the junction.

We now consider the disappearance of the Zero field steps in the $`IV`$ curves. This occurs for intermediate values of $`C_I`$ and $`L_I`$ and gives rise to a large voltage.

Let us recall the equations of the problem

$$\frac{^2\varphi }{t^2}\frac{^2\varphi }{x^2}+\mathrm{sin}\varphi =\gamma \alpha \frac{\varphi }{t}w^{}xw^{}+l,$$
(22)
$$C_I\frac{^2\psi }{t^2}\frac{1}{L_I}\frac{^2\psi }{x^2}=\gamma 0xw^{}\mathrm{and}w^{}+lx,$$
(23)
together with the interface condition at $`x=w^{},w^{}+l`$

$$\varphi =\psi ,\frac{1}{L_I}\frac{\psi }{x}=\frac{\varphi }{x},$$

and homogeneous boundary conditions $`\frac{\psi }{x}|_{x=0}=\frac{\psi }{x}|_{x=}=0`$.

In the limit of large voltages $`V`$

$$\frac{\varphi }{t}V\mathrm{and}\frac{\psi }{t}V$$

so that one can average the above equations which then reduce to

$$\frac{^2\varphi }{x^2}=\gamma \alpha Vw^{}xw^{}+l,$$
(24)
$$\frac{1}{L_I}\frac{^2\psi }{x^2}=\gamma 0xw^{}\mathrm{and}w^{}+lx,$$
(25)
where the voltage $`V=\gamma /(\alpha l)`$ can be derived from the work equation (see Appendix).

The time dependent part of $`\varphi ,\psi `$ is now just $`Vt`$ and the $`x`$ dependent part verifies a boundary value problem with $`\psi _x|_{x=0,}=0`$ and the above given interface conditions. The solution is symmetric with respect to $`x=/2`$ and we obtain for the left and right passive regions and the junction

$$\psi _l=\gamma L_I\frac{x^2}{2}+Vt,$$
(26)
$$\psi _r=\gamma L_I\frac{(x)^2}{2}+Vt$$
(27)
$$\varphi =\gamma \frac{w^{}}{l}x(x)+\gamma \frac{w_{}^{}{}_{}{}^{2}}{2}(2L_I+\frac{2w^{}}{l})+Vt$$
(28)
The above given expressions give a very good approximation of the phase computed numerically as can be seen in Fig. 7 which shows $`\psi (x)`$ and $`\varphi (x)`$ for 11 successive values of time. The electrical parameters are $`C_I=1,L_I=10`$ and the top panel corresponds to a small current density $`\gamma =0.1`$ while the bottom panel is for $`\gamma =0.25`$. The agreement for the latter is excellent, showing complete overlap between the approximation (26,27,28) and the numerical solution. On the top panel the approximations give the overall value but there are some oscillations due to a resonance between the main frequency $`V`$ and a subharmonic.

It is possible to estimate the validity of the approximate solution (26,27,28) by considering perturbations around it. This leads to the linearized equations for the perturbations $`\xi ,\zeta `$.

$$\frac{^2\xi }{t^2}\frac{^2\xi }{x^2}+\xi \mathrm{cos}\varphi +\alpha \frac{\xi }{t}=0w^{}xw^{}+l,$$
(29)
$$C_I\frac{^2\zeta }{t^2}\frac{1}{L_I}\frac{^2\zeta }{x^2}=00xw^{}\mathrm{and}w^{}+lx$$
(30)
with the interface condition at $`x=w^{},w^{}+l`$

$$\xi =\zeta ,\frac{1}{L_I}\frac{\zeta }{x}=\frac{\xi }{x},$$

and boundary conditions $`\frac{\zeta }{x}|_{x=0}=\frac{\zeta }{x}|_{x=}=0`$.

We then assume a uniform time dependence $`e^{\alpha t/2+i\omega t}`$ and obtain a spectral problem which can be solved. We can neglect the damping because it is small.

In the particular case presented in the top panel of Fig. 7 the voltage is $`V=14`$ and the subharmonic observed is $`V/4=3.5`$. This corresponds to the mode $`n_J=11`$ in the junction and $`n_I=7`$ in each lateral passive region of extension 2, taking into account the inductance $`L_I=10`$.

Fig. 8 shows the regions of existence of ZFS in the $`(L_I,C_I)`$ plane. Notice the symmetry. For $`C_I>>1`$ the solution is static while it is dynamic for $`L_I>>1`$, this latter parameter plays a predominant role due to the interface conditions.

## IV 2D window junction

In this section we discuss the dynamics of the kink in the window junction and how to reduce it to a 1D effective problem. It is known that for a pure Josephson junction this reduction is possible. For example in the static case the behavior of a junction of width $`w<2\pi `$ is very well approximated by the equation for the zero order mode Caputo et al. (1996b). For the dynamical case Eilbeck et al. Eilbeck et al. (1985) found that when $`(I/81)`$ the phase is uniform in the $`y`$ direction so that the reduction is possible. For a window junction the situation is different because we have a passive region around the junction which has electrical properties, a capacity and an inductance which fix the velocity $`v_I`$. The comparison between 2D and 1D was done for the two cases of homogeneous electrical properties in different geometries and a fixed geometry with different electrical properties. Before giving these results we discuss how the kink velocity is fixed by the presence of the passive region.

### IV.1 Justification of the effective 1D model

When the velocity in the passive region $`v_I1`$ the formula (21) is completely off for the 2D calculations and the reduction of the 2D window junction problem to a 1D effective problem is impossible.
Here the lateral passive region where waves can only propagate at velocity $`v_I`$ is controlling the kink motion. The velocity -a free parameter for the sine-Gordon equation- can then adjust itself to $`v_I`$ so that the ”dressed” kink can propagate in the device (junction and lateral passive region). The fluxon adapts its velocity to $`v_I`$ so that $`v=v_I`$.

To justify this, we take the time derivative of the first interface condition (15) and obtain

$$\varphi _t=\psi _t\mathrm{on}\mathrm{\Omega }_J,$$
(31)
the second interface condition (15) reads

$$\varphi n=\psi n\mathrm{on}\mathrm{\Omega }_J.$$
(32)
The numerical simulations show that it is legitimate to assume that the solution is a 1D kink travelling in the passive region with the velocity $`v_I`$ and in the junction with velocity $`v`$. Then (32) becomes at the longitudinal boundaries

$$_x\varphi =_x\psi \mathrm{on}\mathrm{\Omega }_J.$$
(33)
The velocity $`v`$ of the kink on the boundary of the junction is defined by

$$v=\frac{\varphi _t}{\varphi _x}\mathrm{on}\mathrm{\Omega }_J.$$
(34)
and the same for the velocity $`v_I`$ of the linear wave on the boundary of the junction

$$v_I=\frac{\psi _t}{\psi _x}\mathrm{on}\mathrm{\Omega }_J.$$
(35)
Thus from (31) and (33) we obtain $`v=v_I`$ so that the estimate for the average voltage is then

$$V_{2\mathrm{D}}=\frac{2\pi }{2w^{}/v_I+l/v_I}$$
(36)
### IV.2  Influence of the geometry on the zero field steps

In this case the electrical properties are homogeneous in the whole device so that ($`L_I=C_I=1`$). The numerical simulations show that for different $`w^{}`$ the solution of the 2D problem are $`y`$ uniform. For $`w^{}=2`$ for example the Fig. 9 represent the contour of the phase in the whole domain at the time $`T=\mathrm{5.024\hspace{0.17em}10}^4`$. We note that the Fig. 9 was normalized by $`\mathrm{max}|\varphi |`$. The levels lines are parallel in the $`y`$ direction so that the solution is very close to a 1D kink. This 1D kink propagates with a velocity $`v_I`$.

The figure (10) shows the IV characteristics for the both 1D and 2D problems. One can see the good agreement between the 1D and 2D calculations. Table (3) reports the limiting values of the zero field step computed numerically and the corresponding estimates from (21) in excellent agreement.
As the numerical results show, in this case the dynamics of the solution in the 2D window problem can be described by a unidimensional fluxon and the reduction to a 1D effective problem is possible.

### IV.3 Influence of the electrical parameters

We have obtained the same average voltage for different values of $`L_IC_I`$ as long as the product $`L_IC_I=v_I^2`$ is kept constant. We observe that the position of the zero field step moves from the right to the left as the velocity in the passive region $`v_I`$ decreases from $`2`$ to $`0.5`$. We expect the ZFS to disappear when $`v_I`$ goes to $`0`$ ie $`L_I`$ or $`C_I`$ go to infinity.

### IV.4  Instabilities

As in 1D case there are two types of instabilities of the kink. The first is due to the geometry of the window, because large widths $`w^{}`$ of the passive region (typically $`w^{}=7`$) lead to a ”stretched” kink in the junction which strongly radiates. It is the impossible to accelerate such an object and have a stable zero field step Benabdallah et al. (2002). This solution which is almost uniform in the junction will oscillate but not rotate and the average voltage will be zero Benabdallah et al. (2002). This situation was observed in the experimental work of Thyssen et al. Parmentier and Pedersen (1995) who never found any zero field step for $`w^{}/w<3`$.

The second type of instability is due to the electrical quantities $`L_I`$ and $`C_I`$. For $`w^{}=2`$, we have investigated the $`(L_I,C_I)`$ plane and found numerically the regions of instability of the kink as shown in Fig. 11. We also plot three hyperbolas corresponding to $`v_I=0.5(),1(\times )`$ and $`2(+)`$. The markers indicate the positions where we have found the ZFS. All the points for a given marker correspond to the same voltage. In Fig. 11 we have isolated the stability region as the interior of the domain bounded by the solid line. We observe that the instability is due essentially to the inductance $`L_I`$ so that the kink motion is stable for $`0.1L_I2`$, independently of $`C_I`$. For example $`C_I=10^4`$ and $`L_I=1`$ gives rise to a ZFS. We can estimate the average voltage using formula (36), ($`V_{2\mathrm{D}}=4.5\times 10^3`$ which is very close to the value found numerically $`V_{\mathrm{num}}2D=4.577\times 10^3`$). Thus even for such a large value of $`C_I`$, the formula (36) gives a good agreement with the numerical calculation.
Small values of $`C_I`$ lead to a static solution as can be seen from equation (13) which reduces to $`\mathrm{\Delta }\psi =0`$ so that $`\psi `$ is constant. Because of the interface conditions $`\varphi `$ is also constant.

The influence of the inductance $`L_I`$ has been analyzed in detail by calculating IV characteristics for $`w^{}=3`$ and $`0.1L_I3`$ which are plotted in Fig. 12. There are three main regions identified by the letters A B and C which correspond to three different dynamical behaviors. In region A $`2L_I3`$, we obtain a static solution and a zero voltage. At this time we do not understand why the kink motion becomes unstable for $`L_I2`$ independently of $`C_I`$. This effect is present only in two dimensions and therefore linked to the presence of a lateral passive region. It may be due to the large jump in the gradient $`\varphi _y=(1/L_I)\psi _y`$ along the lateral interface.
Region B corresponds to a ZFS thus a stable kink motion in the device. The limiting voltage is given to a very good approximation by formula (36) as shown by the inset of Fig. 12. This expression also gives a good estimate even for $`v_I>1`$. Region C $`L_I=0.1,0.2`$ corresponds to much larger values of the voltage and a straight behavior typical of linear resonances in the junction cavity. Small values of $`L_I`$ give rise to a spatially uniform phase in the passive region, which reacts to the phase in the junction. The voltage for $`L_I=0.1`$ corresponds to a higher order cavity mode $`V=2.566\pi /10`$. Fig. 13 shows the corresponding phase at a given instant (top panel) together with the temporal behavior (bottom panel). As expected the phase is almost uniform, the bounds in Fig. 13 are 191.83 and 193.14. The phase increases linearly with time.

Then we conclude that the region of stability of the kink in the window problem is

$$D_{\mathrm{stability}}=\{0.1L_I2\mathrm{and}C_I0.5\}.$$
(37)
## V Conclusion

In this paper we have investigate the electrodynamics of a uniform window junction, using the limiting case termed 1D model where the lateral extension of the passive region is neglected. We use the solution of the 1D model as a guide to discuss the behavior of the full window junction and this enables us to give an estimate of the average voltage for different extensions $`w^{}`$ of the passive region and different electrical parameters which fix the velocity in the passive region $`v_I`$. This estimate is in excellent agreement with the numerical results even when $`v_I1`$. In this situation we show that the fluxon adapts its velocity to $`v_I`$. For both cases $`v_I=1`$ and $`v_I1`$ the numerical results show that the fluxon propagates perpendicularly to the direction of the junction. We have doubled the junction length to $`l=20`$ and find the same estimates for the limiting voltage. In the 2D window junction, the translational symmetry is broken so that we obtain a different result from the case where only a lateral passive region exists.

We have also studied numerically the stability of the zero field step in the window junction and show that when the width of the passive part becomes large the fluxon becomes distorted and gives to another type of solution which is radial. In this case we have not found any ZFS.

This study shows that the kink can travel into the passive region even when the impedance is not adapted as long as the mismatch is in the capacitance and not the inductance. In other words the limit $`ϵ0`$ is not a singular limit. On the other hand if there is an inductance mismatch then this will most likely cause the break-up of kink shuttling leading to the disappearance of the zero field step. The stability of the fluxon depends essentially on the inductance $`L_I`$. This result, we believe, can lead to improving Josephson devices and their coupling to micro-strip lines.

Acknowledgements It is with great pleasure that J. G. C. thanks N. Flytzanis for many useful discussions at the beginning of this work. A. B. thanks the Laboratory of Mathematics of the INSA de Rouen and the Laboratoire de Physique théorique de l’Université de Cergy-Pontoise for support. A. B. and J. G. C. thank Vladislav Kurin and Alexei Ustinov for many useful electronic exchanges. We thank Jérôme Dorignac for useful discussion. This work was supported by the European Union under the RTN project LOCNET HPRNCT-1999-00163.

## Appendix A Derivation of the linear part of the $`IV`$ curve

To determine the Hamiltonian associated with the 1D-window junction equation we multiply Eq. (17) by $`\phi /t`$ and integrate between $`0`$ and $``$ to obtain

$`{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[{\displaystyle _0^{}}{\displaystyle \frac{C(x)}{2}}\phi _t^2+ϵ(x)\left(1\mathrm{cos}\phi \right)\mathrm{d}x\right]`$
$`{\displaystyle _0^{}}_x\left[{\displaystyle \frac{1}{L(x)}}\phi _x\right]\phi _t\mathrm{d}x={\displaystyle _0^{}}\gamma (x)\phi _tdx`$
$`\alpha {\displaystyle _0^{}}ϵ(x)\phi _t^2dx,`$
Integrating by parts the second term in the left-hand side and using the boundary conditions we obtain

$`{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[{\displaystyle _0^{}}{\displaystyle \frac{C(x)}{2}}\phi _t^2+{\displaystyle \frac{1}{2L(x)}}\phi _x^2+ϵ(x)\left(1\mathrm{cos}\phi \right)\mathrm{d}x\right]`$
$`={\displaystyle _0^{}}\gamma (x)\phi _tdx\alpha {\displaystyle _0^{}}ϵ(x)\phi _t^2dx.`$
Recalling the sine-Gordon-wave Hamiltonian

$$H=_0^{}\left[\frac{C(x)}{2}\phi _t^2+\frac{1}{2L(x)}\phi _x^2+\epsilon (x)\left(1\mathrm{cos}\phi \right)\right]dx$$
(38)
The computation of the linear part of the $`IV`$ curve is based on the conservation of the energy of the system that is given by sG-wave Hamiltonian (38) from which it can be shown that $`H`$ satisfies the power balance equation

$`{\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}=\gamma {\displaystyle _0^{}}\phi _tdx\alpha {\displaystyle _w^{}^{w^{}+l}}\phi _t^2dx`$ (39)
In the stationary regime the average $`dH/dt=0`$ , so

$$\gamma _0^{}\phi _tdx\alpha _w^{}^{w^{}+l}\phi _t^2dx=0.$$

Since the voltage $`V`$ is the temporal average of $`\phi _t`$, we obtain in the high voltage limit $`\phi _t^2=\phi _t^2`$, so that the IV curve is given by

$$V=\frac{\gamma }{\alpha l}$$


text:
# BCS-Bose Crossover in Color Superconductivity

## Abstract

It is shown that the onset of the color superconducting phase occurs in the BCS-BE crossover region.

During the last 2-3 years color superconductivity became one of the QCD focal points –see review papers -. To a large extent (but not completely) the basic ideas of the subject are traced back to the BCS theory of superconductivity and its later development. It is known that the discovery of high temperature superconductors (HTSC) gave rise to new ideas and approaches and revealed interest to the problem of the transition from the BCS regime to the Bose-Einstein (BE) condensation. The BCS-BE crossover is important for the physics of HTSC since the underlying distinction of HTSC from ordinary superconductors is that they are characterized by much smaller value of the dimensionless parameter $`\xi n^{1/3}`$, where $`\xi `$ is the coherence length and $`n`$ is the carrier density. In the BCS, or weak coupling regime, $`\xi ^3n(10^810^{10})`$ while in the opposite strong coupling case $`\xi ^3n\begin{array}{c}<\hfill \\ \hfill \end{array}1`$ and we are dealing with the compact pairs of composite bosons, which may undergo BE condensation. It has been suggested (see e.g. ) that the description of the HTSC might require an intermediate approach between the BCS and BE limits. The evolution from weak to strong coupling was theoretically investigated  before the discovery of HTSC. It was shown that the transition proceeds via a smooth crossover though the two limits are physically quite different (see also ).

Having reminded these well known facts we may formulate the question which forms the core of the present note. As model calculations show  the onset of the color superconducting phase in two flavor QCD(the so called 2SC phase) occurs at rather low quark densities $`n`$, namely at $`n`$ only three times larger than the quark density in normal nuclear matter, or even at lower ones .

Thus it is natural to ask in which region regarding the BCS-BE crossover does it happen?

According to  in QCD with two massless flavors transition to the superconducting $`2SC`$ phase occurs at $`n^{1/3}0.2`$GeV (the dimension of $`n`$ is $`1/fm^3`$ or GeV3). As for the corresponding value of $`\xi `$ we may only rely on some estimates since accurate calculations are lacking. One should also keep in mind possible distinction between the correlation length and the pair size. The two quantities coincide in the BCS regime  while in the BE region the pair size is smaller than the coherence length . With these reservations being made we quote the value $`\xi 0.8`$ fm from . Rather close result, namely $`\xi 0.6`$ fm, follows from the BCS estimate $`\xi 1/\pi \mathrm{\Delta }`$ , where $`\mathrm{\Delta }0.1`$ GeV . Thus the diquark pair in the ”newly born” color superconducting phase is rather compact. This is easy to understand from simple physical considerations. In color antitriplet $`\overline{3}`$ state the one-gluon exchange leads to the quark-quark potential which is only a factor of two weaker than the quark-antiquark one. Instanton or NJL models also result in a rather strong $`qq`$ attraction. Consider the NJL ”weak coupling” solution  $`\mathrm{\Delta }=2\omega _D\mathrm{exp}(1/\rho _{NJL}`$), where $`\omega _D`$ is the Debye frequency, $`\rho _{NJL}=8g^2\mu ^2/\pi ^2`$, $`g^22GeV^2`$ is the NJL coupling constant, $`\mu 0.4GeV`$ is the chemical potential corresponding to the onset of the 2SC phase. We immediately see that $`\rho _{NJL}0.3>\rho _{BCS}`$, i.e. the quark-quark interaction is stronger than phonon mediated electron-electron interaction and in this sense the ”newly born” color superconducting phase does not correspond to the standard BCS weak coupling limit.

We conclude that the onset of the color superconducting phase corresponds to $`\xi n^{1/3}(1fm)(0.2GeV)1,\xi \mu (fm)(0.4GeV)2`$.

These values are at least two orders of magnitude smaller than those corresponding to the BCS regime. In order to understand to which region of the BCS-BE ”phase diagram” they correspond one has to resort to model calculations performed for the system of electrons. Most results have been obtained for the system in two dimensions . Crossover in three dimensions has been studied in . Transition between the two regimes occurs in a narrow range of the parameter $`\xi n^{1/3}`$ and for the electron systems with simple model potentials (finite range, separable, Gaussian) the value $`\xi n^{1/3}1`$ corresponds to the lower limit of the BCS-like region. Needless to say that the extrapolation of this result to the system of massless quarks may be considered only as an educated guess and the problem deserves a dedicated study.

The importance of the BCS-BE crossover for the color superconductivity problem was first outlined on  and  -. In  the low density regime was investigated by extrapolating the single-gluon exchange model which is an adequate tool at asymptotically high densities. The smooth transition from $`\xi n^{1/3}1`$ to $`\xi n^{1/3}=10`$ at $`\mu =0.8`$ GeV was observed.

Let us indicate how the crossover problem in color superconductivity theory has to be approached. The starting point in the thermodynamic potential $`\mathrm{\Omega }(\mathrm{\Delta },\mu ,T)`$, where the diquark condensate $`\mathrm{\Delta }`$ is of the form $`\mathrm{\Delta }<q\widehat{0}q>,`$ $`\widehat{0}=\epsilon _{\alpha \beta 3}\epsilon _{ij}C\gamma _5,\alpha ,\beta `$ – are the color indices, $`i,j`$\- flavor ones, $`C`$\- change conjugation operator . The dependence of $`\mathrm{\Omega }`$ on the chiral condensate $`\phi <\overline{q}\widehat{R}q>,R=\epsilon _{\alpha \beta }\epsilon _{ij}`$, in the region where the two condensates $`\mathrm{\Delta }`$ and $`\phi `$ possibly coexist  may be dropped due to the color superconductivity version of the Anderson theorem . For a wide class of models with four-quark interaction the thermodynamic potential $`\mathrm{\Omega }`$ was first analytically calculated in . With $`\mathrm{\Omega }`$ at hands one can write self-constituent set of mean-field equations to determine the gap $`\mathrm{\Delta }`$ and the chemical potential $`\mu `$:

$$\frac{\mathrm{\Omega }}{\mathrm{\Delta }}=0,\frac{\mathrm{\Omega }}{\mu }=n.$$
(1)
Up to now only the first one of equations (1) has been used in color superconductivity theory while the chemical potential has been considered as independent variable. To display the crossover from the BCS to the low density BE regime one also needs the second equation . It enables to consider the region of $`\mu <0`$ values characteristic for the delute gas of tightly bound diquarks $`(\mu =\epsilon _{B/2}`$ in the limit of delute composite bosons with binding energy $`\epsilon _B`$).

In order to find the quantities $`\mathrm{\Delta }`$ and $`\mu `$ vs the dimensionless physical parameter $`\xi n^{1/3}`$ equations (1) have to be complemented by the equation which determines the parameter $`\xi `$

$$\xi ^2=\frac{𝑑𝐫\varphi (r)r^2}{𝑑𝐫\varphi (r)},$$
(2)
where $`\varphi (r)=<q(r)\widehat{0}q(0)>`$. Analytic expressions for $`\mathrm{\Delta }`$ and $`\mu `$ as functions of $`\xi n^{1/3}`$ have been obtained for electron systems with simple model interaction . The necessaty to implement similar program for color superconductivity directly follows from our conclusion that the suggested onset of this phase occurs within the BCS-BE crossover region.

Another question is what are the physical consequences of the fact that the formation of the color superconducting gap is at least partly due to the existence of the preformed Bose pairs of quarks. At present we can again rely only on the corresponding studies of the electron systems . The key point here is that the physical origin of the critical temperature $`T_c`$ is absolutely different in the limits of weak and strong coupling . In the BCS region $`T_c`$ corresponds to the breaking of Cooper pairs while in the BE limit $`T_c`$ corresponds to the pairs center-of mass motion and to the population of zero-momenta state. Transition from the weak to strong coupling regimes results in the decrease of $`T_c`$. comparing to mean-field value. Formally this should also follow from equations (1)-(2).

Finally we note that calculations of the parameters $`\mu `$ and $`T`$ at which the transition into color superconducting phase occurs have been performed neglecting the gluon condensate. General arguments presented in  show that color-magnetic field which is ”frozen” into the quark system in the form of the gluon condensate shifts the transition towards higher densities. (see also ). Therefore equations (1) for the BCS-BE crossover should be embedded into the background gluon field.

The author is thankful to V.I.Shevchenko and Dirk van der Marel for useful discussion and remarks and to Kazunori Itakura and Egor Babaev for bringing attention to Refs.  and . The author acknowledges receipt of financial support from the grants RFFI-00-02-17836, INTAS-110 and NWO-01-250.