format: text

text:
# Phase equivalent chains of Darboux transformations in scattering theory

## I Introduction

The construction of potentials sharing the same scattering data, called in the literature phase equivalent potentials, is one of the most interesting applications of supersymmetric (SUSY) quantum mechanics  to the inverse scattering problem (see e.g. ). As it is known, the SUSY approach, when restricted to the derivation of new exactly solvable quantum problems, is basically equivalent to the Darboux transformations method (see e.g.). Therefore we use SUSY and Darboux transformations as synonyms.

To our knowledge, the simplest way of constructing phase equivalent potentials by using Darboux transformations, is due to Sukumar . It consists in a successive application of two transformations with equal factorization constants. The method has been considerably developed, up to an arbitrary modification of the discrete spectrum . It has been extended to coupled channels  and used to derive deep and shallow phase equivalent potentials describing nucleon-nucleon  and nucleus-nucleus  collisions. In these applications both deep and shallow potentials have been constructed such as to share the same scattering data at all energies. Shallow potentials possess only ”physical” states while the deep potentials, in addition to the physical states also have ”unphysical” bound states which simulate the effect of the Pauli principle (for details see e.g. ).

Based on Darboux transformations, here we propose a new way to derive deep and shallow phase equivalent potentials. Our aim is to improve over the practical procedures of Refs.  and . By replacing a chain of Darboux transformations by an equivalent $`N`$th order transformation and using a Wronskian formula recently obtained in  (see also ) we use a determinant approach to construct phase equivalent potentials. When the transformation functions are sufficiently simple (e.g. they are solutions of the free particle Schrödinger equation) this approach happens to be more efficient than the one of previous works . The advantages of our method will be clearly seen in its application to the nucleon-nucleon scattering potentials considered here. In particular the potential tail will be cured from the unwanted oscillations obtained in Ref. .

As it is known (see e.g. ) there are both deep and shallow potentials corresponding to the same value of the partial wave $`l`$. The value of $`l`$ defines the long distance behaviour of the potential . The short distance behaviour

$$V(x)\frac{\nu (\nu +1)}{x^2},x0$$
(1)
depends on the parameter $`\nu `$ which we shall call the strength of the singularity at short distances or simply the singularity strength. Following Ref.  we generally assume that $`\nu `$ is different from $`l`$. In this work we consider only $`l`$-preserving transformations and restrict our discussion to the $`l=0`$ case, although most of the derived formulae are valid in the general case. In a subsequent paper  we shall consider $`l`$-changing transformations which will allow to introduce higher partial waves.

We assume that no Coulomb interaction is present in the potential i. e. we are dealing with the Schrödinger equation on the half axis, $`x[0,\mathrm{})`$, with the potential $`V(x)`$ satisfying the condition (1) and at long distances going to zero faster than $`x^2`$. For every eigenvalue $`E`$ the Schrödinger equation with such potentials has two linearly independent solutions with different behavior at short distances, one which behaves as $`x^\nu `$, called regular solution, and the other as $`x^{\nu 1}`$, called singular (or irregular) solution. In general Darboux transformations, which use irregular solutions as transformation functions, do not preserve the property of a solution to remain regular at the origin. As a consequence, the change in the Jost function and the phase shift becomes difficult to control since one has to know both the Jost solution of the initial equation and the derivative of this solution at the origin. As far as we know this problem has not been yet solved in the general context. But we have found that there exist particular chains of transformations where the derivative of the Jost solution of the initial equation does not appear in the expression of the Jost function for the transformed Hamiltonian. In this case one can get a closed expression for the transformed Jost function in terms of the initial Jost function and a rational function of the momentum $`k`$ only. Potentials having Jost functions related in such a way are known in the literature as potentials of the Bargmann class (see e.g. ).

Once we know how the Jost function is changed by a chain of Darboux transformations we can select chains that produce potentials with the same $`S`$-matrix and hence with the same phase shift. Such potentials are called in the literature phase equivalent although a more appropriate term could be isophase potentials since they give exactly the same phase shift at all energies. Here we shall use both terms. So, to get phase equivalent potentials we propose to use phase equivalent chains instead of an integral form of Darboux transformations used in . We start from the same exactly solvable reference potential and the action of each chain will result in potentials always having the same $`S`$-matrix. Finally we note that this approach is basically equivalent to the inverse scattering method where the kernel of the Fredholm operator involved in the Gelfand-Levitan-Marchenko equation is degenerate  or to a supersymmetric transformation , but it is more direct and has more computational advantages.

The paper is organized as follows. In the next section we recall the Darboux transformation method in the context of scattering theory. In Sec. III we introduce special chains for which a closed expression for the transformed Jost function is derived. Then we select chains which lead to families of phase equivalent (or isophase) potentials. Sec. IV is devoted to an application of our method to the derivation of deep and shallow phase equivalent potentials for the nucleon-nucleon scattering. Conclusions are drawn in the last section. Details of important analytic calculations are gathered in the Appendix.

## II Darboux transformation method

The method of Darboux transformations, well known in soliton theory , finds now more and more applications in other fields of theoretical and mathematical physics . In this section we recall the definition of Darboux transformations and give their main properties necessary in subsequent sections.

Essentially, the method consists in getting solutions $`\phi `$ of one Schrödinger equation

$$h_1\phi =E\phi ,h_1=\frac{d^2}{dx^2}+V_1(x),$$
(2)
when solutions $`\psi `$ of another equation

$$h_0\psi =E\psi ,h_0=\frac{d^2}{dx^2}+V_0(x),$$
(3)
are known. This is achieved by acting with a differential operator $`L`$ of the form

$$\phi =L\psi ,L=d/dx+w(x),$$
(4)
where the real function $`w(x)`$, called superpotential, is defined as the logarithmic derivative of a known solution of (3) denoted by $`u`$ in the following. One has

$$w=u^{}(x)/u(x),h_0u=\alpha u,$$
(5)
with $`\alpha E_0`$, where $`E_0`$ is the ground state energy of $`h_0`$ if it has a discrete spectrum or the lower bound of the continuum spectrum otherwise. The function $`u`$ is called transformation or factorization function and $`\alpha `$ its factorization constant or factorization energy. The potential $`V_1`$ is defined in terms of the superpotential $`w`$ as

$$V_1(x)=V_0(x)2w^{}(x).$$
(6)
The operators $`L`$ and $`L^+=d/dx+w(x)`$ establish a one-to-one correspondence between two dimensional subspaces of solutions of Eqs. (2) and (3) with a given $`E\alpha `$. When $`E=\alpha `$ the transformed equation (2) has $`v=1/u`$ and $`\stackrel{~}{v}`$ as linearly independent solutions. The latter can be found by using the Wronskian property $`W(v,\stackrel{~}{v})=1`$ which gives

$$\stackrel{~}{v}=v\underset{x_0}{\overset{x}{}}v^2(x)𝑑x,$$
(7)
where $`x_0`$ is an arbitrary constant. Hence, the knowledge of all solutions of the initial equation provides the knowledge of all solutions of the transformed equation and in particular of all physical solutions in the usual quantum mechanical sense.

Since the above procedure does not depend on a particular choice of the potential $`V_0`$ the transformed Hamiltonian $`h_1`$ can play the role of the initial Hamiltonian for the next transformation step. In such a way one gets a chain of exactly solvable Hamiltonians $`h_0,h_1,\mathrm{}h_N`$ with the potentials $`V_0,V_1,\mathrm{}V_N`$. To avoid any confusion we mention that everywhere, except for especially mentioned cases, we shall use subscripts to distinguish between quantities related to different Hamiltonians, $`h_0`$, $`h_1`$,…and shall omit them when discussing general properties regarding all Hamiltonians.

A first order Darboux transformation operator $`L_{j,j+1}`$ as defined by (4) relates two Hamiltonians $`h_j`$ and $`h_{j+1}`$. If one is not interested in the intermediate Hamiltonians $`h_1`$, …, $`h_{N1}`$ and all factorization energies are chosen to be different from each other the whole chain may be replaced by a single transformation given by an $`N`$th order transformation operator denoted by $`L^{(N)}`$, defined as a succession of $`N`$ first order transformations. A compact representation of this operator is given by

$$\psi _N(x,k)=L^{(N)}\psi _0(x,k)=W(u_1,\mathrm{},u_N,\psi _0(x,k))W^1(u_1,\mathrm{},u_N),$$
(8)
where $`\psi _0(x,k)`$ is a solution of Eq. (3) corresponding to the energy $`E=k^2`$ and $`\psi _N(x,k)`$ satisfies

$$h_N\psi _N(x,k)=E\psi _N(x,k),E=k^2.$$
(9)
The transformation functions, although labeled by a subscript, are eigenfunctions of the initial Hamiltonian

$$h_0u_j(x)=\alpha _ju_j(x),\alpha _j=a_j^2.$$
(10)
They should be chosen such as the Wronskian $`W(u_1,\mathrm{},u_N)`$ is nodeless in the interval $`(0,\mathrm{})`$. This condition guarantees the absence of singularities in the potential

$$V_N=V_02\frac{d^2}{dx^2}\mathrm{log}W(u_1,\mathrm{},u_N),$$
(11)
defining the Hamiltonian $`h_N`$ of (9), inside this interval. The formula (8) is valid for any $`E=k^2`$ except for $`k=ia_j`$ ($`j=1,\mathrm{}N`$). For these values of $`k`$ the corresponding solutions are

$$\psi _N(x,ia_j)=W^{(j)}(u_1,\mathrm{},u_N)W^1(u_1,\mathrm{},u_N),j=1,\mathrm{},N,$$
(12)
where $`W^{(j)}(u_1,\mathrm{},u_N)`$ is the $`(N1)`$th order Wronskian constructed from $`u_1,\mathrm{},u_N`$ except for $`u_j`$.

In general, a chain of $`N`$ transformations can have a number of elements corresponding to transformation functions with identical factorization constants. In the simplest $`N=2`$ case where the first transformation is implemented by the function $`u_1=u`$ and the second by $`v=1/u`$ one gets back the initial potential $`V_0`$. This means that the operator $`L^+`$ realizes the transformation in the opposite direction. More generally, the second transformation can be made with the linear combination $`v_1=cv+\stackrel{~}{v}`$ where $`\stackrel{~}{v}`$ is given by (7). The action of such a chain leads to an eigenfunction of $`h_2`$ of the form

$$\psi _2(x,k)=\psi _0(x,k)\frac{u(x)}{W_2(x)}\underset{x_0}{\overset{x}{}}u(t)\psi _0(t,k)𝑑t,$$
(13)
where

$$W_2(x)=c+\underset{x_0}{\overset{x}{}}u^2(t)𝑑t.$$
(14)
Choosing $`x_0=0`$ (or $`x_0=\mathrm{}`$) together with the condition $`c0`$ (or $`c0`$) guarantees a nodeless $`W_2(x)`$ inside the interval $`(0,\mathrm{})`$. The integral representation (13) based on two subsequent transformations with the same factorization constants was used in Refs.  for getting phase equivalent potentials. We now describe an alternative approach which has the advantage that the above $`N=2`$ chain of transformations can be included in a larger chain as a subchain.

The potential $`V_2(x)`$ is a particular case of (11) where $`W(u_1,\mathrm{},u_N)`$ is given by $`W_2(x)`$ of (14). It has been shown  that the function $`W_2`$ is nothing else but the limiting value of the ratio $`W(u_1,u_2)/(\alpha _1\alpha _2)`$ when $`u_2u_1`$. This fact, together with a property of the Wronskians established in , which allows to reduce the order of the determinants involved in calculations, implies that one can use the same formulae (8) and (11) also in the case of pairwise identical factorization constants but the Wronskians should be calculated with the help of another representation. If the Wronskian in (8) is of even order, $`N=2m`$, it becomes proportional to an $`m\times m`$ determinant

