text:
# CHARACTERIZING WIDTH UNIFORMITY BY WAVE PROPAGATION

## I Introduction

Structures defined by the agglomeration of basic elements are often characterized by the presence of objects and voids defined between the latter. Figure 1 illustrates such a situation with respect to neuronal cells grown in vitro on different substrata, namely fibronection (a) and tissue culture plastic (b). Several important properties of such structures can be inferred from the geometry of the objects and voids present in such images, which has motivated the growing use of imaging methods in physics and in biology (e.g. Mandelbrot:1983 ; Hovi:1996 ; Plotnick:1996 ; Einstein:1998 ). For example, the migration of nerve and glial cells and axon elongation in the developing nervous system lead to different patterns of cell agglomeration and separation Dickson:2002 . In material sciences, the geometrical properties of the voids have strong influence over the respective mechanical and electrical properties of the composite and can tell much about the composite formation process. The present work describes a new concept and methodology for characterizing the width uniformity of the objects or voids in such agglomerates.

The first important point is to properly and accurately characterize what is meant by uniformity of objects or voids. A possible interpretation would regard how these elements distribute themselves along the composite space. We could be interested, for example, in quantifying their translational invariance, a problem that could be suitably tackled by using the lacunarity measure described in Mandelbrot:1983 ; Einstein:1998 ; Hovi:1996 . The alternative problem considered in the present work is to quantify the uniformity of the width of objects or voids. Figure 2 illustrates such a situation with respect to two hypothetical structures chosen for didactical purposes. The agglomerate in (a) presents a single void (the black connected region) and objects (in white) that have similar widths, while the case in (b) provides a non uniform collection of objects. The importance of quantifying the kind of width uniformity seen in (a) is particularly relevant because it may provide valuable information about the processes intrinsically related to the composite formation as well as about the interactions between the basic elements of the composite. Ideally, it would be interesting to have a single uniformity index that could be assigned to each agglomerate image. The current work proposes an effective solution to such a problem that can be extended also for analysis of connected components (i.e. objects or voids) involving connected holes such as the voids in the two situations in Figure 2.

Firmly based on concepts from mathematics and physics, wavefront propagation approaches have been extensively used for several purposes in image analysis and applications (e.g. Osher:1988 ; Sethian:1999 ). Typically, such methods involve the numeric propagation, along the object contour normal, of a wavefront, starting from the object of interest, in such a way that shocks between portions or individual fronts define the medial axes of the objects under analysis. A shock is henceforth understood as the collision between two different points along the propagating wavefront(s). Medial axes are often defined as the set of points corresponding to the centers of circles maximally inscribed into the object of interest, which coincide with the shock points. Also known as skeletons, such axes can be informally said to correspond to the ”middle” part of the objects (see Figure 5). Although alternative approaches to medial axes estimation (e.g. Ogniewicz:1995 ) can be used, the current work considers the numeric approach reported in ( Santorini:1999 ; Costa:2000 ; Costa:2001 ), which also allows a spatial scale parameter controlling the degree of detail of the obtained axes. The basic principle underlying the method for quantifying width uniformity described in the current paper is the fact that the distance values along the medial axes provide a natural characterization of the width distribution of objects even when they are not straight or present ramifications and holes. The potential of the use of the distances along the medial axes is illustrated in Figure 3. Any object in this figure is characterized by the peculiar situation that the same distance value is found along the respective medial axes. The concept of width can therefore be extended from straight objects to a more general set of shapes, including those that are curved or present branches. In brief, all shapes in Figure 3(a) can be said to have the same width. This class of objects has been studied in Attali:1994 \] under the name of polyballs. While such objects with constant width represent a somewhat limited class of possible objects, the concept of width can be immediately extended to more general shapes by considering the distance value along the object medial axes. Figure 3(b) illustrates such an approach with respect to a more generic object. The distances d along the respective medial axis represent by the dashed line, provide a natural means to characterize the widths of distribution the object.

Such distributions of distances along the medial axes are henceforth taken as an accurate indication about the object width. While the histogram of the distances along the medial axes of all objects (or voids) in the agglomerate could be taken as a more comprehensive statistical model for characterizing the object width, we restrict our attention to its mean and standard deviation, represented henceforth as $`\overline{\text{w}}`$ and $`\sigma _w`$, as well as the correlation coefficient v between these two measurements. Given its adimensional nature, the latter provides a particularly interesting candidate for quantification of the uniformity of the object’s width. As will be illustrated in this article, these three measures provide meaningful and comprehensive geometrical characterization of the objects or voids in agglomerates not only in terms of width uniformity, but also as far as the size of those objects is concerned. This article starts by presenting the methodology suggested to estimate the medial axes and distances, leading to the three width features $`\overline{\text{w}}`$ , $`\sigma _w`$, and v. Next, the potential of these measures to quantify width uniformity is illustrated with respect to synthetic and real data corresponding to agglomerates of neuronal cells grown over four different types of substrata.

## II Methodology

Given an image of an agglomerate characterized by objects over a background (or any two distinct phases) like that in Figure 2(b), the first step of the proposed methodology involves identifying all connected components for the object or voids. A connected component is understood as the set of all pixels that are connected to its 4 or 8 neighbors Costa:2001 . This work adopts 4-neighborhood. In case we are analyzing voids, a single connected component is obtained for the agglomerates in Figure 2. At the same time, a total of 10 connected components are obtained in case we are interested in the objects. The extraction of the connected components can be done by using standard region growing or inundation algorithms (e.g. Costa:2001 ). Let P be any of the points of the component to be extracted. Such methods search for the neighbors N of P that have the same pixel value as P. The procedure is repeated for the neighbors N until no more neighbors with the same pixel value are found. Once the connected components have been obtained, we need to estimate the multiscale medial axes of each of them. This procedure is illustrated with respect to one of the void connected components from the agglomerate in Figure 2(b), shown in Figure 4(a). First, the border of the connected component is extracted in clockwise or counterclockwise fashion. Traditional algorithms for border extraction (e.g. Costa:2001 ) are used in the present work. The output of such algorithms is a spatially quantified parametric curve defined by the list of ordered coordinates x(i),y(i) , where i=1,2,…,N is the curve parameter corresponding to the order the border elements are visited by the border extraction procedure. Figure 4(c) illustrates a portion of the extracted border of the connected component in (a), parametrized by integer values corresponding to the sequence in which the border was followed. The arc length s(i) along such a curve is then accurately estimated by using the spectral method described in Costa:2001 . This method involves obtaining the Fourier series for the parametric curve defined by the extracted border and calculating the derivatives by using the property $`\dot{x}\mathrm{}^1\{\mathrm{}(s)G_\sigma \}`$, where $`\mathrm{}`$ is the Fourier transform and $`G_\sigma `$ is the Fourier transform of the unit area Gaussian smoothing function $`G_\sigma =1/(\sigma \sqrt{2\pi })\mathrm{exp}\{0.5(\frac{x}{\sigma })^2\}`$ with standard deviation $`\sigma `$. The arc length can be immediately obtained from the calculated derivatives as $`s(i)={\displaystyle \underset{k=1}{\overset{i}{}}}\sqrt{\dot{x}^2(i)+\dot{y}^2(i)}`$ . Additional information about derivative estimation through spectral methods can be found in Costa:2001 . The value of $`\sigma `$ is henceforth fixed at 7 pixels, which minimizes the effect of the spatially quantized nature of the curve - recall that in digital images such curves are represented over the orthogonal lattice - without excessively affecting (smoothing) the curve. The estimated arc length for the object in Figure 4(a) is given in (d). Having obtained the border of the connected component reparametrized in terms of its arc length, a wavefront is propagated from the border, unfolding with constant speed and normal orientation. The shocks implied by such a traveling wavefront as a consequence of curvature variations along the object geometry define the medial axis of the connected component under analysis. The methodology proposed in Santorini:1999 ; Costa:2000 ; Costa:2001 , and reviewed below, is used in this work for such a purpose. The basic idea of this method is to propagate circular waves centered at each of the contour elements, according to the Huyghens principle of wave propagation. Constant radial speed is assumed. Given that the objects in digital images are represented over the orthogonal lattice, which implies specific anisotropies, it is necessary to consider only those distances that are representable in such spaces. For instance, the first such distances are $`0,1,\sqrt{2},2,\mathrm{}`$ . Once such distances, which have been called exact distances, are pre-calculated (see Costa:2001 ) and stored jointly with the respective relative position vectors with respect to the origin, the algorithm consists in visiting each of the border elements for each subsequent exact distance, and updating the empty cells indexed by the respective relative positions with the current distance value or the respectively associated arc length. In the former case, the resulting image corresponds to the distance transform of the object, which associates to each surrounding pixel the smallest Euclidean distance to the object. In case the arc length is updated, we get a labeled image that can be further processed in order to obtain label differences (see Costa:2001 ; Costa:2000 ) corresponding to the intensity of the shocks between different portions of the traveling wavefront. Figures 4(e) and (f) illustrate the distance transform and propagated labels obtained for the object in (a).

