format: text

text:
# Comments on Adiabatic Quantum Algorithms

## 1. Introduction

Recently there has been considerable interest in a proposed scheme for adiabatic quantum computation  and speculation that it may even provide a mechanism for efficient solution of hard problems. Both the validity of the adiabatic theorem and the arguments  for its efficiency depend on the existence of an eigenvalue gap. However, the existence of such a gap has not been proven. It has been conjectured  on the basis of numerical simulations and the so-called non-crossing rule.

The adiabatic quantum algorithm is designed to take the ground state of an initial Hamiltonian $`H_0`$ to that of a final Hamiltonian $`H_1`$ using a linear interpolating Hamiltonian of the form $`H(s)=(1s)H_0+sH_1`$. The quantum adiabatic theorem  can be used to show that the efficiency of the adiabatic approximation is $`O(1/g_{\mathrm{min}}^2)`$ where $`g_{\mathrm{min}}`$ denotes the minimum energy gap $`g(s)=E_1(s)E_0(s)`$ between the ground and first excited states of $`H(s)`$. Thus the efficiency of adiabatic quantum computation depends on how rapidly the eigenvalue gap decreases as the size of the problem increases. The purpose of this note is not to resolve this question, but to discuss some of the issues that arise.

Farhi, et al  (hereafter referred to as FGG) typically use an initial Hamiltonian of the form

(1.1) $`H_0=\frac{1}{2}{\displaystyle \underset{j}{}}a_j[\sigma _x(j)+I]\frac{1}{2}{\displaystyle \underset{j}{}}a_jII\mathrm{}[\sigma _x+I]\mathrm{}I`$
with $`a_j`$ a non-negative integer. The eigenstates of $`H_0`$ are products of eigenstates of $`\sigma _x`$. The so-called “computational basis” $`|k_1\mathrm{}k_n`$ consists of products of eigenstates of the Pauli matrix $`\sigma _z`$. (However, it is customary to identify each $`k_i`$ with an element of $`Z_2=\{0,1\}`$ rather than with the usual eigenvalues of $`\pm 1`$ or, equivalently, replacing $`\sigma _z`$ by $`\frac{1}{2}[\sigma _z+I]`$ as has been done above.) In the computational basis, the eigenstates of $`H_0`$ have the form $`2^{n/2}_{k_1\mathrm{}k_n}\pm 1|k_1\mathrm{}k_n`$ (with all signs $`+1`$ for the ground state). FGG define a non-negative final Hamiltonian $`H_1`$ which is diagonal in the computational basis and has the solution of some problem as its ground state, e.g., they encode the problem known as “3-SAT” in the computational basis. We note here only that this encoding results in a Hamiltonian of the form

(1.2) $`H_1={\displaystyle \underset{k_1\mathrm{}k_n}{}}E_{k_1\mathrm{}k_n}|k_1\mathrm{}k_nk_1\mathrm{}k_n|`$
where the ground state energy is zero and the other $`E_{k_1\mathrm{}k_n}`$ are positive integers with an upper bound that is $`O(n^3)`$.

First, we point out that the non-crossing rule is completely inadequate for the purpose of establishing a gap. There are realistic models of physical systems which do exhibit crossings, despite the absence of a symmetry common to both the initial and final Hamiltonians $`H_0`$ and $`H_1`$. The most well-known such system is $`H_2^+`$, the hydrogen molecule ion . Another example is the Hubbard model for benzene . The non-crossing rule and the limitations on its applicability to adiabatic quantum computation are discussed in Appendix B.

However, the most one could hope to gain from the non-crossing rule is uniqueness of the ground state. Fortunately, a standard argument based on the Perron-Frobenius theorem  suffices for that purpose. This argument is sketched in the Section 2, but does not address the more fundamental issue of the size of the resulting eigenvalue gap.

Some insight into the issues raised by adiabatic computation can be obtained by considering various strategies for replacing FGG’s final Hamiltonian $`H_1`$ by a modification whose ground state is the solution of Grover’s search problem. We argue that if their algorithm is sufficiently robust to solve an NP problem in polynomial time $`O(N^p)`$, then a modification should be able to perform a successful search of an unordered list of $`M`$ items in $`O([\mathrm{log}M]^p)`$ time, violating the conventional wisdom that Grover’s algorithm (which requires $`O(\sqrt{M})`$ time) is optimal. A model problem suggested by FGG shows that this simplistic expectation need not hold; however, their example also suggests that gaps are associated with the presence of a symmetry, not its absence.

The issues raised here involve questions in several subfields fields of physics and computer science. In order to make this note accessible to people with diverse backgrounds ranging from Schrödinger operator theory to computer science, two appendices are included — one on Grover’s algorithm and one on the non-crossing rule.

## 2. Uniqueness of the Ground State

Although proof of an eigenvalue gap is likely to be difficult, proving that systems of the type considered by FGG have a unique ground state is easier. It relies on a standard argument widely used  to prove uniqueness (and positivity) of the ground state in a variety of systems, including quantum lattice models. There is no particular originality in the argument given below. We present it only in the hope of clarifying some issues, particularly the distinction between the uniqueness of the ground state and the existence of a lower bound on the size of the resulting eigenvalue gap.

The ground state of the initial Hamiltonian (1.1) is easily seen to be unique and consists of products of ground states of $`\sigma _x(k)`$. However, we wish to transform $`H_0`$ to the computational basis of products $`|k_1k_2\mathrm{}k_n`$ of eigenstates of $`\sigma _z(k)`$. This is easily achieved using tensor products of the Hadamard transform. Instead of examining $`H_0`$ itself, we consider the operator $`F=e^{H_0}`$. The ground state of $`H_0`$ is the eigenfunction corresponding to the largest eigenvalue of the matrix

(2.1) $`F=e^{H_0}={\displaystyle \underset{k=1}{\overset{n}{}}}e^{a_k[\sigma _x(k)+I]/2}=B_1B_2\mathrm{}B_n`$
where

(2.8) $`B_k`$ $`=`$ $`\frac{1}{2}e^{a_k/2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\left(\begin{array}{cc}e^{a_k/2}& 0\\ 0& e^{a_k/2}\end{array}\right)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)`$
(2.11) $`=`$ $`e^{a_k/2}\left(\begin{array}{cc}\mathrm{cosh}\frac{a_k}{2}& \mathrm{sinh}\frac{a_k}{2}\\ \mathrm{sinh}\frac{a_k}{2}& \mathrm{cosh}\frac{a_k}{2}\end{array}\right)`$
Now, since all elements of each $`B_k`$ are strictly positive, their tensor product $`F=e^{H_0}`$ also has strictly positive elements. Hence, by the Perron-Frobenius theorem  the largest eigenvalue of $`F`$, and the ground state of $`H_0`$, is unique.

We would like to know that the ground state remains unique for a Hamiltonian of the form $`H=(1s)H_0+sH_1`$. For this, it suffices that $`H_1`$ is diagonal (in the computational basis) and has a unique ground state. By the Lie-Trotter formula

(2.12) $`e^H=\underset{m\mathrm{}}{lim}\left(e^{\frac{s}{m}H_1}e^{\frac{1s}{m}H_0}\right)^m`$
The effect of the diagonal matrix $`e^{\frac{s}{m}H_1}`$ is simply to multiply each row of $`e^{\frac{1s}{m}H_0}`$ by a positive number of the form $`e^{\lambda _jt/m}`$. Hence the product $`e^{\frac{s}{m}H_1}e^{\frac{1s}{m}H_0}`$ also has positive elements and so does its $`m`$-th power. (Moreover, because $`e^{\lambda _js/m}1`$ as $`m\mathrm{}`$, these positive elements do not become zero in the limit.) Thus, one can again apply the Perron-Frobenius theorem to conclude that the largest eigenvalue of $`e^H`$ is unique and, hence, $`H`$ has a unique ground state if $`0s<1`$. (The argument breaks down at $`s=1`$ because $`e^{\frac{1s}{n}H_0}=I`$ no longer has strictly positive elements off the diagonal and is completely decomposable.)

Note that, even if $`H_1`$ has a degenerate ground state, the interpolating Hamiltonian $`H=(1s)H_0+tH_1`$ will still have a unique ground state for all $`0s<1`$; however, the difference between the two lowest eigenvalues must $`0`$ as $`s1`$. Thus, uniqueness of the ground state is a very different matter from an eigenvalue gap of minimal size. Indeed, the uniqueness argument above holds for models , such as an adiabatic search for which the gap can be shown to decrease exponentially.

In addition to uniqueness, the Perron-Frobenius theorem implies that the ground state has the form $`_kc_k|k_1k_2\mathrm{}k_n`$ with strictly positive coefficients $`c_k>0`$ in the computational basis for all $`s[0,1)`$. In the limit as $`s1`$ all but one of these $`c_k0`$.

## 3. Adiabatic Search Algorithms

As described in Appendix A, Grover’s  algorithm is designed to efficiently locate an unknown but identifiable target via the use of a unitary operator $`G`$ which can be written as $`G=e^{i\pi A}`$ where

(3.3) $`A|k_1\mathrm{}k_n=\{\begin{array}{cc}0& \text{if}(k_1\mathrm{}k_n)=(t_1\mathrm{}t_n)\hfill \\ |k_1\mathrm{}k_n& \text{otherwise}\hfill \end{array}`$
One can easily implement an adiabatic search by choosing for the final Hamiltonian $`H_1`$ in the FGG algorithm the Grover generator, $`A`$, above. The adiabatic evolution will take the ground state of $`H_0`$, namely, $`|\psi _02^{n/2}_{k_1\mathrm{}k_n}|k_1\mathrm{}k_n`$ to the ground state of $`H_1=A`$, namely $`|t_1\mathrm{}t_n`$. However, FGG have shown that this process takes exponential time. The analysis can be simplified  by modifying the initial Hamiltonian $`H_0`$ to reduce the analysis to a two-dimensional problem.

The reduction to a two-dimensional problem, which plays a critical role in Grover’s algorithm, is associated with a $`(2^n1)`$-fold degeneracy in the adiabatic Hamiltonian $`H_0`$ and $`H_1`$. However, this is not at all essential for the success of an adiabatic search. All that is needed is that the ground state of $`H_1`$ be the target state $`|t_1\mathrm{}t_n`$. This suggests that one try to modify $`H_1`$ so that its ground state is $`|t_1\mathrm{}t_n`$, but the eigenvalue distribution of its excited states is similar to that of a final Hamiltonian known to have a gap.

Suppose that a problem is known to have an efficient solution encoded in the final Hamiltonian $`H_1`$. Then setting $`H_2=GH_1`$, yields a Hamiltonian identical to $`H_1`$ except that the eigenvalue associated with the target state $`|t_1\mathrm{}t_n`$ is multiplied by $`1`$. (Because $`G`$ and $`H_1`$ are both diagonal in the computational basis, they commute and $`H_2`$ is self-adjoint.) Because $`H_1`$ was defined to be non-negative, $`H_2`$ has exactly one negative eigenvalue so that its ground state is now the target state $`|t_1\mathrm{}t_n`$. Applying the adiabatic algorithm to the modified interpolating Hamiltonian $`H_2(s)=(1s)H_0+sH_2`$ should take the ground state of $`H_0`$ to the target state. Moreover, the only effect on the final Hamiltonian is to move one excited state below the previous ground state, without decreasing the final energy gap $`g(1)`$. A similar Hamiltonian which is non-negative could be constructed the replacement

$`E_{k_1\mathrm{}k_n}`$ $`0`$ $`\mathrm{if}k_1\mathrm{}k_n=t_1\mathrm{}t_n`$
$`E_{k_1\mathrm{}k_n}`$ $`E_{k_1\mathrm{}k_n}+1`$ $`\mathrm{otherwise}.`$
in (1.2). In either case, one would not generally expect these modifications to significantly affect $`g_{\mathrm{min}}(s)`$, in which case the adiabatic search would be as efficient as the solution of the problem encoded in $`H_1`$.

Thus, if Farhi, et al’s projection of an efficiency of $`O(n^p)`$ is correct for the problem encoded in $`H_1`$, then one would expect an efficiency of $`O([\mathrm{log}N]^p)`$ for the adiabatic search of a list of $`N=2^n`$ items described above. However, this would imply that a quantum computer could search an unordered list in time $`O([\mathrm{log}N]^2)`$, contradicting the conventional wisdom that a speed-up greater than $`O(\sqrt{N})`$ is not possible . This does not necessarily imply a contradiction. The proofs that $`O(\sqrt{N})`$ is optimal depend on assumptions about the nature of the “oracle query” used in the search. However, van Dam, Mosca and Vazirani  have observed that the encoding of solutions of other problems, such as 3-SAT, in $`H_1`$ implicitly assume the ability of the computer to perform more general queries. Thus, standard complexity query arguments can not rule out the possibility of polynomial time algorithms.

After seeing a preliminary version of this manuscript, FGG  pointed out that the change from $`H_1`$ to $`H_2`$ described above can have a dramatic effect on the gap. We describe their example in the next section.

