text:
# Random background charges and Coulomb blockade in one-dimensional tunnel junction arrays

## I Introduction

Since the discovery of correlated single electron tunneling many electronic devices exploiting this effect have been proposed . A gate is commonly used to control the device: In single electron transistors to modulate the current, in traps to pull electrons on to an island or in pumps and turn-styles to move electrons in a specific manner. The response to the gate is periodic in the induced charge on the gate capacitance, with a period of one electron. These devices are therefore very sensitive to charges trapped in the vicinity of the device, where there is a significant capacitive coupling to the islands. Randomly distributed “background charge” can change the threshold voltage for injecting charge in the array and the gate voltage where the coulomb blockade is lifted. The background charges may be static, causing a shift of the device operation point, or dynamic, causing noise in the device .

The effect of background charges has been observed in experiments with small arrays. The influence of the background charge shows up as an offset in the current versus gate voltage curves . It has also been shown that the threshold voltage could be maximized by tuning gates to each island and thereby compensating for the background charge . In single electron traps it was shown that the discrepancy between experiments and theory was likely to be due to random background charges . Theoretical studies of the effect of random background charges in one dimensional arrays have examined how the threshold voltage is affected  and noise properties when background charges are included .

Here we discuss simulations of the transport properties of one dimensional arrays of series coupled tunnel junctions. The arrays are characterized by nearest neighbor capacitance $`C`$ and gate capacitance $`C_g`$. The simulations include random background charges which are characterized by their density and amplitude. We examine the effect of static background charges on the transport properties, including the gate voltage dependence, for various ratios $`C/C_g`$.

## II Basic Equations

The one dimensional array consists of $`N1`$ small islands connected in series by small capacitance tunnel junctions. An equivalent circuit is shown in fig. 1. A gate voltage, $`V_{g,i}`$, is coupled through the gate capacitance $`C_{g,i}`$ to each island. A transport voltage is applied across the array by the external voltages $`V_L`$ and $`V_R`$. The tunnel junctions are characterized by the capacitance $`C_i`$ and the resistance $`R_i`$. If the tunnel resistance satisfies $`R_ih/e^2=25.81281\mathrm{}\text{k}\mathrm{\Omega }`$, the number of charges on each island becomes sharply defined and the quasi-classical “orthodox” theory of single electron tunneling applies . The theory predicts that electrons tunnel one at a time with a tunneling rate given by

$$\mathrm{\Gamma }_i=\frac{1}{e^2R_i}\frac{\mathrm{\Delta }F}{1\mathrm{exp}(\mathrm{\Delta }F/k_bT)},$$
(1)
where $`\mathrm{\Delta }F`$ is the difference in free energy of the initial and final state. We are interested in the change in free energy when one electron tunnels through one junction in the array, say from island $`k`$ to island $`k\pm 1`$. In this case the free energy difference is given by the electrostatic potential difference of those islands before and after the tunneling event

$$\mathrm{\Delta }F=F^{}F=\frac{e}{2}[(\phi _{k\pm 1}+\phi _{k\pm 1}^{})(\phi _k+\phi _k^{})].$$
(2)
In the case considered here, all junctions and gate capacitances are equal, $`R_i=R`$, $`C_i=C`$ and $`C_{g,i}=C_g`$. Charge conservation gives a relation between the potential on the islands and the charges on the islands:

$`(2+C_g/C)\phi _1\phi _2`$ $`=`$ $`q_1/C+V_L`$ (3)
$`\phi _{i1}+(2+C_g/C)\phi _i\phi _{i+1}`$ $`=`$ $`q_i/C,2iN2`$ (4)
$`\phi _{N2}+(2+C_g/C)\phi _{N1}`$ $`=`$ $`q_{N1}/C+V_R.`$ (5)
where $`q_i=en_i+q_{0,i}+en_g`$ is the charge on island $`i`$, consisting of the induced background charge on that island $`q_{0,i}`$ and the induced gate charge $`en_g`$. The equations (4) must be inverted to find the potential on the islands. This inversion can be done exactly . The result of this inversion can be expressed in matrix form: $`\overline{\phi }=𝐑\overline{q}`$, where the components of the symmetric matrix $`𝐑`$ are

$$𝐑_{ij}=\frac{\mathrm{cosh}(N|ij|)\lambda \mathrm{cosh}(Nij)\lambda }{2\mathrm{sinh}\lambda \mathrm{sinh}N\lambda }$$
(6)
and $`\lambda =\text{arch}(1+C_g/2C)`$ is known as the inverse of the charge soliton length introduced in ref. . The charge soliton length, $`\lambda ^1`$, is the screening length giving the distance over which the potential from one excess electron is screened by the free charge in the ground plane. With equations (1)-(6), all tunneling rates can be calculated. A Monte Carlo scheme  is used for the simulation of the flow of charges. The random numbers used to generate the background charge configuration and for the monte Carlo calculation are chosen to make a good quasi random series . All results are calculated for $`T=0`$.

## III Effect of random background charges

Charges trapped in the substrate or in the tunnel barriers induce an excess charge on the islands, $`q_iq_i+q_{0i}`$. The background charge $`q_{0i}`$, is between -e/2 and e/2 because any value outside this interval would be compensated for by one electron or hole tunneling on to the island.

Starting from a random background charge configuration, under certain conditions tunneling events will lower the free energy by adding or rearranging electrons or holes. Since the electrostatic energy is proportional to the gradient of the potential squared, lowering the energy means smoothing the potential. This possibility to relax the random potential variations depends on the length of the array, the soliton length, and on the actual configuration of the background charges. In fig. 2 the potential before and after relaxation is shown for two different values of $`\lambda `$. One can see by looking at fig. 2a that for short soliton length the potential is nearly unchanged by the relaxation process. However for long soliton length (fig. 2b) we can see how the relaxation reduces variation of the potential. The potentials plotted in figure 2 for different $`C/C_g`$ are normalized by the threshold voltage. The threshold voltage (derived below) is the proper measure of the Coulomb blockade for transport in an array.

The random background charge configuration can be described by two parameters: the number of islands which have an induced charge $`q_{0,i}0`$, or the density of induced charges $`\rho \{0,N1\}`$, and the amplitude $`A`$ of the induced charges $`q_{0,i}\{Ae,+Ae\}`$ with $`0A0.5`$. In fig. 3 the mean value of the “relaxed” potential variation along the array $`\text{std}\phi `$, averaged over 1000 configurations of random background charges, is plotted versus $`\rho `$ and $`A`$. The four plots in fig. 3 correspond to four different values of the soliton length $`\lambda ^1`$. Here we can see that the smoothing of the potential becomes more effective when the soliton length increases. Comparison of figures 3c and 3d show that when $`\lambda ^1\stackrel{>}{_{}}N`$ the value of $`\text{std}\phi `$ is weakly dependent on $`\lambda `$. Further we can see that when $`\lambda ^1\stackrel{>}{_{}}N`$ and $`\rho \stackrel{>}{_{}}N/2`$, $`A\stackrel{>}{_{}}0.25`$, the potential variations of the array are weakly dependent on $`A`$ and $`\rho `$. This weak dependence can also be seen in fig. 4 where the threshold voltage is plotted as a function of the amplitude (right panels) and density (left panels) for various $`\lambda `$. The threshold voltage changes almost linearly with the amplitude and density of the background charges when $`\lambda =1`$ (figs. 4a and 4c), but when $`\lambda =0.01`$ the dependence is very weak for $`A\stackrel{>}{_{}}0.25`$ and/or $`\rho \stackrel{>}{_{}}N/2`$ (figs. 4b and 4d). Thus we see that as far as the electrostatics is concerned, disorder due to random background charges may be effectively reduced by single charge tunneling if the soliton length is larger than the array length and the density and amplitude of random background charges is large enough.