$$W(u_1,\mathrm{}u_{2m})=\underset{lj}{\overset{2m}{}}(\alpha _l\alpha _j)\times \left|\begin{array}{ccccc}W_{\alpha _1,\alpha _2}& W_{\alpha _1,\alpha _4}& W_{\alpha _1,\alpha _6}& \mathrm{}& W_{\alpha _1,\alpha _N}\\ W_{\alpha _3,\alpha _2}& W_{\alpha _3,\alpha _4}& W_{\alpha _3,\alpha _6}& \mathrm{}& W_{\alpha _3,\alpha _N}\\ W_{\alpha _5,\alpha _2}& W_{\alpha _5,\alpha _4}& W_{\alpha _5,\alpha _6}& \mathrm{}& W_{\alpha _5,\alpha _N}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ W_{\alpha _{N1},\alpha _2}& W_{\alpha _{N1},\alpha _4}& W_{\alpha _{N1},\alpha _6}& \mathrm{}& W_{\alpha _{N1},\alpha _N}\end{array}\right|.$$
(15)
For an odd order Wronskian, $`N=2m+1`$, one has

$$\begin{array}{c}W(\psi ,u_1,\mathrm{}u_{2m})=\underset{lj}{\overset{2m}{}}(\alpha _l\alpha _j)\underset{k=1}{\overset{m}{}}(\alpha _{2k}E)\hfill \end{array}$$

$$\begin{array}{c}\\ \times \left|\begin{array}{cccccc}\psi & W_{E,\alpha _2}& W_{E,\alpha _4}& W_{E,\alpha _6}& \mathrm{}& W_{E,\alpha _N}\\ u_1& W_{\alpha _1,\alpha _2}& W_{\alpha _1,\alpha _4}& W_{\alpha _1,\alpha _6}& \mathrm{}& W_{\alpha _1,\alpha _N}\\ u_3& W_{\alpha _3,\alpha _2}& W_{\alpha _3,\alpha _4}& W_{\alpha _3,\alpha _6}& \mathrm{}& W_{\alpha _3,\alpha _N}\\ u_5& W_{\alpha _5,\alpha _2}& W_{\alpha _5,\alpha _4}& W_{\alpha _5,\alpha _6}& \mathrm{}& W_{\alpha _5,\alpha _N}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ u_{N1}& W_{\alpha _{N1},\alpha _2}& W_{\alpha _{N1},\alpha _4}& W_{\alpha _{N1},\alpha _6}& \mathrm{}& W_{\alpha _{N1},\alpha _N}\end{array}\right|,\hfill \end{array}$$
(16)
where

$$W_{\alpha _k,\alpha _l}=\frac{W(u_k,u_l)}{\alpha _k\alpha _l},W_{E,\alpha _l}=\frac{W(\psi ,u_l)}{E\alpha _l}.$$
(17)
The product $`{\displaystyle \underset{lj}{\overset{2m}{}}}(\alpha _l\alpha _j)`$ introduced above contains only the factors $`(\alpha _l\alpha _j)`$ appearing in $`W_{\alpha _k,\alpha _l}`$. For identical factorization constants the function $`W_{\alpha _k,\alpha _l}`$ has to be calculated from the formula (14).

## III PHASE EQUIVALENT CHAINS OF TRANSFORMATIONS

Our method is based on a choice of transformation functions which allows to transform the Jost function through a special chain of Darboux transformations. We follow the usual definitions of the scattering theory which can be found for example in Refs. .

### A Jost function for a special chain of transformations

Let us choose the following set of eigenfunctions of the Hamiltonian $`h_0`$

$$v_1(x),\mathrm{},v_\nu (x),u_{\nu +1}(x),v_{\nu +1}(x),\mathrm{},u_n(x),v_n(x),\nu 0,$$
(18)
$$h_0u_j(x)=a_j^2u_j(x),h_0v_j(x)=b_j^2v_j(x),\mathrm{Im}a_j=\mathrm{Im}b_j=0,$$
(19)
as transformation functions for the Darboux transformation of order $`N=2n+\nu `$. The role of the first $`\nu `$ functions $`v_j(x)`$ which are absent in (18) for $`\nu =0`$ will be clear below. Let the $`a`$’s and the $`b`$’s be different from each other and choose them such as the corresponding factorization energies be smaller than the ground state energy of $`h_0`$ when it has a discrete spectrum or less than the lower bound of the continuous spectrum otherwise.

We distinguish between the functions $`u`$ and $`v`$ from their behaviour at the origin. The functions $`v`$ are regular, $`v_j(0)=0`$, and hence are uniquely defined up to a constant factor which is not essential for our purpose. They form a regular family and are exponentially increasing at large $`x`$ and therefore it is natural to restrict all $`b`$’s to be positive. The functions $`u_j`$ are irregular at the origin, $`u_j(0)0`$, and from a singular family. They are defined as linear combinations of two Jost solutions $`f(x,\pm k)`$:

$$u_j(x)=A_jf(x,ia_j)+B_jf(x,ia_j).$$
(20)
When $`A_j0`$ the functions $`u_j`$ increase exponentially at large $`x`$ provided $`a_j>0`$. At fixed $`a_j`$ they form a one-parameter ($`A_j/B_j`$) family. The case $`A_j=0`$, $`B_j0`$, with $`a_j>0`$ is equivalent to the case $`A_j0`$, $`B_j=0`$, with $`a_j<0`$. They correspond to a decreasing exponential function (20), uniquely defined (up to a constant factor). In this case in order to stress the asymptotical decrease of this function we allow $`a_j`$ to be negative and set, $`B_j=0`$, $`A_j0`$. So, for $`a_j>0`$ the functions (20) have an increasing asymptotic behavior and for $`a_j<0`$ their asymptotic behavior is opposite. It is useful to note that if $`A_j/B_j=F_0(ia_j)/F_0(ia_j)`$ the function (20) becomes regular. In the chain of functions (18) we shall allow to move some functions from the singular family to the regular family. This preserves the number of functions $`N`$ but changes the value of $`\nu `$. The effect of such a modification of a chain will be discussed below. The particular case where the initial potential is zero and $`A_j=B_j`$ has been considered earlier in Ref. .

In the Appendix we show that the chain of transformations with the transformation functions (18) increase the singularity parameter of the initial potential by $`\nu `$ units, $`\nu _0\nu _N=\nu _0+\nu `$, and the transformed Jost function $`F_N(k)`$ is related to the initial function $`F_0(k)`$ by

$$F_N(k)=F_0(k)\underset{j=\mathrm{\hspace{0.17em}1}}{\overset{\nu }{}}\frac{k}{k+ib_j}\underset{j=\nu +1}{\overset{n}{}}\frac{kia_j}{k+ib_j}.$$
(21)
Since a Jost function is analytic in the upper half of the complex $`k`$-plane (see ) all $`b`$’s must be positive whereas the $`a`$’s can have any sign, which agrees with our sign convention for $`a_j`$ and $`b_j`$ and every positive $`a_j`$ corresponds to a discrete level $`E_j=a_j^2`$ of $`V_N`$. A remarkable property of the relation (21) is that it is independent of the ratio $`A_j/B_j`$ with $`A_j`$ and $`B_j`$ from the definition (20). But the potential $`V_N`$ does depend on this ratio. Hence, keeping $`F_N(k)`$ fixed but varying $`A_j/B_j`$ such as $`u_j`$ remains irregular we can obtain a family of isospectral potentials. Nevertheless the potential $`V_N`$ belongs to this isospectral family only if $`A_j0`$ and the function $`u_j(x)`$ increases at large $`x`$. The choice $`A_j=0`$ is equivalent to changing the sign of $`a_j`$ and leads to losing the level $`E=a_j^2`$ in the potential $`V_N`$. This property does not have an analogue in the usual scheme of SUSY transformations .

### B Phase equivalent chains

Now we turn to the construction of phase equivalent or isophase potentials. Such potentials have identical $`S`$-matrices which are defined in terms of the Jost function by

$$S(k)=e^{2i\delta (k)}=\frac{F(k)}{F(k)}.$$
(22)
Therefore if for a real $`k`$ two Jost functions differ from each other by a real factor they correspond to isophase potentials. In particular, these potentials can have equal Jost functions. Taking into consideration the discussion at the end of the previous subsection we see that by changing the ratio $`A_j/B_j`$, which corresponds to a discrete level of the potential $`V_N`$, we can get not only isospectral but also isophase potentials. As it is clear from the Appendix, the potentials are isospectral provided the function $`u_j`$ remains singular. If it moves from the singular to the regular family, the transformed potential does not belong to the family of isospectral potentials any more but as we show below it remains in the family of isophase potentials.

Let us suppose for example that the function $`u_1(x)`$ becomes regular from being singular. Then for a regular $`u_1`$ the Jost function of the potential $`V_N(x)`$ corresponds to (21) with $`\nu =0`$ and for a singular $`u_1`$ it should be calculated by the same formula where $`\nu =2`$ and $`b_1`$ is replaced by $`a_1`$. From here it follows that the ratio of these two Jost functions is equal to $`k^2/(k^2+a_1^2)`$ which is real for a real $`k`$ and hence these two chains of transformations are isophase. It is also clear that if $`\nu `$ functions from the singular family move to the regular family (i.e. $`\nu `$ zeros of the Jost function are transformed into its poles) the potential $`V_N`$ loses $`\nu `$ discrete levels, and acquires an additional singularity strength equal to $`2\nu `$, but remains in the family of isophase potentials. Furthermore, this procedure is completely reversible. If one has a singular potential with $`\nu 2`$ then one of the poles of its Jost function can be transformed into a zero in the upper half of the complex $`k`$-plane which leads to the creation of a discrete level and a decrease of the singularity strength $`\nu \nu 2`$. This procedure can be repeated either up to the point where the potential loses all discrete levels and all free parameters, becoming uniquely defined by its phase shift (see e.g. ) or up to the point where $`\nu `$ becomes equal to zero or one. To avoid additional singularities inside the interval $`(0,\mathrm{})`$ one has to take these poles such as the additional level be located below the existing discrete levels. The latter restriction is well known in the SUSY approach .

The above method of deriving isophase potentials is not unique. As it is known one can create a new energy level independent of the existence of poles of the Jost function, at any desirable position even in the middle of the continuous spectrum, and preserving the $`S`$-matrix. For this purpose one can use a two-step Darboux transformation with the same factorization energy . After a new level has been created the singularity of the potential loses two units $`\nu \nu 2`$ . This procedure can be repeated up to a limit where the singularity parameter becomes equal to either zero or one. For a regular potential or a potential with $`\nu =1`$ no possibility exists to create a level by a phase preserving procedure. In this context the smallest possible value is $`\nu =2`$ when the potential with an additional discrete level becomes regular after transformation. According to Ref.  a chain of two transformations with coinciding factorization constants is equivalent to considering the limit $`a_{j+1}a_j`$ and using the Wronskian formulae (15) and (16). Without loosing generality let us consider $`a_2a_1`$ for our system of transformation functions (18) with $`\nu =0`$. In this limit the Jost function (21) where we set $`\nu =0`$

