text:
# Discrete element simulations of stress distributions in silos: crossover from two to three dimensions

## I Introduction

The stress within a silo packed with granular material has long been of interest in the engineering Janssen1895 ; Nedderman1992  and physics JaegerOct1996  communities. As early as 1895, Janssen, in his pioneering work, constructed a model to describe the vertical stress in a silo Janssen1895 . Treating the granular pack as a continuous medium where a fraction $`\kappa `$ of vertical stress is converted to horizontal stress, Janssen was able to derive a simple functional form for the vertical stress. One main assumption of this model is that the forces of friction between particles and walls are at the Coulomb failure criterion: $`F_t=\mu _wF_n`$, where $`F_t`$ is the magnitude of the tangential friction force, $`F_n`$ is the normal force at the wall, and $`\mu _w`$ is the coefficient of friction for particle-wall contacts. This assumption is also known as incipient failure. This theory qualitatively describes the crossover to a depth-independent vertical stress, although quantitative discrepancies between the Janssen theory and experiment have been observed. The assumption of incipient failure is believed to be the source of the discrepancies, but this has never been conclusively proven due to the difficulty of experimental measurements of internal stress. Numerous improvements have been added to the theory over time, but none of these have gained wide acceptance Nedderman1992 .

Recently, well-controlled experiments have been carried out on granular packs in silos to test the suitability of Janssen’s theory in ideal conditions VanelMay1999 ; VanelFeb2000 . Vanel and Clément VanelFeb2000  constructed their initial packings by tapping a loose poured packing to increase its density. The base of the packing was then slowly moved downward to more fully mobilize the grains, and then the packing was allowed to settle. They measured the apparent mass at the bottom of the silo as a function of the filling mass of the silo. Their experiment did not follow the Janssen form, and they found the best agreement with a phenomenological model, which contains elements of Janssen’s original analysis. We describe this model in more detail in Sec III. Other more recent experimental studies in which the side walls were dragged upward to fully mobilize the packing OvarlezJun2003 ; BerthoApr2003 ; OvarlezSep2003  have found agreement with the Janssen form.

A variety of discrete element and continuous finite element simulation methods have been developed to describe the stresses in a silo due to the difficulty of experimental verification. Most discrete element simulations thus far have been confined to two dimensional ($`2D`$) systems. However, there is wide disagreement within the community on the predictive power of these models and the proper approach to take for accurate simulation MassonApr2000 ; HolstJan1999a ; HolstJan1999b ; SanadOct2001 . Those simulations that are carried out in three dimensions ($`3D`$) usually do so via finite-element methods that yield little information on the internal structure or forces in granular packs ChenOct2001 ; GuinesOct2001 . Some $`3D`$ discrete element simulations MakseMay2000 ; SilbertMar2002  have been performed, but these use periodic boundary conditions in the two directions perpendicular to gravity. Though these studies provide useful information on the internal structure of these packings, they can give no information on vertical stresses or forces induced by the walls of the container. We have recently carried out $`3D`$ discrete element simulations of stress in granular materials in cylinders LandryApr2003 . Here we compares simulations in $`2D`$ and $`3D`$ and study the crossover from two to three dimensions.

We present large-scale discrete element simulations of granular packings in $`2D`$, quasi-$`2D`$, and $`3D`$ containers (silos). Our aim is to understand the vertical stress profiles of these granular packings, particularly the crossover from two to three dimensions, which has largely been unexplored. Both monodisperse and polydisperse particle systems are studied. We explicitly test the suitability of the Janssen theory to the vertical stress profiles produced by our simulations and test the validity of the Janssen theory’s assumptions. We find that the interior of the packing is far from incipient failure, while in general, the forces at the walls are close to incipient failure. In all cases, there is a region at the top of the pile where the particle-wall forces are far from incipient failure and the Janssen form is not observed.

Section II discusses the simulation technique and the method used to generate the packings. Section III presents the vertical stress profiles. We discuss their characteristics and compare our results to the classical theory of Janssen as well as two modified forms of the Janssen analysis. In Section IV we show how forces in the packings are not at incipient failure at the top and in the bulk of the packing and relate this to the failure of the Janssen analysis to fully explain our packing. Finally, we conclude and summarize the work in Section V.

## II Simulation Method

We present discrete element simulations in two and three dimensions of model systems of $`N`$ mono and polydispersed spheres of fixed density $`\rho `$. The mass of a given particle is $`m_i=\frac{\rho }{6}\pi d_i^3`$, where $`d_i`$ is the diameter of a given particle. Particle diameters are uniformly distributed between $`\overline{d}\pm \mathrm{\Delta }`$. For this study, we use $`\mathrm{\Delta }/\overline{d}=0`$, $`0.1`$, and $`0.5`$. $`N`$ varies from 10,000 to 50,000 particles. The system is constrained by either an open rectangular box or an open cylinder with its axis along the vertical $`z`$ direction. The box has length $`L`$, width $`W`$, and height $`H`$, centered at $`x=y=0`$ and bounded below with a flat base at $`z=0`$. The cylinder has radius $`R`$ and is centered on $`x=y=0`$. It is bounded below with a flat base at $`z=0`$. In some cases, a layer of randomly-arranged immobilized particles approximately $`2d`$ high were used to provide a rough base. We vary the length $`L`$ and width $`W`$ of the box and the radius $`R`$ of the cylinder. This work builds on previous 3D simulations of packings in a cylinder LandryApr2003 .

The spheres interact only on contact through a spring-dashpot interaction in the normal and tangential directions to their lines of centers. Contacting spheres $`i`$ and $`j`$ positioned at $`𝐫_i`$ and $`𝐫_j`$ experience a relative normal compression $`\delta =|𝐫_{ij}d|`$, where $`𝐫_{ij}=𝐫_i𝐫_j`$, which results in a force Cundall1979

$$𝐅_{ij}=𝐅_n+𝐅_t.$$
(1)
The normal and tangential contact forces are given by

$$𝐅_n=f(\delta /d)(k_n\delta 𝐧_{ij}m_{eff}\gamma _n𝐯_n)$$
(2)
$$𝐅_t=f(\delta /d)(k_t𝚫𝐬_tm_{eff}\gamma _t𝐯_t)$$
(3)
where $`𝐧_{ij}=𝐫_{ij}/r_{ij}`$, with $`r_{ij}=|𝐫_{ij}|`$. $`𝐯_n`$ and $`𝐯_t`$ are the normal and tangential components of the relative surface velocity, and $`k_{n,t}`$ and $`\gamma _{n,t}`$ are elastic and viscoelastic constants, respectively. $`m_{eff}=\frac{m_im_j}{m_i+m_j}`$ for interactions between particles $`i`$ and $`j`$, and $`m_{eff}=m_i`$ for interactions between particle $`i`$ and a wall, where the mass of the wall is assumed to be infinite. $`f(x)=1`$ for Hookean (linear) contacts, while for Hertzian contacts $`f(x)=x^{D/21}`$ — the two force models differing only in $`3D`$. $`𝚫𝐬_t`$ is the elastic tangential displacement between spheres, obtained by integrating tangential relative velocities during elastic deformation for the lifetime of the contact. The magnitude of $`𝚫𝐬_t`$ is truncated as necessary to satisfy a local Coulomb yield criterion $`F_t\mu F_n`$, where $`F_t|𝐅_t|`$ and $`F_n|𝐅_n|`$ and $`\mu `$ is the particle-particle friction coefficient. Frictionless spheres correspond to $`\mu =0`$. Particle-wall interactions are treated identically, though the particle-wall friction coefficient $`\mu _w`$ is set independently. A more detailed description of the model is available elsewhere SilbertOct2001 . We also tested an alternate form of the $`F_t`$ force due to Haff and Werner Haff1986  that does not have history effects. This model produces a hydrostatic stress profile in all cases. Without the relative transverse displacement in the force model, $`F_t=0`$ in the quasi-static limit and the walls cannot support stress.

We present all the physical quantities in our simulations in terms of $`\overline{d}`$, the average diameter; $`\overline{m}`$ the mass of a particle with average diameter $`\overline{d}`$; and $`g`$, the force of gravity. Most of these simulations were run with $`k_n=2\times 10^5\overline{m}g/\overline{d}`$, $`k_t=\frac{2}{7}k_n`$, and $`\gamma _n=50\sqrt{g/\overline{d}}`$. For Hookean springs we set $`\gamma _t=0`$. In this case, the coefficient of restitution for $`\mathrm{\Delta }=0`$ is $`ϵ=0.8`$. For Hertzian springs, $`\gamma _t=\gamma _n`$. For Hertzian springs in $`D=3`$, $`ϵ0`$ as $`v0`$, so that it takes much longer for energy to drain out of the pack and for the kinetic energy per particle to decrease to sufficiently low levels. For this reason, it is much more computationally expensive to use Hertzian springs to form packings. The convenient time unit is $`\tau =\sqrt{\overline{d}/g}`$, the time it takes a particle to fall its radius from rest under gravity. For this set of parameters, the time step $`\delta t=10^4\tau `$. In most simulations, the particle-particle friction and particle-wall friction are the same, $`\mu =\mu _w=0.5`$.

The two-dimensional ($`2D`$) simulations were performed by fixing $`y=0`$ for all particles and removing any velocity or acceleration components in the $`y`$ direction. Some polydispersity ($`\mathrm{\Delta }/\overline{d}>0`$) is required to prevent crystallization in $`2D`$, so all $`2D`$ simulations were run with a polydispersity of $`\mathrm{\Delta }/\overline{d}=0.1`$.

All of our results will be given in dimensionless units based on $`\overline{m}`$, $`\overline{d}`$, and $`g`$. Physical experiments often use glass spheres of $`\overline{d}=100\mu m`$ with $`\rho =2\times 10^3kg/m^3`$. In this case, the physical elastic constant would be $`k_{glass}10^{10}\overline{m}g/\overline{d}`$. A spring constant this high would be prohibitively computationally expensive, because the time step must have the form $`\delta tk^{\frac{1}{2}}`$ for collisions to be modeled effectively. We have found that increasing $`k`$ in our simulations does not appreciatively change the physical results SilbertOct2001 .

We generate our packings by mimicking the pouring of particles at a fixed height $`Z`$ into the container. For computational efficiency a group of $`M`$ particles is added to the simulation on a single time step. This is done by inserting the $`M`$ particles at non-overlapping positions within a thin cylindrical region of radius $`R\overline{d}`$ with height $`Z\overline{d}<z<Z`$. The $`x`$, $`y`$, and $`z`$ coordinates of the particles are chosen randomly within this insertion region. The height of insertion $`z`$ determines the initial z-velocity $`v_z`$ of the particle — $`v_z`$ is set to the value it would have after falling from a height $`Z`$. After a time $`\sqrt{2}\tau `$, another group of $`M`$ particles is inserted. This methodology generates a steady stream of particles, as if they were poured continuously from a hopper (see Figure 1). The rate of pouring is controlled by setting $`M`$ to correspond to a desired volume fraction of particles within the insertion region. For example, for an initial volume fraction of $`\varphi _i=0.13`$ in a cylinder of $`R=10d`$, the pouring rate is $`45`$ particles/$`\tau `$. This method is similar to the homogeneous “raining” methods used in experiments VanelNov1999 .