It should also be noted that our expectation that the change from $`H_1`$ to $`H_2`$ will not decrease the minimum gap is not based on the presumption that one is a small perturbation of the other. On the contrary, (as FGG  emphasized) the two Hamiltonians differ by a multiple of a projection, which can have a significant effect on the spectrum. Indeed, it is essential to our strategy to effect such a change in the final Hamiltonian. However, unless this also induces a change in the structure of the problem, such as a symmetry-breaking, this need not affect the generic behavior of the spectra of the interpolating Hamiltonian; in particular, it need not lead to an avoided crossing of the two lowest levels.

## 4. A Separable Model

Let $`H_1=\frac{1}{2}_j[\sigma _z(j)+I]`$ and set $`a_j=1`$ in $`H_0`$. Then the interpolating Hamiltonian becomes

(4.1) $`H(s)=H_0+s(H_1H_0)=\frac{1}{2}{\displaystyle \underset{j}{}}\left[(1s)\sigma _x(j)+s\sigma _z(j)+I\right]`$
This system is separable, and exactly solvable. Since $`[(1s)\sigma _x(j)+s\sigma _z(j)+I]`$ has eigenvalues $`\frac{1}{2}(1\pm \sqrt{12s+2s^2})`$, the ground state of $`H(s)`$ has energy $`\frac{n}{2}(1\sqrt{12s+2s^2})`$ and the gap $`g(s)=\sqrt{12s+2s^2}`$ is independent of $`n`$ with $`g_{\mathrm{min}}=\frac{1}{\sqrt{2}}`$. The system also has a high level of symmetry, since $`H(s)`$ commutes with elements of the symmetric group $`S_n`$. In fact, there is an additional accidental degeneracy so that the $`(k+1)`$-st eigenvalue is $`k`$ with a degeneracy of $`\left(\begin{array}{c}n\\ k\end{array}\right)`$ for $`k=0,1,2\mathrm{}n`$.

Now if $`H_1`$ is replaced by $`H_2=GH_1`$, the symmetry is broken. However, there is still some symmetry and a high level of degeneracy. Essentially, only one state from each of the $`(n+1)`$ non-degenerate levels is affected and the problem reduces to an $`(n+1)`$-dimensional one which can be analyzed explicitly and shown to have an exponentially decreasing gap. In effect, an eigenvalue can be associated with the target state and must “cross” the levels lying below it to reach the bottom of the spectrum when $`s=1`$. (A more detailed examination shows that the levels become exponentially close and then bounce away with the target information transmitted to the lower level).

Although this shows that the argument sketched in Section 3 above cannot be made rigorous, this model is not generic. The exponentially decreasing gap is the result of a symmetry breaking which should not occur when the Hamiltonian $`H_1`$ has no symmetry to begin with.

## 5. Discussion

Farhi, et al  analyze several other models for which the gap behavior can be calculated explicitly and shown to decrease slowly (i.e, polynomial in $`n`$). However, as they point out, these models all have a high level or symmetry or structure which would lead to efficient classical algorithms. In some cases symmetry allows a high level of degeneracy which permits one to squeeze $`2^n`$ states into $`[0,O(n^p)]`$ without forcing an exponentially decreasing gap.

In adiabatic computation the typical choices for initial and final Hamiltonians have spectra with high degeneracies and consist of positive integers in a range that is polynomial in $`n`$. In general, the interpolating Hamiltonian $`H(s)`$ breaks these degeneracies and must squeeze $`2^n`$ distinct eigenvalues into a range of the form $`[0,O(n^p)]`$. Thus, most of them must be exponentially close.

It seem that a polynomial gap is more likely to be associated with the presence of symmetry, which allows high degeneracies, than with its absence.

For excited states, it is irrelevant whether or not the observed mergings are “avoided crossings” or true crossings. The simultaneous coalescence of a large number of excited states is essential to the algorithm. The states of the initial and final Hamiltonians are both product states, those of the initial Hamiltonian are products of eigenstates of $`\sigma _x`$, while those of the final Hamiltonian are products of eigenstates of $`\sigma _z`$ (and thus elements of the so-called “computational” basis). Therefore, the final ground state is always an evenly weighted superposition of all eigenstates of the initial Hamiltonian. Unlike standard applications of the adiabatic theorem, in which the main contribution comes from a few low-lying states, all of the excited states must contribute to the first-order correction. To do this, these higher excited states must get close in some sense. Fortunately, a quantum computer can make this first order correction efficiently, mixing in all $`2^n`$ excited states, and this is where the method gets its potential power. However, in order that low order perturbation theory suffice, it is essential that the gap between the ground and first excited state not decrease too rapidly as $`n`$ increases.

The non-crossing rule, which is discussed in Appendix B is based on the belief that “accidents” are extremely rare so that such phenomena as persistent degeneracy, or level crossings do not occur without some underlying physical phenomenon (such as a symmetry) with implications for the associated mathematical model. This viewpoint would suggest that if the lowest gap is to decrease only polynomially when the other eigenvalues are getting exponentially close, there must be some physical mechanism keeping them apart. We are skeptical that such a mechanism can be found for Hamiltonians which encode the solution of problems which do not have enough structure to yield efficient classical solutions.

This raises another question. Is the spectrum of the interpolating Hamiltonian sensitive to the association of particular eigenvalues with particular eigenstates in the final Hamiltonian, or is it primarily dependent on the eigenvalue distribution? For Hamiltonians with a good deal of structure, the first situation clearly holds, and FGG have observed  that the behavior of the eigenvalue gap can depend critically on the choice of initial Hamiltonian. However, one expects Hamiltonians which encode solutions to typical instances of classically intractable problems, to lack the structure needed for this sensitivity.

Because of the difficulty in analyzing the behavior of the gap in problems without structure, it may be worth considering a randomly defined Hamiltonian, i.e., let the energy in (1.2) have the form $`E_{k_1\mathrm{}k_n}=f(k_1\mathrm{}k_n)`$ where $`f`$ is a suitable random process. One can then ask if there is a sense in which the eigenvalue gap is “almost always” exponentially decreasing. If so, this would suggest that an exponentially small gap is generic, and likely to occur for NP-complete problems. There is already an extensive literature  on the spectral behavior of random Schrödinger operators, in which one obtains results about typical Hamiltonians with certain properties rather than one with a fixed potential. Although the model Hamiltonians used here have quite a different structure, and may require the development of new techniques, this approach seems worth considering.

Indeed, one could even define the final Hamiltonian for an adiabatic search by letting $`f`$ be a random variable taking integer values in $`[1,n]`$ and choosing

$`E_{k_1\mathrm{}k_n}=`$ $`0`$ $`\mathrm{if}k_1\mathrm{}k_n=t_1\mathrm{}t_n`$
$`E_{k_1\mathrm{}k_n}=`$ $`f(k_1\mathrm{}k_n)`$ $`\mathrm{otherwise}`$
in (1.2).

Acknowledgment: This note is based on a talk at the Q-Math8 conference in Mexico in December, 2001. It is a pleasure to thank the organizers for the opportunity to present these ideas to an audience of mathematical physicists in a delightful setting; to acknowledge helpful discussions with Professors Charles Bennett, George Hagedorn, Chris King, John Preskill, Barry Simon, and Umesh Vazirani; and to thank Professors Edward Farhi and Sam Gutman for suggesting the model in Section 4. Part of this work was done when the author was visiting the Institute for Theoretical Physics at the University of California, Santa Barbara and thereby also partly supported by the National Science Foundation under Grant PHY-9907949.

## Appendix A Grover’s Algorithm

Grover’s  algorithm is designed to efficiently locate an unknown but identifiable target state $`|t_1\mathrm{}t_n`$. This state may be the key denoting the location of an item in an unsorted list (e.g., the analogue of the name in an alphabetized phone book associated with telephone number one has been given) or a state with a certain verifiable property, e.g., a representation of the factors of a given number or the solution of some NP-complete problem.

Grover showed how to construct a unitary operator $`G`$ whose effect is simply to multiply the unknown target state $`|t_1\mathrm{}t_n`$ by $`1`$ and all others by $`+1`$. This operator can then be used to construct an algorithm which will find the target state with probability greater than $`\frac{1}{2}`$ in $`O(\sqrt{N})`$, i.e., $`O(2^{n/2})`$ time for $`N=2^n`$ states. The operator $`G`$ is a unitary operator whose action on the computational basis is simply

(A.3) $`G|k_1\mathrm{}k_n=\{\begin{array}{cc}\hfill |k_1\mathrm{}k_n& \text{if}(k_1\mathrm{}k_n)=(t_1\mathrm{}t_n)\hfill \\ \hfill |k_1\mathrm{}k_n& \text{otherwise}\hfill \end{array}`$
It is generally assumed that one has an “oracle” which can perform unitary operations to determine whether or not a state has the desired property and outputs a function $`f`$ whose value is $`1`$ if the answer is yes and $`0`$ otherwise. This oracle is described by the unitary operator which takes

(A.4) $`|k_0|k_1\mathrm{}k_n|k_0f(k_1\mathrm{}k_n)|k_1\mathrm{}k_n`$
The action of this oracle when the first bit is in the $`\sigma _x`$ eigenstate $`2^{1/2}(|0|1`$ (for which we use the somewhat unconventional notation $`|1_x`$) is then

(A.5) $`|1_x|k_1\mathrm{}k_ne^{i\pi f(k)}|1_x|k_1\mathrm{}k_n`$
which is exactly the effect of $`G`$ when the ancillary initial bit $`|1_x`$ is omitted. The power of quantum computing is then exploited by applying $`G`$ to superpositions of the form $`_kc_{k_1\mathrm{}k_n}|k_1\mathrm{}k_n`$ rather than to the individual states in the computational basis. The analysis is facilitated by the realization that the problem can be reduced to a two-dimensional one in

$`\mathrm{span}\{|t_1\mathrm{}t_n,{\displaystyle \underset{k_1\mathrm{}k_n}{}}|k_1\mathrm{}k_n\}`$
for which $`|t_1\mathrm{}t_n`$ and $`{\displaystyle \frac{1}{\sqrt{N1}}}{\displaystyle \underset{k_1\mathrm{}k_nt_1\mathrm{}t_n}{}}|k_1\mathrm{}k_n`$ form an orthonormal basis.

## Appendix B The Non-Crossing Non-Rule

The so-called “non-crossing rule” is one of a number of physical principles which arise when symmetry ensures that a critical term in some expression, such as the leading term in a perturbation expansion, is zero. The remaining conditions needed to obtain a crossing, transition, etc. are then more easily satisfied. In the absence of the canonical conditions, such crossings and transitions are not truly “forbidden”, but are either rare events which result from an accidental confluence or the result of another physical circumstance which facilitates the satisfying of certain conditions.

In the case of interest here, the canonical condition for crossing is that the Hamiltonians $`H_1`$ and $`H_0`$ both commute with the operators which generate a symmetry group $`𝒢`$ (for which it is not necessary that $`H_1`$ and $`H_0`$ commute with each other). The irreducible representations of $`𝒢`$ can then be used to classify the eigenspaces of $`H_1`$, $`H_0`$, and the interpolating Hamiltonian $`H(s)=(1s)H_1+sH_1`$. The non-crossing rule asserts that one expects eigenvalues of H(s) to cross only if they belong to different irreducible representations.

A similar situation arises in the Born-Oppenheimer approximation for the hydrogen molecule ion $`H_2^+`$ in which case the role of $`t`$ is replaced by the internuclear distance $`R`$. Standard tables of atomic data, show many instances of crossings of states in the same symmetry class, in apparent violation of the non-crossing rule. In this case, the paradox was rather easily resolved by the discovery of a “hidden symmetry” which can regarded as a generalization of that responsible for the well-known “accidental” degeneracies of states of different angular momentum (but same principle quantum number $`n`$) for hydrogen. However, the generator of this symmetry group is an operator (denoted $`F(R)`$ by Judd ) which depends on the internuclear distance $`R`$. Although this gives a satisfactory physical explanation for the phenomena observed, it points out a difficulty with any attempt to make a mathematically precise theorem out of the “non-crossing rule”. If a crossing exists, the degeneracy in $`H(s)`$ always allows the formal construction of a suitable symmetry group .

A less well-known example of violation of the non-crossing rule occurs in Heilman and Lieb’s study  of the Hubbard model for benzene. Moreover, in this case, a rather detailed analysis  failed to locate a hidden symmetry, even one dependent on a parameter. In addition, their work found persistent degeneracies unexplained by symmetry. This is more serious, as it is far less likely to be an artifact of the numerical methods used. It is worth quoting part of their discussion.

> The \[non-crossing rule\] depends\[s\] crucially on the interpretation of the word symmetry. The conventional meaning is that of a symmetry group independent of \[a parameter\] $`U`$; in this case the “proofs” are false. …\[If\] one allows symmetry groups that are $`U`$-dependent, the “theorems” are mere tautologies, because …one can always invent, post hoc a $`U`$-dependent group to account for any violations.
>
> One may ask what is wrong with the “proofs” quoted above. The fault lies not in the mathematics per se but in the assumptions used to connect the mathematics with the real world: First, in the natural sciences, two real numbers are never equal unless there is a physical reason for it; second, that reason must be the existence of a U-independent symmetry group.
>
> The first assumption has validity, but the second is merely a confession of ignorance …