$$F_N(k)=F_0(k)\frac{k^2+a_1^2}{(k+ib_1)(k+ib_2)}\underset{j=\mathrm{\hspace{0.17em}3}}{\overset{n}{}}\frac{kia_j}{k+ib_j}.$$
(23)
has no singularity. The corresponding potential $`V_N`$ can be considered as a result of the application of the above two-step phase preserving SUSY transformation on a singular potential $`V_{N2}`$ generated by the system of transformation functions $`v_1,v_2,u_3,v_3,\mathrm{},u_n,v_n`$ whose Jost function $`F_{N2}`$ is given by the same formula (21) with $`\nu =2`$. Finally we conclude that $`V_{N2}`$ is isophase with $`V_N`$ since the ratio $`F_N/F_{N2}=(k^2+a_1^2)/k^2`$ is real.

### C The $`S`$-matrix poles

Once the Jost function $`F(k)`$ is known one can calculate the phase shift $`\delta (k)`$ by using the definition (22). The asymptotic behavior of the scattering wave function of the potential $`V_N`$ with $`l=0`$ and singularity strength $`\nu `$ at long distances is

$$\psi _N(x,k)\mathrm{sin}(kx\frac{1}{2}\nu \pi +\delta _N(k)),x\mathrm{}.$$
(24)
In our approach $`F_N(k)`$ differs from $`F_0(k)`$ by a rational function of momentum $`k`$ and therefore the expression for $`\delta _N(k)`$ becomes rather complicated when the number of transformation functions is sufficiently large. An equivalent expression which is more convenient for practical calculations is

$$\delta _N(k)=\delta _0(k)\underset{j=\nu +1}{\overset{n}{}}\mathrm{arctan}(k/a_j)\underset{j=\mathrm{\hspace{0.17em}1}}{\overset{n}{}}\mathrm{arctan}(k/b_j),$$
(25)
in agreement with the result obtained in . From (25) it follows that the $`S`$-matrix acquires additional poles at $`a_j`$ and $`b_j`$ through the application of the transformation $`L^{(N)}`$ consistent with the discussion of the previous section.

For $`l=0`$ the effective range expansion is (see e. g. )

$$\underset{k0}{lim}k\mathrm{cot}\delta (k)=\frac{1}{a}+\frac{1}{2}r_0k^2Pr_0^3k^4+\mathrm{},$$
(26)
where $`a`$ is the scattering length and $`r_0`$ effective range. Expanding the phase shift (25) in a power series one obtains

$$\frac{1}{a^{(N)}}=\frac{1}{a^{(0)}}+\left[\underset{j=\nu +1}{\overset{n}{}}a_j^1+\underset{j=1}{\overset{n}{}}b_j^1\right]^1,$$
(27)
and

$$r^{(N)}=r^{(0)}+\frac{2a^{(N)}}{3}\left[1a_{}^{(N)}{}_{}{}^{3}\left(\underset{j=\nu +1}{\overset{n}{}}a_j^3+\underset{j=1}{\overset{n}{}}b_j^3\right)\right].$$
(28)
These formulae can be used to calculate theoretical values to be compared with experimental results, as for example those obtained in the neutron-proton scattering experiments. But first one has to fit the parameters $`a_j`$ and $`b_j`$ in (25) such as to reproduce the phase shift data as close as possible. As a result of the fit one gets a sequence of $`n_+`$ positive and $`n_{}`$ negative poles of the $`S`$-matrix all having absolute values below the pole associated with the ground state of $`h_0`$ if this has a discrete spectrum or below the lower bound of its continuum spectrum otherwise. (We shall comment below on the cases when this is not so.) Every negative pole is associated with a singular transformation function decreasing at infinity and belonging to the singular family of (18). This is just the Jost solution of the initial equation at $`k_j=ia_j`$ which is uniquely defined. A positive pole corresponds to an increasing transformation function (20), which can be either regular or singular at $`x=0`$. When regular it is uniquely defined, whereas singular functions form a one-parameter family.

According to the choice of transformation functions (18) our method is restricted by the condition $`n_{}n_+`$ and one has to distribute the functions related to positive poles among regular and singular families such as the number of singular transformation functions does not exceed the number of regular ones. This means that at least $`n_{}`$ functions related with positive poles should be regular. For the remaining positive poles there are different possibilities either related to the regular $`\left\{v_j\right\}`$ family or to the singular $`\left\{u_j\right\}`$ family. All these possibilities correspond to the same phase shift and hence will give different phase equivalent chains of transformations. When all the remaining positive poles are associated with the regular family the potential will acquire the largest singularity strength $`\nu =n_+n_{}`$ and according to the formula (21) will have no additional discrete levels with respect to $`V_0`$. In this case every function from the system (18) is uniquely defined and the potential $`V_N`$ has no free parameters. This agrees with the well-known statement of the inverse scattering method (see e.g. ) that in the absence of discrete spectrum a local potential is uniquely defined by its phase shift. Such potentials are called in the literature shallow potentials.

When a maximal number of positive poles (equal to $`\gamma =(n_+n_{})/2`$ for $`N`$ even and $`\gamma =(n_+n_{}1)/2`$ for $`N`$ odd) is included in the singular family every such pole corresponds to an additional discrete level of $`V_N`$ and this potential will acquire the smallest singularity strength $`\nu =0`$ (regular potential) for $`N`$ even and $`\nu =1`$ for $`N`$ odd and $`\gamma `$ additional discrete levels. Every positive pole $`a_j`$ associated with a singular function $`u_j`$ will correspond to a discrete level $`E_j=a_j^2`$ of $`V_N`$. According to (20) the function $`u_j`$ has a free parameter $`A_j/B_j`$. Thus we obtain a $`\gamma `$-parameter family of isophase and isospectral potentials. Such potentials are called in the literature deep potentials. Of course there are also intermediate possibilities. When $`n_+=n_{}`$ or $`n_+=n_{}+1`$ one can construct only a deep (regular if $`n_+=n_{}`$) potential without additional discrete levels. For this potential phase preserving transformations which at the same time do not change the angular momentum $`l`$ are impossible.

When $`\nu =n_+n_{}2`$ one can construct a $`2\gamma `$-parameter family of isophase potentials by introducing $`\gamma `$ energy levels into the shallow potential. For this purpose one has to use a chain of $`N+2\gamma `$ transformations, phase equivalent to the initial chain of $`N`$ transformations, to which $`2\gamma `$ transformation functions with pairwise equal factorization constants $`a_j^2=b_j^2`$ should be added. To calculate the Wronskian of such a system one has to use either the formula (15) or (16). The potentials $`V_N`$ form a $`2\gamma `$-parameter family where $`\gamma `$ energy levels play the role of parameters. The other $`\gamma `$ parameters correspond to the constant $`c`$ of (14).

It may happen that the absolute values of some poles is above the position defined by an excited state energy of $`h_0`$. Then using a phase preserving procedure, we first have to delete all such levels including ground state level from the spectrum of $`h_0`$ and then use the above described technique. This means to include additional pairs of functions with equal factorization constants in the system of transformation functions (18). The corresponding energies have to coincide with the positions of the discrete levels to be deleted, which evidently does not change the fitted phase shift. The Wronskian of the entire system of transformation functions has to be calculated by using either (15) or (16).

When the number of negative poles exceeds the number of positive ones, the problem becomes more sophisticated. In some cases it can be solved by introducing transformations changing the orbital angular momentum $`l`$. We are planning to discuss these problems in a subsequent paper .

In this subsection all potentials have an exponential decrease at infinity, provided the initial potential corresponds to $`l=0`$. Therefore all transformed potentials also correspond to $`l=0`$. An useful remark is that when fitting an experimental phase shift one can consider either real poles of the $`S`$-matrix or pairwise mutually conjugated poles. In either case they should be associated with the regular family of the system (18). The corresponding Wronskian is either real or purely imaginary, since the transformation functions are mutually conjugated. The potential difference $`V_NV_0`$ resulting from (11) is real in either case.

### D The Levinson theorem

Here we show that the formula (25) is in agreement with the Levinson theorem in its generalized form  (see also ). For a potential satisfying the condition (1) and decreasing faster than $`x^2`$ at long distances this theorem states that

$$\delta (0)\delta (\mathrm{})=[𝒩+\frac{1}{2}(\nu l)]\pi ,$$
(29)
where $`𝒩`$ is the number of bound states $`l`$ the partial wave and $`\nu `$ the singularity strength. In our case we have to set $`l=0`$. Suppose that this theorem holds for the initial potential $`V_0`$. We shall show that in this case it also holds for the transformed potential $`V_N`$. Since the right hand side of (29) is a linear form in the variables $`𝒩`$ and $`\nu `$ it remains to show that this formula is valid for the case where $`\delta (0)`$, $`\delta (\mathrm{})`$, $`𝒩`$ and $`\nu `$ represent differences between the same quantities associated with the potentials $`V_0`$ and $`V_N`$. In what follows in this subsection, by these quantities we understand the corresponding increments.

The potential $`V_N`$ is supposed to be obtained from $`V_0`$ by an $`l`$-preserving chain of transformations with transformation functions (18). When all additional positive poles of the $`S`$-matrix are related to the regular family, the potential $`V_N`$ has the same number of bound states as $`V_0`$ but its singularity strength is $`\nu _N=\nu =n_+n_{}`$ . Using (25) one finds that $`\delta (0)=0`$, $`\delta (\mathrm{})=\pi (n_+n_{})/2`$ which is in the agreement with (29) at $`𝒩=0`$. If now we transform $`𝒩(n_+n_{})/2`$ regular functions into singular ones the potential $`V_N`$ acquires $`𝒩`$ additional discrete levels and loses $`2𝒩`$ units of singularity strength which also is in agreement with the relation (29). Hence, we conclude that $`l`$-preserving chains of transformations with transformation functions $`(\text{18})`$ do not violate the Levinson theorem.

## IV Application to neutron-proton elastic scattering

Here we apply our method to derive shallow and deep phase equivalent potentials, both reproducing the neutron-proton ($`np`$) experimental elastic scattering phase shift for the $`{}_{}{}^{1}S_{0}^{}`$ partial wave in the laboratory energy interval 0-350 MeV. The experimental values taken from Ref.  are displayed in Table I. Our choice $`V_0(x)=0`$ which allows to use solutions of the free particle Schrödinger equation simplifies the calculations considerably. The theoretical phase shift also displayed in Table I under $`\delta _6`$ is obtained from the formula (25) with $`\delta _0(k)=0`$ and $`\nu =0`$ where we had to fix $`n`$ optimallly and search for $`a_j`$ and $`b_j`$ ($`j`$ = 1,2,…$`n`$). By using the standard Mathematica package NonlinearRegress  we found that we need six $`S`$-matrix poles in order to fit the experimental phase shift very closely. Thus we have $`N=2n=6`$ and the poles are

$$\begin{array}{c}a_1=0.0401,a_2=0.7540,a_3=4.1650,\hfill \\ b_1=0.6152,b_2=2.0424,b_3=4.6000,\hfill \end{array}$$
(30)
in fm<sup>-1</sup> units. In the notation of Subsection III C we have $`n_{}=2`$ and $`n_+=4`$. From Eqs. (27) and (28) with $`1/a^{(0)}=0`$ and $`r^{(0)}=0`$ we obtain $`a^{(6)}=23.7032`$ fm and $`r^{(6)}=2.6235`$ fm. The comparison with the recommended values $`a_{exp}=23.721\pm 0.020`$ fm and $`r_{exp}=2.658\pm 0.062`$ fm  shows that the phase shift is sufficiently well described by the system of poles (30). The negative poles correspond to transformation functions of type $`u_j=\mathrm{exp}(a_jx)`$. Positive poles can correspond either to regular solutions $`u_j=\mathrm{sinh}(a_jx)`$ or to

$$u_j=A_j\mathrm{exp}(a_jx)+B_j\mathrm{exp}(a_jx),$$
(31)
as implied by (20). According to the previous section the poles $`b_j`$ correspond to regular solutions $`v_j=\mathrm{sinh}(b_jx)`$.