The simulations were run until the kinetic energy per particle was less than $`10^8\overline{m}g\overline{d}`$. The resultant packing is considered quiescent and used for further analysis SilbertMar2002 . In $`3D`$, Hertzian packings have difficulty reaching this criteria due to long-lived pressure waves that traverse the entire pile. In these packings, the kinetic energy per particle tends to stabilize at approximately $`10^3\overline{m}g\overline{d}`$ for very long times. By adding a viscous damping term proportional to $`v`$ after the pack has stabilized at $`10^3\overline{m}g\overline{d}`$, we were able to reduce the kinetic energy per particle of these systems to less than $`10^8\overline{m}g\overline{d}`$. A comparison of the stress profiles before and after the application of this viscous damping term showed no change. Although a change of kinetic energy of this order has a small effect on the coordination number of the packing, it does not have any effect on the vertical stress. All vertical stress profiles for Hertzian packings presented in the following sections were generated with this technique.

These simulations were performed on a parallel cluster computer built with DEC Alpha processors and Myrinet interconnects using a parallel simulation code optimized for short-range interactions SilbertOct2001 ; PlimptonMar1995 . A typical simulation to pour $`50,000`$ particles into a 3D cylindrical silo requires $`5\times 10^6`$ time steps, which requires roughly $`40`$ CPU hours on $`50`$ processors.

We show in Figure 1 how a typical packing is created. The figure shows three snapshots in time as particles are poured into a quasi-$`2D`$ system of length $`30\overline{d}`$ and width $`5\overline{d}`$.

The strong effect of wall friction on the capacity of a silo to support stress as found experimentally is shown in Figure 2 where the particle-wall friction $`\mu _w`$ is varied, while keeping the particle-particle friction constant at $`\mu =0.5`$. The geometry of the system is a three-dimensional U-tube, with both ends open. Particles are poured into one end and the resultant packing is allowed to settle. In the case where $`\mu _w=0`$, no stress can be supported by the side walls, and the granular packing is liquid-like: the particle height is the same in both sides of the U-tube. For $`\mu _w=0.1`$, the right side is considerably lower than the left, though some particles are higher than the bend in the tube. The walls support some of the stress, but the pressure is strong enough to force some particles up the right-hand tube. For $`\mu _w=\mu =0.5`$, there are almost no particles on the right side of the U-tube, and all the weight on the left side is supported by the walls of the tube. In all these cases, the important factor is the particle-wall friction, as the particle-particle friction remains unchanged.

## III Stresses

Of particular interest in the construction of silos is the distribution of stresses Nedderman1992 . In a liquid, hydrostatic pressure increases with depth. Granular materials support shear, so the side walls of a container can support some of this shear, provided $`\mu _w>0`$. The problem of the resultant vertical stress in a silo after filling has a long history, beginning with Janssen in 1895. Janssen’s analysisJanssen1895 ; Rayleigh1906  is still in use today, although it rests on a series of assumptions which have not been fully tested.

For a container with static wall friction $`\mu _w`$ and granular pack of total height $`z_0`$, the Janssen analysis predicts the vertical stress $`\sigma _{zz}(z)`$ at a height $`z`$ is

$$\sigma _{zz}(z)=\rho gl\left[1\mathrm{exp}\left(\frac{z_0z}{l}\right)\right]$$
(4)
where $`\rho `$ is the volumetric density, $`l`$ is the decay length, and $`z_0`$ is the top of the packing. The decay length $`l`$ is determined by the geometry. For a $`2D`$ container of length $`L`$, $`l_{2D}=\frac{L}{2\kappa \mu _w}`$. For a quasi-$`2D`$ container of length $`L`$ and width $`W`$, $`l_{quasi2D}=\frac{LW}{2L+2W}\frac{1}{\kappa \mu _w}`$. For a $`3D`$ cylindrical container of radius $`R`$, $`l_{3D}=\frac{R}{2\kappa \mu _w}`$ Duran2000 .

Another two-parameter fit was proposed by Vanel and Clément VanelMay1999  to reconcile their experimental findings on weakly shaken packings with Janssen theory. The fit assumes a region of perfect hydrostaticity, followed by a region that conforms to the Janssen theory.

$`z_0z<a`$ $`:`$ $`\sigma _{zz}(z)=\rho g(z_0z)`$
$`z_0z>a`$ $`:`$ $`\sigma _{zz}(z)=\rho g\left(a+l\left[1\mathrm{exp}\left({\displaystyle \frac{z_0za}{l}}\right)\right]\right)`$ (5)
The two parameters are $`a`$, the height of the crossover, and $`l`$ the decay length, with the same values for different geometries as in the Janssen theory.

We first present vertical stresses for the two-dimensional ($`2D`$) case. These are shown in Figure 3. Unlike previous studies HolstJan1999b , these piles are extremely deep so that the height-independent region can be easily observed. The Janssen theory does not describe the $`2D`$ stress profiles well. There is a clear linear hydrostatic-like region at the top of the packing that the Janssen theory does not account for. The Vanel-Clément form is a much closer fit to the data. In Figure 3, the deepest pile of length $`30\overline{d}`$ is fit to both forms. The Janssen form predicts $`l/\overline{d}=84.6`$ and thus $`\kappa =0.354`$, while the Vanel-Clement form yields $`a=35`$ and $`l/\overline{d}=49.6`$ and thus $`\kappa =0.605`$. As we will see below, there are a number of differences between $`2D`$ and $`3D`$ packings. One difference is the height of the non-saturated region, which is on the order of $`57L`$ in $`2D`$ and $`6R`$ in $`3D`$. There is about a factor of 2 difference between the diameter or width of the container in $`2D`$ compared to $`3D`$. In addition, the saturated pressure is also much larger in $`2D`$ than the $`3D`$ case. Both of these effects seem to arise from the necessity of all the force being carried by the side walls. A similar effect is seen at the bottom of the packing, where the slight increase in the pressure at the bottom of the pile is over a much larger region than in the $`3D`$ case as shown below. This increase occurs because geometrically the walls can not support forces from particles below this height.

To study the crossover from $`2D`$ to $`3D`$, we studied a series of packings in rectangular containers with increasing width, ranging from $`W=1.2\overline{d}`$ to $`5\overline{d}`$. We show in Figure 4 the vertical stress profiles for these quasi-2D packings. The stress profiles show the same hydrostatic region followed by saturation as the $`2D`$ profiles, though the hydrostatic region at the top of the pile is much smaller in the quasi-$`2D`$ case. In addition, the saturation stress value does not increase monotonically with width. The saturation stress decreases as one goes from a width of $`1.2\overline{d}`$ to $`1.4\overline{d}`$, but increases again as one reaches $`2\overline{d}`$. A width of $`2\overline{d}`$ seems to represent a transition to a more $`3D`$-like behavior, and the saturation stress does not change appreciably for larger widths. This suggests that $`L`$ is the dominant length scale for widths greater than $`2\overline{d}`$, while for smaller widths, the width $`W`$ is the dominant length scale for determination of the vertical stress.

We also investigated quasi-$`2D`$ packings with no friction on the long walls $`\mu _w=0`$, while the short side walls have $`\mu _w=0.5`$. This system was modeled to understand how the long side walls affect the geometry. This system, as shown in Figure 5, produces stress profiles very different from the quasi-$`2D`$ case with friction on all walls. Two widths were used and both show much larger final stress values and much larger hydrostatic regions. This suggests that both sets of walls are important for determining the stress profiles for a wide range of widths. The stress profiles seen in this case resemble much more the $`2D`$ stress profiles presented in Fig 3. These results demonstrate that experiments in thin plates meant to approximate $`2D`$ systems are more representative of $`3D`$ systems than that of pure $`2D`$ systems.

We show in Figure 6 vertical stress profiles for $`3D`$ systems in cylindrical packings. Changing the force model from a Hookean to a Hertzian interaction has a very small effect on the final packing. Changing from a monodisperse to a very highly polydisperse system ($`\mathrm{\Delta }/\overline{d}=0.5`$) also does not change the overall stress profile, although it does increase the height of the packing. This suggests that although polydispersity is required in $`2D`$ systems to avoid crystallization, it has little effect on the stress profile. It is thus reasonable to compare monodisperse $`3D`$ packings to polydisperse $`2D`$ systems.

In this case we see that the important transition is from the $`2D`$ system to the quasi-$`2D`$ system. In most respects the quasi-$`2D`$ system is much closer to the $`3D`$ system than the pure $`2D`$ system. The size of the hydrostatic region changes most markedly from the pure $`2D`$ system to the quasi-$`2D`$ system. The stress values in the depth of the pile are much larger for the $`2D`$ case than the quasi-$`2D`$ case and are comparable to the $`3D`$ case.

## IV Distribution of Forces

The previous section has demonstrated that the stress profiles of our granular packings do not conform to the Janssen theory. The historical assumption has been that the assumption of incipient failure in the Janssen analysis is responsible for the failure of the theory. We thus tested directly the Janssen assumptions of incipient failure in our simulations, i.e. whether the tangential forces at the wall are actually at the Coulomb yield criterion $`F_t=\mu _wF_n`$. We define $`\zeta =F_t/\mu F_n`$ in the bulk of the packing and $`\zeta =F_t/\mu _wF_n`$ for forces at the wall. If a specific force is at the Coulomb failure criterion, $`\zeta =1`$. We analyzed the $`3D`$ systems shown in Figure 6. By examining the distribution of forces in the interior of our packings, we find that no particle-particle contacts are at the Coulomb criterion irrespective of method or level of polydispersity, as shown in Figure 7a. We find similar results in the interior of $`2D`$ packings. However the particle-wall forces in the height-independent stress region are much closer to the Coulomb criterion. In the four systems examined, the bulk of the particle-wall tangent forces are close to incipient failure, in contrast to the particle-particle tangent forces in the bulk, as seen in Figure 7b. The amount of polydispersity and the contact force model have essentially no effect on how close the system is to incipient failure. We have also investigated the role of $`\mu _w>\mu `$ in Ref. LandryApr2003 . In this case, the particle-wall forces are always less than the Coulomb criterion, as the walls cannot support larger tangential forces than those supported in the bulk.

As discussed above, all of our packings exhibit a linear region in the stress profile at the top of the packing. We can also examine the Janssen assumption in this linear region. In contrast to the situation in the height-independent stress region of the packing, near the top of the pile the forces at the wall are far from the Coulomb criterion, as shown in Figure 8. This explains why the Janssen form is not obeyed in this region: the walls in this region do not support stress and thus the stress profile is hydrostatic.

## V Conclusions

We examined the differences between vertical stress in granular packings in $`2D`$, quasi-$`2D`$, and $`3D`$. We focused on the vertical stress of these packings and observed that $`2D`$ packings support much less vertical stress than $`3D`$ packings. In contrast, quasi-$`2D`$ packings (those packings between thin walls) rapidly converge to the behavior of $`3D`$ packings as the wall thickness is increased, though their progression is not monotonic. We also show that polydispersity and different force models do not have a large effect on the vertical stress in $`3D`$ packings. In all cases, we show that the Janssen assumption of incipient failure is not justified everywhere in the packing.