Although the relevance of the non-crossing rule to adiabatic quantum computation is questionable, the principles underlying it are not. If polynomially decreasing gaps are generic for systems in which exponentially many eigenvalues are squeezed into an interval that increases only polynomially, there must be a physical reason for this behavior.


text:
# Triple minima in Free Energy of Semiflexible Polymers

## Abstract

We study the free energy of the worm-like-chain model, in the constant-extension ensemble, as a function of the stiffness $`\lambda `$ for finite chains of length $`L`$. We find that the polymer properties obtained in this ensemble are qualitatively different from those obtained using constant-force ensembles. In particular we find that as we change the stiffness parameter, $`t=L/\lambda `$, the polymer makes a transition from the flexible to the rigid phase and there is an intermediate regime of parameter values where the free energy has three minima and both phases are stable. This leads to interesting features in the force-extension curves.

The simplest model for describing semiflexible polymers without self-avoidance is the so called Worm-Like-Chain (WLC) model . In this model the polymer is modeled as a continuous curve that can be specified by a $`d`$dimensional ($`d>1`$) vector $`\overline{x}(s)`$, $`s`$ being the distance, measured along the length of the curve, from one fixed end. The energy of the WLC model is just the energy due to curvature and is given by

$`{\displaystyle \frac{H}{k_BT}}={\displaystyle \frac{\kappa }{2}}{\displaystyle _0^L}({\displaystyle \frac{\widehat{u}(s)}{s}})^2𝑑s,`$ (1)
where $`\widehat{u}(s)=\overline{x}/s`$ is the tangent vector and satisfies $`\widehat{u}^2=1`$. The parameter $`\kappa `$ specifies the stiffness of the chain and is related to the persistence length $`\lambda `$ defined through $`\widehat{u}(s).\widehat{u}(s^{})=e^{|ss^{}|/\lambda }`$. It can be shown that $`\kappa =(d1)\lambda /2`$ .

The thermodynamic properties of such a chain can be obtained from the free energy which can be either the Helmholtz’s $`(F)`$ free energy or the Gibb’s $`(G)`$ energy. In the former case one considers a polymer whose ends are kept at a fixed distance $`r`$ \[one end fixed at the origin and the other end at $`\overline{r}=(0,\mathrm{}0,r)`$\] by an average force $`f=F(r,L)/r`$, while in the latter case one fixes the force and the average extension is given by $`r=G(f,L)/f`$. It can be shown that in the thermodynamic limit $`L\mathrm{}`$ the two ensembles are equivalent and related by the usual Legendre transform $`G=Ffr`$. For a system with finite $`L/\lambda `$, the equivalence of the two ensembles is not guaranteed, especially when fluctuations become large. We note that real polymers come with a wide range of values of the parameter $`t=L/\lambda `$ \[e.g. $`\lambda 0.1\mu m`$ for DNA while $`\lambda 1\mu m`$ for Actin and their lengths can be varied\] and fluctuations in $`r`$ (or $`f`$) can be very large. Then the choice of the ensemble depends on the experimental conditions. Experiments on stretching polymers are usually performed by fixing one end of the polymer and attaching the other end to a bead which is then pulled by various means (magnetic, optical, mechanical, etc.). In such experiments one can either fix the force on the bead and measure the average polymer extension, or, one could constrain the bead’s position and look at the average force on the polymer. In the former case, the Gibb’s free energy is relevant while it is the Helmholtz in the second case. This point has been carefully analyzed by Kreuzer and Payne in the context of atomic force microscope experiments . Theoretically, the constant-force ensemble is easier to treat, and infact an exact numerical solution has been obtained  (though only for $`t>>1`$). Data on force-extension experiments on DNA  have been explained using this ensemble . The case of constant-extension ensemble turns out to be much harder and no exact solution is available. The $`t0`$ and $`t\mathrm{}`$ cases correspond to the solvable limits of the hard rod and the Gaussian chain. The small and large $`t`$ cases have been treated analytically by perturbation theory about these two limits . Numerical simulations for different values of $`t`$ have been reported by Wilhelm and Frey , who have also obtained series expansions valid in the small $`t`$ limit. A mean-field treatment has also recently been reported .

In this letter we probe the nature of the transition from the Gaussian to the rigid rod with change of stiffness as shown by the form of the Helmholtz free energy of the WLC model (or equivalently the distribution of end-to-end distance). Extensive simulations are performed in two and three dimensions using the equivalence of the WLC model to a random walk with one-step memory. We find the surprising result that, over a range of values of $`t`$, the free energy has three minima. This is verified in a one-dimensional version of the model which is exactly solvable.

We first note that the WLC model describes a particle in $`d`$dimensions moving with a constant speed (set to unity) and with a random acceleration. It is thus described by the propagator

$`Z(\overline{x},\widehat{u},L|\overline{x}^{},\widehat{u}^{},0)={\displaystyle \frac{_{(\overline{x}^{},\widehat{u}^{})}^{(\overline{x},\widehat{u})}𝒟[\overline{x}(s)]e^{H/k_BT}}{𝒟[\overline{x}(s)]e^{H/k_BT}}}`$ (2)
where in the numerator only paths $`\overline{x}(s)`$, satisfying $`\overline{x}(0)=\overline{x}^{},\overline{x}(L)=\overline{x},\widehat{u}(0)=\widehat{u}^{}`$ and $`\widehat{u}(L)=\widehat{u}`$ are considered. It can be shown that the corresponding probability distribution $`W(\overline{x},\widehat{u},L)`$ satisfies the following Fokker-Planck equation :

$`{\displaystyle \frac{W}{L}}+\widehat{u}._{\overline{x}}W{\displaystyle \frac{1}{2\kappa }}_{\widehat{u}}^2W=0`$ (3)
where $`_{\widehat{u}}^2`$ is the diffusion operator on the surface of the unit sphere in $`d`$dimensions. The discretized version of this model is the freely rotating chain model (FRC) of semiflexible polymers . In the FRC one considers a polymer with $`N`$ segments, each of length $`b=L/N`$. Successive segments are constrained to be at a fixed angle, $`\theta `$, with each other. The WLC model is obtained, in the limit $`\theta ,b0`$, $`N\mathrm{}`$ keeping $`\lambda =2b/\theta ^2`$ and $`L=Nb`$ finite.

Here we will consider the situation where the ends are kept at a fixed separation $`r`$ \[with $`\overline{x}^{}`$ at the origin and $`\overline{x}=\overline{r}=(0,\mathrm{}0,r)`$\] but there is no constraint on $`\widehat{u}`$ and $`\widehat{u}^{}`$ and they are taken as uniformly distributed. Thus we will be interested in the distribution $`P(r,L)=\delta (\overline{x}\overline{r})=𝑑\widehat{u}W(\overline{r},\widehat{u},L)`$: this gives the Helmholtz free energy $`F(r,L)=Log[P(r,L)]`$. For the spherically symmetric situation we are considering, $`P(r,L)`$ is simply related to the radial probability distribution $`S(r,L)`$ through $`S(r,L)=Cr^{d1}P(r,L)`$, $`C`$ being a constant equal to the area of the $`d`$dimensional unit sphere. It may be noted that the WLC Hamiltonian is equivalent to spin $`O(d)`$ models in one dimension in the limit of the exchange constant $`J\mathrm{}`$ (with $`Jb=\kappa `$ finite) and all results can be translated into spin language. However, for spin systems, the present free energy is not very relevant since it corresponds to putting unnatural constraints on the magnetization vector.

Numerical simulations: The simulations were performed by generating random configurations of the FRC model and computing the distribution of end-to-end distances. To obtain equivalence with the WLC model the appropriate limits were taken. We note that because these simulations do not require equilibration, they are much faster than simulations on equivalent spin models and give better statistics. The number of configurations generated was around $`10^8`$ for chains of size $`N=10^3`$. We verified that increasing $`N`$ did not change the data significantly. As a check on our numerics we evaluated $`r^2`$ and $`r^4`$. Using Eq. (3) and following  we can compute these (in all dimensions):

$`r^2={\displaystyle \frac{4\kappa L}{d1}}{\displaystyle \frac{8\kappa ^2(1e^{\frac{(d1)L}{2\kappa }})}{(d1)^2}}`$ (4)
$`r^4={\displaystyle \frac{64\kappa ^4(d1)}{d^3(d+1)^2}}e^{\frac{dL}{\kappa }}{\displaystyle \frac{128\kappa ^4(d+5)^2}{(d1)^4(d+1)^2}}e^{\frac{(d1)L}{2\kappa }}+`$ (5)
$`{\displaystyle \frac{64\kappa ^3L(d^28d+7)}{(d1)^4(d+1)}}e^{\frac{(d1)L}{2\kappa }}+{\displaystyle \frac{64\kappa ^4(d^3+23d^27d+1)}{(d1)^4d^3}}`$ (6)
$`{\displaystyle \frac{64\kappa ^3L(d^3+5d^27d+1)}{(d1)^4d^2}}+{\displaystyle \frac{16\kappa ^2L^2(d^33d+2)}{d(d1)^4}}.`$ (7)
Infact it is straightforward to compute all even moments, though it becomes increasingly tedious to get the higher moments. Our numerics agrees with the exact results to around $`0.1\%`$ for $`r^2`$ and $`0.5\%`$ for $`r^4`$.     FIG. 1.:  Monte-Carlo data for $`p(v,t)`$ for the $`2`$dimensional WLC for values of $`t=10(),5,3.33,2`$ and $`1()`$. The inset is a blowup of curves in the transition region ($`t=4,3.33,2.86`$) and clearly shows the presence of the two maxima. Note that because of $`\pm v`$ symmetry, we have plotted data for positive $`v`$ values only. For the fits at large and small $`t`$ see text.       FIG. 2.:  Monte-Carlo data for $`p(v,t)`$ for the $`3`$dimensional WLC for values of $`t=10(),5,3.33,2`$ and $`1()`$. The inset is a blowup of curves in the transition region ($`t=4,3.85,3.7`$) and shows the presence of the two maxima.

The function $`P`$ has the scaling form $`P(r,L)=\frac{1}{L^d}p(r/L,L/\lambda )`$ and we will focus on determining the function $`p(v,t)`$ . In Fig. (1) and Fig. (2), we show the results of our simulations in two and three dimensions. At large values of $`t`$ there is a single maximum at $`v=r/L=0`$ corresponding to a Gaussian distribution while at small $`t`$, the maximum is close to the fully extended value of $`v=\pm 1`$. The transition is first-order-like: as we decrease $`t`$, at some critical value, $`p`$ develops two additional maxima at non-zero values of $`v`$. Further decreasing $`t`$ weakens the maximum at $`v=0`$ until it finally disappears and there are just two maxima which correspond to the rigid chain.

For the limiting cases of small and large values of $`t`$ there are analytic results for the distribution function and as can be seen in Fig. (1,2) our data agrees with them. For large $`t`$ we find that Daniels approximation , which is a perturbation about the Gaussian, fits the data quite well. In the other limit of small $`t`$ the series solutions provided in  fits our data. For intermediate values of $`t`$ neither of the two forms are able to capture, even qualitatively, the features of the free energy. Specifically, we note that all the analytic theories(perturbative, series expansions and mean-field) predict a second-order-like transition and do not give triple minima of the free energy for any parameter value.

It is instructive to study a one-dimensional version of the WLC which shows the same qualitative features (the equivalent spin problem is the Ising model). Consider a $`N`$ step random walk, with step-size $`b`$ which, with probability $`ϵ`$, reverses its direction of motion and with $`1ϵ`$, continues to move in the same direction. The appropriate scaling limit is: $`b0`$, $`ϵ0`$, $`N\mathrm{}`$ keeping $`L=Nb`$, $`t=L/\lambda =2Nϵ`$ finite. Defining $`Z_\pm (x,L)`$ as the probability of the walker to be at $`x`$ with either positive or negative velocity, we have the following Fokker-Planck equation:

$`{\displaystyle \frac{Z_\pm }{L}}={\displaystyle \frac{Z_\pm }{x}}{\displaystyle \frac{1}{2\lambda }}(Z_+Z_{})`$ (8)
This can be solved for $`P(x,L)=Z_++Z_{}=\frac{1}{L}p(x/L,L/\lambda )`$. We get

$`p(v,t)={\displaystyle \frac{te^{t/2}}{4}}[{\displaystyle \frac{I_1(\frac{t}{2}\sqrt{1v^2})}{\sqrt{1v^2}}}+I_0({\displaystyle \frac{t}{2}}\sqrt{1v^2})]`$ (9)
$`+{\displaystyle \frac{e^{t/2}}{2}}[\delta (v1)+\delta (v+1)],`$ (10)
where $`I_0`$ and $`I_1`$ are modified Bessel functions. In Fig. 3 we have plotted $`p(v,t)`$ for different values of stiffness. We find that it always has three peaks. Unlike in $`2`$ and $`3`$ dimensions, the $`\delta `$function peaks at $`v=\pm 1`$ (which corresponds to fully extended chains) persist at all values of stiffness though their weight decays exponentially. Similarly the peak at $`v=0`$ is always present.

