text:
# Theoretical analysis of magnetic Raman scattering in La2CuO4: two-magnon intensity with the inclusion of ring exchange

The magnetic properties of La2CuO4 have been the subject of many detailed investigations over the last decade. Understanding this undoped parent compound of high temperature superconducting cuprates is a precondition for the many theories which describe metallic cuprates by doping carriers into a layered antiferromagnet. The conventional starting point for undoped cuprates is the two-dimensional (2D) spin-1/2 Heisenberg model with nearest-neighbor (nn) exchange interaction $`J`$ . Despite the substantial progress on the theory of the 2D Heisenberg antiferromagnet , some of the experimental facts for La2CuO4 have clearly demonstrated that a complete description of the magnetic excitations requires additional physics not contained in the 2D Heisenberg model with $`J`$ only. Examples include the asymmetric lineshape of the two-magnon Raman intensity  or the infrared optical absorption .

The importance of an additional ring (plaquette) exchange coupling for La2CuO4 recently found direct experimental support from the observed dispersion of the spin-waves along the magnetic Brillouin zone boundary . A fit of the experimental results  using the theoretical spin-wave dispersion, which consistently includes quantum renormalization effects , has provided accurate values for the nn exchange integral $`J=1720K`$ and the ring-exchange coupling $`J_{\mathrm{}}=0.24J`$. The new estimate for $`J`$ corrects previous values which have been used in past years and which were consistently 10 to 15% lower. Also the value for $`J_{\mathrm{}}`$ must be considered surprisingly large, but the spin stiffness and the Néel temperature calculated with the new parameter set obtained in Ref.  were in excellent agreement with the experimental data and thereby confirmed the deduced exchange coupling parameters. The value of the ring exchange $`J_{\mathrm{}}`$ also agreed with the strong-coupling expansion studies of the three-band Hubbard model in the parameter range relevant for CuO2 planes .

The Heisenberg model with $`J`$ only was previously found insufficient to describe the experimentally observed asymmetry and width of the Raman spectra of undoped cuprate compounds . The Raman spectrum of the Heisenberg model was investigated within the ladder approximation for the magnon-magnon scattering vertex many years ago  and gave an almost symmetric narrow peak located at $`E_{\mathrm{max}}=3.3J`$. It was shown  that magnon self-energy and vertex corrections to the ladder approximation as well as 4-magnon scattering contributions are small and negligible. Numerically, the problem was investigated by series expansions around the Ising limit , exact diagonalizations on small clusters , and QMC calculations . The results of these studies led to a peak position and a two-magnon lineshape which were very close to the spin-wave results. Furthermore, the value of the exchange integral $`J=1440`$K for La2CuO4 as extracted from the position of the two-magnon Raman peak appears too small in comparison to early neutron scattering results for the spin-wave spectrum ($`J=1650`$K), see Ref.  and references therein.

The asymmetric lineshape of the two-magnon Raman intensity  has led to proposals that spin-phonon interactions , resonant phenomena , purely fermionic contributions , or cyclic ring exchange  need to be included beyond the nn Heisenberg model. A possible importance of ring exchange for Raman scattering was conjectured also from numerical calculations  which showed that a finite $`J_{\mathrm{}}`$ gives rise to additional high energy contributions in the Raman intensity – yet with little spectral weight.

The new, accurate estimates for the exchange couplings therefore demand a theoretical reanalysis of two-magnon scattering in La2CuO4, which we investigate in the present paper with an emphasis on the effects of ring exchange. The question that we address is to what extent the new values for the exchange parameters – including the finite $`J_{\mathrm{}}`$ – consistently describe the two-magnon intensity in La2CuO4 and whether the Heisenberg model with ring exchange alone is sufficient to explain the experimentally observed lineshape of the Raman intensity. To answer this question we calculate the magnon-magnon interaction vertex in the presence of the ring exchange term and then solve the ladder equation for the Raman scattering vertex.

We start from the Heisenberg model on a square lattice with ring-exchange coupling

$`H`$ $`=`$ $`J{\displaystyle \underset{i,\delta }{}}𝐒_i𝐒_{i+\delta }+J^{}{\displaystyle \underset{i,\delta ^{}}{}}𝐒_i𝐒_{i+\delta ^{}}+J^{\prime \prime }{\displaystyle \underset{i,\delta ^{\prime \prime }}{}}𝐒_i𝐒_{i+\delta ^{\prime \prime }}`$ (3)
$`+J_{\mathrm{}}{\displaystyle \underset{ijkl}{}}[(𝐒_i𝐒_j)(𝐒_k𝐒_l)`$
$`+(𝐒_i𝐒_l)(𝐒_k𝐒_j)(𝐒_i𝐒_k)(𝐒_j𝐒_l)]`$
where $`J,J^{}`$, and $`J^{\prime \prime }`$ are the first ($`\delta `$), second ($`\delta ^{}`$) and third ($`\delta ^{\prime \prime }`$) nearest neighbor in-plane exchanges. We use the Dyson-Maleev representation for the spin operators

$$\begin{array}{c}S_i^+=\sqrt{2S}a_i,S_i^z=Sa_i^{}a_i\hfill \\ S_i^{}=\sqrt{2S}(a_i^{}\frac{1}{2S}a_i^{}a_i^{}a_i)\hfill \end{array}\}iA,$$
(4)
$$\begin{array}{c}S_j^+=\sqrt{2S}b_j^{},S_j^z=S+b_j^{}b_j\hfill \\ S_i^{}=\sqrt{2S}(b_i\frac{1}{2S}b_i^{}b_ib_i)\hfill \end{array}\}jB$$
(5)
where $`A`$ and $`B`$ denote the two sublattices of the antiferromagnet; $`a_i^{},a_i,`$ and $`b_j^{},b_j`$ are Bose operators and $`S=1/2`$. With the Bogoliubov transformation

$`a_𝐤`$ $`=`$ $`u_𝐤\alpha _𝐤+v_𝐤\beta _𝐤^{}`$
$`b_𝐤`$ $`=`$ $`u_𝐤\beta _𝐤+v_𝐤\alpha _𝐤^{}`$
and by introducing the “coherence factors”

$`u_𝐤`$ $`=`$ $`[(A_𝐤+E_𝐤)/(2E_𝐤)]^{1/2}`$
$`v_𝐤`$ $`=`$ $`[(A_𝐤E_𝐤)/(2E_𝐤)]^{1/2}`$
we diagonalize the quadratic part of the Hamiltonian and obtain the spin-wave spectrum (for details see )

$`E_𝐤`$ $`=`$ $`\sqrt{A_𝐤^2B_𝐤^2},`$ (6)
$`A_𝐤`$ $`=`$ $`4S[J\gamma J^{}\gamma ^{}(1\nu _𝐤^\delta ^{})J_{\mathrm{}}S^2(\gamma _0^{\mathrm{}}+\gamma _\delta ^{}^{\mathrm{}}\nu _𝐤^\delta ^{})]`$ (8)
$`4J^{\prime \prime }S\gamma ^{\prime \prime }(1\nu _𝐤^{\delta ^{\prime \prime }}),`$
$`B_𝐤`$ $`=`$ $`4S(J\gamma J_{\mathrm{}}S^2\gamma _\delta ^{\mathrm{}})\nu _𝐤^\delta ,`$ (9)
with the momentum dependent coefficients

$`\nu _𝐤^\delta `$ $`=`$ $`(\mathrm{cos}k_x+\mathrm{cos}k_y)/2;\nu _𝐤^\delta ^{}=\mathrm{cos}k_x\mathrm{cos}k_y;`$ (10)
$`\nu _𝐤^{\delta ^{\prime \prime }}`$ $`=`$ $`(\mathrm{cos}2k_x+\mathrm{cos}2k_y)/2.`$ (11)
Note that the quantum renormalization factors $`\{\gamma _\delta \}`$ in (6) take into account the renormalization of the magnon spectrum due to quartic terms. The resulting Hamiltonian reads

$`H={\displaystyle \underset{𝐤}{}}E_𝐤(\alpha _𝐤^{}\alpha _𝐤+\beta _𝐤^{}\beta _𝐤)`$ (12)
$`+`$ $`{\displaystyle V_{𝐤_1𝐤_2;𝐤_3𝐤_4}^{ijmn}}:R_{𝐤_1,i}^{}R_{𝐤_2,j}^{}R_{𝐤_3,m}R_{𝐤_4,n}:\delta _{𝐤_1+𝐤_2,𝐤_3+𝐤_4}`$ (13)
where we used the vector notation

$$R_𝐤=(\alpha _𝐤,\beta _𝐤^{})$$
(14)
and $`i,j,m,n=1,2.`$ The quartic part in Eq. (12) is normal ordered, since all “Hartree-Fock” renormalizations are already absorbed in the quadratic part of the Hamiltonian, cf. Ref. .

We use the effective Loudon-Fleury  Hamiltonian in $`B_{1g}`$ geometry for the coupling of the incoming and outgoing photons to the localized nn spins which are involved in the two-spin flip Raman process

$$H_R=\mathrm{\Lambda }\underset{j}{}𝐒_j(𝐒_{j+\delta _y}𝐒_{j+\delta _x});$$
(15)
$`\delta _{x(y)}`$ are unit vectors in the directions $`x,y`$ respectively, and $`\mathrm{\Lambda }`$ is a coupling constant which includes the electric field vectors of the two photons. For the calculation of the Raman light scattering intensity in the ladder approximation for repeated magnon-magnon scattering processes only the vertex $`V_{\mathrm{𝐤𝐤}^{};𝐤^{}𝐤}^{(4)}=V_{\mathrm{𝐤𝐤}^{};𝐤^{}𝐤}^{1,2,1,2}`$ is needed . The complicated vertex function $`V_{𝐤_1𝐤2;𝐤_3𝐤_4}^{1,2,1,2}`$ was calculated with the help of computer algebra. After Bogoliubov transformation the ring-exchange term in Eq. (1) gives 1057 (!) contributions with different combinations of $`\alpha `$ and $`\beta `$ operators. The general result for the vertex $`V_{𝐤_1𝐤_2;𝐤_2𝐤_4}^{(4)}`$ is rather involved, however, for equal momenta it simplifies to

$$V_{\mathrm{𝐤𝐤}^{};𝐤^{}𝐤}^{(4)}=\frac{B_{𝐤𝐤^{}}(A_𝐤A_𝐤^{}+E_𝐤E_𝐤^{})A_{𝐤𝐤^{}}B_𝐤B_𝐤^{}R_{\mathrm{}}}{2E_𝐤E_𝐤^{}}$$
(16)
where the effects of the ring-exchange coupling are included in the magnon spectrum (6) and in the additional vertex contribution

$`R_{\mathrm{}}`$ $`=`$ $`8J_{\mathrm{}}[B_𝐤B_𝐤^{}(4+\nu _𝐤^\delta ^{}+\nu _𝐤^{}^\delta ^{}+\nu _{𝐤𝐤^{}}^\delta ^{})+`$ (17)
$``$ $`2A_𝐤^{}B_𝐤(\nu _𝐤^{}^\delta +\nu _{𝐤,𝐤𝐤^{}}^{\mathrm{}})2A_𝐤B_𝐤^{}(\nu _𝐤^\delta +\nu _{𝐤𝐤^{},𝐤^{}}^{\mathrm{}})`$ (18)
$`+`$ $`2A_𝐤A_𝐤^{}(\nu _{𝐤𝐤^{}}^\delta +\nu _{𝐤,𝐤^{}}^{\mathrm{}})+2E_𝐤E_𝐤^{}(\nu _{𝐤𝐤^{}}^\delta +\nu _{𝐤,𝐤^{}}^{\mathrm{}})]`$ (19)
and $`\nu _{𝐤,𝐤^{}}^{\mathrm{}}=(\mathrm{cos}k_x\mathrm{cos}k_y^{}+\mathrm{cos}k_y\mathrm{cos}k_x^{})/2`$.

The first two terms in the numerator of Eq. (16) coincide with the corresponding result for the 2D nn Heisenberg antiferromagnet  but with renormalized coefficients $`A_𝐤,`$ $`B_𝐤`$ , and the renormalized magnon dispersion $`E_𝐤`$ in (6). The time-ordered response function of two-magnon Raman scattering is given by

$`G(\omega )={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑te^{i\omega t}T[H_R(t)H_R(0)]`$where $`H_R(t)`$ is the Heisenberg representation of $`H_R,`$ Eq. (15). The Raman light scattering intensity is then determined by $`I(\omega )=(1/\pi )\mathrm{Im}G(\omega )`$. In the ladder approximation we obtain the equations :

$`G(\omega )=i\mathrm{\Lambda }^2{\displaystyle \frac{d\omega ^{}}{2\pi }\underset{𝐤}{}\frac{A_𝐤f_𝐤}{E_𝐤}}`$ $`G_\alpha ^{(0)}(𝐤,\omega +\omega ^{})`$ (21)
$`G_\beta ^{(0)}(𝐤,\omega ^{})\mathrm{\Gamma }_𝐤(\omega ),`$
$`\mathrm{\Gamma }_𝐤(\omega )={\displaystyle \frac{f_𝐤A_𝐤}{E_𝐤}}i{\displaystyle \frac{d\omega _1}{2\pi }\underset{𝐤_1}{}}`$ $`V_{\mathrm{𝐤𝐤}_1;𝐤_1𝐤}^{(4)}G_\alpha ^{(0)}(𝐤_1,\omega +\omega _1)`$ (23)
$`G_\beta ^{(0)}(𝐤_1,\omega _1)\mathrm{\Gamma }_{𝐤_1}(\omega )`$
where $`G_{\alpha ,\beta }^{(0)}(𝐤,\omega )=1/(E_𝐤\omega i0^+)`$ and $`f_𝐤=\mathrm{cos}k_x\mathrm{cos}k_y.`$ The result reads