We can now derive isophase potentials associated with $`\delta _6`$ by using (11). The system of poles (30) can correspond either to a one level potential (deep potential) with $`E_0=a_3^2`$ or to a shallow potential which has no discrete levels. The latter is the desired case and can be obtained when the function $`u_3`$ is a $`\mathrm{sinh}`$ function. Fig. 1 shows our shallow potential $`V_6`$ obtained in this way together with the Reid68 potential

$$V_{Reid68}(x)=[10.463\mathrm{exp}(0.7x)1650.6\mathrm{exp}(2.8x)+6484.2\mathrm{exp}(4.9x)]/(0.7x),$$
(32)
and the potential of Ref. . The dots represent an updated version  of the Reid soft core potential . In contrast with Ref.  our potential has a correct short range behaviour $`V_6(x)6x^2`$ corresponding to the case $`\nu =n_+n_{}=2`$ (the potential of Ref.  has $`\nu =1`$). In Fig. 2 we present the tail of the absolute value of various potentials in a logarithmic scale. One can see that $`V_6`$ has a smooth behaviour consistent with the one-meson exchange theory as it is the Reid’s potential , in contradistinction to the SUSY potential from Fig. 3 of Ref.  which displays an undesired oscillation around $`r=2`$ fmThe potential drawn in Figs. 1 and 2 under Ref.  represents the output of Eq. (32) of this reference.. Interestingly, our potential is very close to the Reid68 potential (32) in a large interval up to $`5`$ fm.

Actually the choice (31) for $`u_3`$ gives an isospectral and isophase family of $`V_6`$ potentials with the bound state energy $`E_0`$ fixed. In Fig. 3 we display a few deep potentials, members of this family and isophase with the shallow $`V_6`$ potential. These potentials differ from each other only by the value of the coefficient $`A_3^{}=A_3/B_3`$ in (31). One can see that by increasing $`A_3^{}`$ the potential well shifts to the right without changing its depth. The bound state with $`E_0=a_3^2`$ follows the movement, being localized at larger and larger values of $`x`$. A similar behavior has been observed by Bargmann  already in 1949 (see also ). When $`A_3^{}`$ approaches $`1`$ the well concentrates near the origin and becomes very narrow and deep. The ground state becomes well localized near the origin. At very short distances a barrier appears before the well which grows when $`A_3^{}1`$. In the limit $`A_3^{}=1`$ the well becomes infinitely narrow and icompletely collapses, being replaced by an infinitely high barrier which becomes a repulsive core, as for a shallow potential.

One can also construct a deep potential isophase with the shallow potential $`V_6`$ but having a discrete level at a desired place. First one has to delete the bound state $`E_0=a_3^2`$ of $`V_6`$ following the procedure described in Sec III. Next one has to use a $`2`$-parameter family of isophase potentials (Subsection III C). This is possible by adding two new transformation functions, which gives a total of $`8`$ functions. The new functions $`u_4`$ and $`v_4`$ must have $`a_4=b_4`$ ($`a_4<0`$) in order to give $`\delta _8(k)=\delta _6(k)`$. It is enough to take $`u_4=\mathrm{exp}(b_4x)`$ and calculate the Wronskian of the added pair from (14). When $`|a_4|=a_3`$ the potential behaves like a deep $`V_6`$ potential with the energy $`E=a_3^2`$. By varying $`a_4`$ one can change the depth of the well. The family of $`V_8`$ potentials has in fact two parameters: the energy $`E=a_4^2`$ and the parameter $`c`$ of (14). One can select a potential from this family by fixing these two parameters. Here we look for a deep potential which describes the nucleon-nucleon interaction according to microscopic theories based on the quark structure of nucleons (see e. g. Refs.  or ). We wish to obtain a potential close to the deep potential $`V_K`$ of Ref. , given by the formula

$$V_K=1106.21\mathrm{exp}(1.6x^2)10.464\mathrm{exp}(0.7x)/(0.7x)[1\mathrm{exp}(3x)],$$
(33)
where $`x`$ is given in fm. We varied $`a_4`$ and $`c`$ such as to be as close as possible to (33). From a least square fit we obtained $`c=0.155`$ fm and $`a_4=3.7944`$ fm<sup>-1</sup> which gives $`E`$ = 596.42 MeV. In this way we obtained the potential $`V_8`$ which is compared to $`V_K`$ of (33) in Fig. 4. One can see that the two potentials are quite close to each other. The bound state of $`V_K`$ is 442.MeV which is not much different from our value. The form of $`V_K`$ was chosen on qualitative grounds related to the approximations used in microscopic studies of the nucleon-nucleon interaction . The potential $`V_8`$ of Fig. 4 is as legitimate as $`V_K`$ inasmuch as it correctly describes the experimental phase shift. The difference between the two potentials is partly due the difference in the fitted phase shifts and the $`E_{lab}`$ interval chosen for the fit in each case. If instead of fitting the phase shift to the experimental values we were fitting it to the phase shift corresponding to the potential (33), the difference would have been smaller.

## V CONCLUSIONS

In this study we propose an alternative method to derive families of phase equivalent (isophase) local potentials based on the application of an $`N`$th order Darboux transformation or an equivalent chain of $`N`$ first order transformations. In practical applications this method is simpler and more efficient than the conventional SUSY approach of Refs. . First one has to look for imaginary poles of the S-matrix, i.e. one has to find the parameters $`a_j`$ and $`b_j`$ of the phase shift $`\delta _N(k)`$ of (25) which fit a particular experimental phase shift. One can start with the lowest order $`N`$ = 2 in (25) and increase it until a good fit of the phase shift is achieved. The second step consists in finding an appropriate mapping between the $`S`$-matrix poles found in this way and the system of functions (18) which give a nodeless Wronskian. In terms of these functions the formula (11) gives phase equivalent potentials, among which one is shallow and the others form families of deep potentials from which one can select a particular one for a specific physical problem. We have applied our method to derive potentials describing the neutron-proton experimental phase shift for the $`{}_{}{}^{1}S_{0}^{}`$ partial wave. We obtained a shallow potential very close to the Reid’s soft core potential and its deep phase equivalent partner is close to potentials resulting from microscopic derivations of the nucleon-nucleon interaction. In order to obtain phase equivalent potentials describing higher partial waves $`l`$-changing Darboux chains of transformations are necessary. This is the scope of a subsequent study.

## A

Here we prove the formula (21). The proof is split into two steps: 1) Taking an initial potential with $`\nu =\nu _0=0`$ we prove (21) for this case. 2) We introduce a number of regular solutions to the previous case which corresponds to $`\nu >0`$. Finally we comment on singular initial potentials.

Let us consider a chain of $`N=2n`$ Darboux transformations for a potential $`V_0(x)`$ with $`\nu =\nu _0=0`$ and transformation functions given by (18) with $`\nu =0`$. Since all $`a`$’s and $`b`$’s are different from each other this chain is equivalent to an $`N`$th order transformation. The transformed potential $`V_N`$ is given by (11) and the eigenfunctions $`\psi _N`$ by (8) or (12). One can also find the Jost solution of $`h_N`$ by applying the same operator $`L^{(N)}`$ to the Jost solution of $`h_0`$. The result is the function $`\stackrel{~}{f}_N(x,k)=L^{(N)}f_0(x,k)`$, which like $`\psi _N`$ from (8), is proportional to the ratio of two determinants. The long-distance asymptotic value of $`\stackrel{~}{f}_N(x,k)`$ is proportional to the asymptotic value of the Jost solution. Therefore using the usual definition of the Jost function as the value of the Jost solution at the origin  one has

$$F_N(k)=\frac{\stackrel{~}{f}_N(0,k)}{[e^{ikx}\stackrel{~}{f}_N(x,k)]_{x\mathrm{}}^{}{}_{}{}^{}}$$
(A1)
The definition (8) and our special choice of transformation functions (18) implies that

$$\stackrel{~}{f}_N(x,k)=W(u_1,v_1,\mathrm{},u_n,v_n,f_0(x,k))W^1(u_1,v_1,\mathrm{},u_n,v_n)W_fW_N^1$$
(A2)
Therefore both for $`x=0`$ and $`x\mathrm{}`$ we have to calculate the ratio of two determinants of orders $`N+1`$ and $`N`$ respectively, but of a similar structure.

First let us consider the function $`\stackrel{~}{f}_N(x,k)`$ at $`x=0`$. Using the definition $`f_0(0,k)=F_0(k)`$ the determinant appearing in the numerator of (A2) becomes

$$W_f=\left|\begin{array}{cccccccc}C_1\hfill & 0\hfill & C_2\hfill & 0\hfill & \mathrm{}\hfill & C_n\hfill & 0\hfill & F_0(k)\hfill \\ c_1\hfill & D_1\hfill & c_2\hfill & D_2\hfill & \mathrm{}\hfill & c_n\hfill & D_n\hfill & \stackrel{~}{F}(k)\hfill \\ C_1a_1^2\hfill & 0\hfill & C_2a_2^2\hfill & 0\hfill & \mathrm{}\hfill & C_na_n^2\hfill & 0\hfill & F_0(k)(ik)^2\hfill \\ c_1a_1^2\hfill & D_1b_1^2\hfill & c_2a_2^2\hfill & D_2b_2^2\hfill & \mathrm{}\hfill & c_na_n^2\hfill & D_nb_n^2\hfill & \stackrel{~}{F}(k)(ik)^2\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ C_1a_1^{2n2}\hfill & 0\hfill & C_2a_2^{2n2}\hfill & 0\hfill & \mathrm{}\hfill & C_na_n^{2n2}\hfill & 0\hfill & F_0(k)(ik)^{2n2}\hfill \\ c_1a_1^{2n2}\hfill & D_1b_1^{2n2}\hfill & c_2a_2^{2n2}\hfill & D_2b_2^{2n2}\hfill & \mathrm{}\hfill & c_na_n^{2n2}\hfill & D_nb_n^{2n2}\hfill & \stackrel{~}{F}(k)(ik)^{2n2}\hfill \\ C_1a_1^{2n}\hfill & 0\hfill & C_2a_2^{2n}\hfill & 0\hfill & \mathrm{}\hfill & C_na_n^{2n}\hfill & 0\hfill & F_0(k)(ik)^{2n}\hfill \end{array}\right|$$

$$=C_1\mathrm{}C_nD_1\mathrm{}D_nF_0(k)\left|\begin{array}{ccccc}1\hfill & 1\hfill & \mathrm{}\hfill & 1\hfill & 1\hfill \\ a_1^2\hfill & a_2^2\hfill & \mathrm{}\hfill & a_n^2\hfill & (ik)^2\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ a_1^{2n2}\hfill & a_2^{2n2}\hfill & \mathrm{}\hfill & a_n^{2n2}\hfill & (ik)^{2n2}\hfill \\ a_1^{2n}\hfill & a_2^{2n}\hfill & \mathrm{}\hfill & a_n^{2n}\hfill & (ik)^{2n}\hfill \end{array}\right|\times \left|\begin{array}{ccc}1\hfill & \mathrm{}\hfill & 1\hfill \\ b_1^2\hfill & \mathrm{}\hfill & b_n^2\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill \\ b_1^{2n2}\hfill & \mathrm{}\hfill & b_n^{2n2}\hfill \end{array}\right|$$