Discussion: The most interesting result of this paper is the triple minima in the Helmholtz free energy of the WLC. Physically, this results from the competing effects of entropy, which tries to pull in the polymer and the bending energy, which tries to extend it. This form of the free energy leads to a highly counterintuitive force-extension curve, very different from what one obtains from the constant force ensemble or from approximate theories. In Fig. (4) we show the force-extension curve for a two dimensional chain with $`t=3.33`$. We see that there are two stable positions for which the force is zero. In the constant-force ensemble, it is easy to show that $`r/f=r^2r^2`$ and so the force-extension is always monotonic. In the constant-extension ensemble, there is no analogous result (for finite systems), and monotonicity is not guaranteed.     FIG. 3.:  The exact distribution $`p(v,t)`$ of the $`1`$dimensional WLC \[ Eq. (10)\] for different values of $`t`$ ($`10,5,3.33,2,1`$). Even for the most stiff chain considered here ($`t=1`$), the distribution has a peak at the centre (in addition to the $`\delta `$function peaks at ends) though it looks flat.

Most of the recent experiments on stretching DNA have $`t100`$. The distribution is then sharply peaked at zero and one expects the equivalence of different ensembles. Experimentally the value of $`t`$ can be tuned by various means, for example, by changing the length of the polymer or the temperature. Polymer-stretching experiments can thus be performed for intermediate $`t`$ values. Since we consider the tangent vectors at the polymer-ends to be unconstrained an accurate experimental realization of our set-up would be one in which both ends are attached to beads \[see Fig. (5)\]. The beads are put in optical traps and so are free to rotate (this setup is identical to the one used in refn. ). Making the traps stiff corresponds to working in the constant-extension ensemble  and one can measure the average force. Our predictions can then be experimentally verified. We make some estimates on the experimental requirements (for a $`3`$-d polymer with stiffness $`t=3.85`$). Assume that at one end, the origin, the trap is so stiff that the bead can only rotate. We make measurements at the other end. The trap-center is placed at $`\overline{r}_0=(0,0,z_0)`$ and the mean bead displacement $`\mathrm{\Delta }z=(zz_0)`$ gives the mean force $`f`$ on the polymer. We then consider the problem of the polymer in the presence of a trap potential $`V=[k_t(x^2+y^2)+k(zz_0)^2]/2`$. Assume $`k_t>>k`$ so we can neglect fluctuations in the transverse directions. The distribution of the bead’s position in the presence of the potential is given by     FIG. 4.:  The free energy (dotted line) and the corresponding force-extension curve (solid line) for a $`2`$dimensional chain with $`t=3.33`$.       FIG. 5.:  A schematic of the experimental set-up required to realize the constant-extension ensemble discussed in the paper (see refn. ). For a stiff trap the average displacement of the bead $`\mathrm{\Delta }z`$ from the trap center is small and the average force on the polymer is: $`f=k\mathrm{\Delta }z`$.   $`Q(\overline{r})=e^{\beta [F(\overline{r})+V(\overline{r})]}/d^3\overline{r}e^{\beta [F(\overline{r})+V(\overline{r})]}`$. For a stiff trap, we can expand $`F`$ about $`\overline{r}=\overline{r}_0`$ and find that the average displacement of the bead is given by: $`\mathrm{\Delta }z=d^3\overline{r}(zz_0)Q(\overline{r})=f/k^{}`$, where $`k^{}=k+F^{\prime \prime }(z_0)k`$ (valid except when $`z_0L`$) and $`f=F^{}(z_0)`$. The rms fluctuation of the bead about the trap center is given by $`z_{rms}^2=k_BT/k`$. Hence we get $`\mathrm{\Delta }z=fz_{rms}^2/(k_BT)=(fL/k_BT)(z_{rms}^2/L)`$. The scaled force $`fL/(k_BT)`$ is of order $`0.1`$. The different minima are separated by distances $`0.2L`$, hence to see the effect we need to have $`z_{rms}/L0.1`$. Thus finally we find that the typical displacement of the bead $`\mathrm{\Delta }z`$ is about $`0.01z_{rms}`$. This is quite small and means that it is necessary to collect data on the bead position over long periods of time.

As suggested in , a more direct way of measuring the Helmholtz free energy would be to attach marker molecules at the ends of the polymer and determine the distribution of end-to-end distances. Fluorescence microscopy as in  could be another possible method. It is to be remembered of course that real polymers are well-modeled by the WLC model provided we can neglect monomer-monomer interactions (steric, electrolytic etc.). Thus the experiments would really test the relevance of the WLC model in describing real semiflexible polymers in different stiffness regimes.

In conclusion we have presented some new and interesting properties of the WLC model and have pointed out that polymer properties are ensemble-dependent. In this paper we have given one example of qualitative differences in force-extension measurements in different ensembles. Other quantitative differences will occur even in more flexible chains and should be easier to observe experimentally. We hope this work will motivate further experimental and theoretical work on this simplest model for semiflexible polymers.

We thank O. Narayan and J. Samuel for discussions. One of us (D.C.) thanks CSIR, India for support.

Note: After submission of this paper, an exact numerical solution of the WLC model has been obtained, and has reproduced our results .


text:
# 1 Introduction

## 1 Introduction

Supersymmetry (symmetry between the fermionic and bosonic degrees of freedom) has important role  in analyzing of the quantum mechanical systems, since it can study remarkable properties including the degeneracy structure of the energy spectrum, the relations among the energy spectra of the various Hamiltonians and etc.. In particular in this theory, the energy eigenvalues are necessarily non-negative and the energy of non-zero (zero) ground state is related to the broken (unbroken) supersymmetry. For the first time, Rubakov and Spiridonov  extended the supersymmetry which is called parasupersymmetry and it describes an essential symmetry between bosons and parafermions. Then, it was realized that the parasupersymmetry presented in Ref.  is of order $`p=2`$ and Khare  generalized the Rubakov-Spiridonov (R-S) parasupersymmetry to arbitrary order $`p1`$:

$`Q_1^pQ_1^{}+Q_1^{p1}Q_1^{}Q_1+\mathrm{}+Q_1^+Q_1^p=2pQ_1^{p1}H,`$
$`Q_1^{p+1}=0,`$
$`[H,Q_1]=0.`$ (1)
Here, $`Q_1`$ and $`H`$ stand for a parasupercharge and the bosonic Hamiltonian. The relations (1) for $`p=1`$ and $`p=2`$ are reduced to the supersymmetry algebra and the Rubakov-Spiridonov (R-S) parasupersymmetry algebra , respectively. In addition, the relations (1) describe the unitary parasupersymmetry algebra of arbitrary order $`p`$, since $`Q_1^{}`$ is the Hermitian conjugate of the parasupercharge $`Q_1`$. So, it is evident that similar relations are satisfied under interchange of $`Q_1`$ and $`Q_1^{}`$ in (1).

Before the extension of the R-S parasupersymmetry algebra to arbitrary order $`p`$, for the $`p=2`$ case remarkable and interesting discussions had been made \[6-9\]. For instance, the motion of a spin-1 particle along the $`z`$-axis in the presence of a magnetic field can be described by a parasupersymmetric Hamiltonian which is obtained by the simple harmonic oscillator and Morse potentials , and also the problem was generalized  for the spin-$`\frac{p}{2}`$ particles. Meanwhile, in the context of quantum field theory, the R-S parasupersymmetry algebra of order $`p=2`$ leads to the infinite bosonic and parafermionic variables . In a special formulation of the R-S parasupersymmetric quantum mechanics , the Hamiltonian for which the energy spectrum cannot be negative is expressed in terms of an explicit function of the parasupercharge $`Q_1`$. However, in general for the Rubakov-Spiridonov parasupersymmetric theory of arbitrary order the bosonic Hamiltonian cannot be obtained directly in terms of the parasupercharge $`Q_1`$, and the energy eigenvalues are not necessarily non-negative. Moreover, there is no connection between the non-zero (zero) ground-state energy and the broken (unbroken) parasupersymmetry. Also, it has been shown  that there are $`p1`$ other conserved parasupercharges and $`p`$ bosonic constants in the R-S parasupersymmetric theory of arbitrary order $`p`$.

In an early work, Infeld and Hull  studied the factorization and algebraic solutions of the bound state problems and later, Gendenshtein et.al. considered the subject in the framework of the shape invariance symmetry as an important aspect of solvability of wide range of the 1D quantum mechanical models. It must be mentioned that the factorization and the shape invariance symmetry have obtained a very helpful approach in the representation of the supersymmetry theory \[16-24\]. Recently, most of the one dimensional shape invariant solvable quantum mechanical models have been classified into two bunches. The first bunch  includes models for which the shape invariance parameter is the main quantum number $`n`$. On the other hand, in the second bunch  the shape invariance parameter of the models is the secondary quantum number $`m`$. Meanwhile, it has been shown in Ref.  that the R-S parasupersymmetry algebra of arbitrary order $`p`$ is realized by the shape invariant quantum mechanical models so that the algebra can be represented by the quantum mechanical states of the models. For realizing the algebra, it has also been shown that the bosonic Hamiltonian involves $`p+1`$ isospectrum Hamiltonians. This fact has been studied in detail for the second bunch of the shape invariant models in Ref. . In fact, $`p+1`$ isospectrum Hamiltonians are obtained by adding $`p`$ appropriate constants to the factorized Hamiltonians $`\frac{1}{2}A(1)B(1)`$, $`\frac{1}{2}A(2)B(2)`$, …, $`\frac{1}{2}A(p)B(p)`$, and by adding the corresponding constant of the last Hamiltonian to its partner Hamiltonian i.e. $`\frac{1}{2}B(p)A(p)`$ as well. Here, the operators $`B(l)`$ and $`A(l)`$; $`l=1,2,\mathrm{},p`$ are the raising and lowering operators of the quantum states of the shape invariance theory, respectively. One of the other successes of the R-S parasupersymmetry theory of arbitrary order $`p`$ is the fact that it can be realized \[26-28\] by the 2D and 3D solvable quantum mechanical models obtained from the shape invariance symmetry.

In 1990 another formulation of the unitary parasupersymmetry algebra of order $`p=2`$ was introduced by Beckers and Debergh  as follows:

$`[Q_1,[Q_1^{},Q_1]]=2Q_1H`$
$`[Q_1^{},[Q_1,Q_1^{}]]=2Q_1^{}H`$
$`Q_1^3=Q_{1}^{}{}_{}{}^{3}=0`$
$`[H,Q_1]=[H,Q_1^{}]=0.`$ (2)
The appropriate form of the parasupersymmetric Hamiltonian presented in the Eqs. (2) has been constructed in Ref.  by using the generators of the simple harmonic oscillator, and the algebra (2) has been represented by the quantum states of the simple harmonic oscillator. But, a generalization of the Beckers-Debergh (B-D) parasupersymmetry algebra to arbitrary order $`p2`$ which is satisfied by the quantum mechanical systems has not been found yet. In general, the R-S parasupersymmetry algebra has been much more successful than the B-D parasupersymmetry algebra. Nevertheless, for the B-D parasupersymmetry algebra of order $`p=2`$, it has been performed useful discussions. For example, Mostafazadeh has proved that  for both Rubakov-Spiridonov and Beckers-Debergh formulations of the parasupersymmetric quantum mechanics of order $`p=2`$, the degeneracy structure of the energy spectrum can be derived using a thorough analysis of the parasupersymmetry algebra. He showed that the result is independent of the details of the Hamiltonian, for example, the degeneracy structure is not related to the dimension of the coordinate manifolds that the Hamiltonian is defined on it. Moreover, it has been shown that in general the Rubakov-Spiridonov (R-S) and Beckers-Debergh (B-D) systems possess identical degeneracy structures . Also, similar to the Witten index of the supersymmetric quantum mechanics, for the two kinds of ($`p=2`$) parasupersymmetry systems (B-D and R-S) a new set of topological invariants has been obtained .

In this paper it is intended to generalize the unitary parasupersymmetry algebra of B-D. Also, we introduce a special representation of the B-D unitary parasupersymmetry algebra of even arbitrary order by 1D shape invariance solvable models, and some 2D and 3D quantum solvable models as well. Meanwhile, the B-D unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$ is introduced with $`2k`$ independent conserved parasupercharges.

## 2 Towards an appropriate generalization of the B-D parasupersymmetry algebra

In order to generalize the B-D unitary parasupersymmetry algebra, we shall take into account the parafermionic operators $`b`$ and $`b^{}`$ of arbitrary order $`p`$ which have been used by Khare . In fact in moving from statistics to parastatistics, the parafermionic operators are considered as the following $`(p+1)\times (p+1)`$ matrices:

$$(b)_{i,j}=C_j\delta _{i,j+1}(b^{})_{i,j}=C_i\delta _{i+1,j}i,j=1,2,\mathrm{},p+1$$
(3)
where the coefficients $`C_j`$ are given by