We acknowledge discussions with L. Silbert and Jin Ooi, and thank L. Silbert for his critical reading of the manuscript. This work was supported by the Division of Materials Science and Engineering, Basic Energy Sciences, Office of Science, U.S. Department of Energy. This collaboration was performed under the auspices of the DOE Center of Excellence for the Synthesis and Processing of Advanced Materials. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.


format: text

text:
# Current carrying Andreev bound states in a Superconductor-Ferromagnet proximity system

## Abstract

We study the ground state properties of a ferromagnet-superconductor heterostructure on the basis of a quasiclassical theory. We have solved the Eilenberger equations together with Maxwell’s equation fully self-consistently and found that due to the proximity effect a Fulde-Ferrel-Larkin-Ovchinnikov ($`FFLO`$) like state is realized in such system. Moreover this state has oscillations of the pairing amplitude in either one or two directions, depending on the exchange splitting and thickness of the ferromagnet. In particular, using semiclassical arguments (Bohr-Sommerfeld quantization rule) we show that owing to the presence of the Andreev bound states in the ferromagnet, a spontaneous current in the ground state is generated as a hallmark of the $`FFLO`$ state in the direction parallel to the interface. We also discuss the effects of the the elastic disorder and finite transparency of the interface on the properties of the $`FFLO`$ state in the system.

As is well known, ferromagnetism and singlet pairing superconductivity are competing phenomena. Whilst an exchange interaction favors parallel spin aligment, Cooper pairs must be in a singlet (spin) state. But the two phenomena are less mutually exclusive in artificially made ferromagnet-superconductor ($`FM`$/$`SC`$) heterostructures Chien . In such structures the ferromagnetism and the superconductivity can coexist near the $`FM`$/$`SC`$ interface owing to the proximity effect Lambert . In the case when the normal metal is not a ferromagnet, the proximity effect has been studied for a long time Lambert . By now it is rather well understood in terms of the Andreev reflection processes Andreev . By contrast, the proximity effect in $`FM`$/$`SC`$ systems has became a centre of attention only recently. It is not only important from a scientific point of view, as it allows for study of the interplay between magnetism and superconductivity Berk , but also from a technological one, as it may find applications in magnetoelectronics Bauer  and quantum computing Blatter .

A number of new phenomena has been revealed in $`FM`$/$`SC`$ multilayers. The most interesting examples are: non-monotonic behavior of the $`SC`$ transition temperature Wong , oscillations of a pairing amplitude Buzdin ; Radovic ; Demler  and the density of states in the $`FM`$ Kontos ; Zareyan ; Buzdin\_1 ; KGA , paramagnetic Meissner effect Lee , proximity induced very long range triplet superconductivity in $`FM`$ Volkov  or generation of spontaneous currents in the ground state of such systems KGA . These unusual properties, associated with Cooper pairs in an exchange field, can be explained in terms of an phenomenon first identified by Fulde, Ferrell, Larkin and Ovchinnikov ($`FFLO`$) FFLO . Originally it has been studied in a bulk superconductor with the exchange splitting. It turns out that, although of great conceptual interest, the bulk $`FFLO`$ state can be realized only in a very small region of the parameter space near the transition to the normal state FFLO  and usually is destroyed when the exchange splitting $`E_{ex}`$ is larger than $`\sqrt{2}/2\mathrm{\Delta }`$ (Clogston criterion) Clogston , where $`\mathrm{\Delta }`$ is the $`SC`$ energy gap. Moreover the $`FFLO`$ state is very sensitive to both elastic and spin-orbit scattering Takada . The last two effects make the $`FFLO`$ state very difficult to observe experimentally in a bulk samples. The situation is much more favorable in $`FM`$/$`SC`$ heterostructures where, due to the proximity effect, the Cooper pairs can survive even if the exchange field in $`FM`$ is much larger than $`SC`$ gap.

According to our current understanding of the $`FFLO`$ phenomenon in a $`FM`$/$`SC`$ structure, when a Cooper pair enters the ferromagnet it acquires a center of mass momentum $`\mathrm{}Q=2E_{ex}/v_F`$ Demler , where $`v_F`$ is Fermi velocity. Usually the $`SC`$ phase changes linearly with distance $`x`$ from the interface, $`\varphi =Q_xx`$, and this results in a oscillatory behavior of the pairing amplitude in the $`FM`$. It turns out that under certain conditions a $`3D`$-$`FFLO`$ state, featuring a spatial dependence of the pairing amplitude also along the interface, can be realized Izyumov\_1 ; KGA . It is this latter case that we shall deal with here.

The purpose of the present paper is to demonstrate that in a ferromagnetic layer on a superconducting substrate, for particular values of the exchange splitting $`E_{ex}`$ and the layer thickness $`d_F`$, a $`FFLO`$ like state is realized with the pairing amplitude $`f`$ varying both perpendicular and parallel to the $`FM`$/$`SC`$ interface. It will be shown that such a ground state supports a spontaneously generated current flowing in opposite directions in the $`FM`$ and the $`SC`$ regions. The existence of this remarkable state was first predicted on the basis of a simple lattice model KGA . Here we shall deal with the problem by a less model dependent, semiclassical approach and address the issue of the observability of the phenomenon in the presence of disorder within $`FM`$ layer and at the $`FM`$/$`SC`$ interface.

The system we consider is sketched in Fig 1. It consists of thin ferromagnet ($`FM`$) of thickness $`d_F`$ deposited on a semi-infinite superconductor ($`SC`$) and bounded on the other side by an insulator.

In such an $`I`$/$`FM`$/$`SC`$ quantum well there will be bound states corresponding to the closed quasiparticle trajectories deGennes . Each trajectory consists of an electron segment, $`e`$, which includes an Andreev reflection at the $`FM`$/$`SC`$ interface and an ordinary reflection at the $`I`$/$`FM`$ interface plus a hole segment, $`h`$, retracing back the electron trajectory (see Fig 1). The Bohr-Sommerfeld quantization rule Kevin :

$`{\displaystyle _b^a}p_e(\omega )𝑑l{\displaystyle _a^b}p_h(\omega )𝑑l+\delta \varphi \gamma (\omega )=2n\pi `$ (1)
gives the energies of the bound states. The first and second terms represent the total phase accumulated during propagation through the $`FM`$ region from $`b`$ to $`a`$ ($`=4d_F(\omega +\sigma E_{ex})/v_Fcos(\theta )`$), $`\delta \varphi `$ is the phase difference between points $`b`$ and $`a`$, and $`\gamma (\omega )=arccos(\omega /\mathrm{\Delta })`$ is the Andreev reflection phase shift. An example of the $`\theta `$-dependence of the Andreev bound state ($`ABS`$) energies for $`\xi _F/d_F=0.425`$, where $`\xi _F=\mathrm{}v_F/E_{ex}`$ is the $`FM`$ coherence length, and $`\delta \varphi =0`$ is shown in the inset of Fig. 1. Clearly, for any exchange splitting $`E_{ex}0`$ it is possible to find such a $`\theta `$ that the corresponding $`ABS`$ is exactly at zero energy. If there is a large number of such zero-energy $`ABS`$ for some exchange splittings (more precisely for $`\xi _F/d_F`$ ratio) the density of states $`DOS`$ has a large peak at the Fermi energy ($`ϵ_F=0`$). Such a situation turns out to be energetically unfavorable. There is a number of mechanisms which split this peak and thereby lower the energy of the system Sigrist . One of these is a spontaneous current which ’Doppler’ shifts the quasiparticle energies by $`\delta =ev_FA_ycos(\theta )`$ Higashitani , where $`A_y`$ is a vector potential in the $`y`$ direction.

An example of such a density of states is depicted in Fig. 2, where one can see large peak at zero energy (solid line).

This corresponds to zero $`SC`$ phase difference between points $`a`$ and $`b`$ in Fig. 1, namely no spontaneous current. In the inset of Fig. 2 the energy difference between state with spontaneous current $`E_{GS}(\varphi )`$ and state with no current $`E_{GS}(0)`$ is shown. The corresponding $`DOS`$, with the zero energy $`ABS`$ split, is plotted with the dashed line in Fig. 2.

To summarize, we have shown, using a Bohr-Sommerfeld semiclassical argument, that for any exchange splitting, $`E_{ex}`$, there are Andreev bound states at zero energy and for certain $`E_{ex}`$ the number of such states is so large that it produces huge zero energy peak in the density of states. Such a peak in turn is split by a spontaneous current and this lowers the total energy of the system. This current carrying state can be regarded as a realization of the $`FFLO`$ variation of the pairing amplitude in the $`y`$-direction. So, one can say that the system can be switched between $`1D`$ and $`2D`$ $`FFLO`$-like states as the exchange field or thickness of the ferromagnet is changed. In the following we will show that this spontaneous current can be also obtained within a self-consistent quasiclassical theory. This approach allows for a treatment of disorder and finite transparency of the interface, and hence provides further useful insights from the experimental point of view.

The quasiclassical matrix Eilenberger equation Eilenberger  reads

$`𝐯_F\widehat{g}_\sigma (𝐯_F,𝐫)+`$
$`[\stackrel{~}{\omega }_\sigma (𝐫)\widehat{\tau }_3+\widehat{\mathrm{\Delta }}(𝐫)+{\displaystyle \frac{1}{2\tau }}\widehat{g}_\sigma (𝐯_F,𝐫),\widehat{g}_\sigma (𝐯_F,𝐫)]_{}=0`$ (2)
where

$`\widehat{g}_\sigma =\left(\begin{array}{cc}g_\sigma \hfill & \hfill f_\sigma \\ f_\sigma ^+\hfill & \hfill g_\sigma \end{array}\right),\widehat{\mathrm{\Delta }}=\left(\begin{array}{cc}0\hfill & \hfill \mathrm{\Delta }\\ \mathrm{\Delta }^{}\hfill & \hfill 0\end{array}\right)`$ (7)
and

$`\stackrel{~}{\omega }_\sigma (𝐫)=\omega +i\sigma E_{ex}+ie𝐯_F𝐀(𝐫)`$ (8)
Here $`\widehat{\tau }_3`$ is the Pauli matrix, $`\omega =\pi T(2n+1)`$ is the Matsubara frequency, $`\sigma `$ ($`=\pm 1`$) labels the electron spin, $`𝐀(𝐫)`$ is the vector potential and $`\mathrm{}`$ denotes averaging over directions of the Fermi velocity $`𝐯_F`$. The matrix Green’s function has to obey the normalization condition $`\widehat{g}_\sigma ^2(𝐯_F,𝐫)=\widehat{1}`$. The exchange splitting $`E_{ex}`$ is non-zero and constant in the ferromagnet only while $`\mathrm{\Delta }(𝐫)`$ is non-zero in the superconductor and is calculated self-consistently from

$`\mathrm{\Delta }(𝐫)=U\pi \rho _0T{\displaystyle \underset{\omega ,\sigma }{}}f_\sigma (𝐯_F,𝐫),`$ (9)
where we have assumed that the coupling constant $`U<0`$ in the $`SC`$ and $`=0`$ in the $`FM`$. $`\rho _0`$ is the normal state $`DOS`$ and $`T`$ stands for temperature.

The spontaneous current has to be determined self-consistently together with the Maxwell equation (Ampere’s law), which couples the electron current to the magnetic field. The total current in the $`y`$-direction at each point $`x`$ measured from the interface is given by