$$G(\omega )=\frac{\mathrm{\Lambda }^2L_2(\omega )}{D(\omega )}\left\{1+\frac{J}{2}\left[L_0(\omega )\frac{L_1^2(\omega )}{L_2(\omega )}\right]\right\},$$
(24)
$`D(\omega )`$ $`=`$ $`1+(J2J_{\mathrm{}}S^2){\displaystyle \frac{L_2(\omega )}{2}}+(J+6J_{\mathrm{}}S^2){\displaystyle \frac{L_0(\omega )}{2}}`$ (26)
$`+(J2J_{\mathrm{}}S^2)(J+6J_{\mathrm{}}S^2){\displaystyle \frac{L_0(\omega )L_2(\omega )L_1^2(\omega )}{4}}`$
where

$`L_m(\omega )`$ $`=`$ $`i{\displaystyle \frac{d\omega ^{}}{2\pi }\underset{𝐤}{}\frac{f_𝐤^2A_𝐤^m}{E_𝐤^m}G_\alpha ^{(0)}(𝐤,\omega +\omega ^{})G_\beta ^{(0)}(𝐤,\omega ^{})}`$ (27)
$`=`$ $`{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{f_𝐤^2A_𝐤^m}{E_𝐤^m}}{\displaystyle \frac{1}{\omega 2E_k}}.`$ (28)
For $`J_{\mathrm{}}=0`$ the above formulas reduce to the known results of Refs..

For the numerical calculations we use the parameter set determined in Ref.  from the fit to the spin wave dispersion

$$J=151.9\text{meV},J^{}=J^{\prime \prime }=0.025J,J_{\mathrm{}}=0.24J.$$
(29)
For these exchange couplings the renormalization parameters follow as

$`\gamma `$ $`=`$ $`1.158;\gamma ^{}=0.909;\gamma ^{\prime \prime }=0.852`$ (30)
$`\gamma _0^{\mathrm{}}`$ $`=`$ $`2.220;\gamma _\delta ^{\mathrm{}}=1.971;\gamma _\delta ^{}^{\mathrm{}}=1.721`$ (31)
Nevertheless we first explore the $`J_{\mathrm{}}`$ dependence of the two-magnon Raman intensity as shown in Fig. 1. On varying $`J_{\mathrm{}}`$ the renormalization parameters (30) necessarily have to be recalculated each time. For $`J^{}=J^{\prime \prime }=J_{\mathrm{}}=0`$ the result in Fig. 1 coincides with that obtained in Refs. . In a first step we turn on $`J^{}`$ and $`J^{\prime \prime }`$ and observe a shift of the two-magnon peak position to lower frequencies and a further downward shift when also the the ring-exchange coupling is added. The magnitude of the peak shift due to the ring exchange is quite strong. In addition, due to the finite $`J_{\mathrm{}}`$ the spectrum near the upper edge of the two-magnon continuum develops a “foot” structure. For comparison, we plot also the Raman scattering intensity for a twice larger ring exchange $`J_{\mathrm{}}=0.48J`$ where the above effects are even more pronounced. The position of the peak estimated with $`J_{\mathrm{}}=0.24J`$ is $`E_{\mathrm{max}}=2.71J`$ and with the absolute values in (29) is therefore excellent agreement with the experimental result  $`E_{\mathrm{max}}=2.65J=3150\mathrm{c}\mathrm{m}^1.`$ If we turn off the vertex contribution $`R_{\mathrm{}}`$ in Eq. (16) and therefore keep only the influence of the ring exchange in the spin-wave dispersion, we obtain $`E_{\mathrm{max}}=2.86J`$ which is still much smaller than the peak position without ring exchange ($`E_{\mathrm{max}}=3.15J`$). Therefore, the modification of the spectrum due to a finite $`J_{\mathrm{}}`$ is the major origin for the peak shift to lower energies.

In Fig. 2 we compare the Raman lineshape, obtained for the parameter set (29) with the experimental data on La2CuO4 from Ref. . Although the precise peak position is in excellent agreement with the data, its width is only slightly influenced by ring-exchange coupling. The foot structure near the upper edge of the spectrum does result from ring exchange and therefore indeed leads to an asymmetry of the lineshape, its weight, however, is too small to account for the overall linewidth and lineshape. Therefore, other processes must contribute which are not described by the Heisenberg model with ring exchange alone. It was proposed, that the damping due to spin-phonon coupling  is responsible for the peak width which in the current calculation is underestimated by roughly a factor 2. In order to incorporate extrinsic sources for damping beyond the extended Heisenberg Hamiltonian Eq. (1), we have introduced a small damping $`\mathrm{\Gamma }`$ into the denominators of $`L_m(\omega ),`$ Eq. (28). As a reasonable choice $`\mathrm{\Gamma }`$ is assumed constant over the entire Brillouin zone, because it is predominantly the zone boundary magnons which determine the shape of the two-magnon peak. We plot the result of this calculation for $`\mathrm{\Gamma }=0.13J`$ in Fig. 2. One can see that while the width of the experimentally observed spectra can be fitted in such a way, the asymmetry of the spectrum is not accounted for. Since the Raman spectra do depend on the incoming light frequency, the two-magnon scattering is identified as a resonant process and therefore an additional source for the strong asymmetry of the lineshape is expected to arise from the coupling to the charge degrees of freedom .

In summary, we have investigated the effect of ring exchange on two-magnon Raman scattering in the 2D Heisenberg antiferromagnet. Using the previously fixed exchange coupling parameter set we find in a fit-parameter free calculation an excellent agreement with the position of the two-magnon peak in La2CuO4. This reconfirms the magnitude of the ring-exchange coupling $`J_{\mathrm{}}=0.24J`$ in this material. Ring-exchange creates additional high energy spectral weight in the Raman intensity, but $`J_{\mathrm{}}`$ alone is insufficient to explain the overall asymmetric lineshape and linewidth.

It is a pleasure to thank T. Kopp and T. Nunner for insightful discussions. We acknowledge support through Sonderforschungsbereich 484 of the Deutsche Forschungsgemeinschaft.


format: text

text:
# Air Shower Simulations in a Hybrid Approach using Cascade Equations

## I Introduction

The field of highest energy cosmic rays is an exciting subject with many open questions: What is the nature of the primary cosmic ray? What are the highest energies? What are possible sources/acceleration mechanisms? Is there clustering of events? Is there a GZK cutoff due to the microwave background? Ongoing (HIRES, AGASA) and future (Auger, OWL, EUSO) cosmic ray experiments aim to shed light on these mysteries.

At these high energies, direct measurement of the primary cosmic ray is impossible due to the low flux, which is only of order one event per square-kilometer per century at the highest observed energies. But cosmic rays initiate showers in the atmosphere, a cascade of secondary particles from collisions with air molecules, which themselves collide and so on. Experiments measure these air showers and reconstruct from their properties information about the primary ray at the beginning of the reaction.

Air shower models are of crucial importance for the reconstruction of the energy and primary type. The straightforward approach is to model each possible interaction of hadrons, leptons and photons with air molecules, and trace all secondary particles. At high energies this leads quickly to unpractical computation times, since the time grows with the number of particles in the shower and therefore increases rapidly with the primary energy. A shower of even $`10^{19}eV`$ has more than $`10^{10}`$ particles at its maximum and would take months to compute. The thinning algorithm proposed by Hillas Hillas1985  tries to solves this problem: below a fraction $`f_{\mathrm{thin}}`$ of the primary energy only a small sample of the particles is actually followed in detail, attributing them a higher weight. This procedure introduces artificial fluctuations and one must compromise between these and computation time.

People have tried to overcome these difficulties by defining systems of (mostly one-dimensional) transport equations which describe air showers Kalmykov:1986ii ; GaisserBook . The numerical solutions of these equations can then be combined with a Monte-Carlo in order to account for natural fluctuations due to the first interactions and for lateral spread of low-energy particles Dedenko1968 ; Lagutin:1999xh . This is the principle of the hybrid method. Another realization of the hybrid approach is to use shower libraries in which presimulated longitudinal profiles are combined to compute the one-dimensional properties of air showers Gaisser1997 ; Alvarez-Muniz:2002ne .

In a recent paper Bossard:2000jh , a new approach to an old idea was introduced: the method of cascade equations, which allows one to compute longitudinal characteristics of air showers numerically in a very short time.

In this paper we introduce the further development of this approach. Traditional Monte-Carlo methods are combined with cascade equations in a hybrid approach. This allows to construct an efficient model which accounts not only for natural fluctuations due to the first interactions but also for the correct 3-dimensional spreading of low-energy secondary particles. In reasonable computing time it is possible to calculate longitudinal profiles and lateral distribution functions with detailed knowledge about particle momenta and arrival-times, which are reliable on an event-by-event basis.

## II Hybrid approach to air shower modeling

The solution of one-dimensional cascade equations cannot account for natural fluctuations or the lateral spread of particles. The fluctuations can be, as already suggested, solved by doing the first interactions up to a certain fraction $`f`$ of primary energy $`E_0`$, in a classical Monte-Carlo approach, where each collision is treated individually by the chosen hadronic model. All secondary particles below the critical energy $`fE_0`$ are not followed further on, but taken to be initial conditions for the hadronic cascade equations.

The one-dimensional cascade equations are only valid for large enough energies that the emission angles of secondaries can be neglected, and the whole problem can be treated longitudinally. Therefore the lateral spreading of particles cannot be treated in this approach and we return to the Monte Carlo method for the low energy regime. At $`E_{\mathrm{min}}`$, the output of the cascade equations – the number of particles at certain depths and energies – is used as a source function for the Monte Carlo approach, by creating single particles according to the source function and following them individually. This method is able to reproduce the lateral spread of secondary particles even though it is neglected for collisions with $`E>E_{\mathrm{min}}`$.

Fig. 1 illustrates schematically the hybrid approach.

Throughout this paper, which concentrates on establishing the validity of the technique, we use QGSJET QSJET  as high energy hadronic model. Low energy hadrons are treated by GHEISHA GHEISHA . The electromagnetic part is calculated by the EGS4-code system EGS4 . The bremsstrahlung and $`e^\pm `$ pair-production by muons is done using the GEANT3.21 code GEANT . At this stage we neglect the LPM effect, muon-pair production and photo-nuclear reactions. These effects are essentially negligible for hadron primaries Alvarez-Muniz:2002ne ; CORSIKA ; Cillis:1998hf . The program embodying this method which was used for the calculations presented in this paper, SENECA, is available for public use at http://cosmo.nyu.edu/~hjd1/SENECA/.

## III Hadronic cascade equations

In the domain of applicability of the cascade equations, the shower is one-dimensional and relativistic. Therefore it is completely specified by $`h_n`$, where $`h_n(E,X)dE`$ is the number of particles of a given species $`n`$ with energy in the range $`[E,E+dE]`$, at an atmospheric slant depth $`X`$ (with $`X=\rho _{\mathrm{air}}(x)𝑑x`$ ). The reaction probability of a particle in the atmosphere is given by its interaction length and decay length, so

$$\frac{dh_n(E,X)}{dX}=\frac{h_n(E,X)}{\lambda _n(E)}\frac{1}{c\tau _n\gamma \rho _{\mathrm{Air}}}h_n(E,X),$$
(1)
where $`\lambda _n(E)`$ is the mean free path, $`\tau _n`$ the lifetime of the particle and $`\gamma `$ its Lorentz-factor. By writing $`\rho _{\mathrm{Air}}=\rho _0\mathrm{exp}(h/h_0)=X/h_0`$, one can rewrite these equations as

$$\frac{dh_n(E,X)}{dX}=\frac{h_n(E,X)}{\lambda _n(E)}\frac{B_n}{EX}h_n(E,X)$$
(2)
where $`B_n`$ is the decay constant of hadron $`n`$ defined by

$$B_n=mc^2h_0/c\tau _n.$$
(3)
Accounting for particles produced at higher energies gives rise to the following system of hadronic cascade equationsBossard:2000jh

$`{\displaystyle \frac{h_n(E,X)}{X}}`$ $`=`$ $`h_n(E,X)\left[{\displaystyle \frac{1}{\lambda _n(E)}}+{\displaystyle \frac{B_n}{EX}}\right]`$
$`+{\displaystyle \underset{m}{}}{\displaystyle _E^{E_{\mathrm{max}}^{\mathrm{had}}}}h_m(E^{},X)[{\displaystyle \frac{W_{mn}(E^{},E)}{\lambda _m(E^{})}}`$
$`+{\displaystyle \frac{B_mD_{mn}(E^{},E)}{E^{}X}}]dE^{}.`$
Most important are the functions $`W_{mn}(E^{},E)`$, which are the energy-spectra $`\frac{dN}{dE}`$ of secondary particles of type $`n`$ in a collision of hadron $`m`$ with air-molecules. $`D_{mn}(E^{},E)`$ are the corresponding decay-functions. Equation (III) is a typical transport equation with a source term. The first term with the minus-sign accounts for particles disappearing by collisions or decays, whereas the source term accounts for production of secondary particles by collisions or decays of particles at higher energies. The technique for the solution is explained in detail in reference Bossard:2000jh  and we discuss in the following only two major changes.

First, the discretized functions

$$W_{mn}^i(E^j)=_{E_i/\sqrt{c}}^{E_i\sqrt{c}}\frac{E}{E_i}W_{mn}(E^j,E^{})𝑑E^{}$$
(5)
are not calculated by fitting the functions $`W_{mn}`$ to the energy-spectrum given by the hadronic Monte Carlo model, but by direct counting of the number of particles falling in the energy-bin defined by the limits of the integral in equation 5. This gives the desired precision as long as the number of simulated events is high, and avoids introducing systematic errors due to the fitting procedure. The binning of the discrete energies is

