text:
# Spontaneous phase oscillation induced by inertia and time delay

\[

## Abstract

We consider a system of coupled oscillators with finite inertia and time-delayed interaction, and investigate the interplay between inertia and delay both analytically and numerically. The phase velocity of the system is examined; revealed in numerical simulations is emergence of spontaneous phase oscillation without external driving, which turns out to be in good agreement with analytical results derived in the strong-coupling limit. Such self-oscillation is found to suppress synchronization and its frequency is observed to decrease with inertia and delay. We obtain the phase diagram, which displays oscillatory and stationary phases in the appropriate regions of the parameters.

preprint: SNUTP 00-035; KIAS-P00074

\]

Systems of coupled nonlinear oscillators, serving as a prototype model for various oscillatory systems in nature, have been known to exhibit remarkable phenomena of collective synchronization, which have been observed in a variety of physical, biological, and chemical systems . Among them there are systems where finite inertia and time delay are relevant. For example, in a superconducting junction network , the capacitance plays the role of inertia while time delay naturally arises in some physical as well as biological systems, where finite time interval is actually required for the transmission of information . This has motivated recent studies of the effects of finite inertia and of time delay on synchronization: In a system with inertia suppression of synchronization and emergence of hysteresis has been observed , whereas a multitude of coherent states with different synchronization frequencies and suppression of the collective frequency have been reported in systems with delay . Note, however, that the effects of the two have been considered separately and independently of each other. The interplay between the two has not been explored in the most general system with both inertia and time delay. For a single oscillator with inertia, it is known that delayed restoring force can destabilize a fixed point and produce oscillation ; this raises an interesting question as to how such oscillatory behavior affects collective synchronization in a set of globally coupled oscillators with inertia and delay.

This paper is to investigate the interplay between inertia and retarded interaction and their combined effects on synchronization. We consider a system of globally coupled oscillators, each possessing finite inertia and interacting with others via time-delayed interactions. To explore the interplay, we first examine the temporal behavior of the phase velocity, and obtain a nonvanishing ac component in the absence of external periodic driving. Such spontaneous oscillation is found to suppress synchronization, and its frequency observed to diminish with inertia and delay.

We begin with the set of equations of motion governing the dynamics of $`N`$ coupled oscillators, the $`i`$th of which is described by its phase $`\varphi _i`$ $`(i=1,2,\mathrm{},N)`$:

$$\mu \ddot{\varphi _i}(t)+\dot{\varphi _i}(t)+\frac{K}{N}\underset{j=1}{\overset{N}{}}\mathrm{sin}[\varphi _i(t)\varphi _j(t\tau )]=\omega _i,$$
(1)
where $`\mu `$ represents the magnitude of (rotational) inertia relative to damping. The third term on the left-hand side denotes the global coupling of strength $`K/N`$ between oscillators, indicating that each oscillator interacts with others only after the retardation time $`\tau `$. The term $`\omega _i`$ on the right-hand side denotes the intrinsic frequency of the $`i`$th oscillator, and is randomly distributed over the whole oscillators according to $`g(\omega )`$, which is assumed to be symmetric about $`\omega =0`$, and concave at $`\omega =0`$. In the absence of time delay and inertia ($`\tau =\mu =0`$), Eq. (1) exactly reduces to the Kuramoto model , for which analytical results are available. The system with either time delay or inertia (but not both) also exhibits synchronization for the coupling strength strong enough. The phases are stationary or monotonic in time , indicating that the system either relaxes to the minimum-energy configuration or displays constant phase velocity.

Here we investigate the dynamics of phases in the presence of both inertia and delay, first by means of numerical simulations. For convenience, we take the Gaussian distribution with unit variance ($`\sigma ^2=1`$) for $`g(\omega )`$ and choose the coupling strength $`K=3`$, which is larger than the known critical values . We thus probe the interplay between inertia and delay in the synchronized state. Equation (1) has been integrated with discrete time steps of $`\delta t=0.01`$. In computing the order parameter, $`N_t=10^6`$ time steps have been used at each run, with the data from the first $`9.6\times 10^5`$ steps discarded. We have varied both $`\delta t`$ and $`N_t`$ to verify that the steady state has been attained and performed twenty independent runs with different initial configurations, over which averages have been taken. In this manner, we have computed the power spectrum of the phase velocity:

$$S\frac{1}{N}\underset{j=1}{\overset{N}{}}|\psi _j(f)|^2,$$
(2)
where $`\psi _j(f)\dot{\varphi }_j(t)e^{2\pi ift}𝑑t`$ is the Fourier component of the phase velocity at frequency $`f`$.

Figure 1 shows the obtained power spectrum in a system of $`N=100`$ oscillators, for coupling strength $`K=3`$ and time delay $`\tau =2`$ at various inertia values, $`\mu =0.4,\mathrm{\hspace{0.17em}0.6}`$, and $`0.8`$. (We have considered the size $`N`$ up to $`3200`$, and confirmed the absence of appreciable finite-size effects for $`N100`$). For small inertia ($`\mu 0.3`$), the power spectrum possesses only the zero-frequency component. Namely, the system displays stationary or monotonic behavior of the phase, like the one with delay only . When the inertia is raised beyond $`0.3`$, on the other hand, several peaks develop at nonzero fundamental and harmonic frequencies, and the oscillatory phase emerges. As $`\mu `$ is increased further, those peaks shift to the lower frequency sides, eventually disappearing for $`\mu 1.0`$. We have examined such self-oscillatory behavior for various values of $`K`$, $`\tau `$, and $`\mu `$, and observed that the oscillation frequency $`f_0`$ decreases with both inertia and delay, as shown in Fig. 2 for $`\tau =1.8,\mathrm{\hspace{0.17em}2.0}`$, and $`2.2`$.

The resulting phase boundaries in the plane of $`(\mu ,\tau )`$, discriminating the oscillatory state and the stationary one, are displayed in Fig. 3, where the coexistence of the two states is also identified in the region between solid squares and open ones. In the coexistence region, there exist a multitude of solutions, giving stationary or oscillatory behavior depending on the initial conditions. To explore this, we examine the behavior of the order parameter $`\mathrm{\Delta }N^1|_j\mathrm{exp}(i\varphi _j)|`$ with the inertia. Figure 4(a) shows the order parameter computed via numerical simulations, as the inertia is varied between $`0.1`$ and $`1.2`$ with the increment of $`\delta \mu =0.01`$. We have first decreased the inertia from 1.2 to 0.1, which yields the step-like structure consisting of the solid circles in Fig. 4(a). Upon increasing the inertia from given value, however, the order parameter does not reverse the path along the solid circles; instead it exhibits continuous dependence on the inertia, eventually following the open circles. In this manner a multitude of solutions or “bands” have been found, as shown in Fig. 4(a). Note that there are two kinds in shape, flat bands and bent ones. As the behavior of the phase velocity is examined, the system is stationary in the flat bands, existing for small inertia ($`\mu 0.3`$) and for large inertia ($`\mu 1.1`$). In the bent bands (for $`0.3\mu 0.9`$), on the other hand, the system displays phase oscillations. Accordingly, in the range $`0.9\mu 1.1`$, where two kinds of bands coexist, either the oscillatory phase or the stationary one is expected to appear depending on the initial conditions. This coexistence range of inertia (for $`\tau =2`$) is indeed consistent with Fig. 3.

To understand analytically these numerical results, we first divide the population of the oscillators into two groups: the synchronized group $`(S)`$ and the desynchronized one $`(D)`$, and take the ansatz $`\varphi _i=\varphi _i^0+\omega _it`$ for $`D`$ and

$$\varphi _i=\varphi _i^0+A\mathrm{sin}(\mathrm{\Omega }t+\alpha )$$
(3)
for $`S`$. In view of the symmetry of $`g(\omega )`$, we consider the dc component $`\varphi _i^0`$ to be constant. Results of numerical simulations also manifest that near the phase boundary the ac component is much smaller than the dc component, allowing the expansion in terms of the amplitude $`A`$. Upon substituting Eq. (3) into Eq. (1), we thus obtain for the oscillators in $`S`$, to the order of $`A^2`$,