$$=F_0(k)\underset{j=1}{\overset{n}{}}(C_jD_j)\underset{m>j}{\overset{n}{}}(a_j^2a_m^2)(b_j^2b_m^2)\underset{j=1}{\overset{n}{}}(a_j^2+k^2)$$
(A3)
Here $`C_j=u_j(0)`$, $`c_j=u_j^{}(0)`$, $`D_j=v_j^{}(0)`$, $`\stackrel{~}{F}(k)=f_0^{}(0,k)`$. Note that the derivative of the Jost function, $`\stackrel{~}{F}(k)`$, does not appear in the final expression for $`W_f`$. The above result has been obtained by first using the following determinant property

$$\left|\begin{array}{cc}M_{p\times p}& 0\\ \mathrm{}& \stackrel{~}{M}_{q\times q}\end{array}\right|=detM_{p\times p}det\stackrel{~}{M}_{q\times q}$$
(A4)
and next by expanding the two Vandermonde determinants (see e.g. ) from the middle line.

The determinant $`W_N`$ which appears in the denominator of (A2) is of even order. It is obtained from the above determinant by the removal of the last row and column. Its final form is similar to the one above. It has opposite sign and the factors depending on $`k`$ are absent. Finally one finds

$$\stackrel{~}{f}_N(0,k)=(1)^nF_0(k)\underset{j=1}{\overset{n}{}}(k^2+a_j^2)$$
(A5)
Let us now consider the denominator of (A1), i. e. $`\stackrel{~}{f}_N(x,k)^{}`$ for $`x\mathrm{}`$. In this limit all functions entering the Wronskians of (A2) tend to exponentials which cancel out with each other except for $`\mathrm{exp}(ikx)`$ which cancels out with the exponential present in (A1). The remaining factor is a ratio of two Vandermonde determinants of orders $`N+1`$ and $`N`$ respectively containing powers of $`a_j`$ and $`b_j`$. The determinant in the numerator differs from that in the denominator by the presence of a column containing powers of $`ik`$ and an additional row composed of $`(N+1)`$th power of $`a_1,b_1,`$$`a_n,b_n,ik`$. Using once more the expression of Vandermonde’s determinant one obtains

$$\stackrel{~}{f}_N(x,k)|_x\mathrm{}i^N\underset{j=1}{\overset{n}{}}(k+ia_j)(k+ib_j)$$
(A6)
Introducing (A5) and (A6) in (A1) one gets the following relation between the initial and final Jost functions

$$F_N(k)=F_0(k)\underset{j=1}{\overset{n}{}}\frac{kia_j}{k+ib_j}.$$
(A7)
Moreover, since the Wronskian $`W_N`$ is different from zero at $`x=0`$ the singularity strength parameter of $`V_N`$ is zero and it is regular at the origin.

Consider now the case where the first $`\nu `$ functions of (18) belong to the regular family and calculate the singularity strength of $`V_N`$. Note that the system (18) can be divided into two subsystems, one containing an equal number of singular and regular functions, and the other composed of regular functions only. Accordingly the $`N`$th order transformation can be decomposed into a superposition of two transformations. The first transformation of order $`N\nu `$ which is based on functions paired by their regular property will not affect the regularity of the potential at the origin, whereas the second one of order $`\nu `$ composed of regular functions only will increase the singularity parameter of the potential by one unit for each transformation function (see e.g. ). So, the singularity strength parameter of $`V_N(x)`$ is equal to $`\nu `$.

In the following we have to use the definition of the Jost function

$$F(k)=\underset{x\mathrm{\hspace{0.17em}0}}{lim}\frac{(kx)^\nu }{i^\nu (2\nu 1)!!}f(x,k)$$
(A8)
We give detailed calculations for the particular case $`\nu =2`$ from which a general formula can easily be derived. In this case both $`W_f=W(v_1,v_2,u_3,v_3,\mathrm{},u_n,v_n,f_0)`$ and $`W_N=W(v_1,v_2,u_3,v_3,\mathrm{},u_n,v_n)`$ are equal to zero at $`x=0`$. Since $`f_0(0,k)0`$ only one function (either $`v_1`$ or $`v_2`$) in $`W_f`$ is an unpaired regular solution and hence the zero of $`W_f(0)`$ is of the first degree. For $`\nu =2`$ the second logarithmic derivative of $`W_N(x)`$ at $`x=0`$ is of order $`\nu (\nu +1)=6`$. We conclude that the zero of $`W_N(0)`$ is of the third degree. To evaluate the indeterminate form defining the Jost function (see (A8) for $`\nu =2`$ and $`f(x,k)=f_N(x,k)`$) one can use L’Hôspital’s rule to get $`(kx)^2W_f/W_N=k^2W_f/xx^3/W_N3!k^2W_f^{}(0)/W_N^{\prime \prime \prime }(0)`$ when $`x0`$, because the first two derivatives of $`W_N(0)`$ are zero. Since the derivative of an $`N`$th order Wronskian $`W(f_1,\mathrm{},f_N)`$ is a determinant obtained from the same Wronskian by the replacement of the last row composed of the derivatives $`f_j^{(N1)}`$ by the derivatives $`f_j^{(N)}`$ one can use the property (A4) to factorize $`W_f^{}(0)`$ into a product of two Vandermonde determinants $`W_{V1}`$ and $`W_{V2}`$ of orders $`n`$ and $`n+1`$ respectively. When calculating the third derivative of $`W_N`$ at $`x=0`$ one can see that it is equal to a sum of two equal determinants which can be factorized according to (A4) in a product of two Vandermonde determinants $`W_{V3}`$ and $`W_{V4}`$ of orders $`n1`$ and $`n+1`$ respectively. We do not give here the explicit form of these determinants. The factor $`2`$ coming from two equal Vandermonde determinants will cancel out with the corresponding factor coming from $`3!`$ and the factor $`3`$ from this factorial will cancel out with $`(2\nu 1)!!`$ for $`\nu =2`$, the latter being present in the definition of the Jost function (A8). Hence, $`(kx)^2W_f/W_Nk^2W_{V1}W_{V2}/(W_{V3}W_{V4})`$. When evaluating the corresponding Vandermonde determinants one can see that the factors independent of $`k`$ are exactly the same both in the numerator and the denominator. So, one gets the following Jost function for $`V_N(x)`$

$$F_N(k)=F_0(k)\frac{k^2}{(k+ib_1)(k+ib_2)}\underset{j=\mathrm{\hspace{0.17em}3}}{\overset{n}{}}\frac{kia_j}{k+ib_j}$$
(A9)
It is clear from here that every additional regular function in the system (18) with the parameter $`b_j`$ produces a factor $`k/(k+ib_j)`$ in the Jost function which corresponds exactly to the formula (21). Moreover, a regular at the origin potential admits two kinds of transformations, the ones which do not affect its singularity ($`\nu =0`$ in (18)) and others which increase its singularity ($`\nu >0`$). Actually if the initial potential has the singularity parameter $`\nu _00`$, three types of transformations are possible. In addition to two the previous types of transformations existing for the case $`\nu _0=0`$ transformations which decrease the singularity parameter are also possible. All these possibilities become clear if one takes into consideration the fact that all transformations are invertible (this means that from an intermediate element of a chain one can move in both directions) and an intermediate member of a previously considered chain is just a potential with a singularity strength different from zero.

ACKNOWLEDGMENTS.

We are most grateful to Mart Rentmeester for providing details about the Reid93 potential. Illuminating discussions with Daniel Baye and Francis Michel are gratefully acknowledged. One of us (B.F.S) is grateful to the Theoretical Fundamental Physics Laboratory of the University of Liège for warm hospitality during autumn 2001 when this work has been started. He also acknowledges a partial support from RFBR.


text:
# Flatness of the setting Sun

## I Atmospheric refractions

Atmospheric refraction is between the first scientifically approached phenomenon. This effect is responsible for the apparent scintillation of stars, mirages, the puzzling and spectacular green flash, difference between apparent and real position of stars or the asymmetric rim of the Sun near the horizon. The refractive index of dry air is very close to $`1`$. Its small dependence on temperature and pressure leads to a refractive index gradient in the atmosphere. Although this gradient is very small, due to the large distances travelled by the light-rays in the atmosphere, it can result observable bending or dispersion.

Atmospheric refraction is usually divided into three categories: astronomical, terrestrial and geodesic. Terrestrial refraction appears when both the object and observer are within the Earth’s atmosphere. This refraction is responsible for ordinary mirages, and has been extensively studied . The first one reporting and accounting for mirages was Aristotle in Meteorologica . Correctly, he concluded that dense air layers can act as mirrors, and considered this effect responsible for mirages. (A complete historic and bibliographic study of mirages can be found on the splendid web-page of A.T.Young .) Geodesic refraction is a special case of terrestrial refraction where both the object and observer are at low altitudes. A well-known example for this is surveying. We speak about astronomical refraction when a terrestrial observer detects ray-bending effects for the light-rays coming from objects outside the Earth’s atmosphere. This refraction is responsible for the difference between the real and apparent position of stars near the horizon, the green flash or the asymmetric rim of the setting (or rising) Sun. In the present paper we discuss in detail the flatness of the Sun or full-Moon near the horizon. Particularly, we are interested in the extreme values this flatness can take for standard atmospheric profiles, and the dependence of the flatness on the altitude of the observer and meteorological conditions (temperature and pressure). A similar study, considering fixed atmospheric conditions was recently published by Thomas and Joseph . Here we plan a more complete analysis, considering three different theoretical methods, experiments, computer simulations and illustration of the results with a collection of pictures and movies.

The spectacular nature of sunset or sunrise makes the phenomenon attractive to students. After our experience it can effectively be used to exemplify refraction phenomena and exercising the principles of geometrical optics. Our paper is structured as follows. First, we introduce the atmosphere model, and the refractive index profile is given. Then, we present three different methods for computing the ray-path in this atmosphere. Results computed for different meteorological conditions and observation altitude will close the theoretical part. Experimentally, we analyze and compare with theory the measured flatness for some sunsets photographed or video-filmed by us. Finally, we discuss some extreme and unusual conditions which are exemplified by pictures and movies on a web-page supporting this study. From this web-page one can also freely download a computer-program written by us, which visualizes the sunset for arbitrary meteorological conditions and observation altitude.

We emphasize here, that our study is restricted to standard atmospheric conditions with smooth temperature and pressure profile. Non-standard, but often encountered atmospheric profiles leads to non-standard distortions of the solar rim, and are not considered within this study.

## II The optical atmosphere model (refractive index as a function of altitude)

The atmosphere of the Earth is composed mainly of $`N_2`$ ($`79\%`$) and $`O_2`$ ($`20\%`$), and extends up to a few hundreds kilometer height. The refractive index of air is very close to $`1`$, and depends slightly on its pressure and temperature following Edlen’s semi-empirical law , valid for dry air:

$$n=1+10^6(776.2+4.36\times 10^8\nu ^2)\frac{P}{T}.$$
(1)
In the above formula $`\nu `$ is the wave-number of the light in $`cm^1`$, $`P`$ is the pressure in $`kPa`$ and $`T`$ the temperature in $`K`$. Since the pressure and temperature varies within the atmosphere, we get a refractive index gradient which is responsible for atmospheric refractions. Although the variations in $`n`$ are quite small, the large distances travelled by light-rays in the atmosphere makes refraction effects observable and sometimes important.