By thresholding the difference image with different values T, it is possible to obtain a family of medial axes with varying degrees of detail. The higher the value of T, the more the details are filtered out. The application of the above methodology for medial axes estimation to the problem of width uniformity involves selecting a suitable value of T so as to remove unimportant small scale detail such as vertices along the object. Figure 5 illustrates several medial axes obtained for varying values of T. The specific choice of the threshold value should take into account the resolution and specific demands imposed by each application.

If connected component presents holes, as is the case with the composites in Figure 1, each hole has to be identified by using border extraction algorithms and labeled individually with specific labels (i.e. a single distinct label is assigned to each border). By applying the wave propagation method to such borders, the generalized Voronoi diagram of the borders is obtained, and the multiscale medial axes are then calculated inside each of the Voronoi cells by using the same method described above Falcao:2002 . The generalized Voronoi diagram corresponds to a partition of the image space in such a way that all points in each obtained region are closest to the respective object (one of the image connected component) than to any other object. Recall that the distance between a point and an object is defined as the smallest distance between the point and any of the points composing the object. Figure 6 illustrates the generalized Voronoi diagram and the distance values along the borders (fixed-scale medial axes) determined by the Voronoi cells.

At this point, it is possible to derive a distance weighted medial axis by assigning to each point of the selected medial axes the respective distance transform values. The set of distance values found along the medial axis are henceforth called *medial axis distances* or *skeleton distances*. The methodology for width uniformity characterization involves the consideration of the distance values along all the medial axes in the agglomerate. The mean and standard deviation of the distances over every medial axis point of the object are then estimated. Table I gives these values, as well as the correlation coefficient, for the two images in Figure 1 considering objects and voids.

Several interesting facts result, or can be inferred from such measures. To start with, the mean width $`\overline{\text{w}}`$ obtained for the objects (cells) in Figure 2(a) is less than half of the value obtained for the objects in Figure 2(b), which is in full agreement with the visual interpretation of those images. At the same time, the mean widths obtained for the voids are similar in both images. The standard deviations obtained for the objects in those images are also in full agreement with the fact that the width dispersion of the objects in (a) is much smaller than that characterizing (b). Similar values are obtained for the voids. The coefficient of variation v provides an adimensional characterization of the observed width dispersion differences for the composites in Figure 2.

## III Application to Neuronal Agglomerates

During the development of the nervous system the distribution and connections of cells is determined by a variety of factors, including the adhesion of cells and axons to extracellular matrix material such as laminin and fibronectin Dickson:2002  little is know, however, about the determinants of cell shape. Some aspects of cell-matrix adhesions and their relationship to cell shape can be studied in vitro by plating nerve cells derived from a single common precursor cell (clonal cells) on different surfaces and asking how the different surfaces regulate the morphology of the cell. Figure 7 shows that there are great differences in the shape of a clonal nerve cell line called B103 grown on polylysine, laminin, fibronection and tissue culture plastic. 400 X 319 pixel images of the cultured cells were obtained by using a Leitz DMIRB inverted microscope and an Openlab (Leica) image acquisition system connected to a personal computer. Cell shape can be assigned a numerical value and questions asked about how the shape of cells grown on the different substrata differ from each other. The original gray-level images, such as those illustrated in Figure 1, were segmented Costa:2001  in order to separate the cells from the substrate, yielding binary images where the voids are marked as ”0” and the cells (or objects) as ”1”. Attention is concentrated on the analysis of void uniformity. The segmentation was performed by automated segmentation (thresholding) followed by manual editing of the obtained structures in order to remove artifacts such as the presence of spurious particles and shadows. Figure 7 presents examples of binary images obtained by the segmentation methodology described above. Each of the connected voids were isolated by using a standard region growing approach Costa:2001 , and the multiscale medial axes were obtained by using the level set-based approach as described in Section 2. Specific medial axes were obtained by selecting the spatial scale to 10.

Figure 8 shows the phase space (or scatterplot) obtained by considering the mean and standard deviation of the void medial axes distances. The first important feature in this plot is the relatively strong overall correlation coefficient, namely 0.86, between the two considered measures. While an almost linear separation was obtained between classes 3 (laminin, triangles) and 4 (polylysine, circles), the other two classes led to more dispersed, overlapped clusters. Regarding the mean value of $`\overline{\text{w}}`$, calculated for each class, we have $`mean(\overline{\text{w}}_3)<mean(\overline{\text{w}}_2)<mean(\overline{\text{w}}_1)`$ , indicating that the cells grown over these types of substrate tended to be progressively less packed. This is probably a reflection of a competition between cell-cell and cell-substratum adhesion on each surface cells that adhere better to themselves than to the surface on which they are grown tend to aggregate or, in the case of neurons, form fascicles Schubert:1991 . Therefore, cells grown on more adhesive surfaces such as tissue culture plastic would tend to be spread out on the surface, while those on less adhesive substrate would form clumps as observed in Figure 7. The correlation coefficient obtained for each of the four clusters are respectively 0.91, 0.73, 0.84 and 0.87, indicating substantially distinct uniformities characterizing each substrate. Whereas the voids in class 1 are characterized by the highest width uniformity, those in class 2 exhibit the most heterogeneous widths.

## IV Concluding Remarks

This article has reported a novel and comprehensive approach to characterizing width uniformity in objects and voids typically found in images of agglomerates. The methodology is based on the objective characterization of ”object width” in terms of the distances found along its medial axes, which are obtained by using a wave front propagation scheme. In this way, it becomes possible to speak objectively about the width of a wide class of shapes, including those that are curved or present branches. Such an approach to width characterization is further enhanced by the fact that most objects can be approximated or decomposed into varying radius polyballs. An accurate numerical approach has been adopted for the estimation of the medial axis distances that involves object isolation, border extraction, arc length reparametrization through a spectral approach, and the calculation of multiscale medial axes by using wavefront propagation. Several of these methods were only recently described in the literature and are applied in combined fashion for width estimation for the first time in the present work. The potential of the proposed framework was illustrated and corroborated with respect to synthetic agglomerate images and real data concerning neuronal cells grown on different substrata. By examining a clonal nerve cell line grown on different protein substrata it was possible, using this new methodology, to quantify the morphological differences caused by the different substrata. It is shown that when nerve cells are grown on different surfaces, there is a dramatic correlation between the shape of the obtained voids and the tendency of the cells to associate between themselves. The use of the above methodology combined with complementary agglomerate measures such as the fractal dimension and lacunarity provide potential for a rich and comprehensive characterization of the geometry of the shapes and agglomerates of both biological and physical objects.


format: text

text:
# Periodic alternating 0,𝜋-junction structures as realization of 𝜑-Josephson junctions

## Abstract