## IV Gate effects in 1D arrays

When a gate voltage is capacitively coupled to the islands of the array the current voltage characteristics are strongly affected around the threshold. In the case when the gate is uniformly applied ($`V_{g,i}=V_g`$ and $`C_{g,i}=C_g`$) and there are no background charges, the threshold voltage changes with induced charge as shown in fig. (5). For appropriate values of $`N`$ and $`\lambda `$ an array with uniform coupled gate capacitance has ideal transistor characteristics. Such an array would make a good switch because of the initial steep dependence of $`V_t`$ on $`n_g`$, and the wide region of $`n_g`$ with a weaker dependence.

The sensitivity to gate voltage is greatest when $`N`$ is large but $`N<\lambda ^1`$. In this case the threshold voltage is large at zero gate voltage , and the gate drives the threshold almost to zero. The initial linear decrease of $`V_t`$ with $`n_g`$ can be examined analytically because during this decrease the gate voltage is not large enough to bind any charges inside the array. When no charges are localized in the array (electrons, holes or background charges) it was shown in ref.  that the threshold voltage is determined by the tunneling rates through the two end junctions. With eq. (6) we find the relevant potentials to calculate the tunneling rate through the end junction

$$\phi _1=\frac{e}{C}\left[𝐑_{1,1}(1+CV_L/e)+n_g\underset{j=1}{\overset{N1}{}}𝐑_{1,j}\right],$$
(7)
$$\phi _1^{}=\frac{e}{C}\left[𝐑_{1,1}V_LC/e+n_g\underset{j=1}{\overset{N1}{}}𝐑_{1,j}\right],$$
(8)
$$\phi _0=\phi _0^{}=V_L.$$
(9)
Here we have taken $`V_L>0`$, $`V_R=0`$ and $`V_g<0`$ so in this case a hole will tunnel in from the left. Making use of the relations

$$𝐑_{1,j}=\frac{\mathrm{sinh}(Nj)\lambda }{\mathrm{sinh}(N\lambda )},$$
(10)
$$\underset{j=1}{\overset{N1}{}}𝐑_{1,j}=\frac{\mathrm{sinh}[\frac{1}{2}(N1)\lambda ]}{2\mathrm{cosh}(N\lambda /2)\mathrm{sinh}(\lambda /2)}.$$
(11)
We obtain the tunneling rate from eq. (2). Setting the tunneling rate to zero and solving for $`V_L`$ we get the threshold voltage:

$$V_t=\frac{e}{C}\left[\frac{\mathrm{sinh}(N1)\lambda }{\mathrm{sinh}N\lambda \mathrm{sinh}(N1)\lambda }\frac{\mathrm{sinh}[\frac{1}{2}(N1)\lambda ]\mathrm{sinh}N\lambda /2}{\mathrm{sinh}^2\lambda /2\mathrm{cosh}(N\frac{1}{2})\lambda }n_g\right].$$
(12)
If $`n_g`$ is set to zero in this expression we recover the expression for the threshold voltage found in . This equation explains the very steep slope of the decrease in threshold voltage, for large $`\lambda ^1`$ and $`N`$.

Eq. (12) can be used to calculate the voltage gain of the uniform array. If the array is current biased very close to zero and it is assumed that the gate voltage drives the threshold voltage to zero, the voltage gain is given by

$$\eta =\frac{\text{d}V_{out}}{\text{d}V_g}\frac{C_g}{e}\frac{\text{d}V_t}{\text{d}n_g}=4\frac{\mathrm{sinh}[\frac{1}{2}(N1)\lambda ]\mathrm{sinh}N\lambda /2}{\mathrm{cosh}(N\frac{1}{2})\lambda }.$$
(13)
The gain is plotted in fig. 6 for different values of $`\lambda `$. The gain approaches 2 when $`\lambda N\stackrel{>}{_{}}3`$, which resembles the voltage gain in the uniform single electron transistor where one needs $`C_gC`$ to get a maximum gain of 2.

Unfortunately this gain behavior near zero gate voltage is destroyed by random background charges. Melsen et. al.  have shown that the mean value (averaged over background charge configurations) of the threshold voltage was independent of the gate voltage. Our simulations show that it is still possible to modulate the threshold voltage with gate voltage when background charges are present, however the modulation changes. The large peak and high sensitivity around zero gate voltage disappears, and the largest modulation can occur near any value of $`n_g`$. The maximum modulation ($`V_t|_{max}V_t|_{min}`$) also becomes smaller than in the no-background-charge case. In fig. 7 the maximum modulation (averaged over 1000 background charge configurations) is plotted as a function of number of junctions for different values of $`\lambda `$. As in the no-background-charge case, the modulation becomes larger with increasing $`N`$ and $`\lambda ^1`$. Here again, we see weak dependence of the modulation on $`N`$ when $`N\lambda ^1>1`$. The error bars plotted for several points give the standard deviation of the modulation, or the spread in modulation depending on the particular configuration of background charges. Thus, static random background charges cause random shift in the operation point in gate voltage, reduce the maximum modulation with gate voltage, and cause significant variation in the gate modulation of the array.

We have simulated the effects of random background charges on the Coulomb blockade in one dimensional tunnel junction arrays. We find that when $`\lambda ^1\stackrel{>}{_{}}N`$ the static potential along the array is smooth and does not depend very strongly on the density or amplitude of background charges. We derive an analytical expression for the voltage gain of an array transistor, valid when $`n_g`$ is small and no background charges present. However, by simulating the behavior of the threshold voltage as a function of a symmetrically coupled gate voltage with background charges in the array, we conclude that random background charges render arrays useless in transistor applications.

We gratefully acknowledge financial support from the Swedish TFR/NFR and the European Union through MEL ARI project 22953 CHARGE.


format: text

text:
# Fragmentation with a Steady Source

## Abstract

We investigate fragmentation processes with a steady input of fragments. We find that the size distribution approaches a stationary form which exhibits a power law divergence in the small size limit, $`P_{\mathrm{}}(x)x^3`$. This algebraic behavior is robust as it is independent of the details of the input as well as the spatial dimension. The full time dependent behavior is obtained analytically for arbitrary inputs, and is found to exhibit a universal scaling behavior.

PACS numbers: 05.40.+j, 64.60.Ak, 62.20.Mk

Fragmentation underlies numerous natural phenomena. The quantity being “split” can be the mass, momentum, or the area, and typically, fragments continue splitting independently of each other. Examples include polymer degradation , energy cascades in turbulence, breakup of liquid droplets and atomic nuclei, martensitic transformations, shattering of solid objects, and meteor impacts. Fragmentation also arises in several topics of computer science.

The simplest fragmentation models assume that the rate by which fragments are produced is a function of their size only. In this study, we focus on the classic “random scission model” where the cutting is uniform and hence a fragment is cut with a rate proportional to its size. In particular, we are interested in situations where the system is subject to a steady input of fragments. Such “open” systems may be realized physically in impact fragmentation of solid objects where it is possible to constantly supply the system with unfragmented objects.