$`J_y^{tot}(x)=2ie\pi \rho _0T{\displaystyle \underset{\omega ,\sigma }{}}𝐯_Fg_\sigma (𝐯_F,x),`$ (10)
where $`e`$ is the electron charge, while the Maxwell equation (in the Landau gauge) reads

$`{\displaystyle \frac{d^2A_y(x)}{dx^2}}=\mu _0J_y^{tot}(x)`$ (11)
with $`\mu _0`$ being the permeability of free space.

We have solved the Eilenberger equation (2) numerically along each quasiparticle $`2D`$ trajectory using the Riccati parametrization (Schopohl-Maki transformation) Riccati  together with the self-consistency relations (9)-(11).

The most remarkable feature of the self-consistent solution is that the iterations of the Eilenberger equations frequently converge to a solution with a finite value of the current even though there is no external vector potential. The current flows in one direction over the whole ferromagnet and flows back on $`SC`$ side on the scale of $`SC`$ coherence length $`\xi _S`$ (see solid curve in the Fig.3), so the total current is zero, as it should be in the true ground state. Such a distribution of the current can be associated with a quasiparticle current in $`FM`$ and with a supercurrent on the $`SC`$ side mediated by Cooper pairs (see Fig. 1). The fact that the current flows over the whole $`FM`$ is due to the extended nature of the $`ABS`$. There is also magnetic flux associated with such current distribution. Typically the spontaneous magnetic field produced by this current is of order of $`0.1B_{c2}`$, where $`B_{c2}`$ is the upper critical field of the bulk $`SC`$.

Within the present self-consistent calculations we were also able to study the effect of elastic scattering in the $`FM`$. The disorder in the ($`s`$-wave) $`SC`$ can be neglected due to Andreson’s theorem Moradian . The current for a number of mean free paths is shown in Fig. 3.

As we one can see disorder introduces oscillations of the current. The spontaneous current is proportional to the $`DOS`$ at the Fermi energy KGA  and so the oscillations of the current are related to the oscillatory behavior of $`DOS`$ in the disordered sample Zareyan ; Buzdin\_1 . In the clean limit the $`DOS`$ is constant in the whole $`FM`$. As is well known, this is a property of the Eilenberger equations in the clean limit Buzdin\_1 . Moreover disorder also suppresses the current, as expected, since it introduces decoherence of electron-hole pairs in the $`FM`$. Finally if the mean free path $`l`$ is shorter than the $`FM`$ thickness $`d_F`$ the current is completely suppressed. If $`l<d_F`$ the Andreev reflected particles cannot reach the $`I`$/$`FM`$ interface, which is a necessary condition for the formation of the current carrying $`ABS`$, because they are scattered on the impurities and the electron-hole coherence is lost. In this regard the $`FFLO`$ variation of the pairing amplitude in the $`y`$-direction is very sensitive to the elastic disorder. However in the $`x`$-direction $`FFLO`$ state persists until $`l<\xi _F`$ Demler , even if $`l<d_F`$.

To take into account the effect of specular reflections at the $`FM`$/$`SC`$ interface we adapt the approach proposed by M. Zareyan et al. Zareyan , where a certain probability distribution was associated with each semiclassical trajectory (for details see Ref. Zareyan ). In Fig. 4 we show the current for two different transparencies $`0<\eta <1`$.

As we would expect transparency $`\eta <1`$ suppressed the current because it suppresses Andreev reflection processes and at the same time introduces normal (specular) reflections at the $`FM`$/$`SC`$ interface. But even for $`\eta =0.3`$, in the present case, we still get the current. It turns out that $`\eta `$ has similar influence on the properties of the system as $`E_{ex}`$ does. As we can read from the inset of Fig. 4, it changes the period of oscillations of the pairing amplitude. Moreover if we changed $`\eta `$ only, for certain its values we get solutions with a current flowing, whilst for others solutions with no current. For the above set of parameters we have found that current flows in the regions where $`0.7<\eta <1`$ and $`0.2<\eta <0.35`$. So we can also switch between the $`1D`$ and $`2D`$ $`FFLO`$ state by changing the transparency.

Finally we note that so far we took no account of the demagnetization field $`B_d`$ due to the ferromagnet. In the limit of a zero film width this is zero. For finite thickness, $`d_F`$, and magnetization, $`M`$, (which in the Stoner model Mohn  corresponds to $`E_{ex}=IM`$, where $`I`$ is a fenomenological parameter) we estimate $`B_d`$ to be considerably less then $`B_{sp}`$ due to the spontaneous current. Furthemore it should be stressed that in all our calculations the magnetization was constrained to point in the $`z`$ direction. Thus the direction of the spontaneous current was determined by the condition that $`B_{sp}`$ is parallel to $`M`$. In a more general theory where Andreev orbits also occur in the $`x`$-$`z`$ plane and spin orbit coupling is taken into account these issues would need to be reexamined. Also the above calculations were two dimensional, but preliminary studies of a $`3D`$ system indicate that there are no qualitative changes when orbits in the $`x`$-$`z`$ plane are included.

In summary we have demonstrated that under certain conditions the ground state of a $`I`$/$`FM`$/$`SC`$ trilayer features a spontaneous current flowing in opposite directions in the $`FM`$ and $`SC`$ layers. We argued that this state can be viewed as a $`2D`$ $`FFLO`$ proximity state and hence the observation of the above current would be a decisive proof that the surprising behavior of such heterostructures is governed by the $`FFLO`$ phenomenon. We also showed that this state persists only in the clean limit where the mean free path is longer than $`FM`$ thickness and investigated the effect of low transparency of the interface on the observibility of the ground state current.

One of us, M.K., would like to thank Prof. Yuli Nazarov for helpful discussion on the problem of transparency of the interface. This work has been supported by Computational Magnetoelectronics Research Training Network under Contract No. HPRN-CT-2000-00143.


text:
# A hybrid constraint programming and semidefinite programming approach for the stable set problem

## 1 Introduction

This paper describes a hybrid method to solve the maximum weighted stable set problem, or stable set problem111 Alternative names for the stable set problem are vertex packing, coclique or independent set problem. in short. Given a graph with weighted vertices, the stable set problem is to find a subset of vertices of maximum weight, such that no two vertices in this subset are joined by an edge of the graph. In the unweighted case (when all weights are equal to 1), this problem amounts to the maximum cardinality stable set problem, which has been shown to be already NP-hard .

We propose a two-phase approach to solve the stable set problem, either with or without proving optimality. The first phase generates subproblems based upon a semidefinite relaxation, the second phase solves the subproblems using constraint programming. Concerning the first phase, given a model for the stable set problem, we solve its semidefinite relaxation. The solution provides us fractional values for the variables of the model. These fractional values are a good indication for the optimal (discrete) values of the variables. Hence we divide selected variable domains in two parts: a ‘good’ subdomain and a ‘bad’ subdomain. By branching on these subdomains using a limited discrepancy strategy , we obtain first a very promising subproblem, and subsequently less promising subproblems.

The second phase consists of the solution of the subproblems. Since they are much smaller than the original problem, we can easily solve them using a constraint programming solver.

As computational results will show, the semidefinite relaxation is quite informative. In several cases we can simply round the solution of the relaxation and obtain a provable optimal solution already. Otherwise, we are likely to find a good solution in one of the first subproblems. Using a limited number of subproblems to investigate, we yield an incomplete method producing good solutions. In order to obtain a complete search strategy, we need, in principle, to generate and solve all possible subproblems. A good upper bound is necessary to prove optimality earlier. For this reason we investigated the use of ‘discrepancy cuts’ that can be added to the semidefinite program to strengthen the relaxation and thus prune large parts of the search tree. However, computational results will show that they can not be applied efficiently on the instances we considered.

The outline of the paper is as follows. The next section gives a motivation for the approach proposed in this work. Then, in Section 3 some preliminaries on semidefinite programming are given. In Section 4 we introduce the stable set problem, integer optimization formulations and a semidefinite relaxation. A description of our solution framework is given in Section 5. Section 6 presents the computational results. This is followed by an overview of related literature in Section 7. Finally, in Section 8 we conclude and discuss future directions.

## 2 Motivation

Let us first motivate why a semidefinite relaxation is used rather than a linear relaxation. Indeed, one could argue that linear programs are being solved much faster in general. However, for the stable set problem, linear relaxations are not very tight. Therefore one has to identify and add inequalities that strengthen the relaxation. But it is time consuming to identify such inequalities, and by enlarging the model the solution process may slow down.

Several papers on approximation theory following  have shown the tightness of semidefinite relaxations. However, being tighter, semidefinite relaxations are more time-consuming to solve than linear relaxations in practice. Hence one has to trade strength for computation time. For large scale applications, semidefinite relaxations can often be preferred as the relaxation of choice to be used in a branch and bound framework. Moreover, our intention is not to solve a relaxation at every node of the search tree. Typically, we only solve a relaxation at the root node of the first phase (its solution is used to identify the subdomains), and optionally at the root node of each subproblem (in order to strengthen the upper bound). Therefore, we are willing to make the trade-off in favour of the semidefinite relaxation.

Another point of view is the following. Although semidefinite programming has been developing for many years now in the operations research community, no efforts of integration or cooperation with constraint programming have been made to our knowledge. Application of semidefinite programming to problems typical to constraint programming, as was done in the papers on approximation algorithms mentioned in Section 7, is not yet hybrid problem solving. In this paper a first step is being made. The solution of the semidefinite relaxation is used to identify promising subdomains, and also produces a tight upper bound for the constraint programming solver. On the other hand, the solutions found by the constraint programming solver serve as a lower bound inside the semidefinite programming solver.

## 3 Preliminaries on semidefinite programming

Combinatorial optimization problems that are NP-hard are often solved with the use of a relaxation. Typically, the integral restriction on the decision variables is relaxed to obtain a continuous relaxation, which can be polynomially solvable. Here we introduce semidefinite programming as an extension of the more common linear programming. Both paradigms can be used as a continuous relaxation.

In linear programming, combinatorial optimization problems are modeled in the following way:

$$\begin{array}{cc}\hfill \mathrm{max}& c^𝖳x\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& Axb\hfill \\ & x0\hfill \end{array}$$
(1)
Here $`x^n`$ is an $`n`$-dimensional vector of decision variables and $`c^n`$ a cost vector of dimension $`n`$. The ($`m\times n`$) matrix $`A^{m\times n}`$ and the $`m`$-dimensional vector $`b^m`$ define $`m`$ linear constraints on $`x`$. In other words, this approach models problems using nonnegative vectors of variables.

Semidefinite programming makes use of positive semidefinite matrices of variables instead of nonnegative vectors. A matrix $`X^{n\times n}`$ is said to be positive semidefinite (denoted by $`X0`$) when $`y^𝖳Xy0`$ for all vectors $`y^n`$. Semidefinite programs have the form

$$\begin{array}{ccc}\hfill \mathrm{max}& \mathrm{tr}(CX)\hfill & \\ \hfill \mathrm{s}.\mathrm{t}.& \mathrm{tr}(A_jX)=b_j\hfill & (j=1,\mathrm{},m)\hfill \\ & X0.\hfill & \end{array}$$
(2)
Here $`\mathrm{tr}(X)`$ denotes the trace of $`X`$, which is the sum of its diagonal elements, i.e. $`\mathrm{tr}(X)=_{i=1}^nX_{ii}`$. The cost matrix $`C^{n\times n}`$ and the constraint matrices $`A_j^{n\times n}`$ are supposed to be symmetric. The $`m`$ reals $`b_j`$ and the $`m`$ matrices $`A_j`$ define again $`m`$ constraints.