We consider the properties of a periodic structure consisting of small alternating $`0`$\- and $`\pi `$\- Josephson junctions. We show that depending on the relation between the lengths of the individual junctions, this system can be either in the homogeneous or in the phase-modulated state. The modulated phase appears via a second order phase transition when the mismatch between the lengths of the individual junctions exceeds the critical value. The screening length diverges at the transition point. In the modulated state, the equilibrium phase difference in the structure can take any value from $`\pi `$ to $`\pi `$ ($`\phi `$-junction). The current-phase relation in this structure has very unusual shape with two maxima. As a consequence, the field dependence of the critical current in a small structure is very different from the standard Fraunhofer dependence. The Josephson vortex in a long structure carries partial magnetic flux, which is determined by the equilibrium phase.

The minimum energy of Josephson junction (JJ) usually corresponds to the zero phase difference of the superconducting order parameter JJTextbook . However, long time ago it has been predicted that the JJs with ferromagnetic interlayer (S-F-S junctions) may have ground state with phase difference equal to $`\pi `$ (”$`\pi `$-junctions”)SFSpi-junction . Only recently this prediction has been experimentally verified ryazanov ; Kontos2 ; Blum  and now the controllable fabrication of S-F-S $`\pi `$-junctions becomes possible. Another realization of $`\pi `$-junctions became possible due to the d-wave symmetry of the order parameter in high-temperature superconductors Kirtley ; VanHarlingen . Recent experiments on YBa2Cu307-Nb ramp long JJs of zigzag geometry demonstrated that such junctions are composed of alternating facets of $`0`$\- and $`\pi `$\- junctionsSmilde ; Mints . In a pioneering work Kuzii  Bulaevskii *et al.* demonstrated that the spontaneous Josephson vortex (JV) carrying flux $`\mathrm{\Phi }_0/2`$ appears at the boundary between $`0`$\- and $`\pi `$\- JJs. The structure of such semifluxon has been studied in more detail in Refs. Xu, ; Goldobin, .

In the present work we study the properties of the periodic array of $`0`$\- and $`\pi `$\- JJs. We show that depending on the ratio of lengths of $`0`$\- and $`\pi `$\- JJs, such an array can be either in the homogeneous or modulated state. The second order phase transition between these states takes place when the length mismatch between $`0`$\- and $`\pi `$\- JJs is small. At the transition point the screening length of the magnetic field diverges. In the modulated states the average phase difference $`\phi _0`$ can take any value between $`\pi `$ and $`\pi `$. Further, we call such systems $`\phi `$-junctions. Such structures were first predicted by Mints MintsPRB98  in the case of randomly alternating $`0`$\- and $`\pi `$\- JJs. We study the properties of $`\phi `$-junctions: effective Josephson current, Josephson length, and Fraunhofer-like oscillations of the critical current. Magnetic properties of long $`\phi `$-junctions are determined by two types of Josephson vortices carrying partial fluxes $`\mathrm{\Phi }_0\phi _0/\pi `$ and $`\mathrm{\Phi }_0(\pi \phi _0)/\pi `$ (see also MintsPRB98 ; MintsPRB01 ). We find analytically the shapes of these unusual Josephson vortices. We also demonstrate that at the boundary between $`\phi `$\- and usual JJs a Josephson vortex appears carrying a partial flux $`\mathrm{\Phi }_0\phi _0/2\pi `$.

Consider a periodic structure composed of alternating $`0`$\- and $`\pi `$-JJs of lengths $`d_0`$ and $`d_\pi `$ respectively. The energy per period of such a structure is given by

$`F`$ $`={\displaystyle \frac{\mathrm{\Phi }_0j_c}{2\pi c}}{\displaystyle _{d_0}^0}\left({\displaystyle \frac{\lambda _J^2}{2}}\left({\displaystyle \frac{d\varphi }{dx}}\right)^2+1\mathrm{cos}\varphi \right)𝑑x`$
$`+{\displaystyle \frac{\mathrm{\Phi }_0j_c}{2\pi c}}{\displaystyle _0^{d_\pi }}\left({\displaystyle \frac{\lambda _J^2}{2}}\left({\displaystyle \frac{d\varphi }{dx}}\right)^2+1+\mathrm{cos}\varphi \right)𝑑x,`$ (1)
where $`j_c`$ is the Josephson current density, $`\lambda _J`$ is the Josephson length of the individual junctions, $`\lambda _J^2=c\mathrm{\Phi }_0/(8\pi ^2tj_c)`$, and $`t=t_0+\lambda _1+\lambda _2`$ is the effective junction thickness JJTextbook . For simplicity, we focus on the case when the critical current densities $`j_c`$ are the same in $`0`$\- and $`\pi `$-JJs. Generalization for different critical currents is straightforward. Possible realization of such system is the zigzag junction between high-Tc and conventional superconductors as well as S-F-S junctions with periodically modulated thickness of the ferromagnetic interlayer. In the case of the zigzag structure, $`x`$ is the coordinate along the zigzag boundary. The ground state phase distribution is determined by the equation

$$\frac{d^2\varphi }{dx^2}+j(x)\mathrm{sin}\varphi =0,d_0<x<d_\pi ,$$
(2)
where $`j(x)=\lambda _J^2`$ in $`0`$-JJs and $`j(x)=\lambda _J^2`$ in $`\pi `$-JJs. At the boundaries, the phase and its derivative have to be continuous. Periodicity of $`\varphi (x)`$ implies that the solution is symmetric with respect to $`x=d_0/2`$ and $`d_\pi /2`$, i.e., $`d\varphi /dx|_{x=d_0/2}=d\varphi /dx|_{x=d_\pi /2}=0`$. Eq. (2) always has homogeneous solutions $`\varphi =0`$ and $`\varphi =\pi `$. However these solutions give the ground state only in some range of the ratio $`d_0/d_\pi `$. In general, Eq. (2) also allows for the inhomogeneous solution. Consider the case $`d_\pi <d_0`$ and weak modulation around the $`\varphi =0`$ state. In this case Eq. (2) can be linearized and its solution is given by