$$E_i=1GeV\times 10^{\frac{i1}{n_\mathrm{d}^{\mathrm{had}}}},$$

meaning $`n_\mathrm{d}^{\mathrm{had}}`$ logarithmic bins per decade. A typical value is 10, but as we will see below, a higher value can be preferable for some applications.

Second, the equations are modified to account for an arbitrary atmospheric density, since a real atmospheric profile is somewhat more complicated than the simple exponential form. Since the cascade equations are solved in layers $`X_i,X_{i+1}=X_i+\mathrm{\Delta }X`$ ,$`\mathrm{}`$ with, typically, $`\mathrm{\Delta }X=2.5g/cm^2`$, one can approximate the density in each layer as

$$\rho _{\mathrm{Air}}=a_i+b_iX\mathrm{for}X_i<XX_{i+1}.$$
(6)
The parameters $`a_i`$ and $`b_i`$ can easily be calculated from any function of the density. Dropping the indices, function (2) becomes

$$\frac{dh_n(E,X)}{dX}=\frac{h_n(E,X)}{\lambda _n(E)}\frac{B_n/h_0}{E(a+bX)}h_n(E,X)$$
(7)
which has the solution

$$h_n(E,X)=C\mathrm{exp}(X/\lambda )\left(a+bX\right)^{\frac{B_n}{h_0Eb}},$$
(8)
where $`a=a_i`$ and $`b=b_i`$ when $`X`$ is in the range $`X_i<XX_{i+1}`$. Defining the corresponding cascade equations is then straightforward. This generalization allows one not only to implement different atmospheres but also to solve for horizontal showers due to neutrino interactions in the atmosphere.

The initial condition for the cascade equation for a particle of type $`m`$ and energy $`E_m`$ at depth $`X_m`$ is given by:

$$h_n(E,X=X_m)=\delta _{nm}\delta (EE_m).$$
(9)
The initial condition for the cascade equation is thus in general a superposition of many functions like (9), which accounts for the natural fluctuations.

To recapitulate, down to a certain fraction $`f^{\mathrm{had}}=E_{\mathrm{max}}^{\mathrm{had}}/E_0`$ of the primary energy, all particles are followed with Monte Carlo method, meaning that each collision is simulated explicitly by the underlying event generator. Particles falling below $`E_{\mathrm{max}}^{\mathrm{had}}`$ are filled into the initial condition $`h_n(E,X)`$. After all particles above $`E_{\mathrm{max}}^{\mathrm{had}}`$ are processed, one can then proceed by solving the cascade equations.

## IV Electromagnetic cascade modeling

The electromagnetic part of the air-shower in reference Bossard:2000jh  was calculated with the analytic NKG formula. This is advantageous for the speed of the computation but it has some disadvantages: one has no detailed knowledge of particle spectra and it introduces inaccuracies since the NKG formula is only an approximation. The Monte Carlo approach can be used – e.g., the EGS4 EGS4  package provides a detailed Monte Carlo model for electromagnetic showers in any medium – but it is very time consuming for higher energies. We therefore apply the same approach as for hadronic cascades, by defining a system of electromagnetic cascade equations, analogous to (III). Due to the fact that $`e^\pm `$ and photons do not decay, the equations simplify greatly by setting the decay constants B to zero. This basically means that the showering is independent of altitude if one considers path lengths in units of $`g/cm^2`$. This fact allows a further simplification: First one defines energy-bins by

$$E_i=1GeV\times 10^{\frac{i1}{n_{\mathrm{d},\mathrm{em}}}}$$

the limits of each bin being

$$E_i10^{\frac{0.5}{^{n_{\mathrm{d},\mathrm{em}}}}}<EE_i10^{\frac{+0.5}{n_{\mathrm{d},\mathrm{em}}}}.$$

One defines $`V_{ij}^{mn}`$ as the number of particles of type $`n`$ (1=photon, 2=electron/positron) in energy bin $`E_j`$ generated in a electromagnetic shower induced by particle $`m`$ of energy $`E_i`$ , traversing a layer of air with thickness $`\mathrm{\Delta }X`$ of the order of some $`g/cm^2`$. This means $`V_{ij}^{mn}=V_{ij}^{mn}(\mathrm{\Delta }X)`$ is a function of layer thickness $`\mathrm{\Delta }X`$. In our case we choose $`\mathrm{\Delta }X=2.5\text{g}/\text{cm}^2`$.

Let $`g_i^n(X)`$ be the number of particles of type $`n`$ and energy $`E_i`$ at a given depth $`X`$. Then,

$$g_i^n(X+\mathrm{\Delta }X)=\underset{m,j>i}{}g_j^m(X)V_{ji}^{mn}(\mathrm{\Delta }X).$$

The function $`V_{ji}^{mn}(\mathrm{\Delta }X)`$ can be calculated in reasonable time by the showering model EGS4 since $`\mathrm{\Delta }X`$ is quite small. Once calculated, $`V_{ji}^{mn}(\mathrm{\Delta }X)`$ is stored as a table, and can be used to calculate efficiently any electromagnetic shower.

## V Low energy source-functions

We wish to follow particles down to an energy $`E_{cut}`$, below which they produce a negligible signal in the detector. Air showers have a lateral expansion of secondary particles, so the approach of 1-D cascade equations is certainly wrong for calculating particles down to lowest energies. Therefore, we employ the CE only to a certain minimum energy, $`E_{\mathrm{min}}^{\mathrm{had}}`$ and $`E_{\mathrm{min}}^{\mathrm{em}}`$ for hadronic and electromagnetic showers respectively. These are parameters which are determined empirically as described in the next section. The cases of electromagnetic and hadronic showers are analogous, so we describe the hadronic case for definiteness here.

Particles with energies $`E<E_{\mathrm{min}}^{\mathrm{had}}`$, produced in collisions with $`EE_{\mathrm{min}}^{\mathrm{had}}`$ contribute to the source function $`h_n^{\mathrm{source}}(E,X)`$ which is the number of produced particles at depth $`X`$ and energy $`E`$. It obeys the equation

$`{\displaystyle \frac{h_n^{\mathrm{source}}(E,X)}{X}}`$ $`=`$ $`{\displaystyle \underset{m}{}}{\displaystyle _{E_{\mathrm{min}}^{\mathrm{had}}}^{E_{\mathrm{max}}^{\mathrm{had}}}}h_m(E^{},X)[{\displaystyle \frac{W_{mn}(E^{},E)}{\lambda _m(E^{})}}`$ (10)
$`+{\displaystyle \frac{B_mD_{mn}(E^{},E)}{E^{}X}}]dE^{}.`$
The first term from equation (III) is missing because the propagation of these will be described by a Monte Carlo method.

An example of a typical source function of pions for a vertical $`10^{19}`$eV proton induced shower is given in Fig. 2, choosing $`E_{\mathrm{min}}^{\mathrm{had}}=10^4GeV`$. The source function is used to generate particles, which are then traced in the Monte-Carlo part of the air shower simulation code. With unlimited computational speed, the number $`N_n=_{E_{\mathrm{cut}}^{\mathrm{had}}}^{E_{\mathrm{min}}^{\mathrm{had}}}\frac{h_n^{\mathrm{source}}(E,X)}{X}𝑑E𝑑X`$ of particles would be produced for each species $`n`$, however at high energies this is time consuming. Instead, only a certain fraction of the total number of particles is sampled, attributing to each particle a suitable weight larger than 1. A practical way to define the sampling procedure is to specify the total amount of hadronic (or em) energy which is carried by the particles followed in the low energy MC. Because computation time is roughly proportional to energy, this procedure allows one to achieve equally good statistics independent of the energy $`E_0`$ of the primary cosmic ray. To be precise, the procedure is the following. The total energy of hadrons in the low energy regime produced by reactions with energy greater than $`E_{\mathrm{min}}^{\mathrm{had}}`$ is $`E_{\mathrm{low},\mathrm{tot}}^{\mathrm{had}}=_n_{E_{\mathrm{cut}}^{\mathrm{had}}}^{E_{\mathrm{min}}^{\mathrm{had}}}E\frac{h_n^{\mathrm{source}}(E,X)}{X}𝑑E𝑑X`$, with the index $`n`$ summing over the particle types. If low energy particles distributed according to the source function would be generated until their energy totalled $`E_{\mathrm{low},\mathrm{tot}}^{\mathrm{had}}`$, the weight would be 1. Instead we produce particles until their total energy is $`E_{\mathrm{low}}^{\mathrm{had}}<E_{\mathrm{low},\mathrm{tot}}^{\mathrm{had}}`$, so the weight attributed to each particle is $`w=\frac{E_{\mathrm{low},\mathrm{tot}}^{\mathrm{had}}}{E_{\mathrm{low}}^{\mathrm{had}}}>1`$. With this method, a simple adjustment of $`E_{\mathrm{low}}^{\mathrm{had}}`$ controls the final weight of all particles, since the shower of each generated particle is followed in full detail with no further thinning. Thus this method also overcomes a weakness of the normal thinning method, where the weight itself can fluctuate a lot.

It can happen that particles with $`E<E_{\mathrm{min}}^{\mathrm{had}}`$ are produced in the high energy MC stage, in reactions with particles of an energy $`E>E_{\mathrm{max}}^{\mathrm{had}}`$. Following all of these with weight 1 down to $`E_{\mathrm{cut}}^{\mathrm{had}}`$ is time consuming and unnecessary for particles with angle less than about $`5^{}`$ with respect to the shower axis. These are stored in the low energy source function and re-appear in the computation at the stage when low energy particles are created from the source function as discussed above. Low energy particles with larger angles are treated directly in the Monte Carlo part.

Neutral pions have a very short decay length and are therefore treated separately. In the system of cascade equations (III), they appear only as secondary particles and formula (10) is evaluated for $`E<E_{\mathrm{max}}^{\mathrm{had}}`$. The resulting $`\pi ^0`$s can then be fed either into the electromagnetic cascade equations for $`E>E_{\mathrm{min}}^{\mathrm{em}}`$ or into the Monte Carlo part of the code for $`EE_{\mathrm{min}}^{\mathrm{em}}`$. This approach is valid for $`E_{\mathrm{max}}^{\mathrm{had}}<B_{\pi ^0}=3\times 10^{19}eV`$.

The generalization of the source function method to the electromagnetic case is straightforward.

## VI Tests and applications

In the ideal case all the methods described here are just of technical nature, which means that the final result of physical observables of air showers should not depend on whether they are calculated with the traditional Monte Carlo (MC) method or with the hybrid method proposed here, using cascade equations (CE). Therefore a first step is to check the new technique by comparing the results of the two approaches and the influence of the parameters $`E_{\mathrm{min}}^{\mathrm{had}}`$ and $`E_{\mathrm{min}}^{\mathrm{em}}`$ on longitudinal and lateral profiles. In a second step we show some comparisons to the CORSIKA model, which can be configured to use the same external models - QGSJET, GHEISHA and EGS4.

Before doing so we show a comparison of the computation time necessary to simulate $`10^{19}eV`$ proton induced vertical showers in table 1. The results of CE using $`E_{\mathrm{low}}^{\mathrm{had}}=10^6GeV`$ and $`E_{\mathrm{low}}^{\mathrm{em}}=10^5GeV`$ are compared with a pure traditional Monte-Carlo method using various thinning levels $`f_{\mathrm{thin}}`$ ($`f_{\mathrm{thin}}\mathrm{E}_0`$ is the energy below which only one secondary $`i`$ in each reaction is followed with probability $`p_i=E_i/_iE_i`$ having a weight $`w_i=1/p_i`$).

The corresponding lateral distribution functions and the relative fluctuations can be seen in Fig. 3. The relative fluctuations in each lateral bin are defined by

$$\frac{\sigma }{\rho }=\frac{\sqrt{w_i^2}}{w_i},$$
(11)
where $`w_i`$ is the weight of particle $`i`$. One sees in Fig. 3 that the quality of the LDF computed with CE is somewhere between the thinning levels $`10^6`$ and $`10^7`$ (approaching the latter for large distances), whereas the computation time is at least 20 times lower as seen in table 1. As the energy of the primary cosmic ray is increased, the CPU time of the CE stays approximately constant if one keeps the same values for $`E_{\mathrm{low}}^{\mathrm{had}}`$ and $`E_{\mathrm{low}}^{\mathrm{em}}`$. This is because most of the time is used for the low energy Monte Carlo part. In the pure MC method, higher values of the thinning parameter $`f_{\mathrm{thin}}`$ can often be used while maintaining the statistical quality, so the improvement factor using the CE grows only slowly at higher energies. For instance, at $`3\times 10^{20}eV`$, $`f_{\mathrm{thin}}=10^6`$ gives results comparable to the CE approach using the same values for $`E_{\mathrm{low}}^{\mathrm{had}}`$ and $`E_{\mathrm{low}}^{\mathrm{em}}`$ with a time difference of a factor of about 40.

### VI.1 Checks on an average shower basis

If one applies the cascade equations starting from the primary energy, an average shower is calculated. In this case the initial condition consists just of the primary cosmic ray. We can compare such a result to the average of many showers computed by the MC-method. Fig. 4 shows such a comparison for $`10^{19}eV`$ proton induced vertical showers. The lateral and longitudinal profiles agree nicely within a small error. The shower maxima $`X_{max}`$ agree within less than 1%. As for the shower size $`N_{max}`$ (number of particles at shower maximum), we achieve 3% accuracy if we use 10 bins per decade in the numerical solution, but this can be improved to 1% by using 30 bins instead.