Fragmentation in open systems has received much less attention than fragmentation in closed systems. We will show, however, that fragmentation with input is actually simpler than the classical counterpart as the system reaches a stationary state which is remarkably robust. Specifically, fragmentation with a steady source is characterized by an algebraic divergence of the size distribution in the small size limit, and this behavior is independent of the particular form of the input. Additionally, the time dependent behavior, obtained analytically for arbitrary inputs, follows a scaling behavior. These two features are shown to be closely related.

We start with a one-dimensional fragmentation process subject to constant input of segments. Here “one-dimensional” means that the fragments are characterized by a single variable which we shall call “length” (the fragments can be viewed as segments). Let the system be initially empty and intervals whose length is within the range $`(x,x+dx)`$ are added with rate $`f(x)dx`$. Additionally, intervals are cut with a constant spatially homogeneous rate; we set this rate equal to unity without loss of generality. Fragmentation with input has a natural geometric interpretation. Consider the segments as part of an infinite line. The fragmentation process is equivalent to deposition of point “cracks” on the line. The line is initially “immune” to fragmentation, but then segments of length $`x`$ become “susceptible” to fragmentation with rate $`f(x)`$. Hence, fragmentation with input is equivalent to inhomogeneous fragmentation of a growing line.

The density $`P(x,t)`$ of intervals of length $`x`$ at time $`t`$ evolves according to the following rate equation

$$\frac{P(x,t)}{t}=xP(x,t)+2_x^{\mathrm{}}𝑑yP(y,t)+f(x).$$
(1)
The negative term on the right-hand side accounts for loss due to fragmentation with the rate equal to the fragment size since the cutting is uniform. The gain term gives the increase in fragments of size $`x`$ due to cutting of longer fragments. The last term accounts for input of fragments of size $`x`$.

The size distribution can be determined by applying the Mellin transformation. The Mellin transform (or moment) of the distribution, $`M(s,t)=𝑑xx^{s1}P(x,t)`$, satisfies

$$\frac{M(s,t)}{t}=\frac{2s}{s}M(s+1,t)+\widehat{f}(s),$$
(2)
where $`\widehat{f}(s)=𝑑xx^{s1}f(x)`$ is the Mellin transform of the input density $`f(x)`$. Although this hierarchy of equations is infinite, its linear nature makes it tractable, as will be seen below.

We first examine what happens when $`t\mathrm{}`$. In this limit, the length density should approach the stationary distribution, $`P(x,t)P_{\mathrm{}}(x)`$. Setting the time derivative in Eq. (2) to zero gives the corresponding transform $`M_{\mathrm{}}(s)=\left(1+\frac{2}{s3}\right)\widehat{f}(s1)`$. Note that $`\widehat{f}(s+n)`$ is the Mellin transform of $`x^nf(x)`$, and $`(sm)^1\widehat{f}(s)`$ is the transform of $`x^m_x𝑑yy^{m1}f(y)`$. These two facts allow to perform the inverse Mellin transform and yield $`P_{\mathrm{}}(x)`$ explicitly in terms of the input function

$$P_{\mathrm{}}(x)=x^1f(x)+2x^3_x^{\mathrm{}}𝑑yyf(y).$$
(3)
In the small size limit, the integral on the right-hand side of Eq. (3) approaches the average length added per unit time, $`\lambda =\widehat{f}(2)=𝑑xxf(x)`$. Thus, the length density becomes purely algebraic

$$P_{\mathrm{}}(x)2\lambda x^3,\mathrm{when}x0.$$
(4)
This behavior is robust as the first term on the right-hand side of Eq. (3) always diverges slower than $`x^3`$ in the limit $`x0`$ (otherwise, the total length input rate, $`𝑑xxf(x)`$, would be infinite). For a class of input densities, the algebraic behavior may not be necessarily restricted to small sizes. For example, for monodisperse inputs $`f(x)=\lambda \delta (x1)`$, the algebraic behavior extends to all sizes $`x<1`$, $`P_{\mathrm{}}(x)=\lambda x^1\delta (x1)+2\lambda x^3`$.

The general algebraic behavior should be contrasted with the exponential length distribution found generally in the absence of input. Algebraic distributions have been observed experimentally in fragmentation of solid objects such as rods, spheres, and bricks . Although the corresponding exponents measured in these experiments are significantly lower, typically between 1 and 2, it is worth noting that a steady source of fragments can serve as a mechanism for generating algebraic distributions. Curiously, algebraic distributions with an exponent close to $`3`$ were reported recently in social systems which can be viewed as open ones (distributions of citations, of the number of links to sites on the internet, etc.; see e.g. Refs.).

The limiting size distribution is ultimately related to the time dependent behavior. This can be demonstrated using the following heuristic argument. From Eq. (2), the total length $`L(t)=M(2,t)`$ grows linearly with time $`\dot{L}(t)=\lambda `$, and hence, $`L(t)=\lambda t`$. Similarly, the total number of fragments, $`N(t)=M(1,t)`$ satisfies $`\dot{N}(t)=\lambda t+\mu `$, where $`\mu =\widehat{f}(1)=𝑑xf(x)`$ is the number of segments added per unit time, and consequently, $`N(t)=\frac{1}{2}\lambda t^2+\mu t`$. These two time dependent results imply that the average fragment length, $`x=L/N`$, decreases with time according to $`xt^1`$. For the length distribution to follow a scaling form, the corresponding scaling variable must be $`x/x`$. The prefactor is fixed by the total number of fragments, $`N(t)\lambda t^2`$, so the scaling form reads

$$P(x,t)\lambda t^3F(xt).$$
(5)
This scaling form would be consistent with a time independent limiting distribution only when $`F(\xi )\xi ^3`$, thereby implying the algebraic divergence (4).

The full time dependent solution can be found using the Charlesby method . This method starts with a formal expansion of the Mellin transform,

$$M(s,t)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{t^n}{n!}M_n(s),$$
(6)
and proceeds by solving for the functions $`M_n(s)`$ iteratively. Indeed, substituting the expansion (6) into Eq. (2) and equating similar powers of time yields $`M_0(s)=0`$, $`M_1(s)=\widehat{f}(s)`$ and $`M_{n+1}(s)=\frac{s2}{s}M_n(s+1)`$ for $`n2`$. Solving this set of equations recursively gives

$$M_{n+1}(s)=(1)^n\frac{(s1)(s2)}{(s+n1)(s+n2)}\widehat{f}(s+n).$$

To take advantage of the inversion rules used to obtain Eq. (3), we re-write $`M_{n+1}(s)`$ as

$$M_{n+1}(s)=(1)^n\left[1\frac{n(n+1)}{s+n1}+\frac{n(n1)}{s+n2}\right]\widehat{f}(s+n).$$

From Eq. (6), the size distribution can be written as a power series

$$P(x,t)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{t^{n+1}(x)^n}{(n+1)!}P_n(x),$$
(7)
where the inverse transform of $`M_{n+1}(s)`$ has been conveniently written as $`(x)^nP_n(x)`$. The three terms in the above expression for $`M_{n+1}(s)`$ can be inverted using the rules outlined before Eq. (3). The final expression for $`P_n(x)`$ reads