As a function of altitude several layers with different physical properties are distinguishable. The lowest layer extending from sea-level to approximative $`z_t=14km`$ height is called the troposphere. This is the region where the wheatear takes place, i.e. the region of rising and falling packets of air. In this layer the air pressure drops drastically, at the top of the troposphere being only $`10\%`$ of the value measured at sea-level. In the troposphere the temperature decreases almost linearly with the altitude (Fig. 1). A thin buffer zone between the troposphere and the next layer (the stratosphere) is called the tropopause. Within the tropopause the temperature is in good approximation constant. The stratosphere extends from the altitude of approximately $`18km`$ up to $`z_s=50km`$ (Fig. 1). Within the stratosphere the air flow is mostly horizontal. In the upper part of the stratosphere we have the ozone layer. Inside the stratosphere the pressure decreases further with the altitude, but surprisingly the temperature increases with height. Above the stratosphere we find the mesosphere and the ionosphere. In these regions the air is very rare and the temperature profile is shown in Fig. 1.

From the viewpoint of atmospheric refractions only the first two layers, the troposphere and the stratosphere are important. In the upper layers of the atmosphere the air is so rarefied, that the refractive index can be considered $`1`$ within a good approximation.

An accepted and widely used model for the atmosphere is the U.S. Standard Atmosphere, established in 1953 and re-actualized in 1976 . The U.S. Standard Atmosphere consists of single profiles, representing the idealized steady-state atmosphere for moderate solar activity. The listed parameters include temperature, pressure, density, gravitational acceleration, mean particle speed, mean collision frequency, mean free path, etc. as a function of altitude. In our study we consider an atmosphere model with a spherical symmetry, all relevant physical quantities (temperature and pressure) varying only as a function of altitude. We then calculate the refraction index of air as a function of altitude as follows:
(i) The temperature profile in the troposphere is linearly decreasing with a $`\lambda =6.5K/km`$ lapse-rate, as suggested by the U.S. Standard Atmosphere. Following again the U.S. Standard Atmosphere we consider the temperature constant within the tropopause. Although the temperature increases as a function of altitude in the stratosphere, due to the rarefied air (small pressure) the refractive index is close to $`1`$, and for the sake of simplicity we assume the temperature as constant in this region, too. (We checked that this approximation is fully justified.) Above $`z_s=50km`$ height we assume that the refractive index is $`1`$, and do not calculate it anymore by Edlen’s formula. Up to the top of the stratosphere the temperature profile is therefore presumed as

$`T(z)=T_0\lambda zforz<z_t`$ (2)
$`T(z)=T_0\lambda z_tforz_tzz_s,`$ (3)
where $`z`$ is the altitude from sea-level, and $`T_0`$ the temperature at sea-level.
(ii) For the pressure profile we use a barometric formula in which we take into account the variation of temperature with altitude. Considering a vertical slice of air with thickness $`dz`$, the variation of pressure within this slice is due to the hydrostatic pressure

$$dP=\rho (z)g(z)dz,$$
(4)
where $`g(z)`$ is the gravitational acceleration and $`\rho (z)`$ the density of air at height $`z`$:

$$\rho (z)=\frac{M}{V}=\frac{Nm}{V}=\frac{P(z)m}{kT(z)}.$$
(5)
(We denoted by $`m`$ the mass of one molecule, $`N`$ the number of molecules in volume $`V`$, $`T(z)`$ the temperature at height $`z`$ and $`k`$ the Boltzmann constant.) Since we focus on the troposphere and the stratosphere only, the $`z`$ height is small in comparison with the radius of the Earth ($`R_E6378km`$), thus the $`g`$ gravitational acceleration can be considered constant. We can write thus

$$dP=\frac{P(z)mgdz}{kT(z)},$$
(6)
which is a separable differential equation for $`P(z)`$. Integrating this between a height $`z_0`$ where the pressure is $`P_0`$ and an arbitrary height $`z`$, we get the desired barometric formula:

$$P(z)=P_0exp\left(\frac{mg}{k}_{z_0}^z\frac{dz^{}}{T(z^{})}\right).$$
(7)
Considering $`z_0=0`$ (sea-level) and using the temperature profile given by equation (2) and (3) a simple calculus leads us to

$$P(z)=P_0\left[1\frac{\lambda z}{T_0}\right]^{mg/k\lambda }$$
(8)
for $`z<z_t`$, and

$$P(z)=P_0\left[1\frac{\lambda z_t}{T_0}\right]^{mg/k\lambda }exp\left(\frac{mg(zz_t)}{k(T_0\lambda z_t)}\right),$$
(9)
for $`z_tzz_s`$. Plugging this and (2)-(3) in Edlen’s formula (1), we get for the troposphere and stratosphere the variation of the refractive index as a function of altitude. As stated before, for altitudes higher than $`z_s`$ we can merely use $`n=1`$. The refractive index profile as a function of altitude calculated in this manner for the parameters of the U.S. Standard Atmosphere is plotted in Fig. 2.

## III Trajectory of a light-ray in the atmosphere

We present now three different methods for computing the path of a light-ray in our optical atmosphere model. We proved that the results given by these methods are the same, justifying our forthcoming theoretical considerations. None of these methods is purely analytical, they all make use of numerical calculations to derive the light-ray trajectory. In the following we briefly describe these methods and sketch how one can compute the deviation angle due to atmospheric refraction for light-rays coming from distant sources. In order to avoid the phenomenon of dispersion let us first consider monochromatic light.

### A The Integral method

This method closely follows the one discussed by Smart . To illustrate the method we use the geometry from Fig. 3.

An observer in $`M`$ (at altitude $`z_0`$) detects the light-source placed in $`S_r`$ (altitude $`z^{}`$). Let the arc $`S_rM`$ be the presumed path of the light-ray between the $`S_r`$ source and the observer in $`M`$. The observer detects the light ray in the $`MS`$ direction, which is tangent to the ray path in $`M`$. The source will be positioned by the observer in $`S`$ ($`SO`$ and $`S_rO`$ have radial directions, thus the light-rays in these directions would not bend within our optical atmosphere model). We denote by $`\delta `$ the apparent inclination angle, characterizing the direction of the $`S`$ image. Let $`\delta _{dev}`$ be the deviation angle of $`MS`$ relative to the real $`MS_r`$ direction of the source. The meaning of the $`\theta `$ and $`\theta _r`$ angles are obvious from the figure. We are now interested to compute $`\delta _{dev}`$ as a function of $`\delta `$. Presuming that $`z^{}>>z_0`$ (i.e. the source is very far from the Earth), we can approximate $`MS`$, $`MS_r`$ and $`OS`$ by $`z^{}`$. Some elementary geometry will convince us that the following approximations are justified:

$`\theta {\displaystyle \frac{\pi }{2}}\delta arcsin\left[{\displaystyle \frac{z_0+R_E}{z^{}+R_E}}\mathrm{cos}\delta \right]`$ (10)
$`\delta _{dev}\theta _r\theta +{\displaystyle \frac{R_E+z_0}{R_E+z^{}}}(\mathrm{sin}\theta _r\mathrm{sin}\theta )`$ (11)
In order to get the desired $`\delta _{dev}(\delta )`$ dependence we need $`\theta _r`$ as a function of $`\delta `$.

Let us follow now a light-ray approaching the Earth and let us consider the atmosphere stratified in infinitesimally thin layers of thickness $`\mathrm{\Delta }r`$, with slightly different refractive indices (Fig. 4).

A useful relation between the initial incident angle $`\alpha _0`$ (in layer with refractive index $`n_0`$) and a later incident angle $`\alpha _k`$ (for a layer with refractive index $`n_k`$) can be derived. Using Snell’s law and the notations from Fig. 4, we can write:

$`{\displaystyle \frac{\mathrm{sin}\alpha _0}{\mathrm{sin}\alpha _1^{}}}={\displaystyle \frac{n_1}{n_0}}`$ (12)
$`\mathrm{\Delta }\theta _1=\alpha _1\alpha _1^{}.`$ (13)
This leads to:

$$\mathrm{sin}\alpha _1=\frac{n_0}{n_1}\mathrm{sin}\alpha _0\left(\mathrm{cos}\mathrm{\Delta }\theta _1+\mathrm{sin}\mathrm{\Delta }\theta _1\mathrm{cot}\alpha _1^{}\right).$$
(14)
Similarly

$$\mathrm{sin}\alpha _2=\frac{n_1}{n_2}\mathrm{sin}\alpha _1(\mathrm{cos}\mathrm{\Delta }\theta _2+\mathrm{sin}\mathrm{\Delta }\theta _2\mathrm{cot}\alpha _2^{}),$$
(15)
and combining with (14) it leads to:

$`\mathrm{sin}\alpha _2={\displaystyle \frac{n_0}{n_2}}\mathrm{sin}\alpha _0(\mathrm{cos}\mathrm{\Delta }\theta _1+\mathrm{sin}\mathrm{\Delta }\theta _1\mathrm{cot}\alpha _1^{})\times `$ (16)
$`(\mathrm{cos}\mathrm{\Delta }\theta _2+\mathrm{sin}\mathrm{\Delta }\theta _2\mathrm{cot}\alpha _2^{}).`$ (17)
Generalizing the above equation we get:

$`\mathrm{sin}\alpha _k={\displaystyle \frac{n_0}{n_k}}\mathrm{sin}\alpha _0(\mathrm{cos}\mathrm{\Delta }\theta _1+\mathrm{sin}\mathrm{\Delta }\theta _1\mathrm{cot}\alpha _1^{})\times \mathrm{}\times `$ (18)
$`(\mathrm{cos}\mathrm{\Delta }\theta _k+\mathrm{sin}\mathrm{\Delta }\theta _k\mathrm{cot}\alpha _k^{}).`$ (19)
Since the $`\mathrm{\Delta }\theta _i`$ angles are small, the following approximations are justified:

$`\mathrm{sin}\mathrm{\Delta }\theta _i\mathrm{\Delta }\theta _i,`$ (20)
$`\mathrm{cos}\mathrm{\Delta }\theta _i1,`$ (21)
$`\mathrm{cot}\alpha _i^{}{\displaystyle \frac{\mathrm{\Delta }r}{r_i\mathrm{\Delta }\theta _i}}.`$ (22)
Using these approximations and the fact that $`\mathrm{\Delta }r<<r_i`$ we get:

$`\mathrm{sin}\alpha _k={\displaystyle \frac{n_0}{n_k}}\mathrm{sin}\alpha _0{\displaystyle \underset{i=1}{\overset{k}{}}}(1+\mathrm{\Delta }\theta _i{\displaystyle \frac{\mathrm{\Delta }r}{r_i\mathrm{\Delta }\theta _i}})`$ (23)
$`{\displaystyle \frac{n_0}{n_k}}\mathrm{sin}\alpha _0(1+{\displaystyle \underset{i=1}{\overset{k}{}}}{\displaystyle \frac{\mathrm{\Delta }r}{r_i}}){\displaystyle \frac{n_0}{n_k}}\mathrm{sin}\alpha _0(1+{\displaystyle _{r_k}^{r_0}}{\displaystyle \frac{dr}{r}})=`$ (24)
$`{\displaystyle \frac{n_0}{n_k}}\mathrm{sin}\alpha _0(1+\mathrm{ln}{\displaystyle \frac{r_0}{r_k}}).`$ (25)
Presuming now $`r_0/r_k1`$, we get:

$$\mathrm{sin}\alpha _k=\frac{r_0n_0\mathrm{sin}\alpha _0}{r_kn_k}.$$
(26)
Following the geometry from Fig. 4 we get

$$\frac{\mathrm{\Delta }r}{r_i}\mathrm{tan}\alpha _i=\frac{\mathrm{\Delta }r}{r_i}\frac{1}{\sqrt{\frac{1}{\mathrm{sin}^2\alpha _i}1}},$$
(27)
which in the $`\mathrm{\Delta }r0`$ and $`\alpha _i=\alpha (r)`$ continuous limit yields:

$$d\theta =\frac{dr}{r}\frac{1}{\sqrt{\frac{1}{\mathrm{sin}^2\alpha (r)}1}}.$$
(28)
Denoting by $`\alpha _M`$ the final incident angle at the observer and using (26) we get:

$$d\theta =\frac{dr}{r}\frac{1}{\sqrt{\left[\frac{n(r)}{n(r_M)}\right]^2\frac{1}{\mathrm{sin}^2\alpha _M}\frac{r}{r_M}1}}.$$
(29)
By using the fact that $`\mathrm{sin}^2\alpha _M=\mathrm{cos}^2\delta `$ we can finally determine the angle $`\theta _r`$ as a function of $`\delta `$ by integrating (29):

$$\theta _r=_{z_0+R_E}^{z^{}+R_E}\frac{dz}{\sqrt{\left[\frac{n(z)}{n(z_0)}\right]^2\frac{1}{\mathrm{cos}^2\delta }\frac{z+R_E}{z_0+r_E}1}}.$$
(30)
The above integral can be only numerically computed. There is of course a singularity at $`z=z_0`$, which can be eliminated by a Gauss-Chebisev expansion, or by adjusting the step in the numerical integration in the vicinity of the singularity (the method followed by us). Computing numerically $`\theta _r`$ as a function of $`\delta `$ and by getting from (10) $`\theta `$ as a function of $`\delta `$, we are able now to calculate from (11) the desired $`\delta _{dev}(\delta )`$ dependence.

It is important to mention here that applying the method for $`\delta <0`$ is not straightforward. One must first find in this case the closest point $`C`$ of the light-ray trajectory relative to the Earth’s surface. This can be done by using the fact that the trajectory is symmetric in the vicinity of this point. Then, we decompose the trajectory in two parts, the first part is from $`M`$ to $`C`$, and the second from $`C`$ to $`S_r`$. The light-ray can be followed now by the presented methods both on the $`CS_r`$ and $`CM`$ segments. The deviation angle can be also computed.

### B Using the Fermat principle

The Fermat principle states that light travels between two point along that path which requires the least time, as compared to other nearby paths. We can of course reformulate the Fermat principle by using the optical path instead of time. If the trajectory of a light-ray travelling in the $`XOY`$ plane is described by the $`y=y(x)`$ curve (see the geometry in Fig. 5), we have that:

$`s={\displaystyle _{P_{ini}(x_0,y_0)}^{P_{fin}(x_f,y_f)}}n[x,y(x)]𝑑s=`$ (31)
$`{\displaystyle _{x_0}^{x_f}}n[x,y(x)]\sqrt{1+y^{}(x)^2}𝑑x,`$ (32)
should have a local minima. In the above formula $`n[x,y(x)]`$ denotes the refractive index of the medium in point with coordinates $`x,y(x)`$. In our case $`n`$ has a spherical symmetry, depending on the $`z=\sqrt{x^2+y(x)^2}R_E`$ altitude only.

One can write then (32) as:

$$s=_{x_0}^{x_f}f[x,y(x),y^{}(x)]𝑑x,$$
(33)
with:

$$f[x,y(x),y^{}(x)]=n[\sqrt{x^2+y(x)^2}R_E]\sqrt{1+y^{}(x)^2}.$$
(34)
We are looking now for the $`y(x)`$ function that minimizes $`s`$. If we fix the two points $`P_{ini}`$ and $`P_{fin}`$ between which the light travels, the minima of $`s`$ leads to a classical variational problem

$$\delta s=0,$$
(35)
with:

$`\delta [y(x)]_{x=x_0}=0,`$ (36)
$`\delta [y(x)]_{x=x_f}=0.`$ (37)
The solution of this problem is well-known , and given by the second-order differential equation:

$$\frac{f}{y}\frac{d}{dx}\left[\frac{f}{y^{}}\right]=0.$$
(38)
It is straightforward to show that:

$`{\displaystyle \frac{f}{y}}={\displaystyle \frac{n}{z}}[x,y(x)]{\displaystyle \frac{y(x)}{\sqrt{y^2(x)+x^2}}}\sqrt{1+y^{}(x)^2},`$ (39)
$`{\displaystyle \frac{d}{dx}}\left[{\displaystyle \frac{f}{y^{}}}\right]={\displaystyle \frac{n}{z}}[x,y(x)]{\displaystyle \frac{y^{}(x)y(x)+x}{\sqrt{y(x)^2+x^2}}}\times `$ (40)
$`{\displaystyle \frac{y^{}(x)}{\sqrt{1+y^{}(x)^2}}}+n[x,y(x)]y^{\prime \prime }(x){\displaystyle \frac{1}{[1+y^{}(x)]^{3/2}}}`$ (41)
By simple algebra we obtain from this a second-order differential equation for the $`y(x)`$ equation describing the trajectory of the light-ray:

$`y^{\prime \prime }(x)={\displaystyle \frac{n}{z}}[x,y(x)]{\displaystyle \frac{[1+y^{}(x)^2]}{n[x,y(x)]\sqrt{x^2+y(x)^2}}}\times `$ (42)
$`\{y(x)[1+y^{}(x)^2]y^{}(x)[y(x)y^{}(x)+x]\}.`$ (43)
Taking the $`n(z)`$ refractive index profile from our optical atmosphere model, equation (43) can be numerically integrated. We start from a $`P_{ini}(x_0,y_0)`$ point and consider a $`\mathrm{tan}\delta =y^{}(x_0)`$ initial derivative. The angle $`\delta `$ will be the apparent inclination angle of the light-ray in $`P_{ini}`$ (the point where the observer is presumed). The $`y(x)`$ trajectory can be computed by numerically integrating (43) with the $`y(x_0)=y_0`$ and $`y^{}(x_0)=\mathrm{tan}\delta `$ initial conditions. We construct thus the $`y(x)`$ trajectory from point-to-point, up to an altitude $`zz_s`$, where $`n(z)=1`$ can be presumed. The derivative of $`y(x)`$ at this point will determine the $`\delta _r`$ angle of the light-ray at the border of the optical atmosphere. For altitudes higher than this, the trajectory of the light-ray is presumed recti-linear. The deviation angle is then simply approximated as $`\delta _{dev}=\delta \delta _r`$, and the desired $`\delta (\delta _{dev})`$ dependence can be numerically computed.

### C The simulation method

As its name suggests this is simply a computer simulation method in which we follow the light-ray by segments of infinitesimally small and fixed-length lines. The angle between two elementary line is given by estimating the refractive index at their end-points and by using Snell’s law. The trajectory of the light-ray is then constructed starting from the observer with an initial angle $`\delta `$, and computing the ray-path from point-to-point until it leaves the optical atmosphere and the refractive index can be taken as $`1`$. At this point the direction of the light-ray determines the final angle $`\delta _r`$, and the deviation is computed as $`\delta _{dev}=\delta \delta _r`$. Resembling the method based on the Fermat principle, this simulation also works smoothly for arbitrary positive or negative $`\delta `$ values. Due to the fact that for an elementary step the changes in the refractive index is rather small, it is crucial to work with the maximal precision offered by the computing environment.

We mention here that all three methods give identical results, proving their applicability and our theoretical considerations. As an example, for the parameters of the U.S. Standard Atmosphere, results for the $`\delta _{dev}`$ deviation angle versus the apparent inclination angle, $`\delta `$, is plotted on Fig.6 ($`T_0=10^0C`$, $`P_0=1atm=101.3kPa`$, $`\nu =20000cm^1`$ wave-number for the light and $`\lambda =6.5K/km`$ lapse-rate).

These results are in excellent agreement with the one given by Thomas and Joseph  and the report of the U.S. Naval Observatory (1993) . From Fig. 6 we learn that the deviation angle is usually quite small, and becomes important only when viewing objects in the vicinity of the optical horizon. The deviation angle increases sharply for very small inclination angles. This effect is responsible for the flattened shape of the rising (or setting) Sun and Moon, and also for the fact that these heavenly objects appear more flattened at their bottom. It is also interesting to note here that the angular extent of the Sun and full-Moon is 32 arcmin, and the maximum deviation obtained for standard conditions is 34.5 arcmin. This leads us to the observation that the Sun or Moon is visible even when in reality it is below the geometrical horizon.

## IV Computing the flatness

Once we numerically determined the $`\delta _{dev}`$ deviation angle as a function of the $`\delta `$ apparent inclination angle, it is easy to characterize the apparently flat rim of the setting Sun. Atmospheric refraction influences only the $`\mathrm{\Delta }_v`$ vertical angular extent of the Sun (or Moon), which becomes thus smaller than the $`\mathrm{\Delta }_h`$ horizontal extent. The asymmetry ratio (or the flatness) for the rim of the Sun can be described by the

$$\alpha =\frac{\mathrm{\Delta }_h}{\mathrm{\Delta }_v},$$
(44)
ratio (of course $`\alpha [1,\mathrm{})`$). Since both the Sun and the full-Moon are normally visible under a $`\mathrm{\Delta }_0=32arcmin`$ angular extent, we have $`\mathrm{\Delta }_h=\mathrm{\Delta }_0`$. The value of $`\mathrm{\Delta }_v`$ can be derived after numerically computing the $`\delta `$ apparent inclination angle as a function of the $`\delta _r=\delta +\delta _{dev}`$ real inclination, i.e. $`\delta =F(\delta _r)`$. If the apparent inclination angle for the bottom of the Sun is $`\delta _0`$, corresponding to a $`\delta _{r0}`$ real inclination angle, than

$$\mathrm{\Delta }_v(\delta _0)=F(\delta _{r0}+\mathrm{\Delta }_0)\delta _0,$$
(45)
and we get:

$$\alpha (\delta _0)=\frac{\mathrm{\Delta }_0}{F(\delta _{r0}+\mathrm{\Delta }_0)\delta _0}$$
(46)
For fixed observation altitude and meteorological conditions the $`\alpha _c`$ maximal possible flatness corresponds to the situation when the bottom of the Sun touches the horizon. This happens for a $`\delta _0=\delta _c`$ critical inclination angle.

By decreasing the $`\delta `$ angle in small steps, and applying the previously described methods for computing the light-ray trajectory, both the $`\delta =F(\delta _r)`$ function and $`\delta _c`$ is numerically determined. In our calculations we have chosen to decrease $`\delta `$ in steps of $`0.01^0`$.

## V Results

Applying the methods presented in the previous sections and our optical atmosphere model, we systematically computed the $`\alpha _c`$ asymmetry ratio for the rim of the setting Sun as a function of observational altitude and meteorological conditions parameterized by pressure and temperature. We have also shown that the value of $`\alpha `$ is insensitive to details of the considered optical atmosphere model, proving the stability of our results. In the following we detail our findings. For all the calculations we considered monochromatic light with wave-number of $`\nu =20000cm^1`$, corresponding to the green color. If otherwise not specified we considered the height of the troposphere $`z_t=14km`$ and the height of the stratosphere as $`z_s=50km`$.

### A Stability regarding the considered optical atmosphere model

As discussed in section II., for the construction of the refractive index profile we borrowed results from the U.S. Standard Atmosphere model. The main parameters needed by us is the $`z_t`$ height of the troposphere, the $`z_s`$ height of the stratosphere and the $`\lambda `$ lapse-rate.

It is evident, that the exact value of $`z_s`$ does not much influence our results, since in the stratosphere the refractive index is already very close to $`1`$. As a first step we studied thus the influence of $`z_t`$ on the $`\alpha _c`$ asymmetry ratio. We considered normal conditions with $`T_0=0^0C`$, $`P_0=1atm`$, observations at sea-level ($`z_0=0`$) and $`\lambda =6.5K/km`$ standard lapse-rate. As illustrated on Fig. 7, the value of $`z_t`$ (in a reasonable range) has no significant influence (note the scale on the vertical axis). The value of the lapse-rate has already a more noticeable effect on $`\alpha _c`$ (Fig. 8), however this variation is also quite small for the practically important fluctuations around the standard $`\lambda =6.5K/km`$ value. We conclude thus, that our results are quite stable regarding the details of the considered optical atmosphere model.