The other relevant parameters of the CE are $`E_{\mathrm{min}}^{\mathrm{had}}=10^4GeV`$, $`E_{\mathrm{min}}^{\mathrm{em}}=10GeV`$. The performance of the CE depends on these parameters as well as on the binning chosen for the numerical solution. A fine binning takes a long time to compute, whereas a more rough binning might introduce a significant error. Similarly, minimizing computation time argues for a low energy threshold $`E_{\mathrm{min}}`$ for both cascades, but not too low in order to obtain accurate lateral distribution functions. In the following we are going to analyze the influence of these parameters on the performance of the CE.

#### VI.1.1 $`E_{\mathrm{min}}^{\mathrm{em}}`$

The lower threshold $`E_{\mathrm{min}}^{\mathrm{em}}`$ should be chosen in the region where the electromagnetic shower cannot be treated anymore as one-dimensional, and the lateral spread of electrons, positrons and photons becomes important. Fig. 5 shows the longitudinal and lateral profiles for different values of $`E_{\mathrm{min}}^{\mathrm{em}}`$. As of a threshold of 10 GeV, the profiles do not change significantly anymore. A very similar approach was used in reference Lagutin:1999xh . There, a more complicated set of cascade equations was solved which involved also angular deviations from the shower-axis, and secondary particles below 10 GeV were followed in a MC method. Here we see that it is sufficient to treat the problem above 10 GeV in a purely longitudinal way.

#### VI.1.2 $`E_{\mathrm{min}}^{\mathrm{had}}`$

Analogously to the lower threshold in electromagnetic CE, the proper choice of the lower threshold $`E_{\mathrm{min}}^{\mathrm{had}}`$ depends on where the one-dimensional assumption is not valid anymore for the hadronic part of the shower. In order to test this we show in Fig. 6 the longitudinal and lateral and distribution function of muons, which are direct decay products of pions and kaons. A value of $`E_{\mathrm{min}}^{\mathrm{had}}=10^4GeV`$ provides sufficient precision for both profiles.

### VI.2 Tests on a single shower basis

By evolving the high energy part of a shower with the MC-method, we are able to reproduce the natural fluctuations, which are primarily due to the varying depth of the first interactions. All particles which fall below a threshold $`fE_0`$ are used in the initial condition for the CE. In order to show that the CE are solved correctly for an arbitrary initial condition, we compare to the MC method by computing the high energy part in exactly the same way for both approaches. Technically this can be done by choosing the same seed for the pseudo random number generator in the computer program. Fig. 7 shows such a comparison for the longitudinal and lateral profile of electrons and positrons. The thinning level for the MC method is $`10^{}7`$. We see a slight sensitivity of the longitudinal profile to the number of bins used in the numerical solution of the cascade equations. The shower maxima are at 738, 742, 739, and 740 $`g/cm^2`$ for the MC method, and the 10, 30 and 50 bin solution of the CE, respectively. The lateral distribution function is very insensitive to the number of bins.

### VI.3 Statistical properties - fluctuations

Next, we compare the statistical properties of two sets of proton induced $`10^{19}eV`$ showers calculated with CE (500 showers) and MC-method (500 showers, $`10^5`$ thinning). In Fig. 8 one sees the distribution of the shower maxima for the MC and CE methods. The two distributions agree well. The threshold $`f`$, where CE takes the initial condition from the high energy MC part and computes the shower numerically, was chosen to be $`0.001`$. The influence of the parameter $`f`$ on the fluctuations is shown in Fig. 9, by plotting $`\sigma =\sqrt{X_{\mathrm{max}}^2X_{\mathrm{max}}^2}`$ against $`f`$. The smaller the value of $`f`$, the further the initial portions of the shower is followed exactly rather than with the CE. We see that even for $`f`$ approaching unity, the fluctuations seen at small $`f`$ are reproduced to a great extent. This shows that natural fluctuations arise for the most part from the depth of the first interaction of the cosmic ray in the atmosphere.

### VI.4 Comparisons with CORSIKA

Finally, we compare some of our results to CORSIKA simulations. CORSIKA CORSIKA  is a well tested simulation package, which can be configured to use the same external models employed here for our hybrid model: QGSJET for the high energy hadronic part, GHEISHA for low energy hadrons and EGS4 for the electromagnetic part. Fig. 10 shows in the upper panel a comparison of the shower maximum for proton and iron induced 30-degree inclined showers. The shower proton curves agree nicely; the iron curve is slightly higher in our case. This might be due to differences in the computation of nucleus-nucleus cross sections, which are calculated from nucleon-nucleon cross sections using the Glauber method. We use $`f_{\mathrm{max}}^{\mathrm{had}}=0.001<1/56`$ which avoids calculating the functions $`W_{mn}`$ for all 56 possible nuclei. The lower panel compares the average lateral distributions functions of $`5\times 10^{19}eV`$ proton showers. They agree nicely for electrons/positrons and photons; compared to CORSIKA we produce slightly less muons, certainly due to the fact that we neglect photo-nuclear reactions at this stage.

### VI.5 Summary

We introduced a hybrid approach to air shower simulations which uses cascade equations as proposed in reference Bossard:2000jh . This method consists of applying traditional Monte Carlo methods where natural fluctuations or the lateral spread of particles are important. Particles are passed to the cascade equations via the initial condition. The low energy source function obtained from the cascade equations provides the probability distribution of low energy particles, whose further propagation is followed by Monte Carlo. The hybrid approach takes advantage of the fast solutions of the cascade equations and provides detailed knowledge about each low energy particle, such as position, energy and arrival time.

Consistency checks have been made by comparing the hybrid CE approach to traditional Monte Carlo. The two methods agree nicely within a small error. The longitudinal profiles obtained with the CE approach are somewhat sensitive to the binning which enters in the numerical solution of the CE, if less than $`30`$ bins per decade are used. The lateral distribution functions are very stable against these technical parameters.

The hybrid technique introduced here is faster than a traditional MC by at least a factor of 20 at $`10^{19}eV`$.

###### Acknowledgements.

This work was supported by NASA grant NAG-9246 and NSF grants NSF-PHY-9996173, and NSF-PHY-0101738. The computations were made on NYU’s Mafalda: a Linux cluster financed in part by the Major Research Instrumentation grant NSF-PHY-0116590. The authors thank D.Heck for providing $`X_{\mathrm{max}}`$ data from CORSIKA simulations. HJD would like to thank N.N.Kalmykov for introducing him into the subject of air shower simulations and the permission to use his cascade equations code.


text:
# 1 The effective potential as a function of 𝑞 for the parameters 𝜆=0.1,𝑘=7,𝑚=1 (a), and 𝜆=123,𝑘=1,𝑚=1 (b).

Non-perturbative approach for time-dependent

quantum mechanical system

Hyeong-Chan Kim1 and Jae Hyung Yee2

Institute of Physics and Applied Physics, Yonsei University

Seoul 120-749, Korea

1hckim@phya.yonsei.ac.kr

2jhyee@phya.yonsei.ac.kr

Abstract

We present a variational method which uses a quartic exponential function as a trial wave-function to describe time-dependent quantum mechanical systems. We introduce a new physical variable $`y`$ which is appropriate to describe the shape of wave-packet, and calculate the effective action as a function of both the dispersion $`\sqrt{\widehat{q}^2}`$ and $`y`$. The effective potential successfully describes the transition of the system from the false vacuum to the true vacuum. The present method well describes the time evolution of the wave-function of the system for short period for the quantum roll problem and describes the long-time evolution up to 75% accuracy. These are shown in comparison with the direct numerical computations of wave-function. We briefly discuss the large $`N`$ behavior of the present approximation.

PACS number: 11.80.Fv, 73.43.Nq, 11.30.Qc, 11.15.Tk

Phase transition is one of the most important physical phenomena in nature and has wide range of applications to condensed matter physics, particle physics, and cosmology. Most of the studies on this subject have been done in the framework of quasi-static transition or using the Gaussian ansatz developed by Jackiw and Kerman . There have been many attempts  to go beyond the Gaussian approximations. It is our purpose in this paper to go beyond the Gaussian approximation in two respects. First, we need a fully non-perturbative way which links the initial Gaussian packet (GP, false vacuum) to the symmetry broken degenerate vacuum state (true vacuum). Second, we try to find the relevant physical parameters which describes the symmetry breaking effectively.

In this paper, we consider a quantum mechanical model for time-dependent dynamics described by the potential,

$`V(\widehat{q},t)={\displaystyle \frac{\lambda }{24}}[\widehat{q}^2k^2(t)]^2,`$ (1)
where $`k^2(t)`$ increases from a negative value to a positive number $`\kappa ^2`$ asymptotically. The initial GP centered at $`q=0`$ cannot remain as Gaussian during the time-evolution, but evolves to the packet centered around two minima of the potential as $`k^2(t)`$ approaches $`\kappa ^2`$. For $`\kappa ^2\mathrm{}`$, the new ground states are linear sum or difference of two un-correlated GPs centered at each minima. In this case, the two ground states are degenerated.

The dispersion $`\widehat{q}^2`$ of a wavepacket may describe the size of a GP or the distance between two packets of a double Gaussian packet (DGP). To discern the shapes (for example, GP or DGP) of wavepackets of the same dispersion we introduce a dimensionless quantity $`y`$, which we call “shape factor”, in addition to the dispersion ($`q^2`$):

$`q^2(t)\widehat{q}^2,y(t){\displaystyle \frac{\widehat{q}^4}{\widehat{q}^2^2}}.`$ (2)
A similar expectation value as $`y`$ was calculated in Ref.  in relation to the new inflationary scenario. To illustrate the role of $`y`$ variable, consider a wavefunction which is a sum of two GPs of the same size. If $`y=1`$, the density of each GP is a delta function or the two GPs are infinitely far away so that no correlation exist between them, which provides the lower bound of $`y(1)`$. If the two GPs completely overlap, it corresponds to $`y=3`$, a single Gaussian packet. In between these two states, $`1<y<3`$, the two GPs are mixed and interfere with each other. For $`y>3`$, there are no separable packets, and the wavefunctions are better localized than GP .

The effective action in the variational method  is given by

$`\mathrm{\Gamma }={\displaystyle 𝑑t\psi ,t|i_t\widehat{H}|\psi ,t},`$ (3)
where $`\widehat{H}=\frac{\widehat{p}^2}{2m}+V(\widehat{q},t)`$ and we use $`\mathrm{}=1`$. In this paper we use the trial wavefunction,

$`Q|\psi ,t=N^1\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{2\mu ^2}}+i\mathrm{\Pi }\right)Q^4+\left({\displaystyle \frac{x}{\mu }}+ip\right)Q^2\right],`$ (4)
which has both of the DGP ($`x\mathrm{}`$) and the GP ($`x\mathrm{}`$) limits, where we assume $`\mu 0`$. In the static case, the double Gaussian approximation was used in Ref. , where a sum of two Gaussian functions is used as a trial wave-function. However, it is difficult to generalize the double-Gaussian method to the case for time-dependent systems. The normalization factor $`N`$ can be determined by the following integral:

$`N^2={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑Q\mathrm{exp}\left({\displaystyle \frac{Q^4}{2\mu ^2}}+{\displaystyle \frac{2xQ^2}{\mu }}\right)=\sqrt{\mu }f(x),`$ (5)
where $`f(x)`$ is given by

$`f(x)=|x|^{1/2}e^{x^2}{\displaystyle \frac{\pi }{\sqrt{2}}}\left[I_{\frac{1}{4}}(x^2)+sgn(x)I_{\frac{1}{4}}(x^2)\right].`$ (6)
The dispersion and the “shape factor” for this wavefunction are

$`q^2(t)={\displaystyle \frac{\mu f^{}}{2f}},y(t)={\displaystyle \frac{1+2xf^{}/f}{f_{}^{}{}_{}{}^{2}/(2f^2)}}.`$ (7)
$`y(x)`$ is a non-increasing function of $`x`$ from 3 to 1, which makes the inverse function $`x(y)`$ be defined uniquely. We use $`y`$ as a basic variable instead of $`x`$, because its range is bounded below by $`y=1`$ for any kinds of wavepacket  and it has definite physical meaning. The expectation values of other polynomials of $`\widehat{q}^2`$ can be written in terms of these parameters.

With this trial wavefunction the effective action is given by

$`\mathrm{\Gamma }={\displaystyle 𝑑t\left\{\frac{yq^4\dot{\mathrm{\Pi }}}{2}q^2\dot{p}2y\left[y\frac{y3}{Y+1}\right]\frac{q^6\mathrm{\Pi }^2}{m}\frac{2q^2p^2}{m}+\frac{4yq^4\mathrm{\Pi }p}{m}U(q,y)\right\}},`$ (8)
where $`Y(y)=2xf^{}/f=yf_{}^{}{}_{}{}^{2}/(2f^2)1`$ and the effective potential is

$`U(q,y)={\displaystyle \frac{V_F(y)}{8mq^2}}+V,`$ (9)
with the free potential, $`V_F`$, given by

$`V_F(y)=1+{\displaystyle \frac{(3y)(Y+1)}{y}}.`$ (10)
This free potential, coming from the expectation value $`\widehat{p}^2`$, represents the effect of quantum mechanical uncertainty. The expectation value of symmetric potential, $`V(\widehat{q},t)=V_0(t)+\frac{1}{2}k(t)\widehat{q}^2+\frac{\lambda (t)}{4!}\widehat{q}^4+\frac{c(t)}{6!}\widehat{q}^6+\mathrm{}`$, with respect to $`|\psi ,t`$ is

$`V=V_0(t)+{\displaystyle \frac{k(t)}{2}}q^2+{\displaystyle \frac{\lambda (t)}{4!}}yq^4+{\displaystyle \frac{c(t)}{6!}}y\left[y{\displaystyle \frac{y3}{Y+1}}\right]q^6+\mathrm{}.`$ (11)
From the action (8), we notice that $`\mathrm{\Pi }`$ and $`p`$ are the momentum conjugates to $`yq^4/2`$ and $`q^2`$, respectively.