$$P_n(x)=f_1(x)+\frac{n(n+1)}{x}f_2(x)+\frac{n(n1)}{x^2}f_3(x),$$
(8)
with

$`f_1(x)`$ $`=`$ $`f(x),`$ (9)
$`f_2(x)`$ $`=`$ $`{\displaystyle _x^{\mathrm{}}}𝑑yf(y),`$ (10)
$`f_3(x)`$ $`=`$ $`{\displaystyle _x^{\mathrm{}}}𝑑yyf(y).`$ (11)
Summing the three terms separately gives the fragment size distribution

$$P(x,t)=\underset{k=1}{\overset{3}{}}t^kf_k(x)F_k(xt),$$
(12)
with the scaling functions

$`F_1(z)`$ $`=`$ $`z^1(1e^z),`$ (13)
$`F_2(z)`$ $`=`$ $`e^z,`$ (14)
$`F_3(z)`$ $`=`$ $`z^3\left[2(2+2z+z^2)e^z\right].`$ (15)
The function $`F_3(z)`$ has been obtained from the power series $`F_3(z)=_{n0}\frac{(z)^n}{n!(n+3)}`$. One can verify that the previous results for $`N(t)`$, $`L(t)`$, and $`P_{\mathrm{}}(x)`$ agree with this solution. Thus, we have obtained the full time dependent solution for an arbitrary time independent input $`f(x)`$.

The size distribution of Eq. (12) exhibits scaling. Indeed, in the limit $`t\mathrm{}`$, $`x0`$ with the scaling variable $`z=xt`$ kept finite, the third term in the sum on the right-hand side of Eq. (12) dominates, and the anticipated scaling behavior of Eq. (5) is confirmed with $`F(z)=F_3(z)`$. Interestingly, the only parameter relevant asymptotically is the overall length input rate $`\lambda `$.

The limiting behaviors of the scaling distribution are