$$C_j:=\sqrt{j(pj+1}=C_{pj+1}j=1,2,\mathrm{},p+1.$$
(4)
It can be easily shown that

$$C_1C_2\mathrm{}C_p=p!$$
(5)
Using the parafermionic operators $`b`$ and $`b^{}`$, one may obtain a spin-$`\frac{p}{2}`$ representation for the group SU(2). Indeed, by defining the operators $`J_\pm `$ and $`J_3`$ as:

$`J_+:=b^{}J_{}:=b`$
$`J_3:={\displaystyle \frac{1}{2}}[b^{},b]=diag({\displaystyle \frac{p}{2}},{\displaystyle \frac{p}{2}}1,\mathrm{},{\displaystyle \frac{p}{2}}+1,{\displaystyle \frac{p}{2}}),`$ (6)
the following commutation relations corresponding to the Lie algebra SU(2) may be derived

$$[J_+,J_{}]=2J_3[J_3,J_\pm ]=\pm J_\pm .$$
(7)
It is noticed that the commutation relation of the parafermionic operators $`b`$ and $`b^{}`$ is proportional to the third component of the spin-$`\frac{p}{2}`$ representation of the Lie group SU(2).

One of the well-known and important properties of the parafermionic operators is

$$b^{p+1}=b_{}^{}{}_{}{}^{p+1}=0.$$
(8)
Now it can be verified that the parafermionic operators $`b`$ and $`b^{}`$ of order $`p`$ possess the following new multilinear structural relations ($`p2`$):

$$\underset{l=0}{\overset{p}{}}(1)^{l+1}C_l^pb^{pl}b^{}b^l=2\delta _{p,2}(1)^pb^{p1}$$
(9a)
$$\underset{l=0}{\overset{p}{}}(1)^{l+1}C_l^pb_{}^{}{}_{}{}^{pl}bb_{}^{}{}_{}{}^{l}=2\delta _{p,2}b_{}^{}{}_{}{}^{p1}$$
(9b)
where the constants $`C_l^p;l=0,\mathrm{},p`$ are the Newton binomial expansion coefficients. Using the Baker-Hausdorff formula for two arbitrary operators $`A`$ and $`B`$ which is given by

$$e^{\lambda A}Be^{\lambda A}=B+\frac{\lambda }{1!}[A,B]+\frac{\lambda ^2}{2!}[A,[A,B]]+\mathrm{},$$
(10)
the structural relations (9) between the operators $`b`$ and $`b^{}`$ may be written down as

$`\underset{}{[b,[b,\mathrm{},[b},[b^{},b]]\mathrm{}]]=2\delta _{p,2}(1)^pb^{p1}`$
$`(p1)\text{times}`$ (11a)
$`\underset{}{[b^{},[b^{},\mathrm{},[b^{}},[b,b^{}]]\mathrm{}]]=2\delta _{p,2}b^{}^{p1}.`$
$`(p1)\text{times}`$ (11b)
It is seen that the multilinear expressions in the left-hand sides of the relations (11) are described in terms of the parafermionic operators $`b`$ and $`b^{}`$ in the right-hand sides. Actually, the multilinear relations (11) between the parafermionic operators $`b`$ and $`b^{}`$ propose similar relations between the parasupercharges of the parasupersymmetric quantum mechanics of order $`p`$. The relations (8) and (11) indicate that there exist the conserved parasupercharges $`Q_1`$ and $`Q_1^{}`$ of order $`p`$, and the bosonic Hamiltonian so that they generate the following parasupersymmetry algebra:

$`\underset{}{[Q_1,[Q_1,\mathrm{},[Q_1},[Q_1^{},Q_1]]\mathrm{}]]=2(1)^pQ_1^{p1}H`$
$`(p1)\text{times}`$ (12a)
$`\underset{}{[Q_1^{},[Q_1^{},\mathrm{},[Q_1^{}},[Q_1,Q_1^{}]]\mathrm{}]]=2Q_1^{}^{p1}H`$
$`(p1)\text{times}`$ (12b)
$$Q_1^{p+1}=Q_{1}^{}{}_{}{}^{p+1}=0$$
(12c)
$$[H,Q_1]=0$$
(12d)
$$[H,Q_1^{}]=0.$$
(12e)
It can be easily verified that the relations (12) for any arbitrary $`p`$ are closure under Hermitian conjugation. Meanwhile, the algebraic relations (2) are a special case ($`p=2`$) of the relations (12). Therefore, the relations (12) is an appropriate generalization of the B-D unitary parasupersymmetry algebra with the bosonic Hamiltonian $`H`$ and the parafermions of arbitrary order $`p`$ i.e. $`Q_1`$ and $`Q_1^{}`$.

## 3 A special representation for realizing of the even arbitrary order unitary parasupersymmetry algebra of B-D by quantum solvable models

In this section we analyze a special representation for the B-D quantum mechanical unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$ by wide range of the 1D, 2D and 3D solvable models. In fact, the 1D models are the solvable models obtained from the two approaches of the factorization with respect to the main and secondary quantum numbers i.e. $`n`$ and $`m`$. On the other hand, the 2D and 3D models representing the algebraic relations (12) for $`p=2k`$ are some known quantum mechanical models on the homogeneous manifolds $`SL(2,c)/GL(1,c)`$ and the group manifolds $`SL(2,c)`$. Now in order to obtain the mentioned representations, we only introduce two bunches of the shape invariance models which have been classified before .

In master function theory, a function $`A(x)`$ which is at most of second order in terms of $`x`$, and a non-negative weight function $`W(x)`$ defined in an interval $`(a,b)`$ may be chosen so that $`(1/W(x))(d/dx)(A(x)W(x))`$ is a polynomial of at most first order. For a given master function $`A(x)`$ and its corresponding weight function $`W(x)`$, it has been shown that the eigenvalue equations of the one dimensional partner Hamiltonians corresponding to the first bunch of the superpotentials, which are obtained from the factorization with respect to the main quantum number $`n`$, will be

$`B(n)A(n)\psi _n(\theta )=E(n)\psi _n(\theta )`$
$`A(n)B(n)\psi _{n1}(\theta )=E(n)\psi _{n1}(\theta ),`$ (13)
where the variable $`\theta `$ is introduced by means of solving the following first order differential equation

$`d\theta ={\displaystyle \frac{dx}{A(x)}}.`$ (14)
The energy spectrum $`E(n)`$ and the wave function $`\psi _n(\theta )`$ are given by

$`E(n)={\displaystyle \frac{n}{4\left[\left(\frac{A(x)W^{}(x)}{W(x)}\right)^{}+nA^{\prime \prime }(x)\right]^2}}\{4\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)_{}^{}{}_{}{}^{2}`$
$`\times \left(nA_{}^{}{}_{}{}^{2}(0)A(0)\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)^{}\right)`$
$`\left({\displaystyle \frac{AW^{}}{W}}\right)(0)\left(A^{\prime \prime }(x)\left({\displaystyle \frac{AW^{}}{W}}\right)(0)2A^{}(0)\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)^{}\right)`$
$`\times \left(2\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)^{}+nA^{\prime \prime }(x)\right)`$
$`+n^2A^{\prime \prime }(0)\left(A_{}^{}{}_{}{}^{2}(0)2A^{\prime \prime }(0)A(0)\right)\left(nA^{\prime \prime }(0)+4\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)^{}\right)`$
$`10nA(0)A^{\prime \prime }(x)\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)_{}^{}{}_{}{}^{2}\},`$ (15)
$`\psi _n(\theta )=\left[{\displaystyle \frac{a_n}{\sqrt{W(x)}}}\left({\displaystyle \frac{d}{dx}}\right)^n\left(A^n(x)W(x)\right)\right]_{x=x(\theta )}.`$ (16)
Note that the prime symbol indicates the derivative with respect to $`x`$. Moreover, the explicit forms of the raising and lowering operators corresponding to the main quantum number $`n`$ are, respectively:

$`B(n)={\displaystyle \frac{d}{d\theta }}+{\displaystyle \frac{1}{2}}[nA^{}(x)+{\displaystyle \frac{A(x)W^{}(x)}{W(x)}}.`$
$`.+n{\displaystyle \frac{A^{}(0)\left(\frac{A(x)W^{}(x)}{W(x)}\right)^{}A^{\prime \prime }(x)\left(\frac{AW^{}}{W}\right)(0)}{\left(\frac{A(x)W^{}(x)}{W(x)}\right)^{}+nA^{\prime \prime }(x)}}]_{x=x(\theta )}`$
$`A(n)={\displaystyle \frac{d}{d\theta }}+{\displaystyle \frac{1}{2}}[nA^{}(x)+{\displaystyle \frac{A(x)W^{}(x)}{W(x)}}.`$
$`.+n{\displaystyle \frac{A^{}(0)\left(\frac{A(x)W^{}(x)}{W(x)}\right)^{}A^{\prime \prime }(x)\left(\frac{AW^{}}{W}\right)(0)}{\left(\frac{A(x)W^{}(x)}{W(x)}\right)^{}+nA^{\prime \prime }(x)}}]_{x=x(\theta )}`$ (17)
The change of variable $`x=x(\theta )`$ is substituted in the relations (16) and (17) by solving the first order differential equation (14). Now, by choosing the suitable normalization coefficients $`a_n`$ for the wave functions $`\psi _n(\theta )`$, one may write down the shape invariance equations (13) as the raising and lowering relations:

$`B(n)\psi _{n1}(\theta )=\sqrt{E(n)}\psi _n(\theta )`$
$`A(n)\psi _n(\theta )=\sqrt{E(n)}\psi _{n1}(\theta ).`$ (18)
The potentials like Coulomb, Rosen-Morse I, Rosen-Morse II and Eckart are involved in the first bunch of the solvable models.

In master function theory, the eigenvalue equations of the one dimensional partner Hamiltonian corresponding to the second bunch of the superpotentials which are obtained from the factorization with respect to the secondary quantum number $`m`$ are given by

$`B(m)A(m)\psi _{n,m}(\theta )=E(n,m)\psi _{n,m}(\theta )`$
$`A(m)B(m)\psi _{n,m1}(\theta )=E(n,m)\psi _{n,m1}(\theta ).`$ (19)
In the factorization equations (19), the variable $`\theta `$ is introduced by solving the following first order differential equation

$`d\theta ={\displaystyle \frac{dx}{\sqrt{A(x)}}}.`$ (20)
The energy spectrum and the eigenfunctions of the partner Hamiltonians obtained from the factorization with respect to the secondary quantum number $`m`$ are

$`E(n,m)=(nm+1)\left[\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)^{}+{\displaystyle \frac{1}{2}}(n+m)A^{\prime \prime }(x)\right]`$ (21)
$`\psi _{n,m}(\theta )=\left[{\displaystyle \frac{a_{n,m}}{A^{(2m1)/4}(x)W^{1/2}(x)}}\left({\displaystyle \frac{d}{dx}}\right)^{nm}\left(A^n(x)W(x)\right)\right]_{x=x(\theta )},`$ (22)
where $`m=0,1,2,\mathrm{},n`$. The explicit forms of the raising and lowering operators corresponding to the secondary quantum number $`m`$ are given by

$`B(m)={\displaystyle \frac{d}{d\theta }}\left[{\displaystyle \frac{\frac{A(x)W^{}(x)}{2W(x)}+\frac{2m1}{4}A^{}(x)}{\sqrt{A(x)}}}\right]_{x=x(\theta )}`$
$`A(m)={\displaystyle \frac{d}{d\theta }}\left[{\displaystyle \frac{\frac{A(x)W^{}(x)}{2W(x)}+\frac{2m1}{4}A^{}(x)}{\sqrt{A(x)}}}\right]_{x=x(\theta )}.`$ (23)
This time the change of variable $`x=x(\theta )`$ is substituted in the Eqs. (22) and (23) by solving the first order differential equation (20). Similar to the Eqs. (18), one may write down the raising and lowering relations corresponding to the secondary quantum number $`m`$ on the wave functions $`\psi _{n,m}(\theta )`$ by using the factorization equation (19) as:

$`B(m)\psi _{n,m1}(\theta )=\sqrt{E(n,m)}\psi _{n,m}(\theta )`$
$`A(m)\psi _{n,m}(\theta )=\sqrt{E(n,m)}\psi _{n,m1}(\theta ).`$ (24)
The potentials like 3D harmonic oscillator, Scarf I, Scarf II, Natanzon and generalized Pöschl-Teller are included in the second bunch of the solvable models. Moreover, simple harmonic oscillator and Morse potentials belong to both of the solvable models.

Now let us analyze a special realization of the B-D unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$ by means of the one dimensional quantum mechanical solvable models which are obtained from the factorization with respect to the secondary quantum number $`m`$. Similar procedure can be made by means of the one dimensional quantum mechanical models which are obtained from the shape invariance with respect to the main quantum number $`n`$. In order to realize the B-D unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$, one may define the parafermionic generators $`Q_1`$ and $`Q_1^{}`$ of order $`p=2k`$, and the bosonic operator $`H`$ as the following $`(2k+1)\times (2k+1)`$ matrices:

$`\left(Q_1\right)_{ll^{}}:=\delta _{l+1,l^{}}A(m),l=\text{odd}`$
$`\left(Q_1\right)_{ll^{}}:=\delta _{l+1,l^{}}B(m),l=\text{even}`$ (25a)
$`\left(Q_1^{}\right)_{ll^{}}:=\delta _{l,l^{}+1}A(m),l=\text{odd}`$
$`\left(Q_1^{}\right)_{ll^{}}:=\delta _{l,l^{}+1}B(m),l=\text{even}`$ (25b)
$$\left(H\right)_{ll^{}}:=\delta _{l,l^{}}H_l,l,l^{}=1,2,\mathrm{},2k+1$$
(25c)
where $`m=1,2,\mathrm{},n`$. It is evident that by choosing the definitions (25a) and (25b) for $`Q_1`$ and $`Q_1^{}`$, the relations (12c) are satisfied automatically. To satisfy the Eqs. (12a) and (12b), by using the definitions (25) in them we obtain

$`\left({\displaystyle \underset{l=0}{\overset{2k1}{}}}(1)^{l+1}C_l^{2k}\right)\left(A(m)B(m)\right)^kA(m)=2\left(A(m)B(m)\right)^{k1}A(m)H_{2k}`$
$`\left({\displaystyle \underset{l=1}{\overset{2k}{}}}(1)^{l+1}C_l^{2k}\right)\left(B(m)A(m)\right)^kB(m)=2\left(B(m)A(m)\right)^{k1}B(m)H_{2k+1}`$ (26a)
$`\left({\displaystyle \underset{l=1}{\overset{2k}{}}}(1)^{l+1}C_l^{2k}\right)\left(B(m)A(m)\right)^kB(m)=2\left(B(m)A(m)\right)^{k1}B(m)H_1`$
$`\left({\displaystyle \underset{l=0}{\overset{2k1}{}}}(1)^{l+1}C_l^{2k}\right)\left(A(m)B(m)\right)^kA(m)=2\left(A(m)B(m)\right)^{k1}A(m)H_2.`$ (26b)
Considering the following identities

$$\underset{l=0}{\overset{2k1}{}}(1)^{l+1}C_l^{2k}=\underset{l=1}{\overset{2k}{}}(1)^{l+1}C_l^{2k}=1,$$

then the Eqs. (26) give the following results

$`H_1=H_{2k+1}={\displaystyle \frac{1}{2}}A(m)B(m)`$
$`H_2=H_{2k}={\displaystyle \frac{1}{2}}B(m)A(m).`$ (27)
Meanwhile, by using the definitions (25), the Eqs. (12d) and (12e) lead to the following relations, respectively

$`A(m)H_{2l}=H_{2l1}A(m)`$
$`B(m)H_{2l+1}=H_{2l}B(m),l=1,2,\mathrm{},k`$ (28a)
$`A(m)H_{2l}=H_{2l+1}A(m)`$
$`B(m)H_{2l1}=H_{2l}B(m),l=1,2,\mathrm{},k.`$ (28b)
It is noticed that the Eqs. (27) satisfy the relations (28a) and (28b), and additionally, in order to determine the remaining components of the bosonic Hamiltonian $`H`$ it is sufficient to substitute the relations (27) in the recursion Eqs. (28a) and (28b). Then, one may obtain consistent solutions which are the same for the Eqs. (28a) and (28b) as:

$$H_{2l1}=\frac{1}{2}A(m)B(m),l=1,2,\mathrm{},k+1$$
(29a)
$$H_{2l}=\frac{1}{2}B(m)A(m),l=1,2,\mathrm{},k.$$
(29b)
The recent result declares that the isospectrum Hamiltonians of the B-D unitary parasupersymmetry theory of even arbitrary order $`p=2k`$ are the two partner Hamiltonians $`\left(\frac{1}{2}A(m)B(m)((k+1)\text{times})\text{and}\frac{1}{2}B(m)A(m)(k\text{times})\right)`$ of the shape invariance theory with the energy spectrum $`\frac{1}{2}E(n,m)`$.

Clearly, the following $`(2k+1)\times 1`$ columns matrix

$$\mathrm{\Psi }(\theta )=\left(\begin{array}{cc}\psi _{n,m1}(\theta )\hfill & \\ \psi _{n,m}(\theta )\hfill & \\ \psi _{n,m1}(\theta )\hfill & \\ \psi _{n,m}(\theta )\hfill & \\ \mathrm{}\hfill & \\ \psi _{n,m1}(\theta )\hfill & \end{array}\right)_{(2k+1)\times 1}$$
(30)
as the basis represent the B-D unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$. The eigenvalue equation of the bosonic Hamiltonian $`H`$ is written down as

$$H\mathrm{\Psi }(\theta )=E(n,m)\mathrm{\Psi }(\theta ).$$
(31)
The representation of the parafermionic generators $`Q_1`$ and $`Q_1^{}`$ on the basis (30) by using the Eqs. (24) has the following forms

$$Q_1\mathrm{\Psi }(\theta )=\sqrt{E(n,m)}\left(\begin{array}{cc}\psi _{n,m1}(\theta )\hfill & \\ \psi _{n,m}(\theta )\hfill & \\ \psi _{n,m1}(\theta )\hfill & \\ \psi _{n,m}(\theta )\hfill & \\ \mathrm{}\hfill & \\ \psi _{n,m}(\theta )\hfill & \\ \mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}\hfill & \end{array}\right),Q_1^{}\mathrm{\Psi }(\theta )=\sqrt{E(n,m)}\left(\begin{array}{cc}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}\hfill & \\ \psi _{n,m}(\theta )\hfill & \\ \psi _{n,m1}(\theta )\hfill & \\ \psi _{n,m}(\theta )\hfill & \\ \psi _{n,m1}(\theta )\hfill & \\ \mathrm{}\hfill & \\ \psi _{n,m1}(\theta )\hfill & \end{array}\right).$$
(32)
It is easily seen that the states $`Q_1^l\mathrm{\Psi }(\theta )`$ and $`Q_{1}^{}{}_{}{}^{l}\mathrm{\Psi }(\theta )`$ ($`l=1,2,\mathrm{},p=2k`$) are eigenfunctions of the bosonic Hamiltonian $`H`$.