$`\omega _i`$ $`={\displaystyle \frac{K}{N}}{\displaystyle \underset{j}{}^{(S)}}\mathrm{sin}(\varphi _i^0\varphi _j^0)\left(1A^2\mathrm{sin}^2{\displaystyle \frac{\mathrm{\Omega }\tau }{2}}\right)`$ (6)
$`\left[\mu \mathrm{\Omega }^2{\displaystyle \frac{2K}{N}}{\displaystyle \underset{j}{}^{(S)}}\mathrm{cos}(\varphi _i^0\varphi _j^0)\mathrm{sin}^2{\displaystyle \frac{\mathrm{\Omega }\tau }{2}}\right]A\mathrm{sin}(\mathrm{\Omega }t+\alpha )`$
$`+\left[\mathrm{\Omega }+{\displaystyle \frac{K}{N}}{\displaystyle \underset{j}{}^{(S)}}\mathrm{cos}(\varphi _i^0\varphi _j^0)\mathrm{sin}\mathrm{\Omega }\tau \right]A\mathrm{cos}(\mathrm{\Omega }t+\alpha ),`$
where the summations run over $`S`$. Oscillators in $`D`$, the fraction of which over the whole is negligibly small, $`O(\text{e}^{(K\mathrm{\Delta })^2/2\sigma ^2})`$, turn out not to contribute to the harmonics in $`\mathrm{\Omega }`$. A similar expansion for $`D`$ leads to trivial equations.

The zeroth harmonics in Eq. (6) leads to

$$\frac{1}{N}\underset{j}{}^{(S)}\mathrm{sin}(\varphi _i^0\varphi _j^0)=\frac{\omega _i}{K}\left[1+A^2\mathrm{sin}^2(\mathrm{\Omega }\tau /2)\right]$$
(7)
or, in terms of the order parameter,

$$\mathrm{\Delta }^2=1\frac{\sigma ^2}{K^2}\left[1+2A^2\mathrm{sin}^2(\mathrm{\Omega }\tau /2)\right]$$
(8)
where $`\sigma ^2`$ is the variance of the distribution $`g(\omega )`$ and the thermodynamic limit $`(N\mathrm{})`$ together with the strong coupling limit $`(KK_c)`$ has been taken. It is shown in Eq. (8) that the appearance of oscillation $`(A0)`$ reduces $`\mathrm{\Delta }`$, suppressing synchronization.

We next consider the first harmonics in Eq. (6), and obtain

$`K\mathrm{\Delta }^2(1\mathrm{cos}\mathrm{\Omega }\tau )`$ $`=`$ $`\mu \mathrm{\Omega }^2,`$ (9)
$`K\mathrm{\Delta }^2\mathrm{sin}\mathrm{\Omega }\tau `$ $`=`$ $`\mathrm{\Omega },`$ (10)
which yields only the trivial solution $`\mathrm{\Omega }=0`$ in the absence of either inertia or delay. It is thus concluded that the (sustained) spontaneous oscillation requires both inertia and delay. Elimination of $`\mathrm{\Delta }`$ in Eq. (9) leads to the equation for the oscillating frequency $`\mathrm{\Omega }`$:

$$\mathrm{\Omega }=\frac{1}{\mu }\mathrm{tan}\frac{\mathrm{\Omega }\tau }{2},$$
(11)
the nonzero solution $`(\mathrm{\Omega }0)`$ of which may be searched in the range $`\pi <\mathrm{\Omega }\tau <2\pi `$. Equation (9) together with Eq. (11) gives $`\mathrm{\Delta }^2`$ in the form

$$\mathrm{\Delta }^2=\frac{1}{2K\mu }(1+\mu ^2\mathrm{\Omega }^2).$$
(12)
Note in Eq. (8) that $`\mathrm{\Delta }`$ approaches unity as $`\sigma /K`$ is reduced to zero ($`\sigma /K0`$). In particular Eq. (8) gives the range $`0\mathrm{\Delta }^21\sigma ^2/K^2,`$ where the upper bound of $`\mathrm{\Delta }^2`$ determines the boundary between stationary and oscillatory states. Accordingly, on the boundary Eq. (12) reads

$$\mathrm{\Omega }=\sqrt{\frac{2K}{\mu }\left(1\frac{\sigma ^2}{K^2}\right)1},$$
(13)
which in turn yields the boundary

$`\tau `$ $`=`$ $`2\left[{\displaystyle \frac{2K}{\mu }}\left(1{\displaystyle \frac{\sigma ^2}{K^2}}\right){\displaystyle \frac{1}{\mu ^2}}\right]^{1/2}`$ (14)
$`\times `$ $`\left[\pi \mathrm{arctan}\sqrt{2K\mu \left(1{\displaystyle \frac{\sigma ^2}{K^2}}\right)1}\right].`$ (15)
It can be observed that the oscillatory state appears only for $`\mu \mu _c\left[2K(1\sigma ^2/K^2)\right]^1`$. The critical value $`\mu _c`$ decreases, allowing larger regions of the oscillatory state, as the coupling strength $`K`$ is raised and as the variance $`\sigma ^2`$ is reduced. The order parameter $`\mathrm{\Delta }`$ obtained from Eq. (9) for $`K=3`$ and $`\tau =2`$ is plotted by the broken curve in Fig. 4 (b), manifesting good agreement with the numerical result. The tiny discrepancy of $`O(10^2)`$ between the two curves can be attributed to the truncation of higher-order terms in Eq. (6).

Solving Eq. (11) in the appropriate regime, we obtain the oscillation frequency as a function of inertia and delay. Figure 5(a) shows the behavior of the fundamental frequency $`\mathrm{\Omega }/2\pi `$ on the plane of inertia and delay. For comparison, the dependence on inertia for given amount of delay $`\tau =1.8,\mathrm{\hspace{0.17em}2.0}`$, and $`2.2`$ is shown in Fig. 5(b), together with the corresponding numerical data (previously plotted in Fig. 2). Note the good agreement between the analytical results, obtained from Eq. (11), and those from numerical simulations. It is manifested that the oscillation frequency in general decreases with the amount of delay and inertia.

The phase diagram in the three-dimensional space of $`(\mu ,\tau ,K)`$ is displayed in Fig. 6(a), where the boundary surface separates the oscillatory state from the stationary one. Figure 6(b) exhibits the phase boundaries on the plane of inertia and delay for several values of the coupling strength. Also shown are the data obtained from numerical simulations for $`K=3`$ (see Fig. 3), again demonstrating good agreement of the analytical results based on Eq. (14) with the simulation results. It is further observed that the region of the oscillatory state grows with the coupling strength. Note, however, that the analytical approach, focusing on the existence of a nonzero solution for $`\mathrm{\Omega }`$, is not able to discern the coexistence of oscillatory and stationary phases from the oscillatory phase. Thus the numerical data represented by solid squares in Fig. 3 have no counterpart in the analytical results shown in Fig. 6.

In summary, we have studied the synchronization phenomena in a system of coupled oscillators, each possessing finite inertia and interacting with others via time-delayed interactions. The interplay between inertia and time delay has been investigated in the temporal behavior of the phase velocity, which reveals the emergence of (sustained) spontaneous oscillation in the absence of external periodic driving. We have also obtained the phase diagram, which displays oscillatory and stationary states in the appropriate regions of the three-dimensional space consisting of the inertia, delay, and coupling strength.

We acknowledge the hospitality from Korea Institute for Advanced Study, where part of this work was performed, and the partial support from the BK21 Project.


format: text

text:
# Core Strings and Flux Spreading Near Pinning Centers

\[

## Abstract

At the nucleus of a superconducting vortex is a small, circular core region where superconductivity is destroyed. Like an atomic nucleus, this core may become deformed, and such distortions can have important consequences in non-equilibrium situations. Using Ginzburg-Landau theory, we have investigated this phenomenon for vortices in the presence of artificial defects. We show that when a vortex approaches the vicinity of a defect, an abrupt transition occurs in which the vortex core develops a “string” extending to the defect boundary, while simultaneously the supercurrents and associated magnetic flux spread out and engulf the defect. The energetics of stretching the string determines the pinning behavior of the vortex. Experimental consequences of these strings are discussed.

\]

One of the most important properties of superconductors is their ability to carry currents without dissipation, allowing them to generate large magnetic fields. Many superconductors allow fields to penetrate in bundles of quantized magnetic flux, with associated whirlpools of current known as vortices. When these vortices are mobile, they spoil the perfect conductivity that make superconductors so useful. The quest to increase the maximum dissipationless current $`J_c`$ that a superconductor may carry has thus fueled intense study of vortex pinning.

In recent years, pinning environments of artificially fabricated nanoscale defect arrays have been developed in hopes of better understanding and improving the pinning properties of superconductors. Some of the earlier contributions have involved macroscopic measurements (e.g. $`J_c`$, magnetization) performed on periodic “antidot” arrays. The antidot regions contain material which is rendered non-superconducting. Pinning behavior may be studied in these periodic systems using various imaging techniques.