Let us solve $`\mathrm{\Pi }`$ and $`p`$ equations first:

$`{\displaystyle \frac{1}{8}}{\displaystyle \frac{d}{dt}}\mathrm{ln}(yq^4)`$ $`=`$ $`\left[1{\displaystyle \frac{y3}{y(1+Y)}}\right]{\displaystyle \frac{yq^2\mathrm{\Pi }}{m}}+{\displaystyle \frac{p}{m}},`$ (12)
$`{\displaystyle \frac{1}{4}}{\displaystyle \frac{d}{dt}}\mathrm{ln}q^2`$ $`=`$ $`{\displaystyle \frac{yq^2\mathrm{\Pi }}{m}}{\displaystyle \frac{p}{m}}.`$
Removing $`\mathrm{\Pi }`$ and $`p`$ by Eq. (12) is just the Legendré transformation. Introducing new variable $`\eta `$ by $`{\displaystyle \frac{d\eta }{dy}}D={\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{1+Y}{y(3y)}}}`$, we get a quite simple effective action in terms of $`\eta `$ and $`q`$,

$`S={\displaystyle 𝑑t\left\{\frac{mq^2}{2}\dot{\eta }^2+\frac{m}{2}\dot{q}^2U[q,y(\eta )]\right\}}.`$ (13)
The dynamical equations of motion for $`q`$ and $`\eta `$ are given by

$`m{\displaystyle \frac{d}{dt}}\left[q^2\dot{\eta }\right]+{\displaystyle \frac{1}{D}}\left[{\displaystyle \frac{V_F^{}(y)}{8mq^2}}+_yV\right]`$ $`=`$ $`0,`$ (14)
$`m\ddot{q}{\displaystyle \frac{V_F}{4mq^3}}+_qV`$ $`=`$ $`mq\dot{\eta }^2.`$
The free potential, $`\frac{V_F(y)}{8mq^2}`$, has an absolute minimum at $`(y=3,q=\mathrm{})`$ and is positive definite. An interesting point here is that $`y=3`$, the GP, actually corresponds to $`x=\mathrm{}`$. On the other hand, in the effective potential (10), $`y=3`$ is a regular point, which can be extended to larger values. This property of the effective potential implies that the trial wavefunction (4) is insufficient to give a full description of $`y`$-dependence and we need a more general trial wavefunction for complete quantum mechanical description which includes the range $`y>3`$. The generalization of the trial function to $`|Q|\overline{\psi },t|^2=\overline{N}^2(1+|z|Q^2)^1|Q|\psi ,t|^2`$ may allow a full description of $`y`$-dependence, but, is difficult to integrate. Instead, we use a patch for the region $`y>3`$,

$`Q|\psi ,t=N^1\mathrm{exp}\left\{\left({\displaystyle \frac{1}{4G}}+i\mathrm{\Pi }\right)Q^2+(x+ip)|Q|\right\},y>3,`$ (15)
where we assume $`x0`$ and $`G0`$. The wavefunction is singular at the origin $`Q=0`$ due to the $`|Q|`$ term in the exponent. Hence we regard the absolute value as a small $`a`$ limit of $`\sqrt{Q^2+a^2}`$. With this trial wavefunction, we have the same effective action as Eq. (13) with $`D={\displaystyle \frac{_y\alpha }{2\alpha \sqrt{\alpha 1}}},(y>3),`$ and

$`V_F=\left[1+2z{\displaystyle \frac{2}{f}}\sqrt{{\displaystyle \frac{z}{\pi }}}\right]\left(1+{\displaystyle \frac{2}{f}}\sqrt{{\displaystyle \frac{z}{\pi }}}\right),\text{for }y>3,`$ (16)
where $`z=2x^2G`$, $`f(z)=e^z(1\text{erf}\sqrt{z})`$, and $`\alpha ={\displaystyle \frac{1+2z2\sqrt{z}/(\sqrt{\pi }f)}{2z[1/(\sqrt{\pi z}f)1]^2}}`$. The “shape factor” $`y(z)`$ monotonically increases from 3 to 6 as a function of $`z`$,

$`y(z)={\displaystyle \frac{3+12z+4z^22(5+2z)\sqrt{z}/(\sqrt{\pi }f)}{[1+2z2\sqrt{z}/(\sqrt{\pi }f)]^2}}.`$ (17)
As an example, let us consider the time-evolution of an initial wavepacket given by Eq. (4) with $`(y=y_0,q=q_0)`$ in the harmonic potential $`s(t)\widehat{q}^2/2`$. The dynamics of $`y`$ can be evaluated exactly from the elliptic integral

$`{\displaystyle _{\eta _0}^\eta }{\displaystyle \frac{d\eta }{\sqrt{c^2V_F[y(\eta )]}}}=\pm {\displaystyle \frac{1}{2m}}{\displaystyle _0^t}𝑑t^{}{\displaystyle \frac{1}{q^2(t^{})}},`$ (18)
where $`c^2=4m^2[\dot{\eta }(0)q^2(0)]^2+V_F[y(0)]`$ is an integration constant. The equation of motion for $`q`$ becomes

$`\ddot{q}(t)=s(t)q(t)+{\displaystyle \frac{2c^2V_F(y)}{8m^2q^3(t)}}.`$ (19)
If the system is potential free $`[s(t)=0]`$, $`\eta (t)`$ asymptotically approaches to a fixed value $`\eta _f`$ and $`q(t)`$ asymptotically increase with a constant velocity determined by the energy conservation law. The allowed range, for constant $`s(t)=s`$, of $`q`$ and $`y`$ is also determined by the energy conservation law:

$`{\displaystyle \frac{V_F(y)}{8mq^2}}+{\displaystyle \frac{1}{2}}sq^2E_{tot}=U(q_0,y_0).`$ (20)
Let us now consider the effective potential for the classical potential (1). The effective potential (9) naturally determines the true ground state with the condition,

$`_qU(q,y)=0=_yU(q,y).`$ (21)
To see the behavior of the effective potential more clearly, we variationally determine $`y`$, and then write down the effective potential in $`q`$:

$`V_{eff}(q)U(q,y_v(q))={\displaystyle \frac{V_F(y_v(q))}{8mq^2}}+{\displaystyle \frac{\lambda }{24}}\left[q^4y_v(q)2k^2(t)q^2+k^4(t)\right],`$ (22)
where $`y_v`$ means that we determine $`y`$ by minimizing the effective potential with the condition, $`\frac{\lambda }{24}q^4=\frac{V_F^{}(y_v)}{8mq^2}`$. In the static case this $`y_v`$ value corresponds to the minimum position of the potential for a given $`q`$, which takes $`y_v(0)=3`$ and $`y_v(\mathrm{})=1`$. We present figures (Fig. 1) of the effective potential (22) for two typical sets of parameters.

In these figures, $`V_{eff}=U`$, $`V_{G,q}`$, and $`V_{G,q^2}`$ represent the effective potential (22), the Gaussian approximated effective potential for $`\widehat{q}`$, and the Gaussian approximated effective potential for $`\widehat{q}^2`$, respectively. Here, $`V_{G,q}`$ can be calculated from Eqs. (2.9) and (4.6) of Ref.  with slight notational change ($`q=\widehat{Q}`$) and $`V_{G,q^2}`$ from Eq. (2.9) of Ref.  with $`\widehat{Q}=0`$ and $`Gq^2`$. The effective potential, $`V_{eff}`$, is very close to $`V_{G,q^2}`$ for $`q<k/3`$, and it becomes close to $`V_{G,q}`$ for $`q>2k/3`$. This clearly shows that the initial GP is divided into a DGP, with each packet of the DGP moving as if it is a free GP for large $`k`$. The value of $`y_v[1+\left(\frac{3}{\lambda mq^6}\right)^{2/5},`$ for $`y_v1.]`$ is effectively 1 for the most range of $`q`$ if $`k`$ is sufficiently large \[$`(\lambda m)^{1/6}`$\], since the characteristic size of $`q`$ is $`O(k)`$.

Let us explicitly describe the dynamics of an initial GP with $`q=q_0(\kappa )`$ for the time-dependent potential (1). Because of the transition (as $`k(t)`$ increase) $`y`$ eventually goes to 1 for most part of the dynamics. The potential energy difference $`\mathrm{}V=V_{eff}(q_0,y=1)V_{eff}(q_0,y=3)>0`$ and the presence of kinetic energy in $`y`$ prevent $`q`$ from reaching $`q_0`$. The time-dependence of $`k(t)`$ decreases the total energy so that $`q`$ oscillates near the true vacuum. We present, in Fig. 2, a solution of the differential equation (14) and its exact numerical solution for the case of $`k(t)`$ linearly increasing to a finite value for about a half period of $`q`$.

In this example, we do not need the patching process by the wavefunction (15) since states with $`y>3`$ do not appear. To see the long time behavior of the system, we present one more figure(Fig. 3).

The main characteristic feature of the long time behavior is the oscillation of the amplitude of short time oscillation. The error of the oscillation period in Fig. 3 is about 25%. This error comes from the non-exactness of the variational wave-function to the exact time-evolution of the wave-function.

As another example, let us consider the quenched transition potential with $`k(t>0)=\kappa `$ and $`k(t<0)=0`$. The discussions above for the time-dependent $`k(t)`$ transition also applies to the present example. We present a numerical solution of the differential equation (14) and its exact time evolution in Fig. 4.

The state with large $`y(>3)`$ periodically appears. This means that we need the patching process (15) for the time evolution of this system. Comparing the two results in Fig. 4, one may notice the merits and the weakness of the present approach for the quenched transition. The present approach explains the periodic appearance of the large “shape factor” and well presents the period of its occurrence, but the details of the evolution is not exact. This discrepancy is related to the “patched” trial wavefunction (4) and (15) at $`y=3`$. We have chosen this artificial patching method because of its simplicity. A better approach may be to include the excited states of (4) without introducing the patching, (15). One of the excited states, $`Q|\psi _2,t=(y1)^{1/2}\left(\frac{Q^2}{q^2(t)}1\right)Q|\psi ,t`$, is an orthonormal wavefunction to $`Q|\psi ,t`$. Because of the symmetry of the potential (1), the odd function of $`q`$ cannot contributes to the evolution. One may try the variational method by using the following trial wavefunction:

$`Q|\overline{\psi },t=\overline{N}^1\left[Q|\psi ,t+zQ|\psi _2,t\right],`$ (23)
where $`z`$ is a complex valued function of time and the normalization factor is $`\overline{N}^2=1+|z|^2`$. This wavefunction naturally includes the regions with $`y>3`$ due to the contribution of the excited state. In this sense, the appearance of large “shape factor” ($`>3`$) is the signal for the contribution of excited states in the time evolution of the systems. Generally, the accuracy of the approximation (4) increases as the potential varies slowly. We applied the present method to the case of scalar $`\varphi ^4`$ field theory in Ref. , and it would be interesting to apply more realistic quantum mechanical systems in second order phase transition.

Another point we need to speculate is the large $`N`$ limit. It was shown that the large $`N`$ wave-function satisfies

$`i{\displaystyle \frac{\mathrm{\Phi }(z,\tau )}{\tau }}=\left[{\displaystyle \frac{1}{2N^2}}{\displaystyle \frac{^2}{z^2}}+u(z,N)\right]\mathrm{\Phi }(z,\tau ),`$ (24)
where $`\tau =Nt`$, $`r^2{\displaystyle \underset{k=1}{\overset{N}{}}}x_k^2=Nz^2`$, and $`u(z,N)=\frac{(N1)(N3)}{8N^2z^2}+\frac{g}{8}(z^2z_0^2)^2`$. The present approximation for symmetric state is given by

$`\mathrm{\Gamma }={\displaystyle 𝑑\tau \psi ,\tau |i_\tau \frac{1}{2N^2}\mathrm{\Pi }^2u(z,N)|\psi ,\tau }.`$ (25)
This is the same as Eq. (3) with the change of parameters $`t,Q,m,V(Q,t)\tau ,z,N^2,u(z,N)`$. With the use of the trial wave-function (4) in (25) the expectation value of $`\frac{1}{2N^2}\mathrm{\Pi }^2`$ is given by

$`2y\left[y{\displaystyle \frac{y3}{Y+1}}\right]{\displaystyle \frac{q^6\mathrm{\Pi }^2}{N^2}}+{\displaystyle \frac{2q^2p^2}{N^2}}{\displaystyle \frac{4yq^4\mathrm{\Pi }p}{N^2}}+{\displaystyle \frac{V_F(y)}{8N^2}},`$ (26)
where the first three terms in Eq. (26) are $`O(1)`$, and in the large $`N`$ limit the quantum mechanical effects on the potential $`\frac{V_F(y)}{8N^2}`$ effectively vanishes as $`O(1/N^2)`$. In the absence of this quantum mechanical term, the equations of motion in the large $`N`$ limits for the present quartic exponential approximation with $`y=1`$ is the same as that of the Gaussian approximation centered at $`z0`$. Since the Gaussian approximation was proven to be the same as large $`N`$ approximation , the present approximation is equivalent to the large $`N`$ approximation for $`N\mathrm{}`$.

Acknowledgments
This work was supported in part by Korea Research Foundation under Project number KRF-2001-005-D2003 (H.-C.K. and J.H.Y.).


text:
# On two-dimensional finite-gap potential Schrödinger and Dirac operators with singular spectral curves

## 1 Introduction

In the present paper we describe a wide class of two-dimensional potential Schrödinger and Dirac operators which are finite-gap on the zero energy level and whose spectral curves at this level are singular and, in particular, may have $`n`$-multiple points with $`n3`$.