We can view semidefinite programming as an extension of linear programming. Namely, when the matrices $`C`$ and $`A_j(j=1,\mathrm{},m)`$ are all supposed to be diagonal matrices, the resulting semidefinite program is equal to a linear program. In particular, then a semidefinite programming constraint $`A_jX=b_j`$ corresponds to the linear programming constraint defined by the $`j`$-th row of $`A`$ and the $`j`$-th entry of $`b`$.

Applied as a continuous relaxation, semidefinite programming in general produces solutions that are much closer to the integral optimum than linear programming. Untuitively, this can be explained as follows. The problem is ‘lifted’ to a higher dimension, where one can exploit more of the structure of the problem. Unfortunately, it is not a trivial task to catch this structure, and otain a good semidefinite program for a given problem.

Theoretically, semidefinite programs have been proved to be polynomially solvable using the so-called ellipsoid method (see for instance ). In practice, nowadays fast ‘interior point’ methods are being used for this purpose (see  for an overview). Being a special case of semidefinite programming, linear programs are also polynomially solvable using an interior point method. However, they are often solved with a special linear programming solver, the Simplex method. Although this method can have an exponential behaviour in theory, in practice it is the fastest method known for solving linear programs.

## 4 The stable set problem

In this section, the stable set problem is formally defined, and formulated in two different ways as an integer optimization problem. From this, a semidefinite programming relaxation is inferred.

### 4.1 Definition

Consider an undirected weighted graph $`G=(V,E)`$, where $`V`$ is the set of vertices and $`E`$ the set of edges $`\{(i,j)|i,jV,ij\}`$ of $`G`$. Let $`|V|=n`$ and $`|E|=m`$. To each vertex $`iV`$ a weight $`w_i`$ is assigned. A stable set is a set $`SV`$ such that no two vertices in $`S`$ are joined by an edge in $`E`$. The stable set problem is the problem of finding a stable set of maximum total weight in $`G`$. This value is called the stable set number of $`G`$ and is denoted by $`\alpha (G)`$. The maximum cardinality (or unweighted) stable set problem can be obtained by taking all weights equal to 1.

### 4.2 Integer optimization formulation

Let us first consider an integer linear programming formulation. We introduce binary variables to indicate whether or not a vertex belongs to the stable set $`S`$. So, for $`n`$ vertices, we have $`n`$ integer variables $`x_i`$ indexed by $`iV`$, with initial domains $`\{0,1\}`$. In this way, $`x_i=1`$ if vertex $`i`$ is in the stable set $`S`$, and $`x_i=0`$ otherwise. We can now state the objective function, being the sum of the weights of vertices that are in the stable set $`S`$, as $`_{i=1}^nw_ix_i`$. Finally, we define the constraints that restrict two adjacent vertices to be both inside $`S`$ as $`x_i+x_j1`$, for all edges $`(i,j)E`$. Hence the integer linear programming model becomes:

$$\begin{array}{cccc}\hfill \alpha (G)=& \hfill \mathrm{max}& \hfill _{i=1}^nw_ix_i& \\ & \hfill \mathrm{s}.\mathrm{t}.& \hfill x_i+x_j1& (i,j)E\hfill \\ & & \hfill x_i\{0,1\}& iV.\hfill \end{array}$$
(3)
Another way of describing the same solution set is by the following integer quadratic program

$$\begin{array}{cccc}\hfill \alpha (G)=& \hfill \mathrm{max}& \hfill _{i=1}^nw_ix_i& \\ & \hfill \mathrm{s}.\mathrm{t}.& \hfill x_ix_j=0& (i,j)E\hfill \\ & & \hfill x_i^2=x_i& iV.\hfill \end{array}$$
(4)
Note that here the constraint $`x_i\{0,1\}`$ is replaced by $`x_i^2=x_i`$. This quadratic formulation will be used below to infer a semidefinite programming relaxation of the stable set problem.

In fact, both model (3) and (4) can be used as a constraint programming model. We have chosen the first model, since the quadratic constraints take more time to propagate than the linear constraints, while having the same pruning power.

### 4.3 Semidefinite programming relaxation

The integer quadratic program (4) gives rise to a semidefinite relaxation introduced by Lovász , see Grötschel et al.  for a comprehensive treatment. The value of the objective function of this relaxation has been named the theta number of a graph $`G`$, indicated by $`\vartheta (G)`$. Let us start again from model (4). As was indicated in Section 3, we want to transform the current model that uses a nonzero vector into a model that uses a positive semidefinite matrix to represent our variables. In the current case, we can construct a matrix $`X^{n\times n}`$ by defining $`X_{ij}=x_ix_j`$. Let us also construct a $`(n\times n)`$ cost matrix $`W`$ with $`W_{ii}=w_i`$ for $`iV`$ and $`W_{ij}=0`$ for all $`ij`$. Since $`X_{ii}=x_i^2=x_i`$, the objective function becomes tr$`(WX)`$. The edge constraints are easily transformed as $`x_ix_j=0X_{ij}=0`$. The first step in the transformation of model (4) can now be made:

$$\begin{array}{ccc}\hfill \mathrm{max}& \hfill \mathrm{tr}(WX)& \\ \hfill \mathrm{s}.\mathrm{t}.& \hfill X_{ij}=0& (i,j)E\hfill \\ & \hfill x_i^2=x_i& iV.\hfill \end{array}$$
(5)
This model is still a quadratic program, although reformulated. The problem remains how to model the last, very important, constraint. We need a mapping of the diagonal entries $`X_{ii}=x_i^2`$ with the vector entries $`x_i`$. For this, we extend $`X`$ with another row and column that contain vector $`x`$ (both indexed by 0). For this reason, define the $`(n+1)\times (n+1)`$ matrix $`Y`$ as

$$Y=\left(\begin{array}{cccc}1& x_1& \mathrm{}& x_n\\ x_1& & & \\ \mathrm{}& & X& \\ x_n& & & \end{array}\right)$$

where the 1 in the leftmost corner of $`Y`$ is needed to obtain positive semidefiniteness. In this case we can express the required mapping as $`Y_{ii}=\frac{1}{2}Y_{i0}+\frac{1}{2}Y_{0i}`$ (note that $`X`$ and $`Y`$ are symmetric), since then $`x_i^2=Y_{ii}=\frac{1}{2}Y_{i0}+\frac{1}{2}Y_{0i}=x_i`$. The final step in the transformation consists of replacing the constraints on $`X`$ by constraints on $`Y`$. In particular, instead of demanding $`X`$ to be a product of nonnegative vectors, we restrict $`Y`$ to be a positive semidefinite matrix. Namely, if the vector $`x`$ represents a stable set, then the matrix $`Y`$ is positive semidefinite. However, not all positive semidefinite $`Y`$ matrices represent a stable set, in particular its values can take fractional values.

In order to maintain equal dimension to $`Y`$, a row and a column indexed 0 should be added to $`W`$, all entries of which containing value 0. Denote the resulting matrix by $`\stackrel{~}{W}`$. The theta number of a graph $`G`$ can now be described as

$$\begin{array}{cccc}\hfill \vartheta (G)=& \hfill \mathrm{max}& \mathrm{tr}(\stackrel{~}{\mathrm{W}}\mathrm{Y})& \\ & \hfill \mathrm{s}.\mathrm{t}.& Y_{ii}=\frac{1}{2}Y_{i0}+\frac{1}{2}Y_{0i}& (i=1,\mathrm{},n)\hfill \\ & & Y_{ij}=0& (i,j)E\hfill \\ & & Y0.& \end{array}$$
(6)
By construction, the diagonal value $`Y_{ii}`$ serves as an indication for the value of variable $`x_i`$ $`(i\{1,\mathrm{},n\})`$ in a maximum stable set. In particular, this program is a relaxation for the stable set problem, i.e. $`\vartheta (G)\alpha (G)`$. Note that program (6) can easily be rewritten into the general form of program (2). Namely, $`Y_{ii}=\frac{1}{2}Y_{i0}+\frac{1}{2}Y_{0i}`$ is equal to tr$`(AY)`$ where the $`(n+1)\times (n+1)`$ matrix $`A`$ consists of all zeroes, except for $`A_{ii}=1`$, $`A_{i0}=\frac{1}{2}`$ and $`A_{0i}=\frac{1}{2}`$, which makes the corresponding $`b`$ entry equal to 0. Similarly for the edge constraints.

The theta number also arises from other formulations, different from the above, see . In our implementation we have used the formulation that has been shown to be computationally most efficient among those alternatives . Let us introduce that particular formulation. Again, let $`x\{0,1\}^n`$ be a vector of binary variables representing a stable set. Define the $`n\times n`$ matrix $`X`$ as $`X=\frac{1}{x^𝖳x}xx^𝖳`$. Furthermore, let the $`n\times n`$ cost matrix $`W`$ be defined as