Much of the theoretical work on vortex pinning has employed numerical studies to focus on the behavior of large collections of vortices under the influence of a driving force (supercurrent). Such studies usually employ simplified pinning potentials, in part to make possible simulations of large numbers of vortices, but also because information about pinning potentials at the microscopic level is simply unavailable. A few studies have focused on energy scales for pinning varying numbers of flux quanta to the defects, as well as defect-vortex potentials as derived from the London equation. However, the latter approach does not allow for variation of the Cooper pair density, and in particular cannot correctly treat the vortex core. In the London approach, vortex cores are usually assumed to be rigid in shape, and interactions of vortices with their environments are determined by the core position as well as the distribution of currents. Our results demonstrate that the vortex core in fact deforms dramatically near an artificial defect: when the vortex center is sufficiently close to the defect, a string of suppressed order parameter develops from the vortex position to the pinning center edge, while the currents and magnetic flux spread out over a large area (see Fig. 2 below). The pinning potential of the vortex turns out to be dominated by the string when it is present, as we discuss below.

Methods– Our calculations focus on two dimensional arrays of artificial pinning centers in the form of holes in a bulk superconductor. Our goal is to find the lowest energy state of the system for specified locations of a superconducting vortex; from this we can construct a pinning potential. The appropriate description of the superconducting state is in terms of Ginzburg-Landau theory, which focuses on a complex superconducting order parameter $`\psi (\stackrel{}{r})`$, for which $`\left|\psi \right|^2`$ is proportional to the local density of superconducting electrons. Unlike the London theory, Ginzburg-Landau theory is valid at scales as small as the coherence length, $`\xi `$. Written in terms of dimensionless variables, the Ginzburg-Landau energy functional is

$$E_{GL}=\left[\begin{array}{c}\left|\psi ^{}\left(\stackrel{}{}/i\stackrel{}{A}\right)\psi \right|^2\left|\psi \right|^2\\ +\frac{\kappa ^2}{2}\left|\psi \right|^4+B^2\end{array}\right]d^3x.$$
(1)
In Eq. 1, $`\stackrel{}{A}`$ is the vector potential, and the magnetic field $`\stackrel{}{B}(\stackrel{}{r})=\stackrel{}{}\times \stackrel{}{A}`$. $`\kappa \frac{\lambda }{\xi }`$ is the Ginzburg parameter, the ratio of the magnetic penetration depth $`\lambda `$ and the coherence length. In this work, we report on results obtained for $`\kappa =8`$. This choice of $`\kappa `$ is deep enough into the high $`\kappa `$ limit that, apart from scale factors in $`E_{GL}`$ and $`\psi `$, the results vary little as $`\kappa `$ is increased.

To analyze the behavior of the vortex near the pinning center, we employ a mean field approach in which one minimizes $`E_{GL}`$ for a fixed vortex location to find $`\psi `$, $`\stackrel{}{A}`$, and the current $`\stackrel{}{J}`$. Our strategy for calculating $`\psi `$ and $`\stackrel{}{A}`$ self-consistently involves first holding $`\psi `$ fixed at some initial guess, and minimizing $`E_{GL}`$ with respect to $`\stackrel{}{A}`$ and $`\stackrel{}{B}`$. Next, we fix $`\stackrel{}{A}`$ and $`\stackrel{}{B}`$ and minimize with respect to $`\psi `$. These steps are iterated until changes in the variables become negligible. We implement this self-consistent approach numerically by dividing the unit cell into a fine lattice of small unit cells. In this discrete scheme, $`\psi (\stackrel{}{r})`$ is replaced by $`\psi _{ij}`$ with $`ij`$ specifying a grid point on a square lattice, while $`\stackrel{}{A}_{ij}`$ and $`\stackrel{}{J}_{ij}`$ are defined on nearest-neighbor links between the grid points. Derivatives in Eq. 1 are replaced by the corresponding finite differences. The resulting theory is very similar to lattice gauge theories studied in particle physics. To model the defect array as accurately as possible, one desires a fine grid; we find that with a $`128\times 128`$ mesh our results are well-converged with respect to the discretization.

To see how one minimizes $`E_{GL}`$ under the constraint of a specified vortex location, it is useful to write the current $`\stackrel{}{J}`$ in terms of the order parameter $`\psi `$ and the vector potential $`\stackrel{}{A}`$. By minimizing $`E_{GL}`$ with respect to the vector potential and employing a Maxwell equation one has

$`\stackrel{}{J}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\psi ^{}\left(\stackrel{}{}/i\stackrel{}{A}\right)\psi +\psi \left(\stackrel{}{}/i\stackrel{}{A}\right)\psi ^{}\right]`$ (2)
$`=`$ $`\left|\psi \right|^2\left(\stackrel{}{}\varphi \stackrel{}{A}\right).`$ (3)
In Eq. 2, we have used for the order parameter $`\psi =\left|\psi \right|e^{i\varphi }`$. The familiar fluxoid quantization condition arises from the requirement that the order parameter be single valued, i.e. $`\stackrel{}{}\varphi 𝑑\stackrel{}{s}=2\pi n_v`$. Hence, in terms of $`\stackrel{}{J}`$ and $`\stackrel{}{A}`$,

$`2\pi n_v(ij)`$ $`=`$ $`{\displaystyle \left(\stackrel{}{J}/\left|\psi \right|^2\right)𝑑\stackrel{}{s}}+{\displaystyle \stackrel{}{A}𝑑\stackrel{}{s}}`$ (4)
$`=`$ $`{\displaystyle \left(\stackrel{}{J}/\left|\psi \right|^2\right)𝑑\stackrel{}{s}}+\mathrm{\Phi }_B`$ (5)
In the second part of Eq. 4, $`\mathrm{\Phi }_B`$ is the total magnetic flux passing through the area of the contour, which we conveniently choose to be the small unit cell associated with the grid point $`ij`$, while $`n_v(ij)`$ is the total number of “fluxoid quanta” contained in the contour of integration. It is through $`n_v`$ that the vortex location(s) in the full unit cell of the system may be fixed: $`n_v=0`$ except at the grid points where we wish to place a vortex, for which $`n_v=1`$. Armed with knowledge of $`|\psi |`$ and some specified realization of $`n_v(ij)`$, one can solve for $`\stackrel{}{J}`$ and $`\stackrel{}{A}`$ via Eq. 4. Using the expression for the current given in Eq. 2, one obtains $`\stackrel{}{}\varphi `$; inserting $`\stackrel{}{}\varphi `$ and $`\stackrel{}{A}`$ into Eq. 1 yields an expression depending only on $`\left|\psi \right|`$ and $`\kappa `$ which we minimize with respect to $`\left|\psi \right|`$ to obtain the order parameter modulus.

We note that the above method can be generalized to the case of a thin film superconductor. We have performed some calculations for such systems, and find that for antidot systems the results obtained are quite similar to the large $`\kappa `$ results we report here. This may be understood in terms of the effective magnetic penetration depth for a thin film, $`\lambda _{eff}=\frac{\lambda ^2}{d}`$, which is typically much larger than the bulk value $`\lambda `$. This means that the energy stored in the magnetic field generated by the supercurrents is quite small, so that the fact that the field varies as one moves out of the plane has little impact on the state of the system. The resulting energy functional is thus nearly identical to the bulk three dimensional case, with columnar antidots and vortices.

Results– At large vortex-defect separations, the core has the usual compact structure with supercurrents localized about it. Fig. 1 presents a perspective plot of $`\left|\psi \right|`$, as well as a vector plot of the currents. The distances shown are in units of the coherence length, $`\xi `$. Our choice of a unit cell with side spanning 20 coherence lengths is typical of many of the nanoscale arrays studied experimentally. As the flux quantum nears the defect edge, it eventually reaches a critical distance $`d_c`$ (typically several $`\xi `$, with precise value depending on the dot size and shape), where there is a sudden dramatic change. Fig. 2 illustrates the situation after the transition: the vortex core has developed a string extending from the flux quantum position to the defect edge; simultaneously the current now encircles the vortex-defect pair, and the magnetic flux created by these currents spreads over a larger area. The string is energetically favorable because it allows the formerly dense current of the vortex to spread out (engulfing the defect in the process), thereby reducing the kinetic energy of the state.