If we consider the first bunch of the shape invariance models, we will be able to construct the parafermionic generators by means of the raising and lowering operators $`B(n)`$ and $`A(n)`$, therefore we can obtain the bosonic Hamiltonian $`H`$ including two independent partner components $`\frac{1}{2}A(n)B(n)`$ and $`\frac{1}{2}B(n)A(n)`$ with the same energy spectrum $`\frac{1}{2}E(n)`$. In this case, the basis of the representation of the B-D unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$ is constructed by the eigenfunctions $`\psi _n(\theta )`$ and $`\psi _{n1}(\theta )`$.

In Ref.  the following eigenvalue equations have been obtained by using the shape invariance equations (19)

$`L_+L_{}\psi _{n,m}(\theta ,\varphi )=E(n,m)\psi _{n,m}(\theta ,\varphi ),`$
$`L_{}L_+\psi _{n,m1}(\theta ,\varphi )=E(n,m)\psi _{n,m1}(\theta ,\varphi ),`$ (33)
where the explicit differential forms of the operators $`L_+`$ and $`L_{}`$ are given by

$`L_+=e^{i\varphi }({\displaystyle \frac{}{\theta }}+{\displaystyle \frac{i}{2}}{\displaystyle \frac{A^{}(x)}{\sqrt{A(x)}}}|_{x=x(\theta )}{\displaystyle \frac{}{\varphi }}`$
$`[{\displaystyle \frac{1}{2\sqrt{A(x)}}}\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)+{\displaystyle \frac{1}{4\sqrt{A(x)}}}A^{}(x)]_{x=x(\theta )})`$
$`L_{}=e^{i\varphi }({\displaystyle \frac{}{\theta }}+{\displaystyle \frac{i}{2}}{\displaystyle \frac{A^{}(x)}{\sqrt{A(x)}}}|_{x=x(\theta )}{\displaystyle \frac{}{\varphi }}`$
$`[{\displaystyle \frac{1}{2\sqrt{A(x)}}}\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right){\displaystyle \frac{1}{4\sqrt{A(x)}}}A^{}(x)]_{x=x(\theta )}).`$ (34)
Now taking into account the following operators

$`L_3=i{\displaystyle \frac{}{\varphi }}`$
$`I=1`$ (35)
then, the operators $`L_+`$, $`L_{}`$, $`L_3`$ and $`I`$ constitute the Lie algebra $`gl(2,c)`$ i.e.

$`[L_+,L_{}]=A^{\prime \prime }(x)L_3\left({\displaystyle \frac{A(x)W^{}(x)}{W(x)}}\right)^{}I`$
$`[L_3,L_\pm ]=\pm L_\pm `$
$`[\text{L},I]=0.`$ (36)
The change of variable $`x=x(\theta )`$, which is used in the Eqs. (34), is obtained by solving the differential Eq. (20). It has been shown in Ref.  that the Casimir of the generators $`L_+`$, $`L_{}`$, $`L_3`$ and $`I`$ is the corresponding Hamiltonian of the charged particle on the homogeneous manifolds $`SL(2,c)/GL(1,c)`$ in the presence of magnetic monopole with degeneracy group $`GL(2,c)`$. The wave functions $`\psi _{n,m}(\theta ,\varphi )`$ , which represent the Lie algebra $`gl(2,c)`$ as

$`L_+\psi _{n,m1}(\theta ,\varphi )=\sqrt{E(n,m)}\psi _{n,m}(\theta ,\varphi )`$
$`L_{}\psi _{n,m}(\theta ,\varphi )=\sqrt{E(n,m)}\psi _{n,m1}(\theta ,\varphi )`$
$`L_3\psi _{n,m}(\theta ,\varphi )=m\psi _{n,m}(\theta ,\varphi )`$
$`I\psi _{n,m}(\theta ,\varphi )=\psi _{n,m}(\theta ,\varphi ),`$ (37)
describe the two dimensional quantum states of the charged particle on the homogeneous manifolds $`SL(2,c)/GL(1,c)`$ in the presence of the magnetic monopole and they are given by

$$\psi _{n,m}(\theta ,\varphi )=e^{im\varphi }\psi _{n,m}(\theta ).$$
(38)
Now it can be easily shown that the quantum states $`\psi _{n,m}(\theta ,\varphi )`$ also represent the B-D unitary parasupersymmetry algebra of arbitrary order $`p=2k`$. In order to show the mentioned fact it is sufficient to define the parafermionic generators $`Q_1`$ and $`Q_1^{}`$ of order $`p=2k`$, and the bosonic operator $`H`$ as

$`\left(Q_1\right)_{ll^{}}:=\delta _{l+1,l^{}}L_{},l=\text{odd}`$
$`\left(Q_1\right)_{ll^{}}:=\delta _{l+1,l^{}}L_+,l=\text{even}`$ (39a)
$`\left(Q_1^{}\right)_{ll^{}}:=\delta _{l,l^{}+1}L_{},l=\text{odd}`$
$`\left(Q_1^{}\right)_{ll^{}}:=\delta _{l,l^{}+1}L_+,l=\text{even}`$ (39b)
$$\left(H\right)_{ll^{}}:=\delta _{l,l^{}}H_l,l,l^{}=1,2,\mathrm{},2k+1.$$
(39c)
Once again, the relations (39a) and (39b) satisfy the relation (12c) automatically. Using the definitions (39a), (39b) and (39c), the Eqs. (12a) and (12b) lead to the following results

$`H_1=H_{2k+1}={\displaystyle \frac{1}{2}}L_{}L_+`$
$`H_2=H_{2k}={\displaystyle \frac{1}{2}}L_+L_{}.`$ (40)
The Eqs. (12d) and (12e) by using the definitions (39), and the Eqs. (40) lead to

$`H_{2l1}={\displaystyle \frac{1}{2}}L_{}L_+,l=1,2,\mathrm{},k+1`$
$`H_{2l}={\displaystyle \frac{1}{2}}L_+L_{}l=1,2,\mathrm{},k.`$ (41)
Therefore, the operator $`H`$ of the B-D unitary parasupersymmetry algebra of arbitrary order $`p=2k`$ has two Hamiltonian components $`(\frac{1}{2}L_{}L_+((k+1)\text{times})\text{and}\frac{1}{2}L_+L_{}`$ $`(k\text{times}))`$ on the homogeneous manifolds $`SL(2,c)/GL(1,c)`$ with the same energy spectrum $`\frac{1}{2}E(n,m)`$. The representation basis of the parasupersymmetry algebra in terms of the quantum states $`\psi _{n,m}(\theta ,\varphi )`$ which describe the motion of the particle on the homogeneous manifolds $`SL(2,c)/GL(1,c)`$ has the following form

$$\mathrm{\Psi }(\theta ,\varphi )=\left(\begin{array}{cc}\psi _{n,m1}(\theta ,\varphi )\hfill & \\ \psi _{n,m}(\theta ,\varphi )\hfill & \\ \psi _{n,m1}(\theta ,\varphi )\hfill & \\ \psi _{n,m}(\theta ,\varphi )\hfill & \\ \mathrm{}\hfill & \\ \psi _{n,m1}(\theta ,\varphi )\hfill & \end{array}\right)_{(2k+1)\times 1}.$$
(42)
The eigenvalue equation of the parasupersymmetric Hamiltonian is