$$W_{ij}=\{\begin{array}{cc}w_i\hfill & i=j\hfill \\ \sqrt{w_iw_j}\hfill & ij\hfill \end{array}$$

for $`i,jV`$. Then the following semidefinite program

$$\begin{array}{cccc}\hfill \vartheta (G)=& \hfill \mathrm{max}& \hfill \mathrm{tr}(WX)& \\ & \hfill \mathrm{s}.\mathrm{t}.& \hfill \mathrm{tr}(X)=1& \\ & & \hfill X_{ij}=0& (i,j)E\hfill \\ & & \hfill X0& \end{array}$$
(7)
gives exactly the theta number of a graph. Again, it is not difficult to rewrite this program into the general form of program (2).

Program (7) uses matrices of dimension $`n`$ and $`m+1`$ constraints, while program (6) uses matrices of dimension $`n+1`$ and $`m+n`$ constraints. This gives an indication why program (7) is computationally more efficient.

## 5 Solution framework

### 5.1 Overview

The two-phase solution approach proposed here is similar to the one described in . In the first phase subproblems are generated, which are being solved in the second phase. A general overview of the method is presented in Algorithm 1.

Let us first explain how we use the solution of the semidefinite program (7) to partition the domain $`D_i`$ of a variable $`x_i`$ into $`D_i^{\mathrm{good}}`$ and $`D_i^{\mathrm{bad}}`$ (for $`iV`$). As was stated before, the solution of program (7) assigns fractional values between 0 and 1 to the variables. Naturally, if for a variable $`x_i`$ the corresponding fractional value $`Y_{ii}`$ is close to 1, we regard 1 to be a good value for variable $`x_i`$. More specifically, we select the variable $`x_i`$ with the highest fractional value $`Y_{ii}`$, set $`D_i^{\mathrm{good}}=\{1\}`$ and $`D_i^{\mathrm{bad}}=\{0\}`$, and mark it as handled. Then we mark all its neighbours $`j`$ (with $`(i,j)E`$) as being handled, keeping their original domain $`D_j=\{0,1\}`$. This procedure is repeated until all variables are handled. For convenience, we partition $`V`$ into two distinct sets $`V_0=\{iV|D_i^{\mathrm{good}}=\{1\}\}`$ and $`V_1=VV_0`$. Here $`V_1`$ represents the set of neighbours of $`V_0`$, with $`D_j=\{0,1\}`$.

In a similar way, we can use the solution of the semidefinite relaxation to compute a first feasible integer solution. Namely, follow the same procedure, but now instantiate the selected variable $`x_i=1`$ and set its neighbours $`x_j=0`$. The objective value of this feasible integer solution is in many cases already equal to the (downward rounded) solution of the semidefinite relaxation. In that case, we have found an optimal solution and finish. In other cases, we can still use this first solution as a lower bound to be applied during the solution of the subproblem.

Next, we explain how to generate subproblems using these subdomains. The generation of subproblems makes use of a tree structure of depth $`|V_0|`$ in which we branch on $`D_i^{\mathrm{good}}`$ versus $`D_i^{\mathrm{bad}}`$. The tree is traversed using a limited discrepancy strategy (LDS) . LDS visits the nodes of a search tree differently from depth-first search. It tries to follow a given suggestion as good as it can. Branches opposite to the suggestion are regarded as discrepancies and are gradually allowed to be traversed. The first ‘run’ of LDS doesn’t allow any discrepancies, the second allows one, and so on. This means that for a particular discrepancy $`k`$, a path from the root to a leaf is allowed to consist of maximally $`k`$ right branches. Typically this method is applied until a limited number of discrepancies is reached (say 2 or 3), which yields an incomplete search strategy. In order to be complete, one has to visit all nodes, up to discrepancy $`d`$ for a binary tree of depth $`d`$.

In our case, the suggestion that should be followed are branches of the kind $`D_i^{\mathrm{good}}`$, while branches $`D_i^{\mathrm{bad}}`$ are regarded as discrepancies. Hence, our first subproblem is the subproblem defined by program (3), and with $`x_iD_i^{\mathrm{good}}`$ for all $`iV_0`$. The next $`|V_0|`$ subproblems have all $`x_iD_i^{\mathrm{good}}`$, except for one $`x_kD_k^{bad}`$ ($`i,kV_0`$). The next discrepancy generates $`\frac{1}{2}(|V_0|^2+|V_0|)`$ subproblems, each of which contains two variables $`x_{k_1}D_{k_1}^{bad}`$ and $`x_{k_2}D_{k_2}^{bad}`$ ($`k_1,k_2V_0`$), and so on. Since we expect to obtain a very good solution already in the first subproblems, we will only generate subproblems up to a certain maximum discrepancy. In our experiments the maximum discrepancy is chosen 2 and 4 respectively.

Note that the first subproblem, corresponding to discrepancy 0, only contains one solution, namely the one that we obtain in our rounding procedure. By propagation of the edge constraints, all variables $`x_jD_j=\{0,1\}`$ are instantiated automatically to 0. Hence, the subproblem corresponding to discrepancy 0 is obsolete in our current implementation.

### 5.2 Adding discrepancy cuts

In the case one needs to generate and solve subproblems up to a large discrepancy, it is preferable to proof possible suboptimality of a subproblem before entering it. Especially when the subproblems are still relatively large. This can be done in several ways.

First, before entering a subproblem, we can identify variables which have a subdomain of size 1, namely those with $`iV_0`$. For those variables, one can add an extra constraint to the semidefinite program (7), being either $`x_i=1`$ or $`x_i=0`$. Then the semidefinite program can be solved again, and will in general provide a tighter bound, hopefully lower than the current lowerbound, in which case we have proven suboptimality. However, solving the semidefinite program each time before entering a subproblem is very time-consuming and this method will not be very practical.

A better alternative would be to add a specific constraint, a discrepancy cut, that is valid for all subproblems of a given discrepancy. Recall that $`V_0=\{iV|D_i^{\mathrm{good}}=\{1\}\}`$. Hence, all subproblems of discrepancy $`k`$ consist of $`k`$ variables $`x_i`$ with $`x_i=0`$, and $`|V_0|k`$ variables $`x_i`$ with $`x_i=1`$ $`(iV_0)`$. This gives rise to two discrepancy cuts, given discrepancy $`k`$:

$`{\displaystyle \underset{iV_0}{}}x_i=|V_0|k`$ (8)
$`{\displaystyle \underset{iV_0}{}}1x_i=k`$ (9)
We implemented both of them, and cut (9) gives the best results. Stated in terms of semidefinite program (7), the discrepancy cut looks like tr$`(AX)=|V_0|k`$, with $`A_{ij}=1`$ if either $`iV_0`$ or $`jV_0`$ or both, and $`A_{ij}=0`$ otherwise. As mentioned before, solving a semidefinite relaxation is relatively expensive, and one should make a tradeoff between its computation time and the gain in time of not solving the subproblem. For the instances we considered, the time needed to solve a semidefinite program is always larger than the time needed to solve all subdomains we would like to proof suboptimal. However, these cuts might be helpful for larger instances.

## 6 Computational results

Our experiments are being done on a Pentium 1GHz processor, with 256 Mb RAM, using Windows XP. As constraint programming solver we use the ILOG solver library, version 5.1 . As semidefinite programming solver, we use CSDP version 4.1 .

The first instances we consider are randomly generated weighted graphs with $`n`$ vertices and $`m`$ edges. The vertex weights are randomly chosen from a range of 1 up to $`n`$. The edge density is chosen such that the constraint programming solver has difficulties solving them. Namely, the more edge constraints we have, the easier it is solved by constraint programming. The name of the instances represent the number of vertices and the edge density, i.e. g75d015 is a graph on 75 vertices with an approximate edge density of 0.15. For these graphs, we have chosen to generate subproblems up to a maximum discrepancy of 4, based upon earlier experience.

We also considered structured instances, obtained from problems arising in coding theory . These are unweighted graphs, therefore we have set all weights equal to 1. For these graphs, we generate subproblems up to a discrepancy of 2.

The results of our experiments are given in Table 1. It consists of three parts: the first part describes the instances, the next part gives the results of our approach (sdp and cp), the last part concerns the results of a sole constraint programming approach (cp alone).

The columns in this table represent the following. An instance name has $`n`$ vertices and $`m`$ edges. For the part on our approach, the value of the semidefinite relaxation is $`\vartheta `$, the rounded solution of the semidefinite relaxation has value round, and best is the value of the best solution found. This best solution is found in a subproblem generated during discrepancy best discr. Note that we generate subproblems up to discrepancy 4 in all cases, as was mentioned in Section 5. The time spent on solving the semidefinite relaxation is denoted by time sdp. The time spent on solving all generated subproblems is denoted by time subp. These values together form the total time. All times are measured in seconds. The number of all backtracking steps made during the search in our approach is collected in backtracks. Concerning the sole constraint programming approach, we report the best solution found (best), the total time spent during search, and the total number of backtracks. Note that we have set a time limits for the constraint programming solver, to create a fair comparison with our approach. They are 100 seconds for g50d005 up to g150d010, 190 seconds for g150d015 and 324 seconds for the structured instances. Best found solutions that are proven to be optimal are indicated by an asterisk (\*).

For the instances in Table 1 we only solved one semidefinite relaxation per problem, namely at the root node. The reason for this is that the time spent during the subproblem search is less than the time spent on computing another relaxation, as reported in the table. Therefore, we cannot gain time by adding discrepancy cuts and computing another semidefinite relaxation.

In general, our method produces better bounds than the constraint programming approach alone. In many cases, the rounded solution of the semidefininite relaxation is already optimal. However, note that there are two instances for which the constraint programming approach gives better bounds. Note also that the structured instances are handled quite well by our approach, while the constraint programming approach produces very low quality solutions for the larger instances.

## 7 Related literature

Since the stable set problem is NP-hard, no complete (or exact) algorithm is known that solves the stable set problem in polynomial time. Many other techniques have been proposed, including approximation algorithms, heuristics, or branch and bound structured methods. A survey of different formulations, complete methods and heuristics for the maximum clique problem222The maximum weighted stable set problem of a graph is equivalent to the maximum weighted clique problem of its complement graph. is given by Pardalos and Xue  and, more recently, by Bomze et al. .

Although semidefinite programs can be solved in polynomial time theoretically, it lasted until a few years ago until fast solvers for this purpose were implemented. Until then, application inside a branch and bound framework was unrealistic. Still, solving a semidefinite program takes relatively much time, compared to solving a linear program. However, since semidefinite programming solvers are getting faster, semidefinite relaxations become a serious candidate to be used within a branch and bound framework, see for instance the paper by Karisch et al. .

A large number of references to papers concerning semidefinite programming are on the web pages of Helmberg  and Alizadeh . A general introduction on semidefinite programming applied to combinatorial optimization is given by Goemans and Rendl .

Another area that made semidefinite programming useful in practice is that of approximation algorithms. In this field one tries to give a performance guarantee for an algorithm on a particular problem. In particular, the paper  by Goemans and Williamson uses a semidefinite relaxation and randomized rounding to prove such a performance guarantee for the maximum cut problem of a graph and satisfiability problems. Following this result numerous papers appeared, also concerning the approximation of satisfiability problems, including  and .

The solution structure of the current work, namely problem decomposition by branching on promising subdomains, is similar to the method described in . In that paper a linear relaxation is used to identify promising values. Moreover, by exploiting the discrepancy structure of the method combined with reduced costs, suboptimality of subproblems can be proved very fast.

Another hybrid approach using linear programming and constraint programming has been investigated by Ajili et al.  and El Sakkout et al. . A subset of constraints is relaxed as a linear program in such a way that its solution is always integral. The solution to the relaxation serves as a suggestion (a ‘probe’) for solving the complete program using a constraint programming solver. A probe is used to detect infeasibility, to remove inconsistent domain values and to guide the search. During search, many probing steps are being made. This results in a tight cooperation of the linear programming and constraint programming solver.

## 8 Conclusions and further research

We introduced a method that combines semidefinite programming and constraint programming to solve the stable set problem. Our experiments show that constraint programming can indeed benefit greatly from semidefinite programming. On instances that were very difficult to handle for a constraint programming solver, our hybrid method obtained very good results.

The discrepancy cuts we proposed to strengthen the semidefinite relaxation could not be applied efficiently to the instances we considered. However, for larger instances they could be helpful.

Further research in this direction would for instance be to obtain a filtering mechanism similar to the cost-based domain filtering for linear relaxations . In , Helmberg describes such a procedure, called variable fixing, for semidefinite relaxations. It would be interesting to see how his method can be applied in a constraint programming framework.

Also, one could consider a different way of selecting promising values from the solution of the semidefinite relaxation. A strategy that incorporates randomized rounding possibly yields better results. This thought is motivated by the use of randomized rounding of semidefinite relaxations in approximation algorithms, as discussed in Section 7.

Finally, this work has much in common with our previous work . The underlying general principle of decomposing a problem into promising subproblems according to a certain heuristic is currently under research.

## Acknowledgements

Many thanks to Michela Milano for fruitful discussion and helpful comments on this paper.


text:
# Development of an approximation scheme for quasi-exactly solvable double-well potentials

## I Introduction

Quasi-exactly solvable (QES) systems are intermediate to exactly solvable and non-solvable ones. These quantum systems are characterized by the fact that, only a finite number of the eigenstates can be analytically determined . One dimensional examples include anharmonic potentials, some of which have a double-well structure. Interestingly, a QES potential with a centrifugal barrier appears in the treatment of ring lasers, where the Fokker-Planck equation is converted to an equivalent Schr$`\ddot{o}`$dinger eigenvalue problem . Other QES examples have been encountered in similar problems . These type of Hamiltonians, amenable to partial algebraization, have also manifested in various other areas of physics  and have attracted considerable attention in the recent literature . A number of independent methods, from group theoretical to algebraic ones, have been developed for studying the QES systems .

In the conventional approach to quasi-exactly solvable systems, one starts with an operator suitably constructed from the differential realizations of the generators of a given Lie algebra, acting in the finite dimensional space of monomials. A QES Hamiltonian is then arrived at by converting the above problem into a Schr$`\ddot{o}`$dinger eigenvalue equation, through appropriate similarity transformations. Although much effort has gone into studying the structure of the analytically available part of the spectrum, there has not been any attempt in the literature, to the best of the authors’ knowledge, to find the ones, not determined analytically. The fact that, some of the QES potentials exhibit a double-well structure, for which determining the approximate eigenvalues and eigenstates have been quite challenging, makes this problem worth investigating. Further, a number of these states are constrained to lie between certain range of energies; hence, apart from academic interest, developing an approximation scheme to locate these eigenvalues and eigenfunctions, provides an ideal ground for testing the proposed scheme. As mentioned earlier, a number of these type of potentials appear in the Schr$`\ddot{o}`$dinger equation originating from the Fokker-Planck equation, governing the dynamics of non-equilibrium systems. Finding the low-lying states of these problems accurately have physical significance; for example the tunneling rate is related to the energy difference between the ground and first excited state .

The goal of this paper is to develop an accurate approximation scheme for finding the non-exactly determined eigenstates of the QES Hamiltonians. For this purpose, we employ a recently developed method for solving linear differential equations , used earlier for diagonalizing many-body interacting systems . The analytically obtainable part of the spectra of a number of QES systems is first determined, in the following section, to demonstrate the working of the above method. We then proceed, in Sec.III, to compute approximately, those eigenvalues and eigenstates of a prototype QES system, which are not amenable to analytical treatment. The fact that our procedure yields a series solution for a differential equation in terms of energy, allows us to treat the same as a variational parameter. This can be effectively used for the identification of the eigenstates and corresponding eigenvalues, as will be shown in the text. It is particularly useful for the quasi-exactly solvable systems, where the ground state is known and a number of eigenstates are both bounded below and above.

The example of a double-well potential is taken deliberately, since the conventional techniques are not very reliable for the same. The usefulness of the present approach is shown by computing a number of low-lying eigenstates. Our approach yields the approximate eigenfunctions and eigenvalues, whose accuracy can be improved to any desired level, in a controlled manner. We then compare our results, with another convergent numerical scheme ; it is found that, the first few terms in our approximate solutions are enough to yield the excited state eigenvalues, accurate upto the third place of the decimal. We conclude in Sec.IV, after pointing out the advantages and limitations of the present approximation scheme and directions for further investigations.

## II Exact eigenstates of quasi-exactly solvabe systems: a novel approach

In this section, we obtain the analytically solvable part of the eigenspectra, of a class of QES systems, making use of a recently developed method for solving linear differential equations . In this procedure, the solution space of the differential equation is connected with the space of monomials. We concentrate on the QES systems having polynomial potentials, with or without a centrifugal barrier term, although the method can be applied to other systems as well. As will become clear in the subsequent section, the same procedure yields the approximate eigenfunctions and eigenvalues, for the analytically inaccessible states of the QES systems, to the desired accuracy.

A single variable differential equation, after suitable manipulations (which will become clear from the examples in the text), can be written as,

$$\left[F(D)+P(x,d/dx)\right]y(x)=0,$$
(1)
where, $`Dx\frac{d}{dx}`$ is the Euler operator, $`F(D)_{n=\mathrm{}}^{n=\mathrm{}}a_nD^n`$ and $`a_n`$’s are some parameters; $`P(x,d/dx)`$ can be an arbitrary polynomial function of $`x`$, $`\frac{d}{dx}`$ and other operators. It can be straightforwardly shown by direct substitution that, the solution to Eq. (1) can be written in the form,

$`y(x)`$ $`=`$ $`C_\lambda \left\{{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(1)^m\left[{\displaystyle \frac{1}{F(D)}}P(x,d/dx)\right]^m\right\}x^\lambda `$ (2)
provided, $`F(D)x^\lambda =0`$ and the coefficient of $`x^\lambda `$ in $`y(x)C_\lambda x^\lambda `$ is zero (no summation over $`\lambda `$); here, $`C_\lambda `$ is a constant. This straightforward method, not only yields solutions to the familiar differential equations , but also leads to the diagonalization of a number of correlated many-body Hamiltonians .

For the QES case, we first consider the example of the sextic oscillator, whose Hamiltonian (in the units $`\mathrm{}=2m=1`$) is given by

$$H=\frac{d^2}{dx^2}+\alpha x^2+\beta x^4+\gamma x^6.$$
(3)
It is well-known that, this problem is QES, provided a certain relation exists between the parameters $`\alpha ,\beta `$ and $`\gamma `$. Instead of postulating the same, we first illustrate how this condition emerges naturally. Asymptotic analysis indicates a measure of the form $`\widehat{\psi _0}e^{(ax^2+bx^4)}`$, with the unknown parameters, $`a`$ and $`b`$ to be determined from the Hamiltonian parameters $`\alpha `$, $`\beta `$ and $`\gamma `$. A similarity transformation $`\stackrel{~}{H}=\widehat{\psi _0}^1H\widehat{\psi _0}`$, yields;

$`\stackrel{~}{H}`$ $`=`$ $`{\displaystyle \frac{d^2}{dx^2}}+8bx^3{\displaystyle \frac{d}{dx}}+4ax{\displaystyle \frac{d}{dx}}+(\alpha 4a^2+12b)x^2`$ (5)
$`+(\beta 16ab)x^4+(\gamma 16b^2)x^6+2a.`$
Setting the coefficients of $`x^4`$ and $`x^6`$ equal to zero, one obtains,

$$a=\frac{\beta }{4\sqrt{\gamma }}\mathrm{a}ndb=\frac{\sqrt{\gamma }}{4}.$$
(6)
Assuming that the solution of the eigenvalue equation $`\stackrel{~}{H}P_n(x)=EP_n(x)`$ is a polynomial, of degree $`n`$, one observes that, the operators, $`2\sqrt{\gamma }x^3\frac{d}{dx}`$ and $`(\alpha \frac{\beta ^2}{4\gamma }+3\sqrt{\gamma })x^2`$ increase the degree of $`P_n(x)`$ by two. Preserving the degree of the polynomial leads to the above mentioned relationship between the coupling parameters of the sextic oscillator  :

$$\frac{1}{\sqrt{\gamma }}\left(\frac{\beta ^2}{4\gamma }\alpha \right)=2n+3.$$
(7)
It is worth mentioning that, in exactly solvable problems, the differential operators do not increase the degree of the polynomial and hence no additional condition is required there.

The eigenvalue problem can be cast in a form as given in Eq. (1) :

$`\left[D(D1)2\sqrt{\gamma }x^5{\displaystyle \frac{d}{dx}}{\displaystyle \frac{\beta }{\sqrt{\gamma }}}x^3{\displaystyle \frac{d}{dx}}+2n\sqrt{\gamma }x^4+\left(E{\displaystyle \frac{\beta }{2\sqrt{\gamma }}}\right)x^2\right]P_n(x)=0.`$ (8)
The condition $`F(D)x^\lambda =0`$, yields $`\lambda =0`$ and $`1`$. It is easy to see that the two values of $`\lambda `$ separate the Hilbert space into even and odd sectors. The solution, corresponding to the root $`\lambda =0`$, can be expanded as,

$`P_n^0(x)`$ $`=`$ $`C_0{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(1)^m\left[{\displaystyle \frac{1}{D(D1)}}\left(\stackrel{~}{E}x^2{\displaystyle \frac{\beta }{\sqrt{\gamma }}}x^3{\displaystyle \frac{d}{dx}}+2n\sqrt{\gamma }x^42\sqrt{\gamma }x^5{\displaystyle \frac{d}{dx}}\right)\right]^mx^0`$ (9)
$`=`$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}Q_k(\stackrel{~}{E}){\displaystyle \frac{x^{2k}}{2k!}},`$ (10)
where $`\stackrel{~}{E}=E\beta /2\sqrt{\gamma }`$ and $`Q_k(\stackrel{~}{E})`$ is an appropriate polynomial in energy.

Since the degree of the polynomial has been already fixed at $`n`$, an even integer in this case, we impose the condition that the series terminates at the desired point by putting the coefficient of the subsequent term to zero. As will be explicitly seen, the same coefficient appears as a factor in rest of the terms of the series, thereby ensuring that the degree of the polynomial is maintained. Writing $`n=4j`$, where $`j`$ can take semi-integer values, one needs to put $`Q_{2j+1}(\stackrel{~}{E})=0`$ to obtain the above result. This condition leads to $`2j+1`$ independent solutions for energy and the corresponding eigenfunctions. We would like to emphasize that, so far we have not assumed any additional property of the spectral problem, for arriving at these results.

In similar fashion, one can easily incorporate the centrifugal barrier in to QES problems. The Hamiltonian of the sextic oscillator for $`\beta =0`$, with a centrifugal barrier is given by,

$`\widehat{H}={\displaystyle \frac{d^2}{dx^2}}+{\displaystyle \frac{\sigma }{x^2}}+\alpha x^2+\gamma x^6;`$ (11)
the corresponding measure is of the form $`\psi _0=x^{2l}\mathrm{exp}(ax^4/4)`$. Performing an appropriate similarity transformation and setting the coefficients of $`x^6`$ and $`1/x^2`$ to zero, we obtain,

$`\left[{\displaystyle \frac{d^2}{dx^2}}+{\displaystyle \frac{4l}{x}}{\displaystyle \frac{d}{dx}}+2\sqrt{\gamma }x^3{\displaystyle \frac{d}{dx}}+\left(\alpha +3\sqrt{\gamma }+4l\sqrt{\gamma }\right)x^2\right]P_n(x)=EP_n(x),`$ (12)
with $`a=\sqrt{\gamma }`$ and $`l=\frac{1}{4}+\frac{1}{2}\sqrt{\frac{1}{4}+\sigma }`$.

The condition that $`P_n(x)`$ is a polynomial of degree $`n`$, leads to,

$$\frac{\alpha }{4\sqrt{\gamma }}+\frac{1}{2}\sqrt{\frac{1}{4}+\sigma }=n/2+1.$$
(13)
Multiplication of Eq. (12) with $`x^2`$ yields,

$`\left[D(D+4l1)+Ex^2+2n\sqrt{\gamma }x^42\sqrt{\gamma }x^5{\displaystyle \frac{d}{dx}}\right]P_n(x)=0,`$ (14)
which $`\lambda =0`$ and $`14l`$, as the solutions of $`F(D)x^\lambda =0`$. A suitable parameterization in the form of

$`\alpha =4a\left(s+{\displaystyle \frac{1}{2}}+\mu \right)\mathrm{and}\sigma =4\left(s{\displaystyle \frac{1}{4}}\right)\left(s{\displaystyle \frac{3}{4}}\right)`$ (15)
makes the result amenable for comparison with the existing literature . Analogous to the earlier example, one can now find solutions corresponding to various values of $`n`$. For $`n=2`$ one finds,

$`E_{2\pm }=\pm \sqrt{32as}`$ (16)
$`\mathrm{and}P_{2\pm }(x)=\left[ax^2{\displaystyle \frac{E_{2\pm }}{4}}\right].`$ (17)
Even states for other values of $`n`$ can also be similarly obtained.

In the absence of a centrifugal barrier, for the anharmonic oscillator of Eq.(8), under similar condition ($`n=2`$), the energy eigenvalues can be derived from ,

$`Q_2(\stackrel{~}{E})=\stackrel{~}{E}^2{\displaystyle \frac{2\beta \stackrel{~}{E}}{\sqrt{\gamma }}}4n\sqrt{\gamma }=0.`$ (18)
The corresponding two solutions,

$`E_\pm ={\displaystyle \frac{3\beta }{2\sqrt{\gamma }}}\pm \left({\displaystyle \frac{\beta ^2}{\gamma }}+8\sqrt{\gamma }\right)^{1/2},`$ (19)
yield,

$`\psi _{2_\pm }(x)=\mathrm{exp}\left({\displaystyle \frac{\beta x^2}{4\sqrt{\gamma }}}{\displaystyle \frac{\sqrt{\gamma }x^4}{4}}\right)P_{2_\pm }(x);`$ (20)
$`\mathrm{h}ereP_{2_\pm }(x)=1+\left[{\displaystyle \frac{\beta }{\sqrt{\gamma }}}\pm \left({\displaystyle \frac{\beta ^2}{\gamma }}+8\sqrt{\gamma }\right)^{1/2}\right]x^2.`$ (21)
The case of $`n=4`$, with $`\beta =0`$ and for negative values of $`\alpha `$, is a double-well potential. For specificity, we take $`\gamma =1`$, which yields $`\alpha =11`$ and hence

$`\widehat{H}={\displaystyle \frac{d^2}{dx^2}}11x^2+x^6.`$ (22)
The energy eigenvalues for this Hamiltonian are obtained from, $`Q_3(E)=E(E^264)=0`$, which gives $`E=8`$, $`0`$ and $`+8`$. The respective polynomial parts of the wave functions are,

$`P_4_{}(x)`$ $`=`$ $`1+4x^2+2x^4,`$ (23)
$`P_{4_0}(x)`$ $`=`$ $`1(2/3)x^4,`$ (24)
$`\mathrm{a}ndP_{4_+}(x)`$ $`=`$ $`14x^2+2x^4.`$ (25)
It is clear that, the ground state has no nodes on the real line and the subsequent two states have, respectively two and four nodes on the real line, as desired. The above procedure straightforwardly extends to higher values of $`n`$. We note that, polynomial potentials, with the highest degree $`4m+2`$, with $`m`$ integer can also be solved in an analogous manner, provided appropriate conditions are imposed on the couplings. It should also be pointed out that, the solutions obtained so far, correspond to the root $`\lambda =0`$. The other root $`\lambda =1`$, does not lead to a closed form expression for the solution. Hence, the odd states are not amenable to an analytical treatment. This is a manifestation of the QES nature of these quantal problems.

## III Approximation scheme for the analytically inaccessible states

As is clear from the previous example, it is not possible to analytically determine the states, other than the ones allowed by the condition given in Eq.(7). In this section, we present an approximation scheme for finding these analytically inaccessible states, taking the first excited state of the above mentioned double-well potential, as the example. Assuming that, the measure is same for all the states, the residual part, which will be suitably approximated to a finite degree polynomial in the following, should have one zero on the real line. The rest of the zeros, if present, should lie on the complex plane. In developing an approximation scheme, one also needs to take care of the convergence of the series. The above two criteria are utilized for finding out the approximate eigenvalue and eigenfunction of the first excited state, which can then be extended to other states.

Expanding Eq. (8) for $`\lambda =1`$, one obtains,

$`P_4^1(x)=xE{\displaystyle \frac{x^3}{3!}}+(E^236){\displaystyle \frac{x^5}{5!}}+(76EE^3){\displaystyle \frac{x^7}{7!}}+(E^4+8E^23024){\displaystyle \frac{x^9}{9!}}+\mathrm{}`$ (26)
This series does not terminate at any finite order. In developing an approximation scheme, one possible way is to equate the coefficient of an appropriate term in the series to zero, in order to approximate it as a polynomial. This determines the energy eigenvalue and eigenfunction; one also has to show that higher order terms attain sufficiently small values to ensure convergence.

The other possible way is a variational approach. The fact that our procedure of solving the differential equation yields a series solution in terms of energy, allows us to treat the same as a variational parameter. This can be effectively used for the identification of the eigenstates and corresponding eigenvalues, as will be shown below. Since the ground state is known and we want to approximate the excited state as a product of the ground state measure and a finite degree polynomial, we are left to identify only the polynomial part of the excited state. This can be carried out in a convergent manner to the desired degree of accuracy. We first terminate the series given by Eq.(26) at some finite order. Now we check that this approximate polynomial solution $`u(x)`$ has one zero on the the real line, this gives us a range of energy values satisfying the above mentioned condition. Since, we know that $`u(x)`$ is an approximate solution, substitution of the same in the eigenvalue equation $`\stackrel{~}{H}u(x)=Eu(x)`$ will leave a residual term. To get the best approximate value for the energy at this order, we minimize $`\mathrm{\Delta }`$ defined as,

$`\mathrm{\Delta }=|{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{\frac{x^4}{2}}u^{}(x)(\stackrel{~}{H}E)u(x)𝑑x|,`$ (27)
using $`E`$, as a variational parameter. For a better appreciation of the above points, we give below a plot of $`\mathrm{\Delta }`$ versus $`E`$ for a given approximate solution. Minimum of the plot indicates the best approximate excited state for a given degree of the polynomial. The approximate energy eigenvalue obtained by a higher degree polynomial gives better result, as can be easily seen by comparison with the numerically calculated value, given in Table I. This can also be noticed from Fig.1, which clearly shows that, the minimum of $`\mathrm{\Delta }`$ of the higher degree polynomial occurs for a value of $`E`$ closer to the correct result. It should also be noticed that the value of the residual term $`\mathrm{\Delta }`$ is much smaller for the higher order polynomial. Incorporation of higher order terms makes the convergence better in the present scheme. It is found that, the first few terms in our approximate solutions are enough to yield the excited state eigenvalues, accurate upto the third place of the decimal, as shown in Table I.

Table I. Comparison of approximate eigenvalues $`E_1`$ and $`E_2`$ for the solutions of degrees five and nine respectively, with corresponding numerically calculated values $`E_{Num}`$

| Energy State | $`E_{Num}`$ | $`E_1`$ | Deviation | $`E_2`$ | Deviation |
| --- | --- | --- | --- | --- | --- |
| $`1`$ | $`7.917350`$ | $`7.913704`$ | $`0.04\%`$ | $`7.916400`$ | $`0.011\%`$ |
| $`3`$ | $`2.520359`$ | $`2.419229`$ | $`4.01\%`$ | $`2.549348`$ | $`1.15\%`$ |
| $`5^{}`$ | $`14.112964`$ | $`13.403590`$ | $`5.02\%`$ | $``$ | $``$ |

This eigenvalue corresponds to an approximate solution of degree thirteen.

FIG. 1. Plots depicting the minima of $`\mathrm{\Delta }`$ near the numerically evaluated eigenvalue, for the first excited state of the double-well potential. The solid and dashed curves show the variations of $`\mathrm{\Delta }`$, when the approximate polynomial parts of the wave function are of degrees nine and five, respectively.

Apart from finding the states in between the known ones, as has been done above, one can also approximate the even and odd states lying above the analytically determinable part of the spectrum. For this purpose, we plot $`\mathrm{\Delta }`$ for odd and even states in a range of energy values, in Figs. 2 and 3 respectively. This is computed taking approximate solution of degree nine for the odd states; and twelve for the even ones.

FIG. 2. Variations of $`\mathrm{\Delta }`$ with respect to $`E`$, in the odd-parity sector. Here the polynomial part of the wave function has degree nine and the energy range is much bigger, as compared to Fig.1. One clearly notices that several minimum values of $`\mathrm{\Delta }`$ occur closer to the numerically determined eigenvalues.

FIG. 3. Plot of $`\mathrm{\Delta }`$ versus $`E`$, for the even states, when the approximate polynomial solution is of degree twelve. Like the odd sector in Fig.2, here also some of the minima are located close to the corresponding numerical values.

It is extremely interesting to observe that, various minima of the plot show the expected locations of the states in the energy space, with reasonable accuracy, as seen from Table II., containing the numerically calculated eigenvalues of higher excited states. It is worth pointing out that, all the minima of the plot do not correspond to physically acceptable states. One needs to take those values, which fulfill the required conditions on the wave functions. At this stage of approximation, the computed eigenvalues for the higher excited states agree with the numerically obtained ones upto $`5\%`$ accuracy. As has been done earlier, in order to improve the result further, one needs to take still higher degree polynomials.

Below, we give another table containing numerically calculated energy eigenvalues of higher excited states of the double-well potential. It is clear from the spectrum, that higher excited states are approximately equispaced as expected.

Table II. Numerically calculated energy eigenvalues of a number of excited states of the double-well potential

| Energy State | $`E_{Num}`$ | Energy State | $`E_{Num}`$ |
| --- | --- | --- | --- |
| $`6`$ | $`21.1575028`$ | $`13`$ | $`89.0540515`$ |
| $`7`$ | $`28.9747742`$ | $`14`$ | $`100.9447956`$ |
| $`8`$ | $`37.4938424`$ | $`15`$ | $`113.3093212`$ |
| $`9`$ | $`46.6606249`$ | $`16`$ | $`126.1303320`$ |
| $`10`$ | $`56.4324169`$ | $`17`$ | $`139.3922937`$ |
| $`11`$ | $`66.7742029`$ | $`18`$ | $`153.0811528`$ |
| $`12`$ | $`77.6565110`$ | $`19`$ | $`167.1841151`$ |

## IV Conclusions

In conclusion, we have applied a recently developed scheme for solving linear differential equations, to find the exact eigenstates of a wide class of quasi-exactly solvable Hamiltonians and also to develop a convergent approximation scheme to determine the states not obtainable by analytical treatments. The fact that the method used provided an expansion of the wave functions in terms of the corresponding energies, allowed us to treat energy as a variational parameter, in the approximation scheme. An independent numerical scheme was used to check that the low-lying eigenvalues are accurate up to third place of decimal, with only a few terms from the series expansion. The higher excited states need more number of terms, since these have $`5\%`$ accuracy, under similar conditions. It also needs to be pointed out that, the QES systems treated here are equipped with a ground state measure, which facilitates the working of the approximation scheme. It will be interesting to extend the present scheme to other anharmonic potentials.

Acknowledgement : We acknowledge many useful discussions with Prof. V. B. Sheorey, who has also made available to us the algorithm for the numerical solution of the QES eigenvalue problems.