### B Asymmetry as a function of the inclination angle

We computed the observable $`\alpha `$ flatness as function of the apparent inclination angle of the Sun’s bottom. Results for $`T_0=0^0C`$, $`P_0=1atm=101.3kPa`$, $`\lambda =6.5K/km`$ and observations at sea-level are presented on Fig. 9. As emphasized before, the asymmetric rim of the setting (or rising) Sun becomes evident only for very small $`\delta _0`$ values, when the Sun is close to the horizon. For these normal parameters we get that $`\alpha _c`$ is around $`1.2`$.

### C Flatness as a function of observation height

Let us presume now that the observer is at height $`z_0`$ above the sea-level, and there is no obstacle in the direction of the horizon, which is at sea-level (i.e. we are above a vast ocean). It is obvious that from higher altitude the $`\delta _c`$ critical angle will be smaller (becomes negative) and the deviation angle increases more sharply in the neighborhood of $`\delta _c`$. This leads us immediately to the conjecture that the observed flatness should be larger. Taking the $`T_0=0^0C`$ and $`P_0=1atm`$ normal atmospheric conditions at sea-level, the standard $`\lambda =6.5K/km`$, lapse-rate we can compute the $`\alpha _c`$ maximal observable flatness as a function of observation height. Results in this sense are plotted in Fig. 10. As expected, $`\alpha _c`$ increases with $`z_0`$. For these normal conditions we get that from the top of a $`5km`$ height mountain we would observe an $`\alpha _c1.5`$ flatness and from a commercial flight at $`10km`$ height at sunset we would detect an $`\alpha _c1.7`$ maximal flatness. For altitudes above $`30km`$, or observations made from a space-shuttle we can get extreme values for $`\alpha _c`$ up to $`2.5`$.

### D Influence of temperature

The $`T_0`$ temperature measured at sea-level influences in an important manner the observed flatness. Computing the temperature dependence of $`\alpha _c`$ for observations at sea-level, normal $`P_0=1atm`$ atmospheric pressure and standard lapse-rate we get the values presented in Fig. 11. While for a $`T_0=30^0C`$ temperature $`\alpha _c`$ is around $`1.1`$ for $`T_0=40^0C`$ it becomes $`1.3`$, and increases more and more sharply for lower temperatures (arctic conditions).

### E Influence of atmospheric pressure

Increasing the $`P_0`$ pressure at sea-level results the increase of the observable $`\alpha _c`$ flatness. For the reasonable and measurable $`P_0`$ values, $`T_0=0^0C`$, observation at sea-level and standard lapse-rate the variation is almost linear (Fig.12). For an extreme $`P_0=125kPa`$ pressure one can detect an asymmetry ratio of $`1.3`$.

## VI Experiments

By simple experiments it is relatively easy to measure the flatness of the setting (or rising) Sun. We considered photo and video experiments, and analyzed the pictures as a function of the inclination angle of the Sun. With an appropriate filter and calibrated eye-piece, telescope observations were also possible. In order to obtain usable pictures with a nicely visible rim we had to ensure a properly adjusted light intensity. Appropriate filters, specific atmospheric conditions and usually small $`\delta _0`$ inclination angles lead us to usable pictures. When the Sun is close to the horizon, the inclination angle of the Sun can be directly determined by analyzing the taken picture. We can use the fact that the horizontal angular extent of the Sun always corresponds to $`32arcmin`$. The distance of the Sun from the horizon can then be compared with the horizontal extent of the setting Sun and the inclination angle results directly from the picture. Making experiments with the setting or rising Moon is more complicated, since one needs a full-Moon for this, a work during the night and we have to deal with pictures where the horizon is not clearly visible.

First, we have taken several series of pictures in South-Bend (Indiana, altitude $`100m`$) both in winter and late-spring, studying very different temperature conditions. As an immediate confirmation of our theoretical results, from these pictures it was obvious that for the same inclination angle the flatness is bigger in winter, i.e. for lower $`T_0`$ values. Two series of digitized and appropriately enhanced pictures are visible on the web-page dedicated to this study . Results from these pictures, in comparison with the expected (computed) flatness is presented on Fig. 13. The first set of pictures is a sunrise in winter. The mean temperature during the sunrise was $`5^0C`$, and the atmospheric pressure was $`103kPa`$. It is important to note that during this sunrise the temperature remained approximately constant. The second set of picture was made in late-spring for a sunset. During this sunset the mean-temperature was $`12^0C`$, the atmospheric pressure $`98kPa`$ and the temperature dropped detectably during the time the pictures were made. Results for the flatness as a function of the Sun’s inclination angle (bottom of the rim) is plotted with filled circles and triangles for the winter and spring series, respectively.

The theoretical curves were constructed for the mentioned mean temperature, pressure, observation height of $`100m`$ with the horizon at the same altitude and a standard lapse-rate. These results are plotted on Fig. 13 by a continuous and dashed line for the winter and spring conditions, respectively. As observable from Fig. 13 for the winter conditions the measured $`\alpha `$ values are in acceptable agreement with the one given by our theory. However, for the spring series the expected $`\alpha `$ values are higher, and decreasing slower as a function of the $`\delta _0`$ inclination angle, than the measured data. A reason to account for this sharper trend is the decreasing temperature measured during the sunset. Taking this effect into account in the calculations, would definitely result in a trend closer to the observed one. This temperature variation does not account however, for the constantly lower values measured for $`\alpha `$. The only arguments we can give in this sense is that probably the atmosphere at the time of this measurement had a refractive index profile different from the one proposed in our model, behaving in a non-standard manner.

A second set of experiments were realized by video-filming sunsets in Cluj (Romania). The obtained continuous set of picture allowed us to follow-up more precisely the flatness as a function of the inclination angle. The inclination angle was calculated from the images by the same method as in photographs, i.e. by comparing the height between the bottom of the Sun and horizon with the horizontal extent of the setting Sun. Since Cluj is not a flat region like South-Bend, we had to be careful in choosing the observation point, and to determine also the altitude of the optical horizon. Results from video-recording in comparison with theoretical expectations (corresponding to the appropriate atmospheric conditions) are plotted in Fig. 14. For this measurement a quite fair agreement between theoretical and experimental data is achieved. Since the lapse-rate for the theoretical prediction was taken from the standard atmosphere model, the slightly greater $`\alpha `$ values calculated by us are explainable.

## VII Pictures and movies from the Internet

Sunset or sunrise (moon-set or moon-rise) is usually a spectacular phenomena, and thus it is a favorite theme for professional and amateur photographers. The Internet is full of beautiful and useful pictures in this sense. Many of these pictures are taken in extreme conditions (arctic environment, airplanes, space-shuttle or high mountains), offering us an excellent possibility to check our calculations under these conditions, too. Moreover, on the web there are also interesting, freely downloadable or public domain movies exemplifying how the observed flatness increases in the neighborhood of the horizon.

We performed an extensive search on the Internet and collected non-copyrighted materials about sunset, sunrise, moon-set and moon-rise. Since their presentation in the context of this paper is impossible we classified and stored them on the web-page  dedicated to this study. The interested reader can browse this collection and convince himself or herself that these pictures support our theoretical results for the estimated flatness. In agreement with our expectations we found that for a sunrise (or sunset) viewed from a space-shuttle $`\alpha _c`$ should be of the order $`22.5`$, and in arctic environment $`\alpha _c`$ increases up to values of $`1.3`$. For most of the everyday, usual and low-altitude photos one finds $`\alpha _c1.11.2`$. Pictures taken from commercial airplane yield $`\alpha _c1.5`$. We also found a picture-series for a sunset over the ocean, where the rim of the Sun is nicely visible. The scenario presented in this photo is again in agreement with our theoretical predictions.

The movies stored on our web-page will also convince the reader about the sharp variation of $`\alpha `$ as a function of the $`\delta _0`$ deviation angle. A public domain movie, showing a moon-set viewed from the space-shuttle illustrates observations from high altitudes.

## VIII The Green flash

Although the explanation and the study of the green flash is not the aim of the present paper, we briefly discuss here how this phenomenon can be understood and studied through our methods.

From equation (1) it is obvious that atmospheric refraction depends on the frequency (color) of the light-ray, a phenomenon called dispersion. Differently colored light-rays coming from a source emitting a continuous spectra will suffer different deviations. More strongly will bend the light-rays with bigger frequencies, leading to a separation of the colors in the observed Sun. Since the high-frequency visible components (corresponding to violet and blue colors) are strongly scattered by the atmosphere, the green component reaching directly the observer will suffer the strongest bending. When the Sun disappears bellow the horizon this component will be observed lastly, leading to a green-flash on the horizon. Since for standard atmospheric conditions the dispersion from atmospheric refraction is tinny, the effect is hardly observable for usual and everyday conditions. Green flashes become observable exactly under those conditions under which the flatness is more accentuated. This means high observation altitudes, low temperature and high pressure. As an illustration of this, we have just learned that for airplane pilots it is a quite usual phenomenon. Green flashes are often seen in non-standard optical atmosphere, when mirages appear. This is generally the case when layers of air with strongly different temperatures are in contact.

For a nice presentation and discussion on green flashes the interested reader should consult the excellent home-page of A.T. Young .

## IX Sunset simulation program

We also created a computer-code which simulates the sunset on the computer-screen. By using the discussed simulation method, the program computes and visualizes the rim of the Sun during a sunset. After fixing the atmospheric conditions (temperature and pressure) and observation height in the menu, one can follow up how a sunset might look like in our optical atmosphere. The program runs under Windows environment and one can freely download the executable from the web-page accompanying this study .

## X Conclusions

Atmospheric refraction is responsible for the asymmetric rim of the setting (or rising) Sun. Three different methods yielding the same results were presented here to compute the path of a light-ray in an optical atmosphere model, where the refraction index varies continuously as a function of altitude. By determining the deviation angle between the apparent and real inclination of a point-like light-source which is at a large distance from the Earth, we were able to compute the $`\alpha `$ asymmetry ratio for the rim of the setting Sun. We investigated $`\alpha `$ as a function of the inclination angle, observation altitude and atmospheric conditions. We found that the maximal flatness obtained in the vicinity of the horizon increases as a function of observation altitude and pressure, and decreases as the temperature increases. We found that $`\alpha `$ is rather insensitive to the fine-details of the considered optical atmosphere model, which makes our results robust. The maximal flatness observable under normal conditions at sea-level is around $`1.2`$. At extremely low temperatures ($`40^0C`$) one can observe values up to $`1.3`$ and for high altitude observations (space-shuttle) one can get an asymmetry ratio of $`2.5`$. By simple experiments and pictures from the Internet we illustrate and prove our theoretical predictions. The methods presented here can be effectively used to study atmospheric refractions in non-standard atmospheric conditions as well. In this manner one can study mirages or the green flash phenomenon. A freeware computer program created by us and downloadable from the web-page accompanying this study simulates the sunset for arbitrary atmospheric conditions and observation altitude.

## XI Acknowledgments

We thank the Bergen Computational Physics Laboratory in the framework of the European Community - Access to Research Infrastructure of the Improving Human Potential programme for access to their supercomputer facilities. The work of Z. Neda was sponsored by the Sapientia foundation. We thank I. Albert and T. Néda for useful discussions and help in the video and photo experiments.