$$F(z)\{\begin{array}{cc}\frac{1}{3}\frac{1}{4}z\hfill & z1\text{,}\hfill \\ 2z^3\hfill & z1\text{.}\hfill \end{array}$$
(16)
In particular, the large $`z`$ behavior implies the correct asymptotic $`P_{\mathrm{}}(x)x^3`$, in agreement with Eq. (4). Thus, for sufficiently large fragments, $`xt^1`$, the distribution has already reached the final limiting form, while smaller sizes are still created.

The formal solution (12) has an interesting “staircase” structure, a time power series whose terms are products of time independent functions $`f_k(x)`$ and time dependent functions $`F_k(xt)`$. In fact, the solution for the random scission model in the absence of input is also characterized by a similar structure. Indeed, consider the evolution equation

$$\frac{\stackrel{~}{P}(x,t)}{t}=x\stackrel{~}{P}(x,t)+2_x^{\mathrm{}}𝑑y\stackrel{~}{P}(y,t)$$
(17)
corresponding to the above fragmentation process in the absence of input. Given the initial conditions $`\stackrel{~}{P}(x,0)=\stackrel{~}{f}(x)`$, the solution can be obtained following the same steps that led to Eq. (12). Again, the full time dependent solution is a three term expansion:

$$\stackrel{~}{P}(x,t)=\underset{k=1}{\overset{3}{}}t^{k1}\stackrel{~}{f}_k(x)\stackrel{~}{F}_k(xt).$$
(18)
The time independent functions are given by the same expressions (9) as in the input case (with $`f_k`$ replaced by $`\stackrel{~}{f}_k`$), while the time dependent functions are different from (13): $`\stackrel{~}{F}_1(z)=\stackrel{~}{F}_3(z)=e^z`$, and $`\stackrel{~}{F}_2(z)=(2z)e^z`$. In the limit $`t\mathrm{}`$, $`x0`$ with $`z=xt`$ kept finite, the scaling behavior emerges again. Specifically, $`\stackrel{~}{P}(x,t)\stackrel{~}{\lambda }t^2\stackrel{~}{F}(z)`$ with $`\stackrel{~}{\lambda }=𝑑xx\stackrel{~}{f}(x)`$ and the exponential scaling function $`\stackrel{~}{F}(z)=e^z`$.

The above solutions for the input case with empty initial conditions and no input case with arbitrary initial conditions can be used to construct the general solution for Eq. (1). Indeed, the sum of the solutions (12) and (18), $`P(x,t)+\stackrel{~}{P}(x,t)`$, is the solution for a fragmentation process with input $`f(x)`$ starting from an initial distribution $`\stackrel{~}{f}(x)`$. As expected, the initial conditions are “forgotten” in the long time limit as $`P(x,t)`$ given by Eq. (12) dominates over $`\stackrel{~}{P}(x,t)`$ given by Eq. (18). In particular, the scaling solution (5) is recovered, and the $`P_{\mathrm{}}(x)2\lambda x^3`$ divergence of the limiting distribution holds in general.

To examine the robustness of the algebraic behavior above, we consider a natural generalization to $`d`$ spatial dimensions . Given that the most interesting long time behavior is independent of the details of the source term, we focus on the simplest monodisperse inputs, namely, unit hypercubes. For instance, in two dimensions we add unit squares with rate $`\lambda `$. A unit square is divided by choosing a point $`(x_1,x_2)`$ with a uniform probability density, and cutting the original square into four rectangles of sizes $`x_1\times x_2`$, $`x_1\times (1x_2)`$, $`(1x_1)\times x_2`$, and $`(1x_1)\times (1x_2)`$. Similarly, the process is repeated with rectangular fragments.

Let $`P(𝐱,t)`$ with $`𝐱(x_1,\mathrm{},x_d)`$ be the distribution of fragments of size $`x_1\times \mathrm{}\times x_d`$ at time $`t`$. This quantity evolves according to the rate equation

$$\left(\frac{}{t}+|𝐱|\right)P(𝐱,t)=2^d_𝐱𝑑𝐲P(𝐲,t)+\lambda \delta (𝐱\mathrm{𝟏}).$$
(19)
Here, we used the shorthand notations $`\mathrm{𝟏}=(1,\mathrm{},1)`$ and $`|𝐱|=x_1\mathrm{}x_d`$. The $`d`$-dimensional Mellin transform, $`M(𝐬,t)=𝑑𝐱x_1^{s_11}\mathrm{}x_d^{s_d1}P(𝐱,t)`$ with $`𝐬(s_1,\mathrm{},s_d)`$, reduces Eq. (19) to

$$\frac{M(𝐬,t)}{t}=\left(\frac{2^ds_1\mathrm{}s_d}{s_1\mathrm{}s_d}\right)M(𝐬+\mathrm{𝟏},t)+\lambda .$$
(20)
We focus on the limiting size distribution $`P_{\mathrm{}}(𝐱)`$. Its Mellin transform $`M_{\mathrm{}}(𝐬)`$ is found from Eq. (20) by setting the time derivative to zero. One gets $`M_{\mathrm{}}(𝐬)=\lambda \left(1+\frac{2^d}{(s_11)\mathrm{}(s_d1)2^d}\right)`$. Inverting this relation yields

$$P_{\mathrm{}}(𝐱)=\lambda \left[\delta (𝐱\mathrm{𝟏})+2^d|𝐱|^1\mathrm{\Phi }_d(\xi )\right],$$
(21)
with the shorthand notations

$`\mathrm{\Phi }_d(\xi )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\xi ^n}{n!}}\right)^d\mathrm{and}\xi =2\left({\displaystyle \underset{i=1}{\overset{d}{}}}\mathrm{ln}{\displaystyle \frac{1}{x_i}}\right)^{1/d}.`$ (22)
In one dimension, $`\mathrm{\Phi }_1(\xi )=e^\xi =x^2`$, and we recover the one-dimensional result $`P_{\mathrm{}}(x)=2\lambda x^3`$. In two dimensions, $`\mathrm{\Phi }_2(\xi )=I_0(2\xi )`$ where $`I_0`$ is the modified Bessel function, and in general $`\mathrm{\Phi }_d(\xi )`$ can be expressed in terms of hypergeometric functions.

The small size behavior of $`P_{\mathrm{}}(𝐱)`$ can be obtained by using the steepest decent method. The leading tail behavior, $`\mathrm{\Phi }_d(\xi )(2\pi \xi )^{\frac{1d}{2}}e^{\xi d}`$ for $`\xi 1`$, corresponds to the case when at least one of the lengths is small, i.e., $`x_i1`$. Returning to the original variables, we re-write the above asymptotic as

$`P_{\mathrm{}}(𝐱)|𝐱|^1|\mathrm{ln}𝐱|^{\frac{d1}{2d}}\mathrm{exp}\left[2d\left(|\mathrm{ln}𝐱|\right)^{1/d}\right],`$
where $`|\mathrm{ln}𝐱|_{i=1}^d\mathrm{ln}\frac{1}{x_i}`$. Thus, the fragment distribution exhibits a “log-stretched-exponential” behavior.

Let us consider the limiting volume distribution $`P_{\mathrm{}}(V)`$ defined via $`P_{\mathrm{}}(V)=𝑑𝐱P_{\mathrm{}}(𝐱)\delta \left(Vx_1\mathrm{}x_d\right)`$. Its Mellin transform, $`M_{\mathrm{}}(s)=𝑑VV^{s1}P_{\mathrm{}}(V)`$, immediately follows from Eq. (20): $`M_{\mathrm{}}(s)=\lambda \left[1+\frac{2^d}{(s1)^d2^d}\right]`$. Using the identity $`(a^d1)^1=d^1_{k=0}^{d1}\zeta ^k(a\zeta ^k)^1`$, where $`\zeta =e^{2\pi i/d}`$ is the primitive $`d^{\mathrm{th}}`$ root of unity, we can express $`M_{\mathrm{}}(s)`$ as a sum over simple poles at $`1+2\zeta ^k`$. Consequently, the inverse Mellin transform is given by a linear combination of $`d`$ power laws,

$$P_{\mathrm{}}(V)=\lambda \left[\delta (V1)+\frac{2}{d}\underset{k=0}{\overset{d1}{}}\zeta ^kV^{12\zeta ^k}\right].$$
(23)
One can verify that the volume distribution is real since it equals its complex conjugate. The small-volume tail of the distribution can be obtained by noting that the sum in Eq. (23) is dominated by the first term in the series,

$$P_{\mathrm{}}(V)\frac{2\lambda }{d}V^3,\mathrm{for}V0.$$
(24)
Thus, the same $`V^3`$ algebraic behavior occurs in all spatial dimensions. Clearly, this divergence is general. Indeed, the time dependent evolution equation (19) implies that the overall volume grows linearly, and that the overall number of fragments grows quadratically. Therefore, the heuristic scaling argument leading to Eq. (5) extends to higher dimensions, and consequently, the limiting behavior is given by Eq. (24).

In summary, we have studied random fragmentation processes in the presence of a steady source. We have solved for the full time dependent behavior in terms of the input function. In the long time limit, the size distribution exhibits a universal scaling behavior. The limiting distribution diverges algebraically according to $`x^3`$ in the small size limit. This behavior is robust. It applies to arbitrary inputs, and it extends to higher dimensions as well. Interestingly, the only asymptotically relevant parameter is the total volume added per unit time. Additionally, we have shown that the scaling behavior can be used to predict the algebraic nature of the final size distribution. Hence, the scaling behavior and the limiting distribution are closely related.

The solution for the time dependent behavior exhibits an interesting staircase structure. The two progressively weaker corrections to the leading behavior are of the order $`t^1`$ and $`t^2`$, respectively. Such staircase structures may be a useful tool for treating similar integro-differential equations which are expected to exhibit a scaling asymptotic behavior. For instance, one can check whether substituting such an ansatz leads to a closed system of equations for the time dependent and time independent functions.

The above results can be extended in a number of ways. One may try to derive a general solution for $`P(𝐱,t)`$ in higher dimensions for arbitrary input rate $`f(𝐱)`$. We anticipate that geometric features of the fragments will be interesting. Indeed, in the no input case the volume distribution exhibits an ordinary scaling behavior while multiscaling asymptotic behavior underlies the full multivariate size distribution. Additionally, a nontrivial set of conservation laws exists as all moments $`M(𝐬^{},t)`$ with $`_{i=1}^ds_i^{}=2^d`$ are conserved, as seen from Eq. (20). In the presence of input, the same moments grow linearly in time but multiscaling should still hold asymptotically.

We acknowledge support from DOE (W-7405-ENG-36), NSF (DMR9978902), ARO (DAAD19-99-1-0173), and CNRS.


text:
# Ginzburg-Landau equations for superconducting quark matter in neutron stars

## I Introduction

Recently, the possible formation of diquark condensates in QCD at finite density has been reinvestigated in a series of papers following Refs. . It has been shown that in chiral quark models with a nonperturbative 4-point interaction motivated from instantons  or nonperturbative gluon propagators  the anomalous quark pair amplitudes in the color antitriplet channel can be very large: of the order of $`100\mathrm{MeV}`$. Therefore, in two-flavor QCD, one expects this diquark condensate to dominate the physics at densities beyond the deconfinement/chiral restoration transition and below the critical temperature ($`50\mathrm{MeV}`$) for the occurence of this two-flavor color superconductivity (2SC) phase. In a three-flavor theory it has been found  that there can exist a color-flavor locked (CFL) phase for not too large strange quark masses  where color superconductivity is complete in the sense that diquark condensation produces a gap for quarks of all three colors and flavors, which is of the same order of magnitude as that in the two-flavor case.

The high-density phases of QCD at low temperatures are most relevant for the explanation of phenomena in rotating compact stars - pulsars. Conversely, the physical properties of these objects (as far as they are measured) could constrain our hypotheses about the state of matter at the extremes of densities. In contrast to the situation for the cooling behaviour of compact stars where the CFL phase is dramatically different from the 2SC phase , we don’t expect qualitative changes of the magnetic field structure between these two phases. Consequently, we will restrict ourselves here to the discussion of the simpler two-flavor theory first. The 2SC phase window is expected to occur at densities just above the deconfinement transition inside a compact star where the strange quark flavor is either still confined  or the self consistently determined strange quark masses are large enough to entail color-flavor unlocking .

According to Bailin and Love  the magnetic field of pulsars should be expelled from the superconducting interior of the star due to the Meissner effect and decay subsequently within $`10^4`$ years. For their estimate, they used a perturbative gluon propagator which yielded a very small pairing gap and they made the assumption of a homogeneous magnetic field. Since both assumptions seem not to be valid in general, we have performed in Ref.  a reinvestigation of the question with the result that 2SC quark matter is a type-II superconductor which can form vortices in response to the external magnetic field thus behaving similar to protons in the superfluid neutron phase of a neutron star . The magnetic field in the quark core of the neutron star can then exist for periods much longer than the spin-down age of a pulsar.

The authors of Ref.  have considered this question taking into account the rotated electromagnetism. They came to the conclusion that the magnetic field will live in the quark core sufficiently long although it does not form the structure of quantum vortices because it obeys the force-free Maxwell equation.

In the present work we derive the Ginzburg-Landau equations with account of rotated electromagnetism and solve these equations for the example of an isolated quantum vortex in quark matter which demonstrates that there is always an electric and color Meissner current when the quark matter is superconducting. This conclusion holds also for the case of a homogeneous external magnetic field as we have shown in a separate paper .

## II Ginzburg-Landau Equations

In the paper  we have obtained the Ginzburg-Landau equations for relativistic superconducting quarks, supposing that the superconducting quark matter phase with $`ud`$ diquark pairing has $`P_p=0`$, where $`p`$ is the color antitriplet index. We study the consequences of superconducting quark cores in neutron stars for the magnetic field of pulsars. We find that within recent nonperturbative approach to the effective quark interaction the diquark condensate forms a type II superconductor whereas previously quark matter was considered as a type I superconductor . In both cases the magnetic field which is generated in the superconducting hadronic shell of superfluid neutrons and superconducting protons can penetrate into the quark matter core since it is concentrated in proton vortex clusters where the field strength exceeds the critical value. Recently, in the paper , discussing the Meissner effect for color superconducting quark matter, the authors introduce the $`\widehat{Q}`$\- charge generator

$$\widehat{Q}=Q+\eta P_8,$$
(1)
where $`Q`$ is the conventional electromagnetic charge generator and $`P_8`$ is associated with the one of the gluons in the representation of the quarks

$`Q`$ $`=`$ $`\text{diag}({\displaystyle \frac{2}{3}},{\displaystyle \frac{1}{3}},{\displaystyle \frac{1}{3}})\text{in flavor}u,d,s\text{space},`$ (2)
$`P_8`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}\text{diag}(1,1,2)\text{in color}r,g,b\text{space}.`$ (3)
The condition they requested was that the $`\widehat{Q}`$\- charge of all Cooper pairs which form the condensate vanishes

$$\widehat{Q}d_p=0.$$
(4)
Applying this condition (instead of $`P_8=0`$) for the derivation of the Ginzburg-Landau equations, we will see that in the superconducting quark matter two different types of magnetic fields can exist simultaneously. The penetration depths $`\lambda _q`$ for these two magnetic fields are completely different. For one of them the penetration depth is finite and much shorter than has been found in , but for the other one it becomes infinitely large.

The free energy density in the superconducting quark matter with $`ud`$ diquark pairing can be written in the following form

$`f`$ $`=`$ $`f_n+ad_pd_p^{}+{\displaystyle \frac{1}{2}}\beta (d_pd_p^{})^2+\gamma |\widehat{P}d_p|^2+{\displaystyle \frac{(\mathrm{rot}\stackrel{}{A})^2}{8\pi }}+{\displaystyle \frac{(\mathrm{rot}\stackrel{}{G}_8)^2}{8\pi }},`$ (5)
where $`a=dn/dEt`$, $`\beta =dn/dE7\zeta (3)(\pi k_BT_c)^2/8`$, $`\gamma =p_F^2\beta /(6\mu ^2)`$ and $`t=(TT_c)/T_c`$ with $`T_c`$ being the critical temperature, $`p_F`$ the quark Fermi momentum. In zeroth order with respect to the coupling constant, $`dn/dE=\mu p_F/\pi ^2`$. $`|\widehat{P}d_p|^2`$ is the kinetic term of the free energy density (5) which is given by

$$|\widehat{P}d_p|^2=|(\stackrel{}{}+ie\stackrel{}{A}Q+ig\stackrel{}{G}_8P_8)d_p|^2.$$
(6)
If we take into account the condition (4) and use (1) we will have

$$Qd=\eta P_8d.$$
(7)
Inserting (7) into (6) we will get

$$|\widehat{P}d_p|^2=|(\stackrel{}{}+ie\eta \stackrel{}{A}P_8+ig\stackrel{}{G}_8P_8)d_p|^2.$$
(8)
Following the paper , let us introduce instead of the original gauge fields $`\stackrel{}{A}`$ and $`\stackrel{}{G}_8`$ the linear combinations $`\stackrel{}{A}_x`$ and $`\stackrel{}{A}_y`$

$`\stackrel{}{A}_x`$ $`=`$ $`{\displaystyle \frac{\eta e\stackrel{}{A}+g\stackrel{}{G}_8}{\sqrt{\eta ^2e^2+g^2}}}=\mathrm{sin}\alpha \stackrel{}{A}+\mathrm{cos}\alpha \stackrel{}{G}_8,`$ (9)
$`\stackrel{}{A}_y`$ $`=`$ $`{\displaystyle \frac{g\stackrel{}{A}+\eta e\stackrel{}{G}_8}{\sqrt{\eta ^2e^2+g^2}}}=\mathrm{cos}\alpha \stackrel{}{A}+\mathrm{sin}\alpha \stackrel{}{G}_8.`$ (10)
This corresponds to a rotation with the angle $`\alpha `$ in the orthogonal basis of the fields $`\stackrel{}{A}_x`$ and $`\stackrel{}{A}_y`$, where

$`\mathrm{cos}\alpha ={\displaystyle \frac{g}{\sqrt{\eta ^2e^2+g^2}}},`$ (11)
while the normalization is preserved

$$\stackrel{}{A}_x^2+\stackrel{}{A}_y^2=\stackrel{}{A}^2+\stackrel{}{G}_8^2.$$
(12)
At neutron star densities the gluons are strongly coupled ($`g^2/4\pi 1`$) and the photons are of course weakly coupled ($`e^2/4\pi 1/137`$), so that $`\alpha \eta e/g`$ is small. For the diquark condensate, where blue-green and green-blue $`ud`$ quarks are paired, $`\eta =1/\sqrt{3}`$ and therefore $`\alpha 1/20`$. Taking into account (9) and (10), the kinetic term (8) will become

$$|\widehat{P}d_p|^2=(\stackrel{}{}iq\eta \stackrel{}{A}_x)d_p^{}(\stackrel{}{}+iq\eta \stackrel{}{A}_x)d_p,$$
(13)
where

$$q=\sqrt{\eta ^2e^2+g^2}P_8,$$
(14)
and $`P_8=1/\sqrt{3}`$. We see that the kinetic term has the conventional expression but instead of the charge of the $`ud`$ diquark pair equal to $`e/3`$ we have the new charge $`q`$ which is much larger than $`e/3`$ (about 20 times). The Ginzburg-Landau equations are obtained in the usual way. If we demand that the variation of the free energy with respect to the parameters $`d^{}`$, $`\stackrel{}{A}_x`$ and $`\stackrel{}{A}_y`$ have to be equal to zero, we will get

$`0`$ $`=`$ $`ad_p+\beta (d_pd_p^{})d_p+\gamma (i\stackrel{}{}+qA_x)^2d_p,`$ (15)
the equation of motion for the diquark condensate and for the gauge fields

$`0`$ $`=`$ $`\mathrm{sin}\alpha \mathrm{rot}\mathrm{rot}\stackrel{}{A}\mathrm{cos}\alpha \mathrm{rot}\mathrm{rot}\stackrel{}{G}_84\pi iq\gamma [d\stackrel{}{}d^{}d^{}\stackrel{}{}d]8\pi q^2\gamma |d|^2\stackrel{}{A}_x,`$ (16)
$`0`$ $`=`$ $`\mathrm{cos}\alpha \mathrm{rot}\mathrm{rot}\stackrel{}{A}+\mathrm{sin}\alpha \mathrm{rot}\mathrm{rot}\stackrel{}{G}_8.`$ (17)
These equations can be written in the following form

$`\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{A}+\mathrm{sin}^2\alpha \stackrel{}{A}`$ $`=`$ $`i{\displaystyle \frac{\mathrm{sin}\alpha }{2q}}{\displaystyle \frac{(d_p\stackrel{}{}d_p^{}d_p^{}\stackrel{}{}d_p)}{|d|^2}}+\mathrm{sin}\alpha \mathrm{cos}\alpha \stackrel{}{G}_8`$ (18)
$`\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{G}_8+\mathrm{cos}^2\alpha \stackrel{}{G}_8`$ $`=`$ $`i{\displaystyle \frac{\mathrm{cos}\alpha }{2q}}{\displaystyle \frac{(d_p\stackrel{}{}d_p^{}d_p^{}\stackrel{}{}d_p)}{|d|^2}}+\mathrm{sin}\alpha \mathrm{cos}\alpha \stackrel{}{A},`$ (19)
where

$$\lambda _q^1=\sqrt{8\pi \gamma }|d|q=\sqrt{\frac{4q^2tp_F^3}{3\pi \mu }}.$$
(20)
This is the system of Ginzburg-Landau equations for superconducting quark matter which takes into account the “rotated electromagnetism”. When we define in the equations (18) and (19) the order parameter in the following form

$$d_p=|d_p|\mathrm{e}^{i\varphi },$$
(21)
where $`\varphi `$ is the phase of the order parameter, then this system of equations can be written in the following form

$`\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{A}+\mathrm{sin}^2\alpha \stackrel{}{A}`$ $`=`$ $`i{\displaystyle \frac{\mathrm{\Phi }_q\mathrm{sin}\alpha }{2\pi }}\stackrel{}{}\varphi +\mathrm{sin}\alpha \mathrm{cos}\alpha \stackrel{}{G}_8`$ (22)
$`\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{G}_8+\mathrm{cos}^2\alpha \stackrel{}{G}_8`$ $`=`$ $`i{\displaystyle \frac{\mathrm{\Phi }_q\mathrm{cos}\alpha }{2\pi }}\stackrel{}{}\varphi +\mathrm{sin}\alpha \mathrm{cos}\alpha \stackrel{}{A}`$ (23)
Here $`\mathrm{\Phi }_q`$ has the well known form

$$\mathrm{\Phi }_q=\frac{2\pi \mathrm{}c}{q}.$$
(24)
If we take the contour integral on both sides of equation (22) where the contour $`L`$ is chosen as the equatorial ring which limits the quark matter core region and take into account that

$`{\displaystyle _L}\mathrm{rot}\mathrm{rot}\stackrel{}{A}𝑑\stackrel{}{l}`$ $`=`$ $`{\displaystyle _L}\mathrm{rot}\stackrel{}{B}𝑑\stackrel{}{l}=0,`$ (25)
$`{\displaystyle _L}\stackrel{}{G}_8𝑑\stackrel{}{l}`$ $`=`$ $`0,`$ (26)
so that the electrical currents and the vector potential of the gluon field on the surface of the quark core must vanish, then we get

$$_L\stackrel{}{A}𝑑\stackrel{}{l}=_S\mathrm{rot}\stackrel{}{A}𝑑\stackrel{}{S}=_S\stackrel{}{B}𝑑\stackrel{}{S}=\frac{\mathrm{\Phi }_q}{\mathrm{sin}\alpha }N.$$
(27)
This means that the flux of the magnetic field through an arbitrary surface over the contour $`L`$ is equal to the sum of fluxes which are generated by $`N`$ vortices each of which has a flux $`\mathrm{\Phi }_q/\mathrm{sin}\alpha `$ .

When we calculate this flux we obtain

$$\mathrm{\Phi }_q^{}=\frac{6\pi \mathrm{}c}{e}=6\mathrm{\Phi }_0,$$
(28)
where $`\mathrm{\Phi }_0=2\times 10^7`$ G cm2 is the magnetic flux quantum. In order to obtain this result we have used equation (22). It is easy to see that the same result can be obtained using equation (23), as to be expected.

## III Solution of Ginzburg-Landau equations for isolated vortex in quark matter

We now consider the solution of the Ginzburg-Landau equations (22, 23) for an isolated vortex situated in homogeneous superconducting quark matter. To this end it is sufficient to find the solution of these equations without the vortex terms. The contribution of an isolated vortex is taken into account by the condition (27) where $`N=1`$. Consequently we have to solve the following system of equations

$`\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{A}+\mathrm{sin}^2\alpha \stackrel{}{A}`$ $`=`$ $`\mathrm{sin}\alpha \mathrm{cos}\alpha \stackrel{}{G}_8`$ (29)
$`\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{G}_8+\mathrm{cos}^2\alpha \stackrel{}{G}_8`$ $`=`$ $`\mathrm{sin}\alpha \mathrm{cos}\alpha \stackrel{}{A}.`$ (30)
From these equations it is easy to see that

$$\mathrm{rot}\mathrm{rot}\stackrel{}{G}_8=\mathrm{cot}\alpha \mathrm{rot}\mathrm{rot}\stackrel{}{A}.$$
(31)
When we act on both sides of equation (29) with the operator $`\mathrm{rot}\mathrm{rot}`$ and take into account that $`\alpha =\mathrm{const}`$, we obtain

$$\lambda _q^2\mathrm{rot}\mathrm{rot}\mathrm{rot}\mathrm{rot}\stackrel{}{A}+\mathrm{sin}^2\alpha \mathrm{rot}\mathrm{rot}\stackrel{}{A}=\mathrm{sin}\alpha \mathrm{cos}\alpha \mathrm{rot}\mathrm{rot}\stackrel{}{G}_8.$$
(32)
Inserting (31) into (32), we obtain the equation for the determination of the electromagnetic vector potential

$$\lambda _q^2\mathrm{rot}\mathrm{rot}\mathrm{rot}\mathrm{rot}\stackrel{}{A}+\mathrm{rot}\mathrm{rot}\stackrel{}{A}=0.$$
(33)
Let us denote

$$\mathrm{rot}\mathrm{rot}\stackrel{}{A}=\stackrel{}{M}.$$
(34)
Then, equation (31) takes the form

$$\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{M}+\stackrel{}{M}=0.$$
(35)
Consequently, instead of equation (33) we can solve the system of equations (34) and (35).

Let us consider the solution of equation (35). Since our problem has cylindric symmetry, the unknown function shall depend only on the coordinate $`r`$ denoting the distance to the vortex. It is easy to see that the vectors $`\stackrel{}{M}`$ and $`\stackrel{}{A}`$ have only azumutal components $`\stackrel{}{M}_\phi (r)`$ and $`\stackrel{}{A}_\phi (r)`$. Then equation (35) takes the form

$$\frac{d^2M_\phi (r)}{dr^2}+\frac{1}{r}\frac{dM_\phi (r)}{dr}(\frac{1}{r^2}+\frac{1}{\lambda _q^2})M_\phi (r)=0.$$
(36)
The general solution of this equation which tends to zero at infinity is given by

$$M_\phi (r)=cK_1(r/\lambda _q),$$
(37)
where $`K_1`$ is the modified Bessel function of first kind. Next we solve the equation (34). It is easy to see that a special solution of this equation is $`\stackrel{}{A}=\lambda _q^2\stackrel{}{M}`$, and the general solution which vanishes at $`r\mathrm{}`$ will be $`c_2/r`$ . Finally for the vector potential we obtain

$$A_\phi (r)=c_1K_1(r/\lambda _q)+c_2/r.$$
(38)
The nonvanishing component of the magnetic field has $`z`$\- direction and is determined by the formula

$$B_z(r)=\frac{1}{r}\frac{d}{dr}\left(rA_\varphi (r)\right).$$
(39)
Inserting (38) into (39) we finally obtain

$$B_z(r)=c_1K_0(r/\lambda _q),$$
(40)
where $`K_0`$ is the modified Bessel function of zero kind. The constant $`c_1`$ is determined from the equation (27) as

$$c_1=\frac{\mathrm{\Phi }_q^{}}{2\pi \lambda _q^2}.$$
(41)
Then we have

$$B_z(r)=\frac{\mathrm{\Phi }_q^{}}{2\pi \lambda _q^2}K_0\left(\frac{r}{\lambda _q}\right).$$
(42)
If we introduce the gluomagnetic field $`\stackrel{}{B}_8=\mathrm{rot}\stackrel{}{G}_8`$ , then we can determine it using equation (31). The condition that the magnetic fields $`\stackrel{}{B}`$ and $`\stackrel{}{B}_8`$ have to vanish at $`r\mathrm{}`$, immediately determine the solution for $`\stackrel{}{B}_8`$

$$\stackrel{}{B}_8=\stackrel{}{B}\mathrm{cot}\alpha ,$$
(43)
where $`\stackrel{}{B}`$ is given by equation (42).

As can be seen from the solution (42) the field of the quantum vortex in quark matter formally has the same form as the field of a vortex in an ordinary superconductor. Whereas in the case of the ordinary superconductor in the definition of the flux $`\mathrm{\Phi }`$ and the penetration depth $`\lambda `$ enters the same charge of the Cooper pair (i.e. $`2\mathrm{e}`$, where $`\mathrm{e}`$ is the charge of the electron), in quark matter in the definition of $`\mathrm{\Phi }_q`$ enters the charge of the quark Cooper pair $`\mathrm{e}/3`$, and therefore $`\mathrm{\Phi }_q^{}=6\mathrm{\Phi }_0`$, and in the definition of the penetration depth $`\lambda _q`$ enters the charge $`q`$, which is about $`20`$ times larger than $`\mathrm{e}/3`$. Therefore, the penetration depth $`\lambda _q`$ which is proportional to $`q^1`$ is $`20`$ times smaller than that of hadronic matter. It is obvious from equation (43) that the magnetic quantum vortex is accompanied by a gluomagnetic quantum vortex the center of which coincides with that of the former, but its force lines have the opposite direction to the magnetic ones. The flux of the gluomagnetic vortex is greater than that of the magnetic vortex by a factor $`\mathrm{cot}\alpha `$.

Let us consider the equations for the rotated magnetic fields $`\stackrel{}{B}_x=\mathrm{rot}\stackrel{}{A}_x`$ and $`\stackrel{}{B}_y=\mathrm{rot}\stackrel{}{A}_y`$. which we obtain from the equations (16), (17)

$`\lambda _q^2\mathrm{rot}\mathrm{rot}\stackrel{}{B}_x+\stackrel{}{B}_x`$ $`=`$ $`\mathrm{\Phi }_q\widehat{e}\delta (\stackrel{}{r}),`$ (44)
$`\mathrm{rot}\stackrel{}{B}_y`$ $`=`$ $`0.`$ (45)
The energy of isolated quark vortex is given by

$$E=\frac{1}{8\pi }[\stackrel{}{B}_y^2+\stackrel{}{B}_x^2+(\lambda _q\mathrm{rot}\stackrel{}{B}_x)^2]𝑑V$$
(46)
Note that from equations (10) and (43) follows that the magnetic field $`\stackrel{}{B}_y`$ of the isolated quark vortex is equal to zero. The magnetic field $`\stackrel{}{B}_x`$ will be found from the solution of equation (44), which coincides with the well known equation for the magnetic field of an isolated vortex in ordinary superconductors. We can use now these solutions and calculate the energy of an isolated quark vortex in the form

$$E=\left(\frac{\mathrm{\Phi }_q}{8\pi \lambda _q}\right)^2\mathrm{ln}\frac{\lambda _q}{\xi _q}$$
(47)
As we see from this expression the energy of quark vortex depends logarithmically on the renormalized charge $`q`$.

In concluding this section let us add the following: If the quark matter occupies a spherical volume of radius $`a`$, and the vortex passes through the center of this volume then the force lines of the magnetic vortex outside of this volume appear to have the dipole form, whereas the force lines for the gluomagnetic field have to be confined within this volume because of the confinement condition $`\stackrel{}{G}_8(a)=0`$. Outside of the quark matter region we have only the ordinary magnetic field. We also note that we can use these results in the calculation of the distribution of magnetic and gluomagnetic fields in a neutron star if its core consists of superconducting quark matter.

## IV Conclusion

In this article we derived the Ginzburg-Landau equations with account for rotated electromagnetism. From the form of these equations written for the physical fields $`\stackrel{}{B}`$ and $`\stackrel{}{B}_8`$ it is obvious that they are coupled, i.e. the appearance of one of them leads to the immediate appearance of the other. We have solved these equations for an isolated quantum vortex in quark matter and have shown that the centers of electromagnetic and gluomagnetic vortices coincide. The distribution around the continuation of the center of the vortex in quark matter for the fields $`\stackrel{}{B}`$ and $`\stackrel{}{B}_8`$ has the same form as for ordinary superconductors. Nevertheless, the flux of the magnetic vortex is 6 times the magnetic flux quantum $`\mathrm{\Phi }_0`$ and the penetration depth in quark matter is $`20`$ times less than that in hadronic matter. Only the force lines of the magnetic field can leave the quark matter volume whereas the gluomagnetic field is confined therein. As we have demonstrated in this paper for the case of an isolated vortex in quark matter the occurence of electric and color Meissner currents is a consequence of the color-superconducting state of quark matter.

## Acknowledgement

We thank M. Alford, K. Rajagopal, R. Rapp, E. Shuryak, D. Voskresensky and J. Wambach for discussions. We acknowledge financial support from the DAAD for the exchange program between the Universities of Rostock and Yerevan and from the ECT\* Trento for our participation at the workshop: Physics of Neutron Star Interiors. D.S. is grateful to the Observatoire de Paris at Meudon for its hospitality during a research visit; D.B. thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the programs INT-00-1: QCD at Nonzero Baryon Density and INT-01-2: Neutron Stars.