Dirac operators with such spectral curves are important for the Weierstrass representation of tori in $`^3`$ . A study of finite-gap operators with singular spectral curves which usually were not especially considered because of the nongeneric situation is of a special interest for differential geometry where singular curves may serve as the spectral curves of smoothly immersed tori. In particular, the spectral curves of tori in $`^3`$ obtained by a rotation of circles lying in the plane $`y=0`$ around the $`x`$ axis are rational curves with double points.

In the present paper the problem of describing such operators is reduced to a problem which involves only nonsingular curves by using the normalization of spectral curves. That makes this description to be effective.

## 2 Schrödinger and Dirac operators which are finite-gap on the zero energy level

### 2.1 Two-dimensional operators which are finite-gap the zero energy level

The notion of a two-dimensional operator which is finite-gap at one energy level was introduced by Dubrovin, Krichever, and Novikov  for the Schrödinger operator.

First we recall the definition of a Floquet eigenfunction (or a Bloch function) of a differential operator $`L`$ with periodic coefficients. Let $`L`$ act on functions on $`^n`$ and let its coefficients be periodic with respect to a lattice $`\mathrm{\Lambda }`$ isomorphic to $`^n^n`$. A solution of the equation

$$L\psi =\lambda \psi ,\lambda ,$$

is called a Floquet function (or a Bloch function) with the eigenvalue $`\lambda `$, if for any vector $`\gamma \mathrm{\Lambda }`$ we have

$$\psi (x+\gamma )=e^{2\pi ik,\gamma }\psi (x),$$

where $`k,\gamma =_{i=1}^nk_i\gamma _i`$ is the standard scalar product. The components of the vector $`k=(k_1,\mathrm{},k_n)`$ are called the quasimomenta of $`\psi `$. We see that any Floquet function defines a homomorphism

$$\mu :\mathrm{\Lambda }^{}=\{0\},\mu (\gamma )=e^{2\pi ik,\gamma }.$$

By using the Keldysh theorem and assuming that the coefficients of operators are bounded it has been proved that for the Schrödinger operator $`\mathrm{\Delta }+u`$, the heat operator $`_t\mathrm{\Delta }`$ and the two-dimensional Dirac operator, the quasimomenta and the eigenvalues of Floquet functions satisfy analytic relations (i.e. the dispersion laws ) and admissible tuples $`(k_1,\mathrm{},k_n,\lambda )`$ form an analytic subset $`Q`$ in $`^{n+1}`$ . Here it is essential that these operators are hypoelliptic. The set $`Q`$ is invariant under translations by vectors from the dual lattice $`\mathrm{\Lambda }^{}=\{\gamma ^{}:\gamma ^{},\gamma \text{for all}\gamma \mathrm{\Lambda }\}`$. Therefore it is easier to consider the quotient space $`Q/\mathrm{\Lambda }^{}`$.

Let $`n=2`$, i.e. the operators are two-dimensional. Then the intersection of $`Q/\mathrm{\Lambda }^{}`$ with the plane $`\lambda =0`$ is a complex curve (a Riemann surface) $`\mathrm{\Gamma }^{}`$, on which Floquet functions are glued into a function $`\psi (x,P),P\mathrm{\Gamma }^{}`$, which is meromorphic on the surface outside finitely many points. It is said that the operator $`L`$ is finite-gap at the zero energy level $`\lambda =0`$ if the curve $`\mathrm{\Gamma }`$, which is the normalization of the curve $`\mathrm{\Gamma }^{}`$, is a curve of finite genus, i.e. if $`\mathrm{\Gamma }`$ is an algebraic curve. This Riemann surface $`\mathrm{\Gamma }`$ is called the spectral curve of the operator $`L`$ at the zero energy level.

### 2.2 The Schrödinger operator

The two-dimensional Schrödinger operator with a magnetic field has the form

$$L=\overline{}+A(z,\overline{z})\overline{}+u(z,\overline{z}),$$
(1)
where

$$=\frac{}{x}i\frac{}{y},\overline{}=\frac{}{x}+i\frac{}{y},z=x+iy.$$

Following , it is said that such an operator is finite-gap at the zero energy level 111This definition as its analogue for Dirac operators (see §2.3) is given only for operators with nonsingular spectral curves. Operators with singular spectral curves appear from them in the limit under a degeneration which could be rather complicated as we show in this paper. if there exists

a) a nonsingular Riemann surface of finite genus $`g`$ with two marked points $`\mathrm{}_\pm `$ and local parameters $`k_\pm ^1`$ near these points such that $`k_\pm ^1(\mathrm{}_\pm )=0`$;

b) an effective divisor (a formal sum of points on the surface) $`D=P_1+\mathrm{}+P_g`$ of degree $`g`$ and formed by points different from $`\mathrm{}_\pm `$

such that there is a function $`\psi =\psi (x,y,P)`$ on $`\mathrm{\Gamma }`$ meeting the following conditions

1) this function is meromorphic with respect to $`P`$ on $`\mathrm{\Gamma }\{\mathrm{}_\pm \}`$, has poles only at points from $`D`$ and the order of a pole is not greater than the number of appearances of this point in $`D`$: $`(\psi )D`$;

2) $`\psi `$ has the following asymptotics at $`\mathrm{}_\pm `$:

$$\psi (x,y,P)e^{k_+z}(1+\xi (x,y)k_+^1+O(k_+^2))\text{as}P\mathrm{}_+,$$

$$\psi (x,y,P)c(x,y)e^{k_{}\overline{z}}(1+O(k_{}^1))\text{as}P\mathrm{}_{};$$

3) $`\psi `$ satisfies the equation $`L\psi =0`$ at every point $`P\mathrm{\Gamma }\{\mathrm{}_\pm \}`$.

It follows from the theory of Baker–Akhieser functions that

1) for a generic divisor $`D`$ the data $`(\mathrm{\Gamma },\mathrm{}_\pm ,k_\pm ,D)`$ determines a unique function $`\psi `$ meeting conditions 1 and 2;

2) given a function $`\psi `$, a unique Schrödinger operator of the form (1) can be constructed such that $`L\psi =0`$. The explicit formulas take the form

$$A=\frac{\mathrm{log}c}{z},u=\frac{\xi }{\overline{z}}.$$

As it is shown in , if there exists a holomorphic involution of $`\mathrm{\Gamma }`$:

$$\sigma :\mathrm{\Gamma }\mathrm{\Gamma },\sigma ^2=1,$$

such that $`\sigma (\mathrm{}_\pm =\mathrm{}_\pm ,\sigma (k_\pm )=k_\pm `$ and there exists a meromorphic differential (i.e. a $`1`$-form) $`\omega `$ on $`\mathrm{\Gamma }`$ with poles of the first order in the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ and zeroes in the points from $`D+\sigma (D)`$:

$$D+\sigma (D)\mathrm{}_+\mathrm{}_{}K(\mathrm{\Gamma })$$

(the divisor in the left-hand side is equivalent to the canonical divisor of the surface $`\mathrm{\Gamma }`$), then this operator is potential: $`c^21`$ and, therefore, $`A=0`$.

If in addition there exists an antiholomorphic involution

$$\tau :\mathrm{\Gamma }\mathrm{\Gamma },\tau ^2=1,$$

such that

$$\sigma \tau =\tau \sigma ,\tau (D)=D,\tau (\mathrm{}_\pm )=\mathrm{}_{},\tau (k_\pm )=\overline{k}_{},$$

then the potential $`u`$ is real-valued.

It is easy to notice that the form $`\omega `$ is invariant under $`\sigma `$ and hence it descends to a form $`\omega ^{}`$ on the quotient surface $`\mathrm{\Gamma }/\sigma `$ and the form $`\omega ^{}`$ has $`g`$ zeroes and two simple poles. Therefore the genus of $`\mathrm{\Gamma }/\sigma `$ equals $`g/2`$ and the points $`\mathrm{}_\pm `$ are exactly all fixed points of the involution $`\sigma `$ (notice that we assume that the surface $`\mathrm{\Gamma }`$ is nonsingular). The natural covering $`\mathrm{\Gamma }\mathrm{\Gamma }_0=\mathrm{\Gamma }/\sigma `$ is two-sheeted and branched at the points $`\mathrm{}_\pm `$.

Generically the potential $`u`$ is quasi-periodic and if it is periodic then the function $`\psi (x,y,P)`$ is a Floquet function for every $`P\mathrm{\Gamma }\{\mathrm{}_\pm \}`$ and the quasimomenta are locally holomorphic functions on the surface $`\mathrm{\Gamma }`$.

### 2.3 The Dirac operator

The Dirac operator (with potentials) has the form

$$𝒟=\left(\begin{array}{cc}0& \\ \overline{}& 0\end{array}\right)+\left(\begin{array}{cc}U& 0\\ 0& V\end{array}\right).$$

It is said that it is finite-gap on the zero energy level if there exists

a) a nonsingular Riemann surface $`\mathrm{\Gamma }`$ of finite genus $`g`$ with two marked points $`\mathrm{}_\pm `$ and local parameters $`k_\pm ^1`$ such that $`k_\pm ^1(\mathrm{}_\pm )=0`$;

b) an effective divisor $`D=P_1+\mathrm{}+P_{g+1}`$ of degree $`g+1`$ formed by points which differ from $`\mathrm{}_\pm `$,

such that there is a vector function $`\psi =(\psi _1,\psi _2)^{}=\psi (x,y,P)`$ meeting the following conditions:

1) the function $`\psi `$ is meromorphic in $`P`$ on $`\mathrm{\Gamma }\{\mathrm{}_\pm \}`$ and $`(\psi )D`$;

2) there are the following asymptotics:

$$\psi (x,y,P)e^{k_+z}\left[\left(\begin{array}{c}1\\ 0\end{array}\right)+\left(\begin{array}{c}\xi _1^+\\ \xi _2^+\end{array}\right)k_+^1+O(k_+^2)\right]\text{as}P\mathrm{}_+,$$

$$\psi (x,y,P)e^{k_{}\overline{z}}\left[\left(\begin{array}{c}0\\ 1\end{array}\right)+\left(\begin{array}{c}\xi _1^{}\\ \xi _2^{}\end{array}\right)k_{}^1+O(k_{}^2)\right]\text{as}P\mathrm{}_{};$$

3) the equation $`𝒟\psi =0`$ holds on $`\mathrm{\Gamma }\{\mathrm{}_\pm \}`$.

As in the case of the Schrödinger operator for a generic divisor $`D`$ the data $`(\mathrm{\Gamma },\mathrm{}_\pm ,k_\pm ,D)`$ determines a function satisfying the conditions 1 and 2 uniquely and from this function one can construct a unique operator $`𝒟`$ such that $`𝒟\psi =0`$:

$$U=\xi _2^+,V=\xi _1^{}.$$

Again as in the case of the Schrödinger operator generically these potentials are quasi-periodic but when they are periodic the functions $`\psi (x,y,P)`$ are Floquet functions whose quasimomenta locally holomorphically depend on $`P`$.

If there exists a holomorphic involution $`\sigma :\mathrm{\Gamma }\mathrm{\Gamma }`$ such that

$$\sigma (\mathrm{}_\pm )=\mathrm{}_\pm ,\sigma (k_\pm )=k_\pm ,$$

and there exists a meromorphic differential $`\omega `$ with zeroes in $`D+\sigma (D)`$ and two poles in the marked points $`\mathrm{}_\pm `$ with the principal parts $`\pm k_\pm ^2(1+O(k_\pm ^1))dk_\pm ^1`$, then the potentials $`U`$ and $`V`$ coincide: $`U=V`$ ().

If there exists an antiholomorphic involution $`\tau :\mathrm{\Gamma }\mathrm{\Gamma }`$, such that

$$\tau (\mathrm{}_\pm )=\mathrm{}_{},\tau (k_\pm )=\overline{k}_{}$$

and there exists a meromorphic differential $`\omega ^{}`$ with zeroes in $`D+\tau (D)`$ and two poles in $`\mathrm{}_\pm `$ with the principal parts $`k_\pm ^2(1+O(k_\pm ^1))dk_\pm ^1`$, then the potentials $`U`$ and $`V`$ are real-valued: $`U=\overline{U},V=\overline{V}`$ ().

For periodic operators these involutions are easily described on the language of quasimomenta:

$$\sigma (k_1,k_2)=(k_1,k_2),\tau (k_1,k_2)=(\overline{k}_1,\overline{k}_2),$$

and the existence of them immediately follows from from the spectral properties of the Dirac operator . 222Notice that it needs to add the condition that the potentials $`U`$ and $`V`$ are real-valued to the part 2 of Proposition 3 in , where this condition is used in the proof. Obviously in this case these involutions commute.

Remark. In papers  we did some inaccuracy assuming rather strong conditions for the differentials $`\omega `$ and $`\omega ^{}`$ by demanding that they have the following principal parts $`(\pm k_\pm ^2+O(k_\pm ^1))dk_\pm ^1`$ and $`(k_\pm ^2+O(k_\pm ^1))dk_\pm ^1`$ at points $`\mathrm{}_\pm `$. But the exposed proofs work under weaker conditions mentioned above. Indeed,

a) the differential $`\psi _1(P)\psi _2(\sigma (P))\omega `$ has two poles of the first order at the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ and the sum of residues equals $`2\pi i(\xi _2^++\xi _1^{})=0`$ which implies that $`U=V`$;

b) the differentials $`\psi _1(P)\overline{\psi _1(\tau (P))}\omega ^{}`$ and $`\psi _2(P)\overline{\psi _2(\tau (P))}\omega ^{}`$ have first order poles at the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ and the sums of residues are equal to $`2\pi i(\xi _1^{}\overline{\xi }_1^{})=0`$ and $`2\pi i(\xi _2^+\overline{\xi }_2^+)=0`$, respectively, which implies that $`U=\overline{U}`$ and $`V=\overline{V}`$.