When a net supercurrent flows across a superconductor, a Lorentz force is exerted on vortices; if these move in response energy is dissipated. In an artificial defect array, the driving force can be balanced by a pinning force given by the gradient of $`E_{GL}`$ with respect to the vortex position. Hence, it is appropriate to regard $`E_{GL}`$ as the pinning potential. Fig. 3 shows this as a function of distance from the unit cell edge, and it can be seen to have three distinct regions. For large separations, where the vortex core is compact, the pinning potential decreases relatively slowly. As the separation decreases it eventually crosses $`d_c`$ and the core-string structure appears. Hysteresis in the calculations, indicated with arrows, suggests that the transition is first order. Below the transition one observes nearly linear behavior of the pinning potential, suggesting that the string carries an energy proportional to its length. The third region is announced by a discontinuous jump as the vortex is absorbed by the defect, followed by a perfectly flat region inside the defect. This abrupt jump is due to a sudden transformation of the order parameter. Immediately prior to absorption, constant Cooper pair density contours near the flux quantum have a semicircular profile. As the vortex enters the defect, this “semi-vortex” vanishes, removing a finite amount of energy for an infinitesimal change in the vortex position.

Experimental implications– The profile of the pinning potential shown in Fig. 3 suggests a set of measurements one might perform to seek an experimental signature of the unique pinning phenomena discussed above. One possible test would involve an AC driving current of magnitude small enough to leave the vortices pinned by the strings. Such a current would allow the vortices to “rattle around” in the linear part of the pinning well illustrated in Fig. 3 and produce losses, whereas a DC current of equal magnitude would be dissipationless. Upon lowering of the temperature the vortices would be captured by defects, leading to lossless supercurrents for both AC and DC driving forces. An observation of these effects would yield indirect confirmation of the form of the pinning potential we find. Something like this may recently have been observed. Another interesting possibility is that a unique thermal depinning may occur as the temperature is increased in the regime of linear pinning. The presence of a string suggests that this may carry entropy at finite temperature, much as is the case of polymers. This entropy is proportional to the string length and temperature, and at high enough temperatures may overwhelm the energy per unit length found in our mean-field calculations. In analogy with polymer behavior , this leads to unbounded growth of the string and effective depinning of the vortex. However, it is not clear whether the string remains sufficiently well-defined at the temperatures necessary for proliferation that the polymer analogy remains valid up to the transition. Further research into this possibility is currently underway.

Ultimately the best confirmation of our results would involve imaging of the string. An interesting possibility in this context is to fabricate arrays with two defects per unit cell in close proximity. When one magnetic flux quantum passes through each unit cell, we have found solutions in which a string develops between the defects. Fig. 4 is an image of $`|\psi |`$ for an equilibrium configuration in a typical case. The existence of a string between defects opens the possibility of detecting these objects at low temperatures and without driving currents, simplifying the relevant experimental conditions.

Summary– In the context of the Ginzburg-Landau theory, we have given a detailed treatment of the microscopic aspects of pinning phenomena in nanoscale periodic arrays. Strikingly, we see an apparent first order transition involving the creation of a string connecting the vortices to the defect, and an accompanying abrupt transformation of the supercurrent and the magnetic field it generates. The string configuration leads to a region of linear pinning. Absorption of the vortex by the antidot is marked by a jump in the pinning potential. Various aspects of the pinning potential should be observable in experiment.

Acknowledgments– The authors would like to thank L. E. DeLong, S. B. Field, and J. B. Ketterson for useful discussions. This work was supported by NSF Grant No. DMR-9870681 and DMR-0108451.


text:
# A COMPARISON OF TWO APPROACHES: MAXIMUM ENTROPY ON THE MEAN (MEM) AND BAYESIAN ESTIMATION (BAYES) FOR INVERSE PROBLEMS

## 1. Introduction

Inverse problems arises in many areas of science and engineering. In fact, rarely, we can measure directly a quantity $`x`$ and, in general, the unobserved interested $`x`$ is related to the measured quantity $`y`$ via a model. In many area this model can be written in the general form $`g=𝒜(x)+n`$ or in the discrete case:

$$𝒚=𝑨(𝒙)+𝒏,$$
(1)
where $`𝒚`$ stands for the data, $`𝒙`$ for the unknown variables and $`𝒏`$ for the errors (modeling and noise). Since Newton and Gauss, one tries to define a solution to this problem as the optimizer of a criterion, for example the Least Squares (LS):

$$\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{min}}\left\{𝒚𝑨(𝒙)^2\right\}.$$
(2)
But the inverse problems are, in general, ill-posed and the LS criterion may not have a unique optimum or this solution may be very sensitive to noise. Since Tikhonov , the regularization theory became the main approach to give a satisfactory solution by defining it as the optimizer of a compound criterion:

$$\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{min}}\left\{J(𝒙)\right\}\text{with}J(𝒙)=Q(𝒙)+\lambda \mathrm{\Omega }(𝒙)=𝒚𝑨(𝒙)^2+\lambda 𝑫𝒙^2$$
(3)
or in its more general forms :

$$\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{min}}\left\{J(𝒙)\right\}\text{with}J(𝒙)=Q(𝒚𝑨(𝒙))+\lambda \mathrm{\Omega }(𝒙,𝒎).$$
(4)
The questions then raised on how to choose the functionals $`Q`$ and $`\mathrm{\Omega }`$ and the regularization parameter $`\lambda `$ and the default solution $`𝒎`$.

The probabilistic approaches started to give partial answers to this request. In particular in the Bayesian estimation approach and the maximum a posteriori (MAP) estimate:

$$\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{max}}\left\{p(𝒙|𝒚)\right\}=\text{arg}\underset{𝒙}{\mathrm{min}}\left\{\mathrm{log}p(𝒙|𝒚)\right\}=\text{arg}\underset{𝒙}{\mathrm{min}}\left\{\mathrm{log}p(𝒚|𝒙)\mathrm{log}p(𝒙)\right\},$$
(5)
this choice is: $`Q=\mathrm{log}p(𝒚|𝒙)`$ and $`\lambda \mathrm{\Omega }(𝒙)=\mathrm{log}p(𝒙)`$. This approach just pushed a little farther the questions which became how to translate our prior knowledge into a probability law and how to determine their parameters. Even, nowadays, there are many tools for the estimation of the hyperparameters , the main question on how to translate some knowledge about $`𝒙`$ into a probability distribution stays without a complete answer. The maximum entropy (ME) principle gave partial answers . See also  for an extensive discussed bibliography.

At the same time, many authors used the ME principle to find unique solutions to linear inverse problems by considering $`𝒙`$ as a distribution and the data $`𝒚`$ as linear constraints on them. Then, assuming that the data constraints are satisfied by a non empty set of solutions, a unique solution is chosen by maximizing the entropy:

$$\underset{j}{}x_j\mathrm{log}x_j\text{or}\underset{j}{}x_j\mathrm{log}\left[\frac{x_j}{m_j}(x_jm_j)\right],$$
(6)
where $`𝒎`$ is default solution. See for example  and the cited references. However, even if in these methods, thanks to convex analysis and Lagrangian techniques, the constrained optimization of 6 can be replaced by an equivalent unconstrained optimization, the obtained solutions satisfy the uniqueness condition of well-posedness but not the stability one .

Recently, some authors  used the ME principle in a different way by considering $`𝒙`$ not as a distribution but as the mean value of a random vector $`𝑿`$ and the data as the constraints on its distribution $`\text{d}P(𝒙)`$. Then, the ME principle is used to define it uniquely and finally the solution $`\widehat{𝒙}`$ is defined as the expected value of this ME distribution.

Following these authors, some others used, commented and analyzed extensively these ideas ,  and . However, in all these works, the data $`𝒚`$ were considered as exact constraints and the errors on the data were either neglected or partially token account of. (See however new developments in .)

More recently, some authors who were more faced with real applications,  followed the same idea, but by fixing themselves as the objective to use these ideas for describing the solution as the optimizer of a combined convex criteria such as (4) and more on a constructive way to determine these functionals.

The objective of this paper is to make a comparison of the Bayesian approach which we call hereafter BAYES and the maximum entropy in the mean which we refer to as MEM. This comparison is done very pragmatically and is based on the understanding of the author who does not have pretension to know all the details of the both approaches and will be happy to discuss all the following discussions with the pro of the approaches.

## 2. Maximum entropy on the mean approach

### 2.1. Basics

The main references to the basics of this approach are . The original idea and first applications in crystallography are given in . More details and extensions are given in . The mathematical aspects of convex analysis and duality theorems are given in .

The following resumes the different steps of the approach:

* Consider a set $`𝒞`$, assume that $`𝒙𝒞`$ and define a reference measure $`\mu (𝒙)`$:

$$𝒙𝒞,𝒎=_𝒞𝒙\text{d}\mu (𝒙),$$
(7)
where $`𝒎`$ is the mean value of $`𝒙`$ under this reference measure.
* Consider $`𝒙`$ as the mean value of a random vector $`𝑿`$ for which you assume a probability distribution $`P`$:

$$𝒙=\text{E}_P\left\{𝑿\right\}=_𝒞𝒙\text{d}P(𝒙)$$
(8)
and the data $`𝒚`$ as exact equality constraints on it:

$$𝒚=𝑨𝒙=𝑨\text{E}_P\left\{𝑿\right\}=_𝒞𝑨𝒙\text{d}P(𝒙).$$
(9)
* Determine the distribution $`P`$ by:

$$\text{maximize}_𝒞\mathrm{log}\frac{\text{d}P(𝒙)}{\text{d}\mu (𝒙)}\text{d}P(𝒙)\text{s.t.}𝒚=𝑨𝒙=𝑨\text{E}_P\left\{𝑿\right\}.$$
(10)
The solution is calculated via Lagrangian:

$`(𝒙,𝝀)`$ $`=`$ $`{\displaystyle _𝒞\left[\mathrm{log}\frac{\text{d}P(𝒙)}{\text{d}\mu (𝒙)}\underset{i=1}{\overset{M}{}}\lambda _i(y_i[𝑨𝒙]_i)\right]\text{d}P(𝒙)}`$ (11)
$`=`$ $`{\displaystyle _𝒞\left[\mathrm{log}\frac{\text{d}P(𝒙)}{\text{d}\mu (𝒙)}𝝀^t(𝒚𝑨𝒙)\right]\text{d}P(𝒙)}`$
and is given by:

$$\text{d}P(𝒙,𝝀)=\mathrm{exp}\left[𝝀^t[𝑨𝒙]\mathrm{log}Z(𝝀)\right]\text{d}\mu (𝒙),$$
(12)
where

$$Z(𝝀)=_𝒞\mathrm{exp}\left[𝝀^t[𝑨𝒙]\right]\text{d}\mu (𝒙).$$
(13)
The Lagrange parameters are calculated by searching the unique solution (if exists) of the following system of non linear equations:

$$\frac{\mathrm{log}Z(𝝀)}{\lambda _i}=y_i,i=1,\mathrm{},M.$$
(14)
* The solution to the inverse problem is then defined as the expected value of this distribution:

$$\widehat{𝒙}(𝝀)=𝒙\text{d}P(𝒙,𝝀).$$
(15)

These steps are very formal. In fact, it is possible to determine $`\widehat{𝒙}(𝝀)`$ in a more direct manner. Using the following notations:

$$𝒔=𝑨^t𝝀,G^{}(𝒔)=\mathrm{log}Z(𝒔)=\mathrm{log}_𝒞\mathrm{exp}\left[𝒔^t𝒙\right]\text{d}\mu (𝒙),$$
(16)
and

$$H(𝒙)=\underset{𝒔}{\mathrm{max}}\{𝒔^t𝒙G^{}(𝒔)\},D(𝝀)=𝝀^t𝒚G^{}(𝑨^t𝝀)$$
(17)
it is shown that:

$`\widehat{𝝀}=\text{arg}\underset{𝝀}{\mathrm{max}}\left\{D(𝝀)\right\}`$ (Dual criterion) (18)
$`\widehat{𝒙}=\text{arg}\underset{𝒙𝒞}{\mathrm{min}}\left\{H(𝒙)\right\}\text{ s.t. }𝒚=𝑨𝒙`$ $`\text{(Primal criterion)},`$ (19)
$`\widehat{𝒙}(𝒔)={\displaystyle \frac{\text{d}G^{}(𝒔)}{\text{d}𝒔}}`$ $`\text{(Explicit relation)},`$ (20)
where:

* Functions $`G`$ and $`H`$ depend on the reference measure $`\mu (𝒙)`$;
* $`D(𝝀)`$ is the dual criterion which is function of the data and the function G;
* $`H(𝒙)=H(𝒙,𝒎)`$ is the primal criterion which is a distance measure between $`𝒙`$ and $`𝒎`$ which means:
$`H(𝒙,𝒎)0`$, and $`H(𝒙,𝒎)=0\text{iff}𝒙=𝒎`$;
$`H(𝒙,𝒎)`$ is differentiable and convex on $`𝒞`$;
$`H(𝒙,𝒎)=\mathrm{}`$ if $`𝒙𝒞`$.

Now, to be able to go a little more in details, let assume that the reference measure is separable:

$$\mu (𝒙)=\underset{j=1}{\overset{N}{}}\mu _j(x_j)$$
(21)
then, we have:

$$\text{d}P(𝒙,𝝀)=\underset{j=1}{\overset{N}{}}\text{d}P_j(x_j,𝝀)$$
(22)
and

$$G(𝒔)=\underset{j}{}g_j\left(s_j\right),H(𝒙,𝒎)=\underset{j}{}h_j(x_j,m_j),\widehat{x}_j=g_j^{}(s_j).$$
(23)
Replacing $`𝒔=𝑨^t𝝀`$ we obtain:

$$G(𝝀)=\underset{j}{}g_j\left([𝑨^t𝝀]_j\right),H(𝒙,𝒎)=\underset{j}{}h_j(x_j,m_j),\widehat{x}_j=g_j^{}\left([𝑨^t\widehat{𝝀}]_j\right),$$
(24)
where $`h_j`$ and $`g_j`$ depend on the reference measure $`\mu _j`$:

* $`g_j`$ is the log Laplace transform (Cramer transform) of $`\mu _j`$:

$$g(s)=\mathrm{log}\mathrm{exp}\left[sx\right]\text{d}\mu (x);$$
* $`h_j`$ is the convex conjugate of $`g_j`$: $`h(x)=\underset{𝒔}{\mathrm{max}}\{sxg(s)\}`$.

Let give some examples:

$$\begin{array}{cccc}& \mu _j(x)& g_j(s)& h_j(x,m)\\ & & & \\ \text{Gaussian:}& \mathrm{exp}\left[\frac{1}{2}(xm)^2\right]& \frac{1}{2}(sm)^2& \frac{1}{2}(xm)^2\\ \text{Poisson:}& \frac{m^x}{x!}\mathrm{exp}\left[m\right]& \mathrm{exp}\left[ms\right]& x\mathrm{log}\frac{x}{m}+mx\\ \text{Gamma: }& x^{\alpha 1}\mathrm{exp}\left[\frac{x}{m}\right]& \mathrm{log}(sm)& \mathrm{log}\frac{x}{m}+\frac{x}{m}1\end{array}$$

When $`\mu (𝒙)`$ is not separable it is very difficult to do the calculation more in details, excepted the Gaussian case $`\mu (𝒙)=𝒩(𝒎,𝑹_x)`$, where we have:

$$H(𝒙,𝒎)=\frac{1}{2}(𝒙𝒎)^t𝑹_x^1(𝒙𝒎),G(𝝀)=\frac{1}{2}𝝀^2,D(𝝀)=𝝀^t𝒚+\frac{1}{2}𝝀^2.$$
(25)
(See however a new presentation of the method in  trying to extend the method for taking account of the correlations.)

### 2.2. Extensions

How to account for the noise: Two approaches have been developed in :

* Replacing the exact equality constraints $`𝒚=𝑨𝒙`$ by the following inequalities:

$$|y_i[𝑨𝒙]_i|<ϵ,\text{or}𝒚𝑨𝒙^2<\sigma ^2$$
(26)
and using the duality relations they showed:

$$\{\begin{array}{c}\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{min}}\left\{H(𝒙)\right\}\text{ s.t. }|y_i[𝑨𝒙]_i|<ϵ,\text{ with }H(𝒙)=\underset{j}{}h_j(x_j)\hfill \\ \widehat{𝝀}=\text{arg}\underset{𝝀}{\mathrm{max}}\left\{\stackrel{~}{D}(𝝀)\right\}\text{ with }\stackrel{~}{D}(𝝀)=D(𝝀)+\alpha 𝝀^2\hfill \end{array}$$
(27)
where $`\alpha `$ depends on $`ϵ`$ or on $`\sigma ^2`$.
* Replacing $`𝒚=𝑨𝒙`$ by $`𝒚=𝑨𝒙+𝒏`$ and rewriting it as follows:

$$𝒚=[𝑨|𝑰]\left[\begin{array}{c}𝒙\\ 𝒏\end{array}\right]𝒚=\stackrel{~}{𝑨}\stackrel{~}{𝒙}$$
(28)
and assuming $`\mu (\stackrel{~}{𝒙})=\mu _x(𝒙)\mu _n(𝒏)`$ they showed:

$$\widehat{𝒙}=\text{arg}\underset{𝒙𝒞}{\mathrm{min}}\left\{𝒬(𝒚𝑨𝒙)+\alpha H(𝒙)\right\}$$
(29)
$$\text{with}H(𝒙)=\underset{j=1}{\overset{N}{}}h_j(x_j),\text{and}𝒬(𝒛)=\underset{i=1}{\overset{M}{}}q_i(z_j).$$
(30)
Here also $`h_j(x_j)`$ and $`q_i(z_i)`$ depend on the reference measures $`\mu _x(𝒙)`$ and $`\mu _n(𝒙)`$. The determination of $`\alpha `$ is not discussed.

## 3. Bayesian approach

### 3.1. Basics

The different steps of this approach are now well-known:

* From the observation model and the hypothesis (prior knowledge) on the noise derive the likelihood $`p(𝒚|𝒙;𝜷)`$;
* From the hypothesis (prior knowledge) on $`𝒙`$ derive the prior law $`p(𝒙|𝜽)`$;
* Apply the Bayes rule to obtain $`p(𝒙|𝒚;𝜷,𝜽)=p(𝒚|𝒙;𝜷)p(𝒙|𝜽)/p(𝒚;𝜷,𝜽)`$;
* Define an estimation rule via a cost function $`c(𝒙,\widehat{𝒙})`$ by:

$$\widehat{𝒙}(𝒚;𝜷,𝜽)=\underset{𝒛}{\text{arg}\mathrm{min}}\left\{c(𝒙,𝒛)p(𝒙|𝒚;𝜷,𝜽)\text{d}𝒙\right\}.$$
(31)

Different cost functions give different estimators:

* Maximum a posteriori (MAP):

$$C(𝒙,\widehat{𝒙})=1\delta (𝒙\widehat{𝒙})\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{max}}\left\{p(𝒙|𝒚;𝜽,𝜷)\right\}.$$
(32)
* Posterior mean (PM):

$$C(𝒙,\widehat{𝒙})=[𝒙\widehat{𝒙}]^t𝑸[𝒙\widehat{𝒙}]^t\widehat{𝒙}=\text{E}_{𝒙|𝒚}\left\{𝑿\right\}=𝒙p(𝒙|𝒚;𝜽,𝜷)\text{d}𝒙.$$
(33)
* Maximum of the Marginal a posteriori (MMAP):

$$C(𝒙,\widehat{𝒙})=\underset{j}{}1\delta (x_j\widehat{x}_j)\widehat{x}_j=\text{arg}\underset{x_j}{\mathrm{max}}\left\{p(x_j|𝒚;𝜽)\right\},$$
(34)
where

$$p(x_j|𝒚;𝜽)=p(𝒙|𝒚;𝜽)\text{d}x_1\mathrm{}\text{d}x_{j1}\mathrm{}\text{d}x_{j+1}\mathrm{}\text{d}x_n.$$
(35)

To illustrate this, let consider the case of linear inverse problems $`𝒚=𝑨𝒙+𝒏`$ with the following hypothesis:

* $`𝒏`$ is zero-mean, white and Gaussian: $`𝒏𝒩(\mathrm{𝟎},(1/\beta )𝑰)`$ which leads to:

$$p(𝒚|𝒙,\beta )\mathrm{exp}\left[\frac{1}{2}\beta 𝒚𝑨𝒙^2\right].$$
(36)
* $`𝒙`$ is Gaussian: $`𝒙𝒩(𝒙_0,(1/\theta )𝑷_0)`$:

$$p(𝒙|\theta )\mathrm{exp}\left[\frac{1}{2}\theta [𝒙𝒙_0]^t𝑷_0^1[𝒙𝒙_0]\right].$$
(37)

Then, using the Bayes rule it is easy to show that

$$𝒙|𝒚𝒩(\widehat{𝒙},𝑷)\text{with}\widehat{𝒙}=𝑷𝑨^t(𝒚𝑨𝒙_0),𝑷=\left(𝑨^t𝑨+\lambda 𝑷_0^1\right)^1.$$
(38)
The MAP solution is:

$$\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{max}}\left\{p(𝒙|𝒚)\right\}=\text{arg}\underset{𝒙}{\mathrm{min}}\left\{J(𝒙)\right\},\text{ with }J(𝒙)=Q(𝒙)+\lambda \varphi (𝒙),$$
(39)
where

$$Q(𝒙)=𝒚𝑨𝒙^2,\varphi (𝒙)=𝒙^t𝑷_0^1𝒙=𝑫𝒙,\lambda =\frac{\theta }{\beta }$$
(40)
Now, relaxing the second hypothesis, i.e; choosing other prior laws for $`𝒙`$ we obtain other MAP criteria. Let just note some special interesting cases:

* A Generalized Gaussian law for $`𝒙`$:

$$p(x_j)\mathrm{exp}\left[(x_jm_j)^\alpha \right].$$
(41)
The related MAP criterion becomes:

$$J(𝒙)=Q(𝒙)+\varphi (𝒙)\text{ with }\varphi (𝒙)=\underset{j}{}(x_jm_j)^\alpha .$$
(42)
* A Gamma law for $`𝒙`$:

$$x_j𝒢(\alpha ,m_j)p(x_j)(x_j/m_j)^\alpha \mathrm{exp}\left[x_j/m_j\right].$$
(43)
The related MAP criterion becomes:

$$J(𝒙)=Q(𝒙)+\varphi (𝒙)\text{ with }\varphi (𝒙)=\underset{j}{}\alpha \mathrm{log}\frac{x_j}{m_j}+\frac{x_j}{m_j}.$$
(44)
* A Beta law for $`𝒙`$:

$$x_j(\alpha ,\beta )p(x_j)x_j^\alpha (1x_j)^\beta .$$
(45)
The related MAP criterion becomes:

$$J(𝒙)=Q(𝒙)+\varphi (𝒙)\text{ with }\varphi (𝒙)=\alpha \underset{j}{}\mathrm{log}x_j+\beta \underset{j}{}\mathrm{log}(1x_j).$$
(46)
* A Poisson law for $`𝒙`$:

$$p(x_j)\frac{m_j^{x_j}}{x_j!}\mathrm{exp}\left[m_j\right].$$
(47)
The related MAP criterion becomes:

$$J(𝒙)=Q(𝒙)+\varphi (𝒙)\text{ with }\varphi (𝒙)=\underset{j}{}x_j\mathrm{log}\frac{x_j}{m_j}+(x_jm_j).$$
(48)
* Markovian models for $`𝒙`$:

$$J(𝒙)=Q(𝒙)+\varphi (𝒙)\text{ with }\varphi (𝒙)=\alpha \underset{j}{}\underset{iN_j}{}V(x_j,x_i).$$
(49)

### 3.2. Extensions

The Bayesian approach can be exactly applied when all the direct (prior) probability laws ($`p(𝒚|𝒙,𝜷)`$ and $`p(𝒙|𝜽)`$) are assigned. Even, choosing an appropriate law is done in general by hand, another difficulty is to determine their parameters $`(𝜷,𝜽)`$. This problem has been addressed by many authors and the subject is an active area in statistics. See ,  and also .

All these methods can mainly be divided in three main families:

* Generalized MAP: In this approach one tries to estimate both the hyperparameters and the unknown variables $`𝒙`$ directly from the data by defining:

$$(\widehat{𝒙},\widehat{𝜽},\widehat{𝜷})=\text{arg}\underset{(𝒙,𝜽,𝜷)}{\mathrm{max}}\left\{p(𝒙,𝜽,𝜷|𝒚)\right\}$$
(50)
where

$$p(𝒙,𝜽,𝜷|𝒚)p(𝒚|𝒙,𝜷)p(𝒙|𝜽)p(𝜽)p(𝜷)$$
(51)
and where $`p(𝜽)`$ and $`p(𝜷)`$ are appropriate prior laws. Many authors used the non informative prior law for them.
* Marginalization: In this approach one tries to estimate first the hyperparameters by marginalizing over the unknown variables $`𝒙`$:

$$p(𝜽,𝜷|𝒚)p(𝜷)p(𝜽)p(𝒚|𝒙,𝜷)p(𝒙|𝜽)\text{d}𝒙$$
(52)
and then, using them in the estimation of the unknown variables $`𝒙`$:

$$(\widehat{𝜽},\widehat{𝜷})=\text{arg}\underset{(𝜽,𝜷)}{\mathrm{max}}\left\{p(𝜽,𝜷|𝒚)\right\}$$
(53)
* Nuisance parameters: In this approach the hyperparameters are considered as the nuisance parameters, so, marginalized:

$$p(𝒙|𝒚)=p(𝒚,𝒙,𝜽,𝜷)\text{d}𝜽\text{d}𝜷$$
(54)
and the unknown variables $`𝒙`$ are estimated by:

$$\widehat{𝒙}=\text{arg}\underset{𝒙}{\mathrm{max}}\left\{p(𝒙|𝒚)\right\}$$
(55)
To see some more discussions and different possible implementations of these approaches see .

## 4. Discussed points in each approach

### 4.1. Mem:

* Choice of $`𝒞`$ and $`\mu (𝒙)`$:

$`𝒞`$ must be a convex set, such as: $`\text{IR}^N`$, $`\text{IR}_+^N`$, $`[a,b]^N`$

– Up to now, the whole analysis can be done for separable measures.

– The only reference measures $`\mu _j`$ which permits to go through all the steps are those for which we have analytical expression for the Laplace transform of their logarithms.
* Accounting for the noise:
In the first approach only the support and the energy of the noise is used. In the second approach we have more choices via the reference measures $`\mu _n(𝒏)`$, but determination of $`\alpha `$ stays adhoc. In fact, in general, when the reference measures $`\mu _n(𝒏)`$ and $`\mu _x(𝒙)`$ depend on any parameters, this approach lacks any tool to determine them.
* Effective calculation of the solution:
No problem, and more, this is probably the main interest of the approach which defines the solution, by construction, as the optimizer of a convex criterion.
* Characterization of the solution:
A sensitivity analysis has been proposed by , but, in my opinion, this is not enough to characterize a solution to an inverse problem.
It is not easy to use the notions of variance or covariance of the solution, because this approach does not define a posterior distribution for the solution.
* Possibility of the extension of the approach:
I have not yet seen any extension of this approach to non linear inverse problems, or linear inverse problems in which the operator depend on unknown parameters, such as blind deconvolution or antenna array processing;
The fact that we have to choose a convex set $`𝒞`$ on which the solution is defined excludes the inverse problems in which we know a priori that the solution is discrete-valued (binary, $`n`$ary images for example). This excludes the use of this approach in image segmentation or communication inverse problems (canal equalization, blind deconvolution, etc.).

### 4.2. Bayes:

* Choice of $`p(𝒙|𝜷)`$:
$`p(𝒙|𝜷)`$ can be chosen separable or not. Evidently, separable $`p(𝒙|𝜷)`$ (Entropic prior laws) simplifies the calculations. Accounting for correlations is easily done via Markovian models; In both cases (Entropic or Markovian prior laws), there are some tools for choosing them either by physical considerations, or by scale invariance arguments .
* Choice of the cost function or equivalently of an estimator MAP, PM, MMAP:
This choice is done more on the basis of cost calculation. MAP calculation needs, in general, global optimization, but does not need any integration.
MP or MMAP needs multidimensional integration, so in general, greater cost. However, there are approximate calculation techniques based on Monte Carlo methods and Gibbs sampling.
* Effective calculation of different solutions:
For MAP estimate, when the posterior law is unimodal, we can use any gradient descent based method, but if this is not the case, there are two categories of methods: Simulated Annealing or Deterministic relaxation (GNC). For more discussions on Bayesian calculations see  in this volume.

## 5. Comparisons and discussions

The following main items are discussed:

* In MEM, the unknowns $`𝒙`$ are considered as the mean values of a random vector $`𝑿`$ for which a prior probability measure $`\text{d}\mu (𝒙)`$ is defined.
* In BAYES, the unknowns $`𝒙`$ are considered as a sample of a random vector $`𝑿`$ for which a prior probability measure $`p(𝒙)`$ is defined.
* In MEM, a probability distribution $`p(𝒙)`$ is defined as the minimizer of the Kullback distance $`K(p,\mu )`$ subject to the data constraints, and the solution is defined to be $`\text{E}_p(𝑿)`$. What is interesting here is that this solution can equivalently be obtained as the minimizer of a convex criterion $`J(𝒙)`$ subject to the data constraints, and what is more attractive is that, thanks to the convex analysis, this solution can also be obtained as the stationary point of a dual criterion which can easily be calculated numerically.
* In BAYES, the posterior law $`p(𝒙|𝒚)`$ is calculated using the Bayes’ rule. In fact, the data $`𝒚`$ are considered as a sample of a random vector $`𝒀`$ for which we can define a conditional probability law $`p(𝒚|𝒙)`$ which, when used in conjunction with the prior $`p(𝒙)`$ in the Bayes’ rule will give us the posterior law, from which we can define an estimator. One of these estimators is the posterior mean $`\text{E}_p(𝑿)`$, but others can also be defined. This posterior law is used not only to define an estimate (a solution), but also to calculate any other probabilistic information about that solution.
* In MEM, in its original version, the data are considered as the exact linear constraints. The uncertainty on the data are not considered, and consequently, the uncertainty on the solution is not handled. However, some extensions are recently presented to take account of the errors on the data and to calculate the sensitivity of the solution to these errors.
* In BAYES, the errors are naturally considered through $`p(𝒚|𝒙)`$ and the uncertainty of the solution through the posterior probability $`p(𝒙|𝒚)`$. Naturally then, we can compare the information content of the data and the prior model using their entropies. We can also measure the relative information content of the posterior to prior model by $`K(p(𝒙|𝒚),p(𝒙))`$.
* In MEM, even in their extended versions, it is not easy to handle with the hyperparameters. In BAYES, there are the necessary tools to handle them.
* In MEM, one can not yet handle with non linear problems. This is not the case of the BAYES.

As a final conclusion, we have to mention that, even if the two approaches are different, they can, in some cases result to the same definition of the solution as the minimizer of the same criterion, and consequently, to give exactly the same numerical solutions to a given inverse problem. However, we can give different interpretations to the obtained result depending on the approach used to reach it. The main objective of this paper was to give a succinct presentation of the two approaches for the resolution of the inverse problems.

It is important to note that the two approaches give different views and interpretations which can be used advantageously for any application. Also, even the Bayesian approach is now really mature, the MEM approach is more recent. So, many of the conclusions I made today may be altered in future. In particular, new presentation of the method by Heinrich et al  in this volume will probably give new possibilities to the MEM method and will push the limits of the method to greatest generality.


text:
# References

Spontaneous CPT asymmetry of the Universe***Talk at V International Conference on Cosmoparticle Physics COSMION-2001, dedicated to the 80-th Anniversary of Andrei D. Sakharov (Moscow-St. Petersburg, May 21-30, 2001)

A.A.Andrianov1,2, P.Giacconi1, R.Soldati1

1 Dip. di Fisica & INFN, Univ. di Bologna, via Irnerio 46, 40126, Bologna

2 Dept. of Theor. Phys., V. A. Fock Inst. of Physics, 198504, St.Petersburg

The extended QED with renormalizable interactions breaking the Lorentz and CPT symmetry is considered and the phenomenological consequences of such a symmetry breaking are illuminated in view of recent discussion of large scale anisotropy of the Universe. Other physical effects in QED with CPT violation are examined, in particular, mass splitting between electrons of different helicities and decay of very high energetic electrons into lower energy electrons and positrons.

So far the Lorentz symmetry has been proven to hold with a very high accuracy. Nevertheless, one can inquire about whether the Special Relativity Theory is only approximate. The modern, quantum field-theoretical approach admits that the spontaneous Lorentz symmetry breaking is not excluded. At this expense the CPT symmetry can be also broken in a local field theory.

The occurrence of a very small deviation from the Lorentz invariance has been discussed recently \[1-3\], within the context of the Standard model of electroweak interactions. There, some ”background” or ”cosmological” fields are implied , leading to deviations in the propagation of certain particles, within the present experimental limits.

As the photon is a test particle of the Special Relativity the most crucial probes concern Lorentz- and CPT-symmetry breaking modifications of electrodynamics.

When one retains its fundamental character provided by renormalizability and the gauge invariance of its action, then it might be induced in the ordinary 3+1 dimensional Minkowski space-time by supplementing QED with the additional CPT-odd Chern-Simons (CS) coupling of photons to the vacuum  mediated by a constant CS vector $`\eta _\mu `$. This modification of Electrodynamics does not break the gauge symmetry of the action but changes the dispersion relations for different photon polarizations . As well a constant axial-vector field $`b_\mu `$ may couple to the fermions  splitting their masses and breaking the CPT and Lorentz symmetry.

All together the minimal modification of QED is given by the following lagrangian,