$$\varphi =\{\genfrac{}{}{0.0pt}{}{A\mathrm{cosh}\frac{x+d_0/2}{\lambda _J},d_0<x<0,}{B\mathrm{cos}\frac{xd_\pi /2}{\lambda _J},0<x<d_\pi .}.$$
(3)
Matching the logarithmic derivative $`d\mathrm{ln}\varphi /dx`$ of these solutions at $`x=0`$, we obtain the following condition

$$\mathrm{tanh}(d_0/2\lambda _J)=\mathrm{tan}(d_\pi /2\lambda _J),$$
(4)
for the onset of the modulated solution. The energy analysis shows that the uniform $`\varphi =0`$ solution is favorable at $`\mathrm{tanh}\left(d_0/2\lambda _J\right)>\mathrm{tan}\left(d_\pi /2\lambda _J\right)`$. In the opposite case, $`d_\pi >d_0`$, the uniform $`\varphi =\pi `$ solution is favorable at $`\mathrm{tanh}\left(d_\pi /2\lambda _J\right)>\mathrm{tan}\left(d_0/2\lambda _J\right)`$. In the case of a finite junction split into $`0`$\- and $`\pi `$\- pieces, condition (4) was first derived in Ref. Kuzii, . These results are summarized in the phase diagram shown in Fig. 1. In the limit $`d_0\lambda _J`$ the condition for the modulated solution is given by $`d_\pi >(\pi /2)\lambda _J`$. In the case $`d_0,d_\pi \lambda _J`$ the region of the modulated solution is given by

$$\frac{\left|d_0d_\pi \right|}{d_0}<\frac{d_0^2}{6\lambda _J^2}1,$$
(5)
i.e., the modulated solution exists only for a very small length mismatch. Since $`\lambda _J`$ is temperature dependent, the transition into the modulated state may occur with decreasing temperature.

In the case $`d_0=d_\pi d`$, the phase distribution is always modulated. Lets find this distribution. From symmetry $`\varphi (x)=\pi \varphi (x)`$, in particular, $`\varphi (0)=\pi /2`$. The first integral of Eq. (2 ) for $`x<0`$ is given by

$$\frac{\lambda _J^2}{2}\left(\frac{d\varphi }{dx}\right)^2+\mathrm{cos}\varphi =\mathrm{cos}\varphi _0$$
(6)
with $`\varphi _0=\varphi (d/2)`$. The solution is given by

$$x(\varphi )=\frac{\lambda _J}{\sqrt{2}}_\varphi ^{\pi /2}\frac{d\varphi ^{}}{\sqrt{\mathrm{cos}\varphi _0\mathrm{cos}\varphi ^{}}},$$
(7)
where $`\varphi _0`$ is determined by the condition $`x(\varphi _0)=d/2`$. Dependence $`\varphi _0(d/\lambda _J)`$ and shapes of the phase variation $`\varphi (x)`$ for different values of $`d/\lambda _J`$ are presented in Fig. 2. For small length $`d\lambda _J`$, it is easy to find that the phase distribution for $`d<x<0`$ is $`\varphi (x)=\pi /2+x(x+d)/\left(8\lambda _J^2\right)`$, i.e., it varies weakly around $`\pi /2`$.

Consider a structure composed of short alternating $`0`$ and $`\pi `$\- JJs with a small length mismatch, $`|d_0d_\pi |d_0`$. In this case the function $`j(x)`$ in Eq. (2) oscillates rapidly and we can present $`\varphi (x)`$ as $`\varphi (x)=\phi (x)+\xi (x)`$, where the function $`\phi (x)`$ varies slowly over the distance $`d_0`$ and $`\xi (x)`$ oscillates rapidly with $`\overline{\xi (x)}=0`$. In the limit $`d_{0,\pi }<\lambda _J`$, the oscillating term is small (as one can see from the right panel of Fig. 2 even at $`d_0,d_\pi =\lambda _J`$, (i)$`|\xi (x)|0.131`$ and (ii) the oscillating part is smaller than the average phase). In this case one can use a coarse-grained approximation (see, e.g., Ref. LanLifMech, ), meaning that we can average out the rapidly oscillating function $`\xi (x)`$ and derive the equation for the slowly varying function $`\phi (x)`$ only:

$$\frac{d^2\phi }{dx^2}=\frac{d_0d_\pi }{24\lambda _J^4}\mathrm{sin}\left(2\phi \right)+\frac{d_0d_\pi }{(d_0+d_\pi )\lambda _J^2}\mathrm{sin}\left(\phi \right),$$
(8)
Here the first term in the right hand side appears due to current modulations and the second term is proportional to the average Josephson current. A similar equation was derived in Ref. MintsPRB98 . The neglected terms in Eq. (8) are smaller by a factor of $`(d_0/\lambda _J)^2`$. Without the external magnetic field, when condition (5) is fulfilled, the equilibrium phase difference $`\phi _0`$ is:

$$\mathrm{cos}\phi _0=f\frac{12\lambda _J^2\left(d_0d_\pi \right)}{d_0d_\pi (d_0+d_\pi )},$$
(9)
and for $`|f|<1`$, the average phase difference lies in the region $`0<\left|\phi _0\right|<\pi `$,i. e. *we have the realization of a JJ with arbitrary ground state phase difference (*$`\phi `$*-junction)*. For $`|f|>1`$ homogeneous $`0`$\- or $`\pi `$\- phase is realized. An example of the dependence of the equilibrium phase $`\phi _0`$ on length mismatch is shown in the inset of Fig. 1.

The energy in terms of the coarse-grained phase is given by

$$F=\frac{\mathrm{\Phi }_0j_c\lambda _J^2}{2\pi c}\left(\frac{1}{2}\left(\frac{d\phi }{dx}\right)^2+\frac{\left(\mathrm{cos}\phi f\right)^2}{2\lambda _\phi ^2}\right)𝑑x,$$
(10)
with $`\lambda _\phi ^2=12\lambda _J^4/d_0^2\lambda _J^2`$. The coarse-grained Josephson current density flowing through the array is given by

$$j(\phi )=\frac{d_0^2j_c}{12\lambda _J^2}\mathrm{sin}\left(\phi \right)\left(f\mathrm{cos}\phi \right).$$
(11)
This current-phase relation is quite peculiar (see left panel in Fig. 3). The current has two maxima and two minima, which are achieved at $`\mathrm{cos}\phi _\pm =\left(f\pm \sqrt{f^2+8}\right)/4`$. The absolute maximum

$$j_{\mathrm{max}}=\frac{d_0^2j_c}{48\lambda _J^2}\sqrt{\frac{1}{2}\frac{f}{8}\left(f\sqrt{f^2+8}\right)}\left(3f+\sqrt{f^2+8}\right)$$

is achieved at $`\phi =\phi _{}`$.

Consider screening of the weak magnetic field by the periodic structure. In the modulated phase linearizing Eq. (8) near $`\phi _0`$, $`\phi =\phi _0+\stackrel{~}{\phi }`$, we obtain a linear equation for the phase variation $`\stackrel{~}{\phi }`$

$$\frac{d^2\stackrel{~}{\phi }}{dx^2}=\frac{\mathrm{sin}^2\phi _0}{\lambda _\phi ^2}\stackrel{~}{\phi },$$
(12)
One can see that the effective screening length

$$\lambda _{\mathrm{eff}}\lambda _\phi /\mathrm{sin}\phi _0=\lambda _\phi /\sqrt{1f^2}\text{, for }|f|<1$$
(13)
diverges at the transition point. In the homogeneous phase, $`|f|>1`$, similar derivation gives

$$\lambda _{\mathrm{eff}}=\lambda _\phi /\sqrt{|f|1}\text{, for }|f|>1.$$
(14)
Because $`\lambda _\phi (T)`$ diverges at $`TT_c`$ faster than $`f(T)`$, $`\lambda _{\mathrm{eff}}(T)`$ is nonmonotonic above the transition point. The minimum value of $`\lambda _{\mathrm{eff}}`$, $`\lambda _{\mathrm{eff}}=2d_0^2/\sqrt{3}\left(d_0d_\pi \right)`$ is reached at $`f=2`$

Exactly at the transition point, $`f=1`$, screening is nonlinear and the phase variation obeys equation

$$2\lambda _\phi ^2d^2\stackrel{~}{\phi }/dx^2=\stackrel{~}{\phi }^3.$$
(15)
The general solution of this equation is given by $`\stackrel{~}{\phi }=\pm 2\lambda _\phi /\left(x+C\right)`$. Using the phase-field relation $`H=(\mathrm{\Phi }_0/2\pi t)d\stackrel{~}{\phi }/dx`$, we derive from this solution that the magnetic field decays inside the structure as

$$H(x)=H_0\left(1+x/\sqrt{\frac{\lambda _\phi \mathrm{\Phi }_0}{\pi tH_0}}\right)^2$$
(16)
where $`H_0`$ is the external magnetic field. It is interesting to note that the field tail at large $`x`$ does not depend on the applied field

$$H(x)=\frac{\lambda _\phi \mathrm{\Phi }_0}{\pi tx^2}\text{, for }x\sqrt{\frac{\lambda _\phi \mathrm{\Phi }_0}{\pi tH_0}}.$$
(17)
In addition to the usual $`2\pi `$-degeneracy, in the case $`f<1`$, the ground state energy is also degenerate with respect to sign change of the equilibrium phase $`\phi _0\phi _0`$. This degeneracy leads to the appearance of two new kinds of solitons in which the phase sweeps either from $`\phi _0`$ to $`\phi _0`$ or from $`\phi _0`$ to $`2\pi \phi _0`$. Rewriting Eq. (8) in a more convenient form

$$\frac{d^2\phi }{dx^2}=\frac{1}{\lambda _\phi ^2}\mathrm{sin}\left(\phi \right)\left(\mathrm{cos}\phi _0\mathrm{cos}\phi \right),$$

we derive the shapes of the two solitons

$`\phi _1(x)`$ $`=2\mathrm{arctan}\left(\mathrm{tanh}\left({\displaystyle \frac{x\mathrm{sin}\phi _0}{2\lambda _\phi }}\right)\mathrm{tan}\left({\displaystyle \frac{\phi _0}{2}}\right)\right),`$ (18a)
$`\phi _2(x)`$ $`=\pi +2\mathrm{arctan}\left({\displaystyle \frac{\mathrm{tanh}\left(x\mathrm{sin}\phi _0/2\lambda _\phi \right)}{\mathrm{tan}\left(\phi _0/2\right)}}\right).`$ (18b)
The soliton energies per unit length are given by
$`\epsilon _1`$ $`=\epsilon _{J0}(\lambda _J/4\lambda _\varphi )\left(\left|\mathrm{sin}\phi _0\right|\phi _0\mathrm{cos}\phi _0\right),`$ (19a)
$`\epsilon _2`$ $`=\epsilon _{J0}(\lambda _J/4\lambda _\varphi )\left(\left|\mathrm{sin}\phi _0\right|+\left(\pi \phi _0\right)\mathrm{cos}\phi _0\right),`$ (19b)
where $`\epsilon _{J0}=4\mathrm{\Phi }_0j_c\lambda _J/(\pi c)`$ is the energy of a single soliton in a uniform JJ. At the transition point, $`f=1`$, the first soliton vanishes and the second one acquires the following shape, $`\phi _2(x)=\pi +2\mathrm{arctan}\left(x/\lambda _\phi \right)`$. Even though for $`d_0>d_\pi `$ the first energy is smaller than the second one, both solitons play a role in the magnetic properties of the $`\phi `$-junction. This is because without the second solitons the system can not contain two or more first solitons of the same sign due to topological constrains. In the external magnetic field $`H`$, the energy of the system with a small number of solitons of both kinds, $`N_1`$ and $`N_2`$, can be written as
$$E=\epsilon _1N_1+\epsilon _2N_2\left(N_1\frac{\phi _0}{\pi }+N_2\frac{\pi \phi _0}{\pi }\right)\frac{H\mathrm{\Phi }_0}{4\pi },$$
(20)
and has to be supplemented by the topological constrain $`\left|N_1N_2\right|1`$. The penetration scenario is determined by the ratio of the soliton energy $`\epsilon _\alpha `$ to its flux $`\mathrm{\Phi }_\alpha `$. Simple analysis shows that for $`d_0>d_\pi `$ these ratios are arranged in the order $`\epsilon _1/\mathrm{\Phi }_1<(\epsilon _1+\epsilon _2)/\mathrm{\Phi }_0<\epsilon _2/\mathrm{\Phi }_2`$. These inequalities imply that a single first soliton appears at the field $`H_1=4\pi ^2\epsilon _1/\phi _0\mathrm{\Phi }_0`$. The soliton lattice, composed of alternating solitons of two types, penetrates into the structure at higher field $`H_{c1}=4\pi (\epsilon _1+\epsilon _2)/\mathrm{\Phi }_0`$. In the case $`d_0=d_\pi `$, the fields $`H_1`$ and $`H_{c1}`$ coincide.

An interesting possibility to have an equilibrium $`\phi _0`$ Josephson vortex is realized at the boundary between an alternating structure of $`0`$\- and $`\pi `$\- JJs and usual $`0`$-JJ. We assume that the alternating structure occupies the region $`x<0`$, while the $`0`$-JJ with the Josephson penetration length $`\mathrm{\Lambda }`$ occupies a positive semi-axis $`x>0`$. In this case the phase difference must vary continuously from $`\phi _0`$ for $`x\mathrm{}`$ to zero for $`x+\mathrm{}`$. Due to the condition $`\lambda _\phi \mathrm{\Lambda }\lambda _J`$, the value $`\phi (0)1`$, which means that the boundary vortex is almost completely localized in region $`x<0`$ and it is equivalent to half of the first soliton (18a). The metastable $`\pi 2\phi _0`$ boundary vortex is also possible and in the limit $`\lambda _\phi \mathrm{\Lambda }`$ its shape corresponds to the second soliton (18b).

Consider a structure with the total length $`L`$ smaller that the screening length $`\lambda _{\mathrm{eff}}`$ in an external magnetic field $`H`$, smaller than $`\mathrm{\Phi }_0/(d_0+d_\pi )t`$. In this case the coarse-grained approximation for the phase is justified and we can use the current-phase relation (11) where the coarse-grained phase depends on the coordinate as $`\phi (x)=hx+\beta `$ with $`h=2\pi tH/\mathrm{\Phi }_0`$ and $`\beta `$ is an arbitrary phase shift. The total current per unit thickness flowing through the structure is given by

$$I(\beta )=\frac{j_cd_0^2}{12\lambda _J^2}\underset{L/2}{\overset{L/2}{}}𝑑x\left(f\mathrm{sin}\left(hx+\beta \right)\frac{\mathrm{sin}\left(2(hx+\beta )\right)}{2}\right).$$

Calculating the integral, we reduce it to the form

$$I(\beta )=\frac{I_0}{\eta }\mathrm{sin}\left(\beta \right)\mathrm{sin}\eta \left(f\mathrm{cos}\left(\beta \right)\mathrm{cos}\eta \right),$$
(21)
where $`I_0=j_cLd_0^2/(12\lambda _J^2)`$ and $`\eta =hL/2=\pi \mathrm{\Phi }/\mathrm{\Phi }_0`$ and $`\mathrm{\Phi }=tLH`$ is the total flux through the structure. The current reaches a maximum at

$$\mathrm{cos}\left(\beta \right)=\frac{\mathrm{cos}\phi _0\sqrt{\mathrm{cos}^2\phi _0+8\mathrm{cos}^2\left(\eta \right)}}{4\mathrm{cos}\left(\eta \right)}.$$
(22)
Shapes of the field dependencies of the critical current at different values of $`\phi _0`$ are plotted in the right panel of Fig. 3. These dependencies differ significantly from the Fraunhofer dependence in usual JJs. They have large a component with half of the main period. At $`\phi _0=\pi /2`$ the dependence has the Fraunhofer shape but the period is two times smaller than in usual case. Note that the observed in the alternating $`0`$-$`\pi `$ junction structure Smilde  finite critical current at zero magnetic field is naturally obtained in the framework of our analysis. The theoretical analysis in Ref. Smilde, , giving zero current at $`H=0`$, is incomplete, because it does not take into account the current term $`\mathrm{sin}(2\phi )`$ coming from the rapidly oscillating phase.

In conclusion, the possibility to realize the transition into the $`\phi `$junction state with decreasing the temperature from $`T_c`$, may be very helpful for experimental verification of the predicted effects. In particular, the observation of the striking non-monotonous variation of the screening length with temperature would provide an unambiguous proof of such a transition. Also, the studies of the fine structure of the critical current dependence vs magnetic field and the observation of a periodicity two times smaller than that expected for standard JJ, could be of considerable interest. Finally, note that the SQUID microscope may directly probe the specific shapes of the partial flux vortices. We would like to thank A. Lopatin and W. K. Kwok for critical reading of the manuscript. This work was supported by the U. S. DOE, Office of Science, under contract # W-31-109-ENG-38.


text:
# On the vacuum entropy and the cosmological constant1footnote 11footnote 1This essay received an “honorable mention” in the 2003 Essay Competition of the Gravity Research Foundation.

## Abstract

It is generally accepted that the entropy of an asymptotically de Sitter universe is bounded by the area, in Planck units, of the de Sitter horizon. Based on an analysis of the entropy associated to the vacuum quantum fluctuations, we suggest that the existence of such a holographic bound constitutes a possible explanation for the observed value of the cosmological constant, theoretically justifying a relation proposed $`35`$ years ago by Zel’dovich.

As extensively discussed in the literature 9 , there is a fundamental problem related to the existence of a positive cosmological constant. If $`\mathrm{\Lambda }`$ originates from vacuum energy, its expected value, on the basis of quantum field calculations with a cutoff given by the Planck energy, has the order of $`l_{Planck}^210^{70}`$m-2. This is $`122`$ orders of magnitude greater than the observed value $`\mathrm{\Lambda }10^{52}`$m-2 8 . Even considering smaller cutoffs, as the energy scales of electroweak or QCD phase transitions, the expected value is still over 40 orders of magnitude too high. This huge discrepancy is known as the cosmological constant problem.

We have shown elsewhere GS ; preprint  that this problem can be related to other open issues in cosmology, as the large numbers coincidence and the cosmic coincidence problem, with the help of the holographic principle 11 . In this essay, we will try to show that the application of this principle to an asymptotically de Sitter universe leads to a value for $`\mathrm{\Lambda }`$ in accordance with observation.

In a simplified form, the holographic principle can be described as the extension, to any gravitating system, of the Bekenstein-Hawking formula for black-hole entropy. More precisely, it establishes that the number of degrees of freedom of the system is not bounded by its volume in Planck units, as expected from quantum theories of space-time, but by the area, in Planck units, of its delimiting surface.

The use of this principle to the universe as a whole depends on the definition of such a surface at a cosmic scale. The existence of a positive cosmological constant naturally introduces a characteristic surface of radius $`\mathrm{\Lambda }^{1/2}`$ 18 . Another possibility is to use the Hubble horizon, with radius $`H^1`$, for it defines the scale of causal connections for any observer 17 . It is clear that in a homogeneous and isotropic, infinite universe, filled with dust and a positive cosmological constant, this last version implies the former, for the Hubble radius tends asymptotically to $`\sqrt{3/\mathrm{\Lambda }}`$. Therefore, in the context of an asymptotically de Sitter universe, the holographic bound establishes that the maximum number of available degrees of freedom is given by

$$N_{max}\mathrm{\Lambda }^1.$$
(1)
How much entropy has our universe? We know that the number of baryons is of the order of $`10^{80}`$, and that the cosmic background radiation contains about $`10^8`$ photons per baryon. It is reasonable to believe that dark matter contributes with a similar figure. But the major contribution to the entropy of matter seems to come from massive black-holes present in galactic nuclei, which represent an entropy of the order of $`10^{101}`$ Penrose .

But what about the entropy associated to the vacuum fluctuations? To put this question in a proper way and to clarify its role in solving the cosmological constant problem, let us inquire more carefully on the origin and meaning of this problem.

The difficulties appear when we calculate the vacuum energy density with the help of quantum field theories in flat space-time. The vacuum energy comes from two kinds of contributions. The first one is the energy associated to the vacuum expectation value of self-interacting fields, as the Higgs field, the quark and gluon condensates of QCD, or any other field associated to vacuum phase transitions. The second kind comes from the zero-point fluctuations of the fields, which lead to infinite results. To regulate them, it is used to impose some energy cutoffs, but, as said before, this leads to results still many orders of magnitude too high compared to observation.

Nevertheless, as discussed by some authors (see for example Ralf ), any contribution, infinite or not, to the vacuum energy density predicted by quantum field theories in flat space-time must be exactly canceled by a bare cosmological constant in the Einstein equations, because in the flat space-time the right-hand side of those equations is identically zero. Therefore, to properly pose the problem, we have to calculate the vacuum energy density in a curved background. As before, we find a divergent result as well. But now a physically meaningful (renormalized) value for $`\mathrm{\Lambda }`$ should be obtained by subtracting the Minkowskian result. Since the space-time of our universe is not strongly curved at cosmic scale, we expect to obtain a small value for $`\mathrm{\Lambda }`$, in accordance with observations.

The situation is analog to what occurs in the Casimir effect. There, the zero-point fluctuations of the electromagnetic field give rise to an infinite contribution to the vacuum energy density, inside and outside the region between the Casimir plates. But what is physically meaningful, leading to observable effects, is the difference between the values in the two regions, which is shown to be finite. In our case, the role of Casimir plates is played by gravity.

Calculating the vacuum energy density in a curved background is a difficult task. An example of a rough estimation in the line of the above reasoning was recently given by Schützhold Ralf . He estimates the contribution for $`\mathrm{\Lambda }`$ from the chiral anomaly of QCD in a curved, expanding space-time, obtaining $`\mathrm{\Lambda }H\mathrm{\Lambda }_{QCD}^3`$, where $`\mathrm{\Lambda }_{QCD}`$ is the energy scale of the chiral phase transition. In the limiting case of a de Sitter universe, $`H\sqrt{\mathrm{\Lambda }}`$, and his result leads to

$$\mathrm{\Lambda }m^6$$
(2)
(where we have made $`m=\mathrm{\Lambda }_{QCD}`$).

This expression was derived $`35`$ years ago by Zel’dovich, from empirical arguments Zel'dovich . Although it is sensible to the parameter $`m`$, (2) leads to the correct order of magnitude: using $`\mathrm{\Lambda }_{QCD}150`$ MeV, we obtain $`\mathrm{\Lambda }10^{51}`$ m-2, in good agreement with observation, considering that it was not taken into account numerical factors.

An alternative to circumvent the difficulties involved in quantum field calculations in a curved background is to use a thermodynamic approach, which does not depend on the details of the field dynamics. The idea is to obtain a superior limit for the vacuum entropy, instead of its energy density, and to compare the result with the holographic bound (1). The reader may argue that it is not trivial to define the number of virtual particles in curved backgrounds. Let us remind, however, that our universe has a quasi-flat space-time. Therefore, the estimation given below can be considered a good approximation.

It is clear that non-trivial vacuum configurations of classical fields (as the vacuum expectation value of the Higgs field or the QCD condensates) do not contribute to the vacuum entropy. In what concerns the zero-point fluctuations, they have, properly speaking, an infinite entropy density, because (if we do not impose any energy cutoff) the number of modes is infinite. But if we regulate their energy, by introducing an ultraviolet cutoff $`m`$, we also regulate their entropy. A simple estimation of the resulting entropy bound can be derived as follows.

Limiting the energy-momentum space associated to the zero-point fluctuations leads to the quantization of their configuration space, with a minimum size given by $`lm^1`$. This results in a superior bound to the number of available degrees of freedom in a given volume, say, the volume inside the Hubble horizon. The maximum number $`N`$ of observable degrees of freedom will be of the order of $`V/l^3`$, where $`V`$ is the Hubble volume. That is,

$$N\left(\frac{m}{H}\right)^3.$$
(3)
Now, if we take for $`H`$ the de Sitter asymptotic value $`H\sqrt{\mathrm{\Lambda }}`$, we can identify (3) with the holographic bound (1). It is easy to verify that this leads to Zel’dovich’s relation (2).

But why $`m`$ coincides to the energy scale of the QCD phase transition? The common belief is that a natural cutoff should be given by the Planck energy, for at the Planck scale the classical picture of space-time breaks down. Nevertheless, it is not difficult to see that, equating (1) to (3), with $`H\sqrt{\mathrm{\Lambda }}`$ and $`m=m_{Planck}`$, one obtains a Hubble radius of the order of $`l_{Planck}`$, which is not consistent with our universe.

One can also argue that the zero-point fluctuations of other fields than quarks and gluons contribute to the entropy as well. It is then intriguing that just $`\mathrm{\Lambda }_{QCD}`$ enters in Zel’dovich’s relation. Note, however, that in a curved space-time the different sectors of the standard model of particles interactions are coupled by gravity. On the other hand, the de Sitter universe is a stationary space-time, and, therefore, all the (interacting) vacuum fields should tend to a state of thermodynamic equilibrium, at the temperature of the last vacuum phase transition. But the last of such transitions was the chiral transition of QCD, at a temperature given by $`\mathrm{\Lambda }_{QCD}`$.

Finally, let us note that the vacuum entropy already dominates the universe entropy. Indeed, taking for $`H`$ the value observed nowadays, $`H65`$ km/(sMpc) 8 , we obtain from (3) (with $`m\mathrm{\Lambda }_{QCD}`$) $`N10^{122}`$, a figure that predominates over the matter entropy referred above, of order $`10^{101}`$.

I am thankful to F.S. Navarra, M.R. Robilotta, A. Saa, A.E. Santana, and R. Schützhold, for useful discussions. I am also grateful to Prof. P.F. González-Díaz, G.A. Mena Marugán and L. Garay, from CSIC (Madrid), where I have started to work with this subject.


text:
# Self-consistent random phase approximation in the Schütte-Da-Providencia fermion-boson model

## Abstract

The Self-consistent random phase approximation (SCRPA) is applied to the exactly solvable model with fermion boson coupling proposed by Schütte and Da-Providencia. Very encouraging results in comparison with the exact solution of the model for various observables are obtained. The transition from the normal phase to the phase with a spontaneously broken symmetry is carefully investigated. The strong reduction of the variance in SCRPA vs HF is pointed out.

During the last decade the so-called Self-Consistent version of the Random Phase Approximation (SCRPA) has seen very encouraging successes in a number of non trivial model cases (see for example b1  for a detailed description of the method and b2  for the application of SCRPA to the many level pairing model). In spite of these performances of the theory, there are remaining problems. In first place this concerns situations with spontaneously broken symmetries. Such situations were treated in b3 ; b4 ; b1 . Whereas in the Lipkin model b3  the symmetry broken (’deformed’) phase caused no problem because the broken symmetry is discrete (parity), in the other two cases b4 ; b1 , with a continuous broken symmetry, problems appeared with the low-lying mode known to be exactly at zero energy in the standard HF-RPA approach (the spurious or Goldstone mode). In the two cases cited b4 ; b1  the low-lying mode does not appear at zero energy in SCRPA because the RPA operator does not contain the symmetry operator as a limit case. Indeed, e.g. the number operator in quasiparticle (BCS) representation contains a purely hermitian piece $`\alpha _k^{}\alpha _k`$ which can not be incorporated in the RPA operator which by definition is non hermitian. The same situation is present in the Schütte-Da-Providenica boson-fermion model b5  where the symmetry operator contains the boson and fermion number operators. The violation of the Goldstone theorem signifies that the Ward identities and conservations laws are not respected. Though this violation seems relatively mild and to go away in macroscopic systems (the hermitian pieces becoming of zero weight), the situation remains annoying for finite systems.

In this note which can be considered as a sequel of b4  we want again to investigate the Schütte-Da-Providenica model

$$H=\overline{n}+\alpha b^{}b+G(\tau ^+b^{}+\tau ^{}b)$$
(1)
with $`b^{}`$, $`b`$ ideal boson operators and, $`n=_{i=1}^Na_{0i}^{}a_{0i}`$, $`\overline{n}=_{i=1}^Na_{1i}^{}a_{1i}`$, $`\tau ^+=_{i=1}^Na_{1i}^{}a_{0i}`$, $`\tau ^{}=_{i=1}^Na_{0i}^{}a_{1i}`$, $`\tau ^0=\frac{1}{2}(\overline{n}n)`$, where the $`a^{}`$, $`a`$ are fermion operators. In analogy to the work in b4  we will introduce the more general RPA operator

$$Q_\nu ^{}=X_\nu t^+Y_\nu t^{}+\lambda _\nu B^{}\mu _\nu B+U_\nu \beta ^{}\beta ^{}V_\nu \beta \beta ,\nu =1,2,3$$
(2)
where

$$t^\pm =\frac{T^\pm }{\sqrt{2T^0}}\text{and}\beta ^{}\beta ^{}=\frac{B^{}B^{}}{\sqrt{2(1+2B^{}B)}}.$$
(3)
The operators $`T^\pm `$$`T^0`$ are obtained from $`\tau ^\pm `$$`\tau ^0`$ by writing the latter ones in the ”deformed” basis

$$\left(\begin{array}{c}\alpha _{1k}^{}\\ \alpha _{0k}\end{array}\right)=\left(\begin{array}{cc}u& v\\ v& u\end{array}\right)\left(\begin{array}{c}a_{1k}^{}\\ a_{0k}^{}\end{array}\right),u^2+v^2=1.$$
(4)
The bosons operators $`B^{}`$ and $`B`$ are obtained from the original ones by a shift transformation $`Bb\sigma `$, where $`\sigma `$ is a $`c`$-number characterizing the appearance of the Bose condensate. The introduction of the boson pair operators $`\beta ^{}\beta ^{}`$ is motivated by the fact that otherwise there exists a certain dissymmetry between fermions and bosons, the fermions being in any case bilinear whereas the bosons are otherwise only contained to linear order in (2). Also the symmetry operator $`P=b^{}b\overline{n}`$ contains the bosons quadratically and the extended ansatz (2) may therefore show improved behavior with respect to the Goldstone mode. The formalism goes exactly in the same way as in b1 ; b2 ; b4  using the equation of motion method

$$[\delta Q,[H^{},Q_\nu ^{}]]=\mathrm{\Omega }_\nu [\delta Q,Q_\nu ^{}]$$
(5)
to determine the amplitudes in (2). As in b4 , in order to fix the value $`L=P`$, we use in (5) the cranked Hamiltonian $`H^{}=H\mu P`$ in the symmetry broken phase, otherwise $`H^{}=H`$. The mean-field amplitudes $`u`$, $`v`$ and $`\sigma `$ are readily obtained from a minimization of the ground-state energy, leading to $`[H^{},t^+]=[H^{},B^{}]=0`$. The amplitudes in (2) form a complete orthonormal set when calculated from (5). Then (2) can be inverted and with the usual condition for the RPA ground-state

$$Q_\nu |\text{RPA}=0,\nu =1,2,3.$$
(6)
all expectation values appearing in (5) like for example $`t^+B^{}`$, $`t^+B`$ and $`B^{}B`$ can directly be expressed in terms of the RPA amplitudes. The only unknown quantity at this point preventing a fully self-consistent solution of the SCRPA equations Eq. (5) is the expectation value $`T^0`$. However, in analogy to our previous study for the two-level pairing model b1  this quantity can be expressed as an expression in the operators $`T^+`$ and $`T^{}`$ up to any order in a fast converging series according to

$$T^0=\frac{N}{2}+\frac{1}{N}T^+T^{}+\frac{1}{N^2(N1)}T_{}^{+}{}_{}{}^{2}T_{}^{}{}_{}{}^{2}+\mathrm{}.$$
(7)
With this relation the SCRPA equations are completely closed and we can proceed to the numerical solution. We notice that with respect to b4  the SCRPA Eq. (5) is a $`6\times 6`$ dimensional problem whereas before it was $`4\times 4`$.

We now come to the presentation and discussion of the results. In what follows, SCRPA(6) refers to the SCRPA method with RPA excitation operator quadratic in the bosons operators i.e. Eq. (2), while SCRPA(4) refers to the same method but with RPA excitation operator linear in the bosons operators i.e. $`U_\nu =V_\nu =0`$. In the following we also use the set of parameters $`\alpha =3`$, $`N=30`$ as in b4 , for which in this model the phase transition point is localized at $`x=1.0`$, where $`x=G\sqrt{N/\alpha }`$. SCRPA always shows a clear superiority over standard RPA, though, besides for some quantities, the differences are not very pronounced. Concerning the ground-state energy we do not give results but we only notice that we arrive practically at the same interpretations as in b4 . However, in order to test the accuracy of our approach it is instructive to calculate the differences of energies of the ground-state band with $`L`$ values just one unit away from the absolute ground-state. One such quantity is the ”chemical potential” which should be identified with the Lagrange parameter used for restoring the symmetry

$$\mu =\frac{1}{2}\left(E_{L+1}^0E_{L1}^0\right).$$
(8)
In Table 1 we show $`\mu `$ when we calculate separately $`E_{L\pm 1}^0`$ (in standard RPA and SCRPA) and then take the difference. We also give in Table 1 the $`L`$ values which correspond for a given $`x`$ value to the absolute ground-state. In Table 1 we see a strong improvement of SCRPA(6) and SCRPA(4) over standard RPA and the high quality of the results in comparison with the exact values in the region around the phase transition point. Also SCRPA(6) is still improved over SCRPA(4). We could also have taken the $`\mu `$ values found from adjusting the correct $`L=P`$ values in the standard RPA and SCRPA calculations; we have checked numerically that the results are always practically identical.

Two other interesting quantities to be calculated within the SCRPA formalism and closely related to the chemical potential are the energy differences of the absolute ground-state with its ”left” and ”right” neighbors just one unit in $`L`$ away

$$\mathrm{\Delta }E_{\pm 1}=\pm \left(E_{L\pm 1}^0E_L^0\right).$$
(9)
These quantities are interesting because, as we will explain below, they should be closely related to the lowest RPA eigenvalue $`\mathrm{\Omega }_1`$ in the symmetry broken phase. Because we obtain similar interpretation and conclusions for both $`\mathrm{\Delta }E_{\pm 1}`$, we will present and discuss only the result for $`\mathrm{\Delta }E_1`$. In Fig. 1 we see a very good agreement with the exact results of both SCRPA(6) and SCRPA(4). However, we note that in this quantity no clear superiority of SCRPA(6) over SCRPA(4) can be detected, the results being at times in favor of the one or the other. The good quality of the results for $`\mathrm{\Delta }E_{\pm 1}`$ shows that the SCRPA method is able to reproduce the full spectrum. We also should notice that the smallness of $`\mathrm{\Delta }E_{\pm 1}`$ means that two neighboring ground-states with $`L`$ and $`L\pm 1`$, respectively, are almost degenerate which indicates that the system is in the phase of spontaneously broken symmetry. Furthermore, one can check that in the large $`N`$ limit $`\mathrm{\Delta }E_{\pm 1}`$ tends to zero. The zero eigenvalue (Goldstone mode) which is one of the solutions of the standard RPA in the deformed region corresponds to this degeneracy of the neighboring ground-state energies in the large $`N`$ limit.

Let us now discuss the eigenvalues of RPA and SCRPA matrices. As it is well known b5 ; b4 , in standard RPA the lowest eigenvalue corresponds in the symmetry broken (deformed) region to the spurious mode $`\mathrm{\Omega }_1=0`$, whereas the second eigenvalue gives the excitation of the intrinsic system. Before coming to this point we should mention again that the RPA eigenvalues in the deformed region are calculated with the ”intrinsic” Hamiltonian $`H^{}=H\mu P`$. Therefore, when the symmetry is restored due to the appearance of the Goldstone mode, the RPA eigenvalues give the excitation energies of the system. The results for the mode $`\mathrm{\Omega }_2`$ are not shown in a figure because we obtain the same interpretation as already given in b4 , however, with still improved results from SCRPA(6).

Concerning the low-lying eigenvalue of SCRPA which in standard RPA corresponds to the zero-energy eigenvalue (Goldstone or spurious mode) in the ”deformed” region, we present in Fig. 2 a comparison between standard RPA, SCRPA(4) and SCRPA(6) with exact results. In the ”spherical” phase we notice that the eigenvalue $`\mathrm{\Omega }_1`$ is identified with the exact ”intraband” excitation $`\mathrm{\Delta }E_1`$. Furthermore, we see the important improvement of the SCRPA results in both cases SCRPA(4) and SCRPA(6) with respect to standard RPA result. After the phase transition $`\mathrm{\Delta }E_1`$ remains finite but very small, slowly decreasing for increasing $`x`$, while the lowest eigenvalue in standard RPA corresponds to the spurious mode $`\mathrm{\Omega }_1=0`$. Concerning the low-lying eigenvalue in SCRPA calculation we see that SCRPA(6) improves the result with respect to SCRPA(4) but it is still quite far from the exact result. We can notice that the Goldstone theorem is not correctly fulfilled in this case. Therefore, the problem of the identification of the low-lying eigenvalue is not yet solved in SCRPA in spite of the introduction of the quadratic boson terms in (2).

From the presence of the quadratic boson terms in the RPA excitation operator the SCRPA method produces an supplementary eigenvalue which is noted $`\mathrm{\Omega }_3`$. Let us now discuss the results for this mode, i.e. the third eigenvalue RPA of the SCRPA(6) which is presented in Fig. 3. In the spherical region we notice that this mode is identified with very good accuracy to the exact ”intraband” excitation $`\mathrm{\Delta }E_{+2}`$. However, in the deformed region (not shown), the result obtained for this mode cannot be identified with one of the exact excitations of the system. At present we do not have an explanation of this fact. It is likely related to the failure of the Goldstone theorem mentioned above. See also further discussion of this point at the end of this note. The difficulty may be of the same origin as with the $`\mathrm{\Omega }_1`$ mode.

A quantity which is particularly sensitive to the correct treatment of correlations in the ground-state is the mean boson number (not shown in Ref. b4 ). This expectation value can be obtained in terms of the RPA amplitudes according to

$$N_b=b^{}b=B^{}B+\sigma ^2,$$
(10)
where $`B^{}B`$ is given, in SCRPA(6), by $`B^{}B=\mu _1^2+\mu _2^2+\mu _3^2`$. In Fig. 4 we show the results of the SCRPA(6), SCRPA(4), standard RPA and mean field methods for this quantity. Again with SCRPA(6) one notices a significant improvement over standard RPA and HF method for which the agreement with the exact result is not satisfying. Also, we note in Fig. 4 that the SCRPA(6) improves slightly the result over the SCRPA(4) specially in the deformed region. Furthermore, we note that in standard RPA method, because we have a Goldstone mode in the deformed region, we cannot calculate this quantity.

Another quantity which is very interesting to investigate in the SCRPA method is the variance of the symmetry operator $`P`$ given by

$$\mathrm{\Delta }P^2=P^2P^2.$$
(11)
In Fig. 5 we present the results corresponding to this quantity calculated with SCRPA and HF methods. This is a new result which was not elaborated in b4 . We see that the variance is strongly reduced compared to HF values. We, however, see that $`\mathrm{\Delta }P`$ even in the SCRPA acquires sizable non vanishing values. This simply means that the symmetry $`P`$, broken at the level of the mean field theory, is not completely restored. Furthermore, we do not present the standard RPA results concerning this quantity in Fig. 5 because the RPA amplitudes originating from the Goldstone mode are divergent. This constitutes the same kind of situation as that for the boson number in the ground-state calculation for which we also have not given the standard RPA results in deformed region.

In conclusion we reconsidered the work of Bertrand et al. b4  who treated the schematic Schütte-Da-Providencia model for interacting bosons and fermions within the SCRPA scheme. In b4  the RPA operator consisted only out of one boson and fermions pairs. Here we extended this configuration space and included in addition bosons pairs. One of the motivations to do this was to see whether the Goldstone theorem which was quite strongly violated in b4  is improved. It was found that the low-lying mode in the ”deformed” zone, i.e. symmetry broken region, indeed is lowered by $`30\%`$ when boson pair terms are added to the RPA operator. However, with respect to the first physical state, the position of the spurious mode is still too high and one therefore can not say that it decouples to a good approximation from the physical spectrum. However, in spite of this somewhat disappointing result, the introduction of the extra terms allowed to reproduce very well a further excited state of the spectrum in the symmetry conserved phase and additionally the quantities which had already been calculated in b4  without the boson pair operators are still improved. We also calculated further quantities as the number of bosons in the ground-state and the fluctuation of the symmetry operator. For instance the latter becomes strongly reduced with respect to its mean field value. However, also for this quantity a substantial non vanishing value remains. All in all we can say that the inclusion of the extra two boson piece to the RPA operator allowed to calculate one more state in the spherical region, improve existing results from calculations without these terms, and lower the spurious state. However, with respect to the latter feature no real breakthrough could be observed and further ideas are needed to substantially improve the situation in the symmetry broken region whereas the ”spherical” region seems to be well under control. In this respect this may be similar with other approaches treating correlations beyond mean field like coupled cluster theory, Jastrow, correlated basis function etc.

###### Acknowledgements.

We are specially grateful to Guy Chanfray and Peter Schuck for many interesting discussions and stimulating comments. We thank Raouf Bennaceur and Jorge Dukelsky for criticism and helpful comments. We gratefully acknowledge stimulating discussions with M. Oertel, D. Davesne and H. Hansen.