Let us demonstrate the involutions $`\sigma `$ and $`\tau `$ by the following simple example.

Let the potential $`U=V=c`$ equals a real nonzero constant $`c`$. Then the spectral curve is the complex projective line $`\mathrm{\Gamma }=P^1`$ realized as the $`\lambda `$-plane completed by the infinity point $`\lambda =\mathrm{}`$. There two marked points $`\mathrm{}_\pm `$ on $`\mathrm{\Gamma }`$ such that $`\lambda =\mathrm{}`$ at $`\mathrm{}_+`$ and $`\lambda =0`$ at $`\mathrm{}_{}`$. Let us define near these points local parameters $`k_\pm ^1`$ by the formulas

$$k_+=\lambda ,k_{}=\frac{c^2}{\lambda }.$$

The function $`\psi `$ takes the form

$$\psi =\left(\begin{array}{c}\frac{\lambda }{\lambda c}\mathrm{exp}\left(\lambda z\frac{c^2}{\lambda }\overline{z}\right)\\ \frac{c}{c\lambda }\mathrm{exp}\left(\lambda z\frac{c^2}{\lambda }\overline{z}\right)\end{array}\right),$$

the divisor $`D`$ is just the point $`\lambda =c`$:

$$D=c,$$

and the involutions $`\sigma `$ and $`\tau `$ are defined by the formulas

$$\sigma (\lambda )=\lambda ,\tau (\lambda )=\frac{c^2}{\overline{\lambda }}.$$

The differentials $`\omega `$ and $`\omega ^{}`$ have the form

$$\omega =\left(1\frac{c^2}{\lambda ^2}\right)d\lambda ,\omega ^{}=\frac{(\lambda c)^2}{\lambda ^2}d\lambda .$$

## 3 Some facts on singular algebraic curves

We expose some necessary facts on singular algebraic curves following mostly to the book by Serre .

In the following we refer to a (complex) algebraic curve as a curve. Assuming that algebraic varieties are embedded into $`P^n`$ we say that a mapping between them is regular if it is defined by polynomials in homogeneous coordinates.

If a curve $`\mathrm{\Gamma }^{}`$ has singularities, then there is a normalization

$$\pi :\mathrm{\Gamma }\mathrm{\Gamma }^{},$$

where

1) $`\mathrm{\Gamma }`$ is a nonsingular curve with a finite set $`S`$ of marked points on it and given an equivalence relation $``$ between these points;

2) the mapping $`\pi `$ maps the set $`S`$ exactly onto the singular locus $`S^{}`$ of the curve $`\mathrm{\Gamma }^{}`$, and the preimage of every point from $`S^{}`$ consists of a class of all equivalent points;

3) the mapping $`\pi :\mathrm{\Gamma }S\mathrm{\Gamma }^{}S^{}`$ is a smooth one-to-one projection;

4) any regular mapping $`F:X\mathrm{\Gamma }^{}`$ of a nonsingular variety $`X`$ with an everywhere dense image $`F(X)\mathrm{\Gamma }^{}`$ descends through $`\mathrm{\Gamma }`$: $`F=\pi G`$ for some regular mapping $`G:X\mathrm{\Gamma }`$.

We recall that for any point $`P`$ from an algebraic variety there is a corresponding local ring $`𝒪_P`$ defined as ring of functions on the variety which are induced by rational functions $`f/g`$ where $`f`$ and $`g`$ are homogeneous polynomials of the same degree and $`g(P)0`$ (here we assume that the variety is embedded into $`P^n`$). A point is nonsingular exactly when its local ring is integrally closed.

For a point $`P\mathrm{\Gamma }^{}S^{}`$ its local ring $`𝒪_P^{}`$ is $`𝒪_{\pi ^1(P)}=𝒪_P`$. If $`PS^{}\mathrm{\Gamma }^{}`$, then $`𝒪_P^{}`$ is a subring of the ring

$$𝒪_P=\underset{QP}{}𝒪_Q$$

and moreover $`𝒪_P^{}`$ differs from $`𝒪_P`$ and for some integer $`n`$ we have the following inclusions

$$+R_P^n𝒪_P^{}+R_P𝒪_P,$$

where $`R_P`$ is an ideal, of the ring $`𝒪_P`$, consisting of all functions vanishing at $`\pi ^1(P)`$.

There is a particular case of constructing a singular curve $`\mathrm{\Gamma }_D`$ from a nonsingular curve $`\mathrm{\Gamma }`$ and an effective divisor $`D=n_PP`$ with degree $`\mathrm{deg}D=n_P2`$ on the curve $`\mathrm{\Gamma }`$. Let us denote by $`S`$ the set of points from $`\mathrm{\Gamma }`$ with $`n_P>0`$ (the support of the divisor) and put $`\mathrm{\Gamma }_D=(\mathrm{\Gamma }S)\{\text{pt}\}`$, i.e. contract all points from $`S`$ into one point which we denote by $`Q`$. Let us denote by $`C_Q`$ an ideal consisting of all functions $`f`$ which have at points $`PS`$ zeroes of order not less than $`n_P`$. Now put $`𝒪_Q^{}=+C_Q`$. This is the set consisting of all functions which have the same value at all points $`PS`$ and whose first $`(n_P1)`$ derivatives vanish at such a point. The natural projection $`\mathrm{\Gamma }\mathrm{\Gamma }_D`$ is the normalization. For $`D=P_1+\mathrm{}+P_n`$ where all points $`P_i`$ are pair-wise different the curve $`\mathrm{\Gamma }`$ has an $`n`$-multiple point $`Q`$ with different tangents.

To every singular point $`P\mathrm{\Gamma }^{}`$ there corresponds an integer-valued invariant

$$\delta _P=dim_{}𝒪_P/𝒪_P^{}<\mathrm{}.$$

It is obvious that for a singular point $`Q`$ of a curve $`\mathrm{\Gamma }_D`$ we have

$$\delta _Q=dim𝒪_Q/(+C_Q)=dim𝒪_Q/C_Q1=\mathrm{deg}D1.$$

The genus of the nonsingular curve $`\mathrm{\Gamma }`$ which is the normalized curve is called the geometric genus of the curve $`\mathrm{\Gamma }^{}`$ and it is denoted by $`p_g(\mathrm{\Gamma }^{})`$, and the quantity

$$p_a(\mathrm{\Gamma }^{})=p_g(\mathrm{\Gamma }^{})+\underset{PS}{}\delta _P$$

is called the arithmetic genus of the curve $`\mathrm{\Gamma }^{}`$.

We notice that a meromorphic $`1`$-form (a differential) $`\omega `$ on the curve $`\mathrm{\Gamma }`$ is called a differential regular at a point $`P\mathrm{\Gamma }^{}`$, if the equality

$$\underset{QP}{}\text{Res}(f\omega )=0$$

holds for all $`f𝒪_Q^{}`$. It is evident that there are more regular differentials on $`\mathrm{\Gamma }^{}`$ than regular differentials on $`\mathrm{\Gamma }`$, since regular differentials on $`\mathrm{\Gamma }^{}`$ may have have poles in the preimages of singular points. For instance, for a curve $`\mathrm{\Gamma }_D`$ forms which are regular at a singular point $`Q`$ are distinguished by the following conditions: a form $`\omega `$ may have poles only at $`PD`$ with their orders not greater than $`n_P`$ and

$$\underset{QP}{}\text{Res}\omega =0.$$

It is easy to see that the dimension of the space of regular differentials equals $`p_a(\mathrm{\Gamma }^{})`$.

Let the support of an effective divisor $`D`$ on the curve $`\mathrm{\Gamma }^{}`$ pairwise not intersect with the support of the divisor $`S^{}`$. We denote by $`dimL(D)`$ the space of meromorphic functions on $`\mathrm{\Gamma }^{}`$ with poles only at points from $`D=n_PP`$ with orders not greater than $`n_P`$, and denote by $`\mathrm{\Omega }^{}(D)`$ the space of regular differentials on $`\mathrm{\Gamma }^{}`$ which have at every point $`PS`$ a zero whose order is not less than $`n_P`$. The Riemann–Roch theorem reads that

$$dimL(D)dim\mathrm{\Omega }^{}(D)=\mathrm{deg}D+1p_a(\mathrm{\Gamma }^{}).$$

For a generic divisor $`D`$ we have $`dim\mathrm{\Omega }^{}(D)=0`$ and the Riemann–Roch theorem takes the form

$$dimL(D)=\mathrm{deg}D+1p_a(\mathrm{\Gamma }^{}).$$

## 4 Schrödinger and Dirac operators corresponding to singular spectral curves

We consider curves of the form $`\mathrm{\Gamma }_{B_1,\mathrm{},B_n}`$ which are successively constructed from effective divisors $`B_1,\mathrm{},B_n`$ and a curve $`\mathrm{\Gamma }`$ by the same procedure which constructs the curve $`\mathrm{\Gamma }_D`$ from $`\mathrm{\Gamma }`$ and $`D`$. Of course we assume that the supports of divisors $`B_i,i=1,\mathrm{},n`$, are pair-wise nonintersecting.

###### Theorem 1

1) Let $`\mathrm{\Gamma }^{}=\mathrm{\Gamma }_{B_1,\mathrm{},B_n}`$ be a singular curve, let $`\pi :\mathrm{\Gamma }\mathrm{\Gamma }^{}`$ be its normalization, let $`S_1,\mathrm{},S_n`$ be the supports of divisors $`B_1,\mathrm{},B_n`$ and let $`Q_1,\mathrm{},Q_n`$ singular points of the curve $`\mathrm{\Gamma }^{}`$: $`\pi (S_i)=B_i,i=1,\mathrm{},n`$.

Let $`\mathrm{}_+,\mathrm{}_{}`$ be a pair of different points from $`\mathrm{\Gamma }^{}\{Q_1,\mathrm{},Q_n\}`$ with local parameters $`k_\pm ^1`$ near these points such that $`k_\pm ^1(\mathrm{}_\pm )=0`$, and let $`D=P_1+\mathrm{}+P_g`$ be a generic effective divisor on $`\mathrm{\Gamma }\{Q_1,\mathrm{},Q_n,\mathrm{}_+,\mathrm{}_{}\}`$ of degree $`g=p_a(\mathrm{\Gamma }^{})`$.

Then there exists an unique function $`\psi (x,y,P),P\mathrm{\Gamma }^{}`$, such that it is meromorphic everywhere outside the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ where it has the following asymptotics

$$\psi (x,y,P)e^{k_+z}(1+\xi (x,y)k_+^1+O(k_+^2))\text{as}P\mathrm{}_+,$$

$$\psi (x,y,P)c(x,y)e^{k_{}\overline{z}}(1+O(k_{}^1))\text{as}P\mathrm{}_{};$$

and it has poles only in points from $`D=P_1+\mathrm{}+P_g=n_PP`$ of order not greater than $`n_P`$.

The function $`\psi `$ satisfies the equation $`L\psi =0`$, where the operator $`L=\overline{}+A(z,\overline{z})\overline{}+u(z,\overline{z})`$ is uniquely reconstructed from $`\psi `$:

$$A=\frac{\mathrm{log}c}{z},u=\frac{\xi }{\overline{z}}.$$

2) Let there be a holomorphic involution $`\sigma `$ on $`\mathrm{\Gamma }^{}`$ which preserves the marked points $`\mathrm{}_\pm `$ and inverts the local parameters near them: $`\sigma (\mathrm{}_\pm )=\mathrm{}_\pm `$, $`\sigma (k_\pm ^1)=k_\pm ^1`$. Let all singular points of the curve be fixed points of the involution and let the involution preserve the branches of the curve in these points (i.e. the pullback of the involution on $`\mathrm{\Gamma }`$ preserves all points from $`S`$).

If there is a differential $`\omega `$ on $`\mathrm{\Gamma }_{2B_1,\mathrm{},2B_n}`$ such that it is regular everywhere outside the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ in which it has first order poles with residues $`\pm 1`$ and it has zeroes exactly in the points of $`D+\sigma (D)`$, then $`L`$ is a potential operator:

$$L=\overline{}+u.$$

3) If in addition there is an antiholomorphic involution $`\tau :\mathrm{\Gamma }^{}\mathrm{\Gamma }^{}`$ such that $`\sigma \tau =\tau \sigma ,\tau (D)=D,\tau (\mathrm{}_\pm )=\mathrm{}_{},\tau (k_\pm )=\overline{k}_{}`$, then the potential $`u`$ is a real-valued function.

Proof. First we remark that for a nonsingular (smooth) curve the statement 1 is the theorem by Dubrovin, Krichever, and Novikov , and statements 2 and 3 were proved by Veselov and Novikov . We already exposed them in § 2.

1) For a smooth curve of genus $`g`$ a construction of a function $`\psi `$ with given properties was done in . Let us apply this construction to the normalized curve $`\mathrm{\Gamma }`$ and the divisors $`D_1=P_1+\mathrm{}+P_l,D_2=P_1+\mathrm{}+P_{l1}+P_{l+1},\mathrm{}`$, $`D_{gl+1}=P_1+\mathrm{}+P_{l1}+P_g`$, where $`l=p_g(\mathrm{\Gamma }^{})`$ is the genus of the curve $`\mathrm{\Gamma }`$. We obtain the functions $`\psi _1,\mathrm{},\psi _{gl+1}`$. The desired function $`\psi `$ has the form

$$\psi =c_1\psi _1+\mathrm{}+c_{gl+1}\psi _{gl+1},$$

where the coefficients $`c_i`$ are found by using two conditions:

1) the function $`\psi `$ descends to a function on $`\mathrm{\Gamma }^{}`$ (these are $`(gl)=(p_a(\mathrm{\Gamma }^{})p_g(\mathrm{\Gamma }^{}))`$ equations);

2) the asymptotics $`e^{k_+z}`$ holds at $`\mathrm{}_+`$.

For instance, given $`B=Q_1+\mathrm{}+Q_{gl+1}`$, the first condition is written as

$$\psi (Q_1)=\psi (Q_2),\psi (Q_1)=\psi (Q_3),\mathrm{},\psi (Q_1)=\psi (Q_{gl+1}),$$

and for $`B=mQ`$ it takes the form

$$\frac{\psi (Q)}{w}=\mathrm{}=\frac{^{m1}\psi (Q)}{^{m1}w}=0,$$

where $`w`$ is a local parameter on $`\mathrm{\Gamma }`$ near $`Q`$. The second condition looks the same in both cases:

$$c_1+\mathrm{}+c_{gl+1}=1.$$

The uniqueness of the functions $`\psi _i`$ follows from the Riemann–Roch theorem . Together with conditions for $`c_i`$ this implies the uniqueness of the function $`\psi `$ for a generic divisor $`D`$. For given functions $`A`$ and $`u`$ the function $`L\psi `$ is proportional to $`\psi `$, but its asymptotics as $`P\mathrm{}_+`$ are $`\alpha (z,\overline{z})e^{k_+z}k_+^1`$, and therefore $`L\psi `$ vanishes everywhere .

2) Let us consider the form $`\psi (P)\psi (\sigma (P))\omega (P)`$. It is meromorphic and the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ have residues $`1`$ and $`c^2`$, respectively. Let $`Q_i`$ be a singular point, of $`\mathrm{\Gamma }^{}`$, corresponding to the divisor $`B_i=n_QQ`$. Since the function $`\psi (P)\psi (\sigma (P))`$ is invariant under the involution, near a fixed point $`Q`$ this function expands in a series in even degrees of a local parameter $`k`$ where $`k(Q)=0,\sigma (k)=k`$:

$$\psi (P)\psi (\sigma (P))=\psi (Q)^2+a_1k^2+\mathrm{}+a_nk^{2n}+\mathrm{},k(P)=k,$$

and moreover, since $`\psi `$ descends to a function on $`\mathrm{\Gamma }^{}`$, we have $`a_j=0`$ for $`j<n_P`$. Since the form $`\omega `$ is regular on $`\mathrm{\Gamma }_{2B_1,\mathrm{},2B_n}`$, it has a pole of order not greater than $`2n_P`$ at the point $`P`$:

$$\omega =b_{2n_P}\frac{dk}{k^{2n_P}}+\mathrm{}+b_1\frac{dk}{k}+\text{the regular terms}.$$

Therefore, there is the following formula for the residue:

$$\text{Res}[\psi (P)\psi (\sigma (P))\omega ]|_{P=Q}=b_1(Q)\psi (Q_i)^2.$$

The regularity of $`\omega `$ at $`Q_i\mathrm{\Gamma }^{}`$ implies that the sum of the residues of $`\omega `$ over the preimage of this point vanishes:

$$\psi (Q_i)^2\underset{QS_i}{}b_1(Q)=0.$$

Therefore every singular point $`Q_i`$ does not contribute to the sum of the residues of the differential $`\psi (P)\psi (\sigma (P))\omega `$ and, since this sum equals

$$1c^2=0,$$

we have $`c^2=0`$ and $`A=\frac{\mathrm{log}c}{z}=0`$.

3) For nonsingular curves this statement was proved by Novikov and Veselov and in the nonsingular case the proof works without changes as follows. Let us consider the expansion for the function $`\psi `$ at $`\mathrm{}_{}`$:

$$\psi ce^{k_{}\overline{z}}(1+\eta k_{}^1+O(k_{}^2)).$$

Since $`L=\overline{}+u`$ and $`c^2=1`$, we obtain

$$L\psi =ce^{k_{}\overline{z}}((u+\eta )+O(k_{}^1))=0,$$

which implies the formula for the potential in terms of the asymptotics of $`\psi `$ near $`\mathrm{}_{}`$:

$$u=\eta .$$

Since the function $`\psi `$ is uniquely reconstructed from the data $`\mathrm{\Gamma }^{},\mathrm{}_\pm ,k_\pm ,D`$, this uniqueness theorem implies the equality

$$\overline{\psi (\tau (P))}=\overline{c}\psi (P)=c\psi (P).$$

In particular, $`\xi =\overline{\eta }`$, and comparing the two formulas for $`u`$ we have

$$u=\overline{}\xi =\eta ,$$

We conclude that the potential $`u`$ is real-valued: $`u=\overline{u}`$.

This proves the theorem.

The essential part of this theorem is the condition for the differential $`\omega `$ which has to be regular not on $`\mathrm{\Gamma }_{B_1,\mathrm{},B_n}`$ but on $`\mathrm{\Gamma }_{2B_1,\mathrm{},2B_n}`$. Together with that the proof of the following theorem on Dirac operators is obtained by a modification of the proof for a nonsingular case (see ). 333Although it is not difficult to formulate the following theorem using Theorem 1 and the results exposed in §2.3 we do that for the completness of exposition.

###### Theorem 2

1) Let $`\mathrm{\Gamma }^{}=\mathrm{\Gamma }_{B_1,\mathrm{},B_n}`$ be a singular curve, let $`\pi :\mathrm{\Gamma }\mathrm{\Gamma }^{}`$ be its normalization, let $`S_1,\mathrm{},S_n`$ be the supports of the divisors $`B_1,\mathrm{},B_n`$ and let $`Q_1,\mathrm{},Q_n`$ be singular points of $`\mathrm{\Gamma }^{}`$: $`\pi (S_i)=B_i,i=1,\mathrm{},n`$.

Let $`\mathrm{}_+,\mathrm{}_{}`$ be a pair of different points from $`\mathrm{\Gamma }^{}\{Q_1,\mathrm{},Q_n\}`$ with local parameters $`k_\pm ^1`$ near these points such that $`k_\pm ^1(\mathrm{}_\pm )=0`$, and let $`D=P_1+\mathrm{}+P_g+P_{g+1}`$ be a generic effective divisor on $`\mathrm{\Gamma }\{Q_1,\mathrm{},Q_n,\mathrm{}_+,\mathrm{}_{}\}`$ of degree $`g+1=p_a(\mathrm{\Gamma }^{})+1`$.

Then there exists a unique vector function $`\psi (x,y,P),P\mathrm{\Gamma }^{}`$, such that it is meromorphic everywhere outside the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ where it has the following asymptotics

$$\psi (x,y,P)e^{k_+z}\left[\left(\begin{array}{c}1\\ 0\end{array}\right)+\left(\begin{array}{c}\xi _1^+\\ \xi _2^+\end{array}\right)k_+^1+O(k_+^2)\right]\text{as}P\mathrm{}_+,$$

$$\psi (x,y,P)e^{k_{}\overline{z}}\left[\left(\begin{array}{c}0\\ 1\end{array}\right)+\left(\begin{array}{c}\xi _1^{}\\ \xi _2^{}\end{array}\right)k_{}^1+O(k_{}^2)\right]\text{as}P\mathrm{}_{};$$

and it has poles only in points from $`D=P_1+\mathrm{}+P_{g+1}=n_PP`$ of order not greater than $`n_P`$.

The function $`\psi `$ satisfies the equation $`𝒟\psi =0`$, where the operator

$$𝒟=\left(\begin{array}{cc}0& \\ \overline{}& 0\end{array}\right)+\left(\begin{array}{cc}U& 0\\ 0& V\end{array}\right)$$

is uniquely reconstructed from $`\psi `$:

$$U=\xi _2^+,V=\xi _1^{}.$$

2) Let there be a holomorphic involution $`\sigma `$ on $`\mathrm{\Gamma }^{}`$ which preserves the marked points $`\mathrm{}_\pm `$ and inverts the local parameters near them: $`\sigma (\mathrm{}_\pm )=\mathrm{}_\pm `$, $`\sigma (k_\pm ^1)=k_\pm ^1`$. Let all singular points of the curve be fixed points of the involution and the involution preserves the branches of the curve at these points (i.e. the pullback of the involution on $`\mathrm{\Gamma }`$ preserves all points from $`S`$).

If there is a differential $`\omega `$ on $`\mathrm{\Gamma }_{2B_1,\mathrm{},2B_n}`$ such that it is regular everywhere outside the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ in which it has second order poles with the principal parts $`\pm k_\pm ^2(1+O(k_\pm ^1))dk_\pm ^1`$ and its has zeroes exactly at points from $`D+\sigma (D)`$, then the potentials $`U`$ and $`V`$ coincide:

$$U=V.$$

3) Let there be an antiholomorphic involution $`\tau :\mathrm{\Gamma }^{}\mathrm{\Gamma }^{}`$ such that it interchanges the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$:

$$\tau (\mathrm{}_\pm )=\mathrm{}_{},\tau (k_\pm )=\overline{k}_{},$$

and its pullback onto $`\mathrm{\Gamma }`$ preserves all points from $`S_1\mathrm{}S_n`$, changing local parameters $`k`$ by the formula $`\tau (k)=\overline{k}`$.

Let there exist a differential $`\omega ^{}`$, on $`\mathrm{\Gamma }_{2B_1,\mathrm{},2B_n}`$, which is regular everywhere outside the points $`\mathrm{}_+`$ and $`\mathrm{}_{}`$ in which it has second order poles with the principal parts $`k_\pm ^2(1+O(k_\pm ^1))dk_\pm ^1`$ and which has zeroes exactly at points from $`D+\tau (D)`$. Then the potentials $`U`$ and $`V`$ are real-valued:

$$U=\overline{U},V=\overline{V}.$$

Remarks.

1) Potential Schrödinger operators whose spectral curves have only double points were described in the initial paper by Novikov and Veselov  as obtained by contracting invariant cycles on nonsingular curves into points (see pic. 1 where it is demonstrated by a deformation of an elliptic curve. Here by an involution we mean a rotation by $`\pi `$ around the horizontal line which lies in the plane of the picture and the right arrow denotes the normalization mapping).

These potentials are described in terms of the Prym theta functions of double-sheeted coverings of singular curves. The curve $`\mathrm{\Gamma }`$ has double points which are fixed points of the involution and moreover the branches are not permuted by the involution (i.e. on the normalized curve the preimages of such fixed points are fixed by the involution, see pic. 1). In this case a complete Prym variety is defined as a principally-polarized Abelian variety in a limit under the degeneration of nonsingular curves $`\mathrm{\Gamma }`$. In this limit the Jacobian variety of the curve $`\mathrm{\Gamma }/\sigma `$ is a non-complete Abelian variety, i.e. it has a form $`^g/Z`$ where the rank of a lattice $`Z`$ is less than $`2g`$.

If the branches in a double point are permuted by the involution (see pic. 2), then the limiting Prym variety is not complete but the limiting Jacobian variety of the curve $`\mathrm{\Gamma }/\sigma `$ is complete.

2) As we already pointed out above (in § 2), for a nonsingular curve corresponding to a finite-gap potential Schrödinger operator the following equality holds

$$p_a(\mathrm{\Gamma })=2p_a(\mathrm{\Gamma }/\sigma ),$$

which relates the arithmetic genera of the curve and its quotient under the involution $`\sigma `$. It also holds for curves with double points. Let $`\mathrm{\Gamma }_B`$ be a curve such that $`\mathrm{\Gamma }`$ is a hyperelliptic curve of genus two, $`\sigma `$ is the hyperelliptic involution, and $`\mathrm{}_\pm B`$ are six fixed points of $`\sigma `$: ($`B=Q_1+\mathrm{}+Q_4`$). We have $`p_a(\mathrm{\Gamma }_B)=p_g(\mathrm{\Gamma })+3=5`$ and $`\mathrm{\Gamma }_B/\sigma `$ is a rational curve (a sphere) with a quadruple point. Therefore, for this example we have

$$p_a(\mathrm{\Gamma }^{})=5,p_a(\mathrm{\Gamma }^{}/\sigma )=3.$$

In this case the potential of the Schrödinger operator is written in terms of the Prym variety of the double-sheeted covering $`\mathrm{\Gamma }^{}\mathrm{\Gamma }^{}/\sigma `$. This Prym variety is isomorphic to the Jacobian variety of the curve $`\mathrm{\Gamma }`$.

3) Notice, that, by the Krichever theorem , all smooth real potentials of the Schrödinger operator are approximated by the potentials which are finite-gap on the zero energy level to arbitrary precision. Moreover this approximation is generated by approximations of their Floquet spectra by nonsingular spectral curves of finite genera. For the Dirac operator the analogue of this theorem is not proved but it is clear that its proof can be obtained by some modification of Krichever’s reasonings.

4) One-dimensional Schrödinger operators $`L=_x^2+u`$ have hyperelliptic spectral curves which parameterize Bloch functions for all values of the energy . Degenerations of such curves (inside the class of hyperelliptic curves) lead to soliton potentials on the background of finite-gap potentials. .

For instance, on the language of § 3 the rational potential $`u(x)=2x^2`$ is constructed from the curve $`\mathrm{\Gamma }:w^2=E`$ and the point $`P`$ with $`E=w=0`$. Its spectral curve is $`\mathrm{\Gamma }_{2P}`$, and the spectral curve of the potential $`u(x)=2/\mathrm{cosh}^2(x)`$ has the form $`\mathrm{\Gamma }_{Q+\sigma (Q)}`$, where $`\sigma `$ is the hyperelliptic involution and $`Q\sigma (Q)`$.

5) Notice that the operators $`L=i_y_x^2+u`$ whose spectral curves on the zero energy level are singular and normalized by a rational curve were described in the paper .