$$H\mathrm{\Psi }(\theta ,\varphi )=E(n,m)\mathrm{\Psi }(\theta ,\varphi ).$$
(43)
The representation of the parafermionic generators $`Q_1`$ and $`Q_1^{}`$ by using the representation of the Lie algebra $`gl(2,c)`$ given in the relations (37) is the relations (32) and the only difference is the fact that we must substitute the quantum states $`\psi _{n,m1}(\theta ,\varphi )`$ and $`\psi _{n,m}(\theta ,\varphi )`$ instead of $`\psi _{n,m1}(\theta )`$ and $`\psi _{n,m}(\theta )`$, respectively.

By choosing the generators of the Lie algebra $`gl(2,c)`$ in terms of three variables, given in Ref. , we have taken into account the solvable quantum models on the group manifolds $`SL(2,c)`$. In this case, like the two dimensional models, it can be also shown that the three dimensional solvable quantum models on the group manifolds $`SL(2,c)`$ described in Ref.  represent the B-D unitary parasupersymmetry algebra of arbitrary order $`p=2k`$.

## 4 B-D unitary parasupersymmetry algebra with $`p=2k`$ conserved parasupercharges

For the B-D unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$ introduced by the relations (12), it can be shown that there are $`2k1`$ independent conserved parasupercharges in addition to $`Q_1`$ . We denote these parasupercharges by $`Q_2,Q_3,\mathrm{},Q_{2k}`$. They and their Hermitian conjugates are defined as (In this section, we only follow the discussion for the second bunch of the 1D shape invariance models):

$`\left(Q_r\right)_{ll^{}}`$ $`:=`$ $`\delta _{l+1,l^{}}A(m)rll=\text{odd}`$
$`:=`$ $`\delta _{l+1,l^{}}A(m)r=l`$
$`\left(Q_r\right)_{ll^{}}`$ $`:=`$ $`\delta _{l+1,l^{}}B(m)rll=\text{even}`$ (44a)
$`:=`$ $`\delta _{l+1,l^{}}B(m)r=l`$
$`\left(Q_r^{}\right)_{ll^{}}`$ $`:=`$ $`\delta _{l,l^{}+1}A(m)rl^{}l=\text{odd}`$
$`:=`$ $`\delta _{l,l^{}+1}A(m)r=l^{}`$
$`\left(Q_r^{}\right)_{ll^{}}`$ $`:=`$ $`\delta _{l,l^{}+1}B(m)rl^{}l=\text{even}`$ (44b)
$`:=`$ $`\delta _{l,l^{}+1}B(m)r=l^{}.`$
where $`r=2,3,\mathrm{},2k`$. In this section, $`Q_1`$ is considered as the relations (25). Once more, the bosonic Hamiltonian $`H`$ is defined as the relation (25c). The algebraic relations (12c) are satisfied automatically for the definitions (44). The equations (26a) and (26b) are also obtained from the equations (12a) and (12b) for the parasupercharges (44a) and (44b), respectively. Meanwhile, the algebraic relations (12d) and (12e) lead to the equations (28) again. Therefore, the solutions (29) are obtained for the bosonic Hamiltonian $`H`$ once again. Each of $`2k`$ conserved parasupercharges along with a bosonic Hamiltonian $`H`$ satisfy separately the parasupersymmetric algebraic relations. So, we have ($`r=1,2,\mathrm{},2k`$):

$`C_0^{2k}Q_r^{2k}Q_r^{}+C_1^{2k}Q_r^{2k1}Q_r^{}Q_r\mathrm{}C_{2k}^{2k}Q_r^{}Q_r^{2k}=2Q_r^{2k1}H`$ (45a)
$`C_0^{2k}Q_{r}^{}{}_{}{}^{2k}Q_r+C_1^{2k}Q_{r}^{}{}_{}{}^{2k1}Q_rQ_r^{}\mathrm{}C_{2k}^{2k}Q_rQ_{r}^{}{}_{}{}^{2k}=2Q_{r}^{}{}_{}{}^{2k1}H`$ (45b)
$`Q_r^{2k+1}=Q_{r}^{}{}_{}{}^{2k+1}=0`$ (45c)
$`[H,Q_r]=0`$ (45d)
$`[H,Q_r^{}]=0.`$ (45e)
Similar situation is occurred in the R-S arbitrary order parasupersymmetry algebra.

In addition to the relations (45c), the parasupercharges $`Q_1,Q_2,\mathrm{},Q_{2k}`$ and their Hermitian conjugates satisfy a generalised form of the relations (45c) which are given by

$`Q_1^{a_1}Q_2^{a_2}\mathrm{}Q_{2k}^{a_{2k}}=0`$ (46a)
$`Q_{1}^{}{}_{}{}^{a_1}Q_{2}^{}{}_{}{}^{a_2}\mathrm{}Q_{2k}^{}{}_{}{}^{a_{2k}}=0,`$ (46b)
where $`a_1+a_2+\mathrm{}+a_{2k}=2k+1`$. Furthermore, there exist $`(2k+1)`$ independent bosonic constants as

$$\begin{array}{cccc}& & \left(I_1\right)_{ll^{}}:=\delta _{l,l^{}}\hfill & \\ & & \begin{array}{cccc}& & \left(I_s\right)_{ll^{}}:=\delta _{l,l^{}}ls\hfill & \\ & & :=\delta _{l,l^{}}l=s\hfill & \end{array}\}s=2,3,\mathrm{}2k+1\hfill & \end{array}$$

which commute with the bosonic Hamiltonian $`H`$:

$`[H,I_s]=0s=1,2,3,\mathrm{},2k+1.`$ (47)
It is also noticed that the commutation relations of the parasupercharges and the bosonic constants are closure, that is,

$`[I_s,Q_r]={\displaystyle \underset{l=1}{\overset{l=2k}{}}}d_lQ_ls=1,2,3,\mathrm{},2k+1\text{and}r=1,2,3,\mathrm{},2k`$ (48)
where the coefficients $`d_l`$ are constants. Similar relations exist for the Hermitian conjugate of the parasupercharges which are obtained by taking Hermitian conjugate of the relations (49). The bosonic constants $`I_s`$ and the parasupercharges $`Q_r`$ satisfy the mixed multilinear relations which are a generalisation of the relations (45a) and (45b). For example, one may introduce the B-D unitary parasupersymmetry algebra of order $`p=2`$ with two parasupercharges and three bosonic constants:

$`Q_1=\left(\begin{array}{ccc}0& A(m)& 0\\ 0& 0& B(m)\\ 0& 0& 0\end{array}\right)Q_1^{}=\left(\begin{array}{ccc}0& 0& 0\\ B(m)& 0& 0\\ 0& A(m)& 0\end{array}\right)`$ (55)
$`Q_2=\left(\begin{array}{ccc}0& A(m)& 0\\ 0& 0& B(m)\\ 0& 0& 0\end{array}\right)Q_2^{}=\left(\begin{array}{ccc}0& 0& 0\\ B(m)& 0& 0\\ 0& A(m)& 0\end{array}\right)`$ (62)
$`I_1=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)I_2=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)I_3=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)`$ (72)
$`H={\displaystyle \frac{1}{2}}\left(\begin{array}{ccc}A(m)B(m)& 0& 0\\ 0& B(m)A(m)& 0\\ 0& 0& A(m)B(m)\end{array}\right).`$ (76)
Clearly, the three generators $`H,Q_1,Q_1^{}`$ and also the three generators $`H,Q_2,Q_2^{}`$ satisfy separately the algebraic relations (45) with $`k=1`$. Moreover, we have

$`Q_1^2Q_2=Q_1Q_2^2=0[I_1,Q_1]=[I_1,Q_2]=0`$
$`[I_2,Q_1]=2Q_2[I_2,Q_2]=2Q_1`$
$`[I_3,Q_1]=Q_1Q_2[I_3,Q_2]=Q_2Q_1`$
$`[I_s,I_s^{}]=[I_s,H]=0s,s^{}=1,2,3.`$ (77)
The mixed multilinear relations which are a generalisation of the relation (45a) with $`k=1`$, are given by

$`C_0^2Q_r^2Q_r^{}^{}I_2+C_1^2Q_rQ_r^{}^{}Q_rC_2^2Q_r^{}^{}Q_r^2I_3=2Q_r^{}Hrr^{}`$
$`C_0^2Q_r^{}Q_rQ_r^{}+C_1^2Q_rQ_r^{}Q_r^{}C_2^2Q_r^{}Q_r^{}Q_rI_3=2Q_r^{}Hrr^{}`$
$`C_0^2Q_rQ_r^{}Q_r^{}I_2+C_1^2Q_r^{}Q_r^{}Q_rC_2^2Q_r^{}Q_rQ_r^{}=2Q_r^{}Hrr^{}.`$ (78)
Note that when $`r=r^{}`$, the multilinear relations are the same as the relations (45a) with $`k=1`$. The generalised mixed multilinear relations of the relations (45b) with $`k=1`$ are obtained by taking Hermitian conjugate of the relations (52).

## 5 Conclusions

In this paper an appropriate generalization of the B-D unitary parasupersymmetry algebra of arbitrary order $`p`$ is presented. It is shown that in a special approach the generalization for even arbitrary order $`p=2k`$ can be represented by the one dimensional shape invariance quantum models. In this approach, the partner Hamiltonians of the shape invariance theory are the isospectrum components of the bosonic Hamiltonian in the B-D unitary parasupersymmetric theory. Also, as we mentioned before, the 2D and 3D quantum models on the homogeneous manifolds $`SL(2,c)/GL(1,c)`$ and the group manifolds $`SL(2,c)`$ obtained from the shape invariance approach with respect to the secondary quantum number $`m`$ realize the B-D unitary parasupersymmetry algebra of arbitrary order $`p=2k`$. At the same time, the bosonic Hamiltonian $`H`$ for the 2D and 3D models has two isospectrum components. Meanwhile, the B-D unitary parasupersymmetry algebra of even arbitrary order $`p=2k`$ is generalised to a B-D unitary parasupersymmetry algebra in the presence of the $`2k`$ conserved parasupercharges with the mixed multilinear relations.

## Acknowledgment

We would like to thank the referee for fruitful comments which lead to a new section (section 4) in our paper, and the research council of Sahand University of Technology for financial support.


text:
# Corotation: its influence on the chemical abundance pattern of the Galaxy

## 1 Introduction

The radial distribution of heavy elements in the Galaxy is a signature of the Galactic history and structure. Especially important information is brought by various features in the abundance pattern.

In early papers the radial variation of abundance was approximated by a function close to the linear one in logarithmic scale (see e.g. Harris, 1981; Harris & Pilachowski, 1984 and references therein). However, in the last decade some features of the chemical distribution were revealed that do not support the simple linear model. Indeed, Simpson et al (1995) and Vílchez & Esteban (1996) derived a plateau-like (or step-like) abundance structure of the H II regions. The same pattern was found by Twarog et al (1997) for open clusters, by Smart & Rollestone (1997) for B stars, by Caputo et al. (2001) for Cepheids. In recent papers Andrievsky et al. (2002, a,b) on the basis of precise spectroscopic studies of Cepheids derived a ”bimodal” abundance gradient with a plateau beyond the solar distance and a sufficiently steep gradient inside the solar circle. Finally, Maciel & Quireza (2000) found a parabolic metallicity distribution for the planetary nebulae, with a minimum near the solar distance.

What is the mechanism of such features formation? In this paper we will inspect the idea that the bimodal abundance gradient can be a consequence of the corotation resonance.The underlying idea is that at corotation the star formation rate presents a minimum, and consequently, the metallicity accumulated after a few billion years at this radius can also have minimum. Bearing in mind that since Marochnik et al (1972) and Crézé & Mennessier (1973) a lot of studies showed that the corotation is situated very close to the Sun (Nelson & Matsuda, 1976; Mishurov et al, 1979; Mishurov et al, 1997, Amaral & Lépine, 1997; von Linden et al., 1998; Leicht & Vasisht, 1998; Mishurov & Zenina, 1999a,b; Lépine et al, 2001; Fernandez et al., 2001 etc) it would be interesting to analyze the influence of the corotation on the abundance pattern.

On the other hand, according to Weinberg (1994), Englmaier & Gerhard (1999), Dehnen (2000) and others the corotation is situated very close to the galactic center at a distance $``$ 4 kpc. So, investigating the chemical abundance structure of the Galactic disk, we can hope to make a choice between the above models.

## 2 Basic idea and equations

Standard equations governing the chemical evolution of the galactic disk take into account both the creation of heavy elements and mass transfer between gaseous and stellar components (Tinsley, 1980). They are essentially the equations of continuity and convey the mass conservation law for a mixture: gas + stars \+ heavy elements (e.g. Lacey & Fall, 1989, hereafter LF; Portinari & Chiosi, 1999, hereafter PC). By means of these equations researchers usually try to simultaneously explain both the radial behavior of metal abundance and gaseous and stellar distributions. However the corresponding equations must be completed by the dynamical ones for instance by the Euler gasdynamic equation, the Boltzman equation for the stellar component, the Poisson equation for gravitation field etc. Only such an approach would enable us to derive all radial profiles in a self-consistent manner. Of course, this way is very complicated and most authors prefer to exclude the dynamical part of the problem (a rare exception is the paper of Martinet & Friedly, 1997) and consider the task in the framework of the standard scheme, specifying some model laws for the star formation rate (SFR), the infall rate of matter onto the Galactic disk, and the radial flows within the disk, so as to explain the observed radial distribution for chemical abundance, star density and interstellar gas density.