$$=\frac{1}{4}F_{\mu \nu }F^{\mu \nu }Bn_\mu A^\mu \frac{1}{2}\eta _\mu A_\nu \stackrel{~}{F}^{\mu \nu }+\overline{\psi }\left(i\gamma ^\mu _\mu +e\gamma ^\mu A_\mu m\gamma ^\mu \gamma _5b_\mu \right)\psi $$
(1)
in the axial gauge. The Lorentz and CPT symmetry breaking is parameterized by two vectors $`\eta _\mu `$ and $`b_\mu `$ which are not necessarily collinear. Their dynamical origin represents an interesting problem and one of the possible ways to induce LSB by a dynamical mechanism has been suggested recently . Namely, the spontaneous breaking of the Lorentz symmetry via the Coleman-Weinberg mechanism  has been proved for a class of models with the Wess-Zumino interaction between abelian gauge fields and pseudoscalar massless axion $`\theta (x)`$(AWZ models). Then the condensation of the gradient of axion field $`<_\mu \theta >\eta _\mu b_\mu `$ may be a natural implementation for LSB  just relating its origin to the existing “quintessence” fields . On the other hand, a part of the background vector $`b_\mu `$ may represent the constant torsion $`b_\mu ϵ_{\mu \nu \rho \sigma }T^{\nu \rho \sigma }`$ in the large scale Universe \[8-10\].

When having the modified QED (1) one searches for the phenomenological consequences of the symmetry breaking which can be registered in the experiments. As well one should examine the consistency of this QFT (stability  and microcausality) once the latter one pretends to describe the particle physics in a wide range of energies.

Distorted photons

In the CPT odd QED the linearly polarized photons created in the Universe may reveal the birefringence , i.e. the rotation of the polarization direction depending on the distance, when they propagate in the vacuum. The consistent QED exists only for space-like CS vectors $`\eta ^2<0`$ (in the presence of fermions one has also to take into account the radiatively induced CS term ) and it can be canonically quantized in the frames where this vector is entirely space-like, $`\eta _0=0`$ (see). In this case there are always four real solutions of the equation for energy spectrum

$$k_{0\pm }^2(\stackrel{}{k},\stackrel{}{\eta })=\stackrel{}{k}^2+\frac{1}{2}\stackrel{}{\eta }^2\pm \left|\stackrel{}{\eta }\right|\sqrt{\left(\stackrel{}{k}\stackrel{}{\eta }\right)^2+\frac{1}{4}\stackrel{}{\eta }^2}.$$
(2)
As it can be seen for very small wave vector $`\stackrel{}{k}`$ two photon modes (the “left” ones) remain massless and other two (the “right” ones) become massive. Therefore monochromatic plane wave solutions are possible only with a definite chiral polarization whereas the linearly polarized photons reveal the rotation of the polarization direction in time = with a distance. We stress that in the case of CPT odd QED the notion of handedness or chirality which is conserved only approximately corresponds to that one in the Maxwell QED, i.e. the left-(right-)handed photons do not have exactly circular polarizations . Meanwhile the group velocities for these modes turn to be always smaller than $`c`$. This supports a birefringence phenomenon, because the group velocities of the linearly polarized wave-packets made out of chirally polarized CS-photons are smaller than $`c`$.

The earlier analysis of radiowaves from distant galaxies had been performed  in the assumption of purely time-like LSB when $`\stackrel{}{\eta }=0`$ which however is not consistent at the quantum level. The space-like scenario was recently reexamined . The present situation can be conservatively characterized by an upper bound, $`|\stackrel{}{\eta }|10^{32}eV`$.

Fermions in the constant axial field

The consistent Dirac-type quantization of massive fermions in the constant axial-vector background is frame-dependent and available for time-like $`b^2>0`$ if the space part is small $`\stackrel{}{b}^2<m^2`$. In this case one has four real solutions: two of them with positive energies and other two with negative energies . In the rest frame of $`b=(b_0,0,0,0)`$ their dispersion law is

$$p_0=\pm \sqrt{\left(\left|\stackrel{}{p}\right|\pm \left|b_0\right|\right)^2+m^2};p^2=b_0^2+m^2\pm 2\left|b_0\right|\left|\stackrel{}{p}\right|.$$
(3)
One can see that a pair of solutions of a definite helicity always describes the massive particle but another pair of opposite helicity corresponds to a massive particle only for $`|\stackrel{}{p}|<(b_0^2+m^2)/2|b_0|`$. Above this bound the phase velocity exceeds the conventional speed of light $`c`$. Nevertheless if $`b^2>0`$ then the group velocity turns out to be $`<c`$.

Thus in such a theory the two characteristic scales arise: one is supposedly small $`|b_\mu |m`$ and another one is extremely high $`Mm^2/2|b_0|m`$. Respectively two phenomena occur at low and high energies. The first one is the mass splitting between electrons of different helicities ($`\pm `$ in (3)). if the value of $`b_0`$ is universal for all fermions, then the lightest massive ones - electrons and positrons - should give the best precision in its determination. As the precision of a measurement of the electron mass  is of the order $`10^8`$, the upper bound on the time component of CPT breaking vector is not very stringent: namely, $`|b_0|<10^2\mathrm{eV}`$ (or even weaker, being controlled by the energy resolution in accelerator beams). On the other hand, some much more stringent bounds were obtained for space components of $`b_\mu `$, from the experiments with atomic systems using hydrogen masers  - i.e. $`|\stackrel{}{b}|<10^{18}\mathrm{eV}`$ \- and with a spin-polarized torsion pendulum  - i.e. $`|\stackrel{}{b}|<10^{20}\mathrm{eV}`$.

The second phenomenon takes place at very high energies, $`M>m^2/|b_0|10^2TeV`$ when some of electron states obtain the space-like energy-momentum. Then conventional wisdom of relativity theory fails and a high energetic electron may well decay into an electron of the same helicity and a pair of positron and electron of the positive and negative helicities respectively. This process starts for $`|\stackrel{}{p}|>2m^2/|b_0|`$ and the final momenta of particles are roughly three times less than the initial momentum. Such a process will “wash out” from the asymptotic states very high energetic electrons and positrons and therefore the natural momentum cutoff arises being stuck to the rest frame of the vector $`b_\mu `$. The ultraviolet divergences are then cured by this physical cutoff and the effective CS vector radiatively induced by one fermion loop is defined uniquely , independently on renormalization scheme (compare to ),

$$\mathrm{\Delta }\eta _\mu =\frac{2\alpha }{\pi }\underset{a=1}{\overset{N}{}}b_\mu ^a,$$
(4)
where the summation has to be performed over all the $`N`$ internal charged fermions degrees of freedom and the possibility to have different axial charges for different fermions is taken into account.From the recent experimental bounds obtained in Refs.  it is possible to estimate the magnitude of the spatial components of the induced CS vector to be of the order $`|\mathrm{\Delta }\stackrel{}{\eta }|<10^{19}\mathrm{eV}`$, under the assumption that all the vectors $`\stackrel{}{b}^a`$ are of the same order. If the vectors $`b_\mu ^a`$ are related to some background torsion, then one expects them to be identical. If, however, they are generated by e.g. vacuum expectation values of gradients of axion fields, then it is conceivable that $`b_\mu ^a`$ might have different values.

Consistency between photons and fermions

The consistent quantization of fermions in a constant axial-vector field asks for the vectors $`b_\mu ^a`$ to be time-like  which, nevertheless, does not mean that the sum in eq. (4) is also time-like, because some of the fermions may have the opposite axial hypercharges.

On the other hand, the consistent quantization of photons can be achieved when the full dressed CS vector $`\stackrel{~}{\eta }_\mu =\eta _\mu +\mathrm{\Delta }\eta _\mu `$ turns out to be essentially purely space-like .

As we suppose that there is some dynamical mechanism to generate - with the help of axion condensation - the purely bosonic part $`\eta _\mu `$ of the full CS four-vector, the compatibility of the consistent quantization of both fermions and bosons is believed to be quite possible, contrary to the claim in Refs. .

However, from the practical point of view, the cancellation between time components of $`\eta _\mu `$ and $`\mathrm{\Delta }\eta _\mu `$ should be extremely precise, in order to satisfy the experimental bounds  and to fulfill the microcausality requirement of the photo-dynamics. As well the severe experimental bounds on $`\stackrel{}{b}^a`$ for electrons, and protons , together with the estimation of birefringence of radio-waves from remote galaxies and quasars, do not leave too much room to eventually discover the CPT and Lorentz symmetry breaking in the quantum spinor and photon dynamics.

A. A. A. is partially supported by Grant RFBR 01-02-17152, Russian Ministry of Education Grant E00-33-208 and by The Program Universities of Russia: Fundamental Investigations (Grant 992612); P. G. is supported by Grant MURST-Cofin99.