To make more prominent the effects of the spiral arms and of the corotation on the abundance pattern which are not masked by other factors, in the present paper we shall only concentrate our attention on the chemical enrichment. For this task we adopt the stationary model of Haywood et al. (1997) for Galactic evolution. According to these authors, gas density and infall rate are constant in time. This enables us to consider the chemical part of the problem separately for various laws of infall rate, radial flow and star formation rate (SFR) 444Of course, we abstain from explaining the stellar density distribution. In any case it is impossible to adequately derive the complicated structure of the stellar disk in the framework of models which do not take into account dynamical processes like Galactic contraction as a whole, influence of the bar and of the spiral structure etc.. The corresponding equation we write in a form close to that of LF:

$`\mu _g{\displaystyle \frac{Z}{t}}`$ $`=`$ $`P_Z\psi +f(Z_fZ)\mu _gu{\displaystyle \frac{Z}{R}}+`$ (1)
$`{\displaystyle \frac{1}{R}}{\displaystyle \frac{}{R}}(\mu _gD{\displaystyle \frac{Z}{R}}),`$
where $`\mu _g`$ is the surface density of interstellar gas, $`Z`$ is the abundance of heavy elements, $`\psi `$ is the star formation rate, $`P_Z`$ is the mass fraction of a star ejected as newly synthesized heavy elements (Tinsley, 1980), $`u`$ is the macroscopic radial velocity of the gas, $`f`$ is the infall rate of matter onto the Galactic disk per unit area, $`Z_f`$ is the abundance of infalling gas, $`t`$ is time, $`R`$ is galactocentric distance. In addition to the corresponding equation of LF we include the effect of diffusion (in the last term in the right hand side of Eq (1), $`D`$ is the coefficient of heavy elements diffusion). The matter is that, due to stellar winds, supernovae explosions, Galactic fountains, etc., the interstellar gas is highly turbulent. This leads to gas mixing, which tends to smooth out the abundance gradients. For a phenomenological description of this process we introduced the diffusion term into Eq. (1).

From previous studies it is known that the influence of the SFR radial law happens to be very important on the Galactic abundance structure. Numerous approximation formulas were proposed for it. As a rule the SFR $`\psi `$ is assumed to be proportional to some power of gas surface density $`\mu _g`$. However this cannot be recognized as completely satisfactory. Indeed, it is well known that supernovae of all types - sources of heavy elements - are strongly concentrated in spiral arms (Bartunov et al, 1994). Hence synthesis of heavy elements in a given volume of matter takes only place when the volume occurs to be inside an arm. To take into account the effects of spiral arms and corotation on chemical enrichment we represent the SFR in a manner similar to Oort (1974), Wyse & Silk (1989) and PC:

$$\psi =\beta \mu _g^k|\mathrm{\Omega }(R)\mathrm{\Omega }_P|\mathrm{\Theta }(\nu )$$
(2)
where $`\mathrm{\Omega }(R)`$ is the angular rotation velocity of the Galactic disk, $`\mathrm{\Omega }_P`$ is the rotation velocity of the Galactic density waves responsible for spiral arms (recall that the spiral pattern rotates as a solid body, i.e. $`\mathrm{\Omega }_P=const`$, whereas the Galactic disk rotates differentially, i.e. $`\mathrm{\Omega }`$ is a function of $`R`$, the distance $`R_C`$ where both the velocities coincide, i.e. $`\mathrm{\Omega }(R_C)=\mathrm{\Omega }_P`$, is called the corotation radius), $`\beta `$ and $`k`$ are some constants (we select the product $`P_Z\beta `$ so that $`Z(R_{})=Z_{},`$555According to various authors (see e.g. PC and references therein) the solar metallicity is greater than that of its surroundings. However Andrievsky et al. (2002a) showed that the Sun has a metallicity which is typical of its vicinity. and consider two cases $`k=1`$ and $`k=1.4`$ Kennicutt, 1998), $`\mathrm{\Theta }`$ is a cutoff factor (see below).

The factor $`|\mathrm{\Omega }\mathrm{\Omega }_P|`$ in Eq (2) means a frequency of crossing spiral arms by some element of matter of Galactic disk. It appears due to the fact that the enrichment by heavy elements takes only place when the element of Galactic matter occurs to be within a spiral arms. An equivalent interpretation is that the SFR depends on the rate at which the gas penetrates the spiral arms.

Wyse & Silk (1989) and PC considered the case when the corotation resonance is situated at the very end of the Galaxy (model of Lin et al, 1969). So, they simply neglected by the value $`\mathrm{\Omega }_P`$ in Eq (2).

In our paper we do not neglect $`\mathrm{\Omega }_P`$ and consider several cases: with the corotation close to the position of the Sun ($`R_CR_{}=8.5kpc`$); corotation near the galactic center ($`R_C=4kpc`$) and far corotation ($`R_C=14kpc`$).

The cutoff factor $`\mathrm{\Theta }`$ is required because the self-sustained Galactic density waves are known to exist only in the so-called wave zone between inner Lindblad resonance $`R_{in}`$ where $`\nu (R_{in})=1`$ and the outer Lindblad resonance $`R_{out}`$ where $`\nu (R_{out})=+1`$ (Lin et al., 1969). So, we define $`\mathrm{\Theta }=1`$ for $`|\nu |<1`$ and $`\mathrm{\Theta }=0`$ for $`|\nu |>1`$. Here, $`\nu =m(\mathrm{\Omega }_P\mathrm{\Omega })/\kappa `$ is the dimensionless wave frequency, $`m`$ is the number of arms, $`\kappa =2\mathrm{\Omega }(1+Rd\mathrm{\Omega }/dR/2\mathrm{\Omega })^{1/2}`$ is the epicyclic frequency.

The expression (2) for SFR means that we only consider an ordered star formation process triggered by spiral arms. Besides this, in the Galaxy, some sporadic (or spontaneous) star formation may take place. Hence, in general $`\psi `$ should include the corresponding term. But we restrict ourselves to the representation (2) for the SFR; in other words, we consider that the spontaneous SFR is negligible compared to that induced by the spiral arms.

From the above equation, for a given set of input functions, we can compute the time history of chemical evolution of the Galactic disk. However, in this paper we only consider the simplest asymptotic solution. Indeed, due to diffusion term in Eq (1), we expect that $`Z/t0`$ for $`t\mathrm{}`$. A typical diffusion time for the achievement of the asymptotic state is $`t_DL^2/D`$ where $`L`$ is a typical length, say $`L1kpc`$. Let us estimate $`t_D`$ and show that we may use the asymptotic regime. For our phenomenological theory we consider interstellar gas clouds as particles. So, we use the gas kinetic theory for the estimation of the diffusion coefficient:

$$\alpha \mu _gD=\frac{1}{3\pi \sqrt{2}}\frac{m_C}{a^2}HV_T$$

where $`m_C`$ is the mass of a cloud, $`a`$ is its typical radius, $`H`$ is the thickness of gas layer ($`H=0.13kpc`$) and $`V_T`$ is the thermal velocity of a cloud ($`V_T=6.6kms^1`$, Stark, 1984). Following Falgarone & Puget (1986; see also Elmegreen, 1987) we adopt $`m_C/a^2=100M_{}/pc^2`$. Notice that for constant $`H`$ and $`V_T`$ the value $`\alpha `$ is constant as well. For a typical $`\mu _g10M_{}/pc^2`$ $`t_D1.6Gyr`$. So, in the major part of the Galaxy we may use the above asymptotic approximation although in a more refined theory departures from it should be investigated.

As boundary conditions we adopt: at the ”end” of the galactic disk, say at $`R=R_G=25kpc`$ $`Z=Z_f`$ ($`Z_f=0.17Z_{}`$, $`Z_{}=0.02`$ Tinsley, 1980)666Results depend weakly on the exact value $`R_G`$ in the range 20 - 30 kpc and at the galactic center the abundance must be finite, therefore $`Z/R=0`$ for $`R=0`$ (in what follows we shall not be interested in the abundance structure of the central part of the Galaxy since chemical composition there is masked by very complicated processes due to the bar, the bulge etc.).

Let us now discuss the functions entering Eq (1). For the infall rate we adopt:

$$f=f_{}exp(\frac{RR_{}}{\mathrm{\Delta }}),$$

where $`f_{}=3.5M_{}pc^2Gyr^1`$ (Haywood et al, 1997).

For the surface gas density we use a simple representation

$$\mu _g=\mu _0exp(\frac{R}{d}\frac{d}{R})$$

which imitates the hole in the central part of the disk with maximum $`\mu _g`$ at $`R=d`$ and is close to observed radial distribution of interstellar gas (Dame, 1993) for $`d=4kpc`$, the dependence being normalized as to give $`\mu _g(R_{})=10M_{}pc^2`$ (Haywood et al, 1997).

In all computations the rotation curve of Allen & Santillan (1994) with $`R_{}=8.5kpc`$ and $`V_{rot}(R_{})=220kms^1`$ was used.

## 3 Results and Discussion

First of all we consider a model with corotation close to the Sun ($`R_C=R_{})`$ and radial inflow $`u=0`$. Our experiments with $`\mathrm{\Delta }`$ in the range 4 - 5 kpc (close to usually adopted) give similar results. In Fig 1 by solid lines are shown the abundance patterns for $`k=1.4`$ (thick line) and $`k=1`$ (thin line) for $`\mathrm{\Delta }=5kpc`$. To feel the effects of corotation, diffusion and cuttoff $`\mathrm{\Theta }`$-term on chemical enrichment we give the stationary solution of Eq (1) neglecting the diffusion ($`k=1.4`$ \- dotted line; $`k=1`$ \- dashed line). Note abrupt decrease of abundance in the last two cases beyond the Lindblad resonances due to cuttoff $`\mathrm{\Theta }`$-term and the gap in the vicinity of the corotation radius. As it was expected those sharp features are smoothed out due to diffusion.

In Fig 2 thick line is the same as in Fig 1 but the dashed line represent the solution of Eq (1) for the corotation slightly closer to the galactic center, namely $`R_C=8kpc`$ and $`k=1.4`$. Both the curves are superimposed on the observed abundance derived from Cepheids by Andrievsky et al. (2002 b). From this Figure it can be seen that the model with the corotation close to the Sun describes rather well the bimodal abundance gradient derived by Andrievsky et al (2002 b), with a plateau beyond the solar circle and a steep gradient inside of it.

In Fig 3 is shown the influence of radial gas flow: i) with constant velocity $`u`$: $`u=+0.5kms^1`$ and $`u=0.5kms^1`$ and ii) with outflow from the corotation with the same amplitude $`0.5kms^1`$ (Lépine et al, 2001). Thick line is the same as in Fig 1.

From the above Fig 1-3 one can see that it is possible to explain the abundance features like bimodal abundance structure of Andrievsky et al (2002) by means of combined influence of spiral arms, closeness of corotation to the solar galactocentric distance and radial gas flow on chemical enrichment.

In Figs 4,5 by solid lines the results are given for models with the corotation near the galactic center ($`R_C=4kpc`$, Fig 4) and far corotation ($`R_C=14kpc`$, Fig 5). In both pictures dotted lines represent stationary solutions obtained by neglecting the diffusion term. From these picture it is obvious that these models for Galactic spiral waves are not consistent with the observed abundance pattern.

We are grateful to S.M.Andrievsky for very useful discussion of the problem.

Figure captions

Fig. 1 - The theoretical radial abundance distribution for the corotation radius close to the Sun position, i.e. $`R_C=R_{}=8.5kpc`$. Solid lines: thick - $`k=1.4`$, thin - $`k=1`$. Dotted and dashed lines describe the stationary solution of Eq (1) neglecting the diffusion term for $`k=1.4`$ and $`k=1`$ correspondingly. Radial inflow $`u=0`$.

Fig. 2 - The theoretical abundance gradient (solid and dashed lines) superimposed on the observational data derived by Andrievsky et al (2002 b) for Cepheids (crosses). The thick line is the same as in Fig 1, the dashed one corresponds to the corotation radius slightly closer to the Galactic center then in the previous case, namely $`R_C=8kpc`$. Radial inflow $`u=0`$.

Fig. 3 - The theoretical abundance gradient taking into account the radial inflow: i) $`u=+0.5kms^1`$ \- dashed line; ii) $`u=0.5kms^1`$ \- dotted line; iii) outflow from the corotation (Lépine et al, 2001) with the amplitude $`0.5kms^1`$ \- dashed - dotted line, the corotation radius being situated at the Sun galactocentric distance. The thick solid line is the same as in Fig. 1.

Fig. 4 - The theoretical abundance gradient for the model with the corotation near the galactic center, $`R_C=4kpc`$ and $`k=1.4`$ \- solid line. Dashed line represents the stationary solution neglecting the diffusion.

Fig. 5 - The same as Fig. 4 but for the far corotation model, $`R_C=14kpc`$.