text:
# The commensurate-disordered phase transition in 2D classical ATNNI model studied by DMRG

## Abstract

The classical two-dimensional anisotropic triangular nearest-neighbor Ising (ATNNI) model is studied by the density matrix renormalization group (DMRG) technique when periodic boundary conditions are imposed. Applying the finite-size scaling to the DMRG results a commensurate-disordered (C-D) phase transition line as well as temperature and magnetic critical exponents are calculated. We conclude that the (C-D) phase transition in the ATNNI model belongs to the same universality class as the ordered-disordered phase transition of the Ising model.

Analysis of semi-finite systems of small size in one or more directions has been used as a powerful tool in extracting of critical properties of two-dimensional classical models and corresponding one-dimensional quantum models. Although finite or 1D systems themselves do not display any critical behavior, it is, however, possible to extract critical parameter values as well as critical exponents. Temperature, ordering magnetic field, and finite-size deviations from criticality are all described by the same set of the critical exponents . This paper is focused on infinite strips of finite width where the relevant numerical data are obtained from the transfer matrix methods, in particular, the Density Matrix Renormalization Group (DMRG) method.

In 1992 the DMRG technique has been invented by S. R. White  in real space for one–dimensional (1D) quantum spin Hamiltonians. Three years later T. Nishino  applied this numerical technique to classical spin 2D models that is based on the renormalization group transformation for the transfer matrix for the open boundary conditions. DMRG treatment of 2D classical systems exceeds the classical Monte Carlo approach in accuracy, speed, and size of the systems .

Recently, we have modified the DMRG method for the 2D classical models, imposing periodic boundary conditions (PBC) on strip boundaries, and found a relation that helped to determine an optimal strip width $`L^{\mathrm{opt}}`$ in order to obtain correct values of critical temperature and exponents  using the finite-size scaling (FSS). We have obtained results of very high accuracy exceeding the DMRG method with standard open boundary conditions. Our method does not require any extrapolation analysis of the data.

The use of DMRG for 2D classical models may follow one of two different approaches:

(i) DMRG method is applied to strips of finite width and from two largest transfer-matrix eigenvalues or the free energy estimated with high precision, the critical properties of the system are calculated by the FSS analysis (here, we use this approach).

(ii) The strip width is enlarged until a steady state is reached (in the thermodynamic limit) when the output from the DMRG does not depend on the lattice size. Then, the DMRG yields properties of the 2D infinite system with spontaneously broken symmetry and mean-field-like behavior close to the criticality. This approach was used recently to study the high-field part of the ATNNI-model phase diagram , where approach (i) ran into convergence problems. We were able to show that the phase transition between the commensurate phase and the disordered phase proceeds via a narrow strip of an incommensurate phase. This approach gives also accurately the low-field part of the phase diagram, but it is not convenient for determination of the critical properties of the system by FSS. In distinction to the finite-width approach (i), it explicitly undergoes the phase transition, but its critical behavior is mean-field-like and the speed of calculation suffers from critical slowing-down at the phase transition line. Therefore, we use here approach (i) to find the low-field critical behavior of the ATNNI model.

The FSS approach should give the correct critical properties of the system in the limit of infinite strip width. Nevertheless, it was shown in  that in approximate DMRG treatment for given size of the transfer matrix (limited by computer capacity), it is not useful to enlarge the strip width to too large values, because here the the DMRG results do not satisfy the scaling laws assumed by the FSS. Thus, an optimal width, up to which the results systematically improve, must exist. It was also shown that the estimation of critical properties of the Ising and Potts models by DMRG with the periodic boundary conditions are much better than those with the open ones, despite the latter yields better results for the finite-width strips .

Below the optimal strip width $`L^{\mathrm{opt}}`$ the ratio

$$R\frac{\frac{}{L}T_\mathrm{C}^{}(L)}{\frac{^2}{L^2}T_\mathrm{C}^{}(L)}$$
(1)
is almost linear function of $`L`$ while above it, it is not. ($`L`$ in (1) is the width of the strip and $`T_\mathrm{C}^{}(L)`$ is the critical temperature for given $`L`$.)

The deviation of $`R`$ from linearity above the optimal strip width is very fast and the ratio $`R`$ becomes zero or infinity within enlargement of the strip by one lattice constant. Thus, if $`R=0`$ or $`R\mathrm{}`$ (i.e. if the numerator or the denominator tends to zero or changes its sign), we accept that $`L`$ as the strip width for further calculations and call it the optimal width $`L^{\mathrm{opt}}`$ of the strip. The critical temperature for the optimal width $`T_\mathrm{C}^{}(L^{\mathrm{opt}})`$ is taken as the best approximation of the critical temperature of the 2D system studied, and at this temperature the critical exponents of the system are calculated. The values of the critical exponents are sensitive to $`T_\mathrm{C}^{}`$ and must be determined with the due care.

In the FSS approach, the critical exponents are derived from the scaling behavior of the correlation length and free energy at critical point, where they depend on strip width $`L`$ in the following way :

$$K_L^\mathrm{h}L^{2y_\mathrm{h}^{(\beta )}}K_L^\mathrm{T}L^{y_\mathrm{T}^{(\nu )}}c_\mathrm{L}L^{2y_\mathrm{T}^{(\alpha )}d}$$
(2)
where $`K_L^\mathrm{T}`$ and $`K_L^\mathrm{h}`$ are the derivatives of inverse correlation length $`K`$ with respect to temperature $`T`$ and second derivative with respect to ordering (magnetic) field $`h`$, respectively, and $`c_\mathrm{L}`$ is the specific heat, i.e. the second derivative of the free energy with respect to temperature. The two temperature exponents $`y_\mathrm{T}^{(\alpha )}`$ and $`y_\mathrm{T}^{(\nu )}`$ should be equal to each other. The exponents $`y_\mathrm{T}`$ and $`y_\mathrm{h}`$ determine the critical behavior of all statistical quantities characterizing the system. The critical exponents of specific heat, magnetization and correlation length can be calculated from $`y_\mathrm{T}`$ and $`y_\mathrm{h}`$ as follows: $`\alpha =2\frac{2}{y_\mathrm{T}}`$, $`\beta =\frac{2y_\mathrm{h}}{y_\mathrm{T}}`$, $`\nu =y_\mathrm{T}^1`$. Other critical exponents can be obtained from the scaling equations $`y_\mathrm{h}=\beta +\gamma `$, $`\gamma =\beta (\delta 1)`$, and $`\eta =2\gamma y_\mathrm{T}`$ .

Further, we demonstrate the capabilities of our approach to find the critical properties of 2D spin lattice model on Ising model with different symmetries of the lattice, where critical temperatures and critical indices are known from exact solutions, and ATNNI model where the phase diagram is generally unknown and the critical indices are predicted from symmetry considerations.

The 2D classical anisotropic triangular nearest-neighbor Ising (ATNNI) model is given by the Hamiltonian

$$=\underset{i}{}J\left(\underset{\widehat{\delta }=\widehat{1},\widehat{2}}{}\sigma _i\sigma _{i+\widehat{\delta }}+a\sigma _i\sigma _{i+\widehat{3}}\right)H\underset{i}{}\sigma _i$$
(3)
with the antiferromagnetic coupling $`J<0`$ and spins $`\sigma _i=\pm 1`$. The numbers $`\widehat{1}`$, $`\widehat{2}`$, and $`\widehat{3}`$ are lattice directions in the ATNNI model. The coupling J is multiplied by the parameter $`a`$ ($`0<a<1`$) along the direction $`\widehat{3}`$ as depicted in Figure 1 (a).

This model was studied by Domany and Schaub  and in , and it was shown that its phase diagram, as a plot of temperature $`T`$ and external magnetic field $`H`$ (for $`a=0.4`$), exhibits four different phases: two commensurate phases $`I`$ and $`II`$, a disordered phase, and an incommensurate phase, see Figure 1 (b). Commensurate phase $`I`$ occurs at magnetic field $`H<2.4`$. This structure satisfies the Lifshitz condition, and it is characterized by a one-dimensional representation of the lattice symmetry group, i.e. its phase transition is predicted to belong to the Ising universality class . Domany and Schaub tried to confirm this prediction by numerical calculation of the exponent $`y_\mathrm{T}`$, but due to the low-order approximation it differed from the expected value by more than 10% and the magnetic exponent was not calculated at all.

We have calculated critical properties of the ATNNI model at the phase transition line between the commensurate $`I`$ and disordered phase. To illustrate the accuracy of the method, we have calculated critical properties of the exactly solvable models: ferromagnetic and antiferromagnetic Ising model on square and triangular lattices at zero magnetic field as well as the zero-magnetic-field ATNNI model which critical temperature is given by the equation

$$\mathrm{sinh}^2\left(\frac{2J}{T_\mathrm{C}}\right)=\mathrm{exp}\left(\frac{4aJ}{T_\mathrm{C}}\right).$$
(4)
We have used the FSS analysis of DMRG results with superblock consisting of 8 Ising spins and 4 multi-spin variable, each acquiring 85 values $`(m=85)`$. The computational effort at this approximation is less than for the classical transfer matrix method of strip width equal to 17 lattice constants. However, the DMRG enables to treat wider strip (of tens of lattice constants) up to the optimal width further improving the values of the critical parameters.

The first, important step of the calculations is determination of the critical temperature $`T_\mathrm{C}^{}`$, see Table 1, of which the best estimation for given $`m`$ is $`T_\mathrm{C}^{}(L^{\mathrm{opt}})`$ calculated from FSS approach . At this temperature the values of the critical exponents are derived from the scaling laws (2).

As the quantities appearing in (2) are first and second derivatives of the free energy and correlation length, the effect of approximation starts to manifest at lower strip width than $`L^{\mathrm{opt}}`$. The criterion determining strip width at which the value of the critical exponent may be still acceptable, was taken completely analogous to that for critical temperature, Eq. (1). The accepted values of the critical exponents are denoted by filled symbols in Figs 2 (a) and (b).

The critical exponent $`y_\mathrm{T}`$ is determined more precisely from the free energy $`y_\mathrm{T}^{(\alpha )}`$ than from the correlation length $`y_\mathrm{T}^{(\nu )}`$, as for the evaluation of the former only the largest eigenvalue of the superblock matrix is needed in distinction to the correlation length, to whose calculation the ratio of the largest and the second largest eigenvalue is necessary. This point is irrelevant for the models with a symmetric transfer matrix (Ising models in Table 1), but significant for the ATNNI model with a non-symmetric transfer matrix . The plot of thermal critical exponent $`y_\mathrm{T}^{(\alpha )}`$ vs. strip width is shown in Figure 2 (a). For increasing lattice size they both tend to the Ising value 1. The convergence also depends on the magnetic field. It gets worse for magnetic field close to the multi-critical point $`H=2.4`$. Here the reliability of the DMRG breaks down at rather small strip width, as well. The accepted values depicted by black symbols are listed in Table 2.

The critical exponent $`y_\mathrm{h}^{(\beta )}`$ describes the decay of the order parameter at the phase transition line from the commensurate phase $`I`$ to the disordered phase. The structure $`I`$ consists of two ferromagnetically ordered sublattices each with different magnetization. As the external magnetic field $`H`$ is generally non-zero in ATNNI model, the total magnetization (sum of both sublattice magnetizations) is non-zero, as well. The difference between the two magnetization is taken as the order parameter in this case. The small ordering field $`h`$ used for calculation of the derivative $`K_L^h`$ acquires opposite sign at each of the two sublattices. The accuracy of the calculations of the magnetic exponent is smaller than that of the thermal exponent in case of exactly solvable models listed in Table 2. Thus, we can expect a lower accuracy also for ATNNI model. All the exponents depicted in Fig. 2 (b) are below 1.871. Extrapolations to $`L\mathrm{}`$ for $`H=0.5`$–1.5 give values of $`y_\mathrm{h}^{(\beta )}`$ about 1.872, i.e. $`\beta =\frac{1}{7.81}`$, which still differs from the Ising value $`y_\mathrm{h}^{(\beta )}=1.875`$ and corresponding $`\beta =\frac{1}{8}`$. Note that the value of $`y_\mathrm{h}^{(\beta )}`$ is extremely sensitive to the correct determination of the critical temperature. A very small decrease of its value would shift $`y_\mathrm{h}^{(\beta )}`$ to the expected Ising value. At modest magnetic field, where our calculations are assumed to be more accurate, the plots of $`y_\mathrm{h}^{(\beta )}`$ for different magnetic field lie on the same curve what suggests that not only $`y_\mathrm{h}^{(\beta )}`$ is a universal quantity independent of $`H`$, but the corrections to it for finite $`L`$ are universal, as well.

In conclusion, it can be stated that the DMRG method with periodic boundary conditions reproduces with a high accuracy the critical properties of exactly solvable models and confirms the prediction that the C-D phase transition for magnetic fields $`H=0`$–2.4 belongs to the universality class of the Ising model.

A. G. thanks to T. Nishino for the fruitful correspondence and help. Most of this work has been done by Compaq Fortran on the HPC Alpha 21264 in Kobe University in Japan. This work has been supported by Slovak Grant Agency, VEGA n. 2/7174/20.


format: text

text:
# Near threshold Λ and Σ production in 𝑝𝑝 collisions

## Abstract

A model calculation for the reactions $`ppp\mathrm{\Lambda }K^+`$ and $`ppN\mathrm{\Sigma }K`$ near their thresholds is presented. It is argued that the experimentally observed strong suppression of $`\mathrm{\Sigma }^0`$ production compared to $`\mathrm{\Lambda }`$ production at the same excess energy could be due to a destructive interference between the $`\pi `$ and $`K`$ exchange contributions in the reaction $`ppp\mathrm{\Sigma }^0K^+`$. Predictions for $`ppp\mathrm{\Sigma }^+K^0`$ and $`ppn\mathrm{\Sigma }^+K^+`$ are given.

In a recent measurement of the reactions $`ppp\mathrm{\Lambda }K^+`$ and $`ppp\mathrm{\Sigma }^0K^+`$ near their thresholds it was found that the cross section for $`\mathrm{\Sigma }^0`$ production is about a factor of 30 smaller than the one for $`\mathrm{\Lambda }`$ production . We want to report on an exploratory investigation of the origin of this strong suppression of the near-threshold $`\mathrm{\Sigma }^0`$ production . In particular we want to examine a possible explanation that was suggested in Ref. , namely effects from the strong $`\mathrm{\Sigma }N`$ final state interaction (FSI) leading to a $`\mathrm{\Sigma }N\mathrm{\Lambda }N`$ conversion. We treat the associated strangeness production in the standard distorted wave Born approximation. We assume that the strangeness production process is governed by the $`\pi `$\- and $`K`$ exchange mechanisms. In order to have a solid basis for our study of possible conversion effects we employ a microscopic $`YN`$ interaction model developed by the Jülich group (specifically model A of Ref. ). This model is derived in the meson-exchange picture and takes into account the coupling between the $`\mathrm{\Lambda }N`$ and $`\mathrm{\Sigma }N`$ channels.

The vertex parameters (coupling constants, form factors) appearing at the $`\pi NN`$ and $`KNY`$ vertices in the production diagrams are taken over from the Jülich $`YN`$ interaction. The elementary amplitudes $`T_{KN}`$ and $`T_{\pi NKY}`$ are taken from corresponding microscopic models  that were developed by our group. However, for simplicity reasons we use the scattering length and on-shell threshold amplitudes, respectively, instead of the full (off-shell) $`KN`$ and $`\pi NKY`$ t-matrices. The off-shell extrapolation of the amplitudes is done by multiplying those quantities with the same form factor that is used at the vertex where the exchanged meson is emitted. Only s-waves are considered.

We do not take into account the initial state interaction (ISI) between the protons. Therefore we expect an overestimation of the cross sections by a factor of around 3 in our calculation . But since the thresholds for the $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^0`$ production are relatively close together and the energy dependence of the $`NN`$ interaction is relatively weak in this energy region the ISI effects should be very similar for the two strangeness production channels and therefore should roughly drop out when ratios of the cross sections are taken.

The cross section ratio for $`K`$ exchange alone and based on the Born diagram is 16, cf. Table 1. Including the $`YN`$ FSI, i.e. possible conversion effects $`\mathrm{\Sigma }N\mathrm{\Lambda }N`$, leads to a strong enhancement of the cross section in the $`\mathrm{\Lambda }`$ channel but only to a moderate enhancement in the $`\mathrm{\Sigma }^0`$ channel. As a consequence, the resulting cross section ratio becomes significantly larger than the value obtained from the Born term and, in fact, exceeds the experimental value. In case of pion exchange the Born diagram yields a cross section ratio of 0.9. Adding the FSI increases the cross section ratio somewhat, but it remains far below the experiment.

Thus, it’s clear that, in principle, $`K`$ exchange alone could explain the cross section ratio - especially after inclusion of FSI effects. However, we also see from Table 1 that $`\pi `$ exchange is possibly the dominant production mechanism for the $`\mathrm{\Sigma }^0`$ channel and therefore it cannot be neglected. Indeed, the two production mechanisms play quite different roles in the two reactions under consideration. $`K`$ exchange yields by far the dominant contribution for $`ppp\mathrm{\Lambda }K^+`$. The influence from $`\pi `$ exchange is very small. In case of the reaction $`ppp\mathrm{\Sigma }^0K^+`$, however, $`\pi `$\- and $`K`$ exchange give rise to contributions of comparable magnitude. This feature becomes very important when we now add the two contributions coherently and consider different choices for the relative sign between the $`\pi `$ and $`K`$ exchange amplitudes. In one case (indicated by “$`K+\pi `$” in Table 1) the $`\pi `$ and $`K`$ exchange contributions add up constructively for $`ppp\mathrm{\Sigma }^0K^+`$ and the resulting total cross section is significantly larger than the individual results. For the other choice (indicated by “$`K\pi `$”) we get a destructive interference between the amplitudes yielding a total cross section that is much smaller. Consequently, in the latter case the cross section ratio is much larger and, as a matter of facts, in rough agreement with the experiment (cf. Table 1) - suggesting a destructive interference between the $`\pi `$ and $`K`$ exchange contributions as a possible explanation for the observed suppression of near-threshold $`\mathrm{\Sigma }^0`$ production.

Table 1: Total cross section of the reactions $`ppp\mathrm{\Lambda }K^+`$ ($`\sigma _\mathrm{\Lambda }`$) and $`ppp\mathrm{\Sigma }^0K^+`$ ($`\sigma _{\mathrm{\Sigma }^0}`$) at the excess energies $`Q`$ = 13.2 MeV ($`\sigma _\mathrm{\Lambda }`$) and $`Q`$ = 13.0 MeV ($`\sigma _{\mathrm{\Sigma }^0}`$).

| diagrams | $`\sigma _\mathrm{\Lambda }[nb]`$ | $`\sigma _{\mathrm{\Sigma }^0}[nb]`$ | $`\frac{\sigma _\mathrm{\Lambda }}{\sigma _{\mathrm{\Sigma }^0}}`$ |
| --- | --- | --- | --- |
| K (Born) | 739 | 46 | 16 |
| K (FSI) | 2426 | 57 | 43 |
| $`\pi `$ (Born) | 71 | 77 | 0.9 |
| $`\pi `$ (FSI) | 113 | 105 | 1.1 |
| $`K+\pi `$” (FSI) | 2471 | 251 | 9.9 |
| $`K\pi `$” (FSI) | 2607 | 73 | 36 |
| exp.  | 505$`\pm `$ 33 | 20.1$`\pm `$3.0 | 25$`\pm `$6 |

It is now interesting to look also at corresponding results for other $`\mathrm{\Sigma }`$ production channels. E.g., for the reaction $`ppn\mathrm{\Sigma }^+K^+`$ the predicted cross sections at the excess energy of 13 MeV are 86 (“$`K+\pi `$”) and 229 $`nb`$ (“$`K\pi `$”), respectively. Thus, the interference pattern is just the opposite as for $`ppp\mathrm{\Sigma }^0K^+`$, cf. Table 1. For the “$`K\pi `$” case favoured by the experimental $`\sigma _\mathrm{\Lambda }/\sigma _{\mathrm{\Sigma }^0}`$ ratio our calculation yields a cross section for $`ppn\mathrm{\Sigma }^+K^+`$ that is about 3 times larger then the one for $`ppp\mathrm{\Sigma }^0K^+`$. Such a ratio is in fair agreement with data and model calculations at higher energies, see, e.g. Ref. . The other choice, “$`K+\pi `$”, leads to a $`\sigma _{ppn\mathrm{\Sigma }^+K^+}`$ that is a factor of about 3 smaller than $`\sigma _{ppp\mathrm{\Sigma }^0K^+}`$ \- a result which is rather difficult to reconcile with the present knowledge about these reactions at higher energies. These features are displayed graphically in Fig. 1 as a function of the excess energy. Note that all curves are multiplied by a common reduction factor of 0.3  to compensate for ISI effects .

For the reaction $`ppp\mathrm{\Sigma }^+K^0`$ the predicted cross sections at the excess energy of 13 MeV are 725 (“$`K+\pi `$”) and 423 $`nb`$ (“$`K\pi `$”), respectively. Thus, in this case the interference pattern is the same as for $`ppp\mathrm{\Sigma }^0K^+`$. Furthermore, for either choice (“$`K\pm \pi `$”) the cross section for $`ppp\mathrm{\Sigma }^+K^0`$ is about 3.3 times larger then the one for $`ppp\mathrm{\Sigma }^0K^+`$. Note that the experimental evidence at higher energies suggests a ratio of around 1 for those channels.


text:
# Possible composite-fermion liquid as a crossover from Wigner crystal to bubble phase in higher Landau level

\[

## Abstract

The ground state cohesive energy per electron of the composite fermion (CF) Fermi sea, the Laughlin state and the charge density wave (CDW) at higher Landau levels (LLs) are computed. It is shown that whereas for $`n2`$ LL, the CDW state is generally more energetically preferable than those of the CF liquid and the Laughlin liquid, the $`\nu =4+1/6`$ CF liquid state unexpectedly has lower ground state energy than that of the CDW state. We suggest this CF liquid between the Wigner crystal and the bubble phase may lead to the crossover from the normal integer quantum Hall liquid to the novel re-entrant integer quantum Hall state observed in the recent magneto-transport experiments.

\]

While fractional quantum Hall effects (FQHE) are impressive for their odd-denominator fillings, the remarkable phenomena for even-denominator fillings have caused great interests in the last decade . By using the concept of CF , Halperin et al. suggested that the $`\nu =1/2`$ system can be viewed as a spin-polarized Fermi liquid of CF . The only incompressible state at even-denominator filling $`(\nu =5/2)`$ in the single-layer 2-dimensional electron gas is now widely accepted as the pairing of CF in a vanishing effective magnetic field . For higher LLs, recent magneto-transport experiments on high mobility samples in GaAs/AlGaAs heterostructures revealed new classes of correlated many-electron states . The most prominent findings are the discoveries of the giant anisotropy in the resistivity near half filling of the topmost LL and the observation of re-entrant integral quantum Hall (RIQH) states in the flanks of these same levels  as well as the substantial non-linearity of the resistivity. It is considered that the highly anisotropic transport is related to the formation of the unidirectional charge density wave (UCDW) state, i.e., the stripe phase . Specifically, the state of this system may be classified by their symmetries, which are highly analogous to those of liquid crystals . The possible states include stripe crystals , smectic, and nematic phases . The RIQH effect was thought to be the depinning and sliding of the Wigner crystal (WC) and reformation of the bubble phase. But no further discussions were given on the details of this transition.
In the present work, we carry out systematic calculations of the cohesive energies of various competitive ground states in the lowest as well as higher LLs. We find that while in the lowest and second LLs the CF and the Laughlin liquids prevail for $`\nu _n<1/7`$ ($`\nu _n`$ is the filling factor at the $`n`$-th LL), the CDW generally has a lower ground energy than those of the corresponding CF or Laughlin liquid for LL indices $`n2`$. In particular, the stripe phase dominates around $`\nu _n=1/2(n2)`$, which is in agreement with previous calculations . However, we find an interesting exception occurs at $`\nu _2=1/6`$, where the CF liquid has a lower cohesive energy than that of the CDW phase. This result may lead to an important experimental phenomenon. Our calculation shows that the CDW has one electron in each bubble (the WC phase) for $`\nu _2<1/6`$ while two electrons for $`\nu _2>1/6`$ (the bubble phase). We relate this result to the recent observed phenomenon of the RIQH, i.e., the metallic CF liquid phase may appear as an intermediate state between the WC phase and the bubble phase. We suggest that the WC first melts into the CF liquid and then the system recrystalizes into the bubble phase as the filling at the topmost LL increases.
The fractional quantum Hall states at higher LLs were first suggested by MacDonald et al by raising the Laughlin wavefunction to higher LLs. The state wave function at the $`n`$-th LL is defined as follows:

$$|\mathrm{\Psi }_L^n=\underset{i}{}\frac{(a_i^{})^n}{\sqrt{n!}}|\mathrm{\Psi }_L^0$$
(1)
Here $`a_i^{}`$ is the inter-LL ladder operator, promoting the $`i`$-th elestron to the next LL, and $`|\mathrm{\Psi }_L^0`$ is the Laughlin state in the lowest LL. It should be noted that this state in the higher LL is in fact not the realistic electronic state because the lower LLs are empty. The effect of electrons in the lower filled LLs was treated by Aleiner and Glazman  and other authors  by integrating out the electron degrees of freedom in the lower LLs, which leads to a renormalization of the Coulomb interaction between electrons at the topmost LL.
The correlation energy per electron can be calculated using the density-density correlation function in the higher LL:

$$h_n(𝐫)\frac{\rho _n(𝐫)\rho _n(0)\rho _n^2}{\rho _n}$$
(2)
where $`\rho _n(𝐫)`$ is the projection of the density operator onto the $`n`$-th LL. This can be most effectively done in the Fourier space because $`h_n(q)`$ is very simply related to $`h_0(q)`$ (the correlation for $`n=0`$) :

$$h_n(q)=h_0(q)\left[L_n\left(\frac{q^2l_0^2}{2}\right)\right]^2$$
(3)
where $`L_n(\frac{q^2l_0^2}{2})`$ is the Laguerre polynomial. The correlation energy can be written as:

$$E_{cor}=\frac{1}{2}\frac{d^2q}{(2\pi )^2}\frac{v(q)}{\epsilon (q)}h_n(q)$$
(4)
The cohesive energy per electron is then, according to the definition

$$E_{coh}=E_{cor}E^{UEL}$$
(5)
where

$$E^{UEL}=\frac{\nu _n}{2}𝑑q\frac{e^2}{\kappa l_0}\frac{F_{nn}^2(q)}{\epsilon (q)}$$
(6)
is the interaction energy per particle in the uncorrelated electron liquid formed at high temperature. Here $`F_{nn}(q)=L_n(\frac{q^2l_0^2}{2})e^{q^2l_0^2/4}`$ and $`ϵ(q)`$ is the dielectric function which accounts for the screenings of the Coulomb interaction between electrons at the topmost LL by electrons in the lower LLs. Make use of the analytical expressions of the pair-distribution functions for $`\nu =1/3`$ and $`1/5`$ in the lowest LL obtained earlier by MacDonald et al. We reproduce the pair-distribution functions of the Laughlin states for $`\nu _n=1/3`$ and $`1/5`$ in the $`n=0,1,2`$ LLs (Fig.1). It shows that the electrons have a larger probability to approach each other in higher LLs than in lower LLs.
The trial wave function for the CF liquid ground state at higher LL can also be obtained in the same way described above  if the Rezayi-Read wave functions for the CF at the lowest LL are considered . We consider a system of $`N`$ electrons in the external magnetic field. Through attaching $`2p`$ ($`p`$ integer) flux quanta to each electron, one constructs the CF which experiences a vanishing effective magnetic field at even-denominator fillings. The interaction energy of electrons is transmuted into the kinetic energy of CFs. To compute the pair-distribution function of the CF liquid, we work on a spherical geometry in which the $`N`$ electrons move on the two-dimensional surface of a sphere under the influence of a radial magnetic field $`B`$ originating from a magnetic monopole of strength $`Q`$ at the center, which corresponds to a total flux of $`2Q\varphi _0`$, where $`\varphi _0=hc/e`$ is the flux quantum. We consider a filled shell system with $`N=n^2`$, because the ground state is a single Slater determinant, which makes computations easier. Employing the projection scheme proposed by Jain and Kamilla , we compute the pair-distribution function of CFs for $`N=25`$ electrons in the lowest LL by using the standard Metropolis Monte Carlo method. We carry out $`5\times 10^6`$ steps and the results are shown in Fig.2. Although the CFs hardly interact with each other in the Fermi sea state, the electrons carrying more flux avoid each other more steadly. The cohesive energies of CF liquids $`E_{coh}^{CF}`$ and Laughlin states $`E_{coh}^L`$ for $`\nu _n=1/2p`$ ($`p`$=1,2,3,4) are listed in Tables 1-4.
Due to the screenings by the lower LLs, the repulsive interaction between any two electrons at the topmost LL as a function of the separation between the guiding centers of their orbits abruptly drops at the distance of two cyclotron radii. Such a ”box-like” component in the interaction potential makes the uniform distribution of the electron density at the topmost LL unstable. The translational symmetry is spontaneously broken and crystalline domains with filling factor equal to one and zero are formed . The cohesive energy of the CDW state in the Hartree-Fock approximation is expressed as follows:

$$E_{coh}=\frac{n_L}{2\nu _n}\underset{𝐪0}{}u_{HF}(q)|\mathrm{\Delta }(𝐪)|^2$$
(7)
where $`n_L=1/2\pi l_0^2`$ is the density of one completely filled LL and $`\mathrm{\Delta }(q)`$ is the order parameter of the CDW.
For the bubble state, in the limit of a weak magnetic field a simple quasiclassical picture can be given. In this case electrons can be viewed as classical particles rotating in cyclotron orbits. It is shown that the optimum number of electrons in a bubble is $`\stackrel{~}{M}3n\nu _n`$, which corresponds to the separation $`3R_c`$ between nearest bubbles. However, to calculate the cohesive energy accurately, one can not use the quasiclassical approach and the CDW should be defined more precisely. According to Fogler and Koulakov,  it can be shown that for the bubble states,

$$\mathrm{\Delta }(q)2\nu _nA\frac{J_1(ql_0\sqrt{2M})}{ql_0\sqrt{2M}}$$
(8)
which is just the Fourier transform of a uniform disk with a radius $`l_0\sqrt{2M}`$. Here the asymptotic formula for the Laguerre polynomials for $`q\sqrt{M}/l_0`$ is used. The cohesive energy of the bubble state can be calculated in the same way as it has been done for the Wigner crystal. The results are different for each $`M`$. Therefore, one has to find the most optimum $`\stackrel{~}{M}`$ corresponding to the lowest energy. The final results for the stripe phase $`E_{coh}^{UCDW}`$ and the bubble phase $`E_{coh}^B`$ are also listed in Tables 1-4 for comparison.
The cohesive energies of the ground states for the Laughlin liquid, the CF liquid and the CDW are plotted in Fig.3. It can be seen that in the lowest LL the Laughlin liquids and the CF liquids have lower ground energies than those of corresponding CDW states for $`\nu <1/7`$. In the second LL, an incompressible pairing state of CF is considered to be more energetically preferable at $`\nu _1=1/2`$. The energy gap disappears when an in-plane magnetic field is applied to the system . It was ever interpreted as indicating the existence of substantial spin reversal in the ground state . The early tilt field experiments, however, missed an important point. In addition to suppressing the $`\nu =5/2`$ FQHE state, the tilted field leaves the transport in the $`n=1`$ LL highly anisotropic. Recently, the present authors demonstrated that the pairing of CFs is destroyed by the in-plane field and finally transforms to the UCDW . On the other hand, our calculations show that a CF liquid may exist at $`\nu _1=1/4`$ and $`1/6`$, in agreement with Morf and d’Ambrumenil .

For the $`n2`$ LL, the CDW states are generally favorable. The Laughlin liquid and the CF liquid are ruled out in the region of $`\nu _n<1/6`$, and specifically, the stripe phase dominates around $`\nu _n=1/2`$. This result agrees with earlier works using the Hatree-Fock method and recent numerical studies carried out by Rezayi, Haldane and Yang . However, for $`\nu _2=1/6`$, the CF liquid state is unexpectedly lower in energy than the CDW state. There is one electron in each bubble for $`\nu _2<1/6`$ (The WC phase) and two electrons for $`\nu _2>1/6`$ (the bubble phase). This means that a possible metallic phase may exist in the crossover from the WC to the bubble phase. The RIQH effect was previously explained as depinning and sliding of the WC and crystalizing of the bubbles . Unlike at the lowest LL, however, the lattice constant at higher LLs does not change as one increases the LL filling but remains of the order of $`R_cl_0`$. This implies the quantum fluctuations are hard to depin the WC in the higher LLs . In our picture, the WC melts into the CF liquid and subsequently the bubbles containing $`\stackrel{~}{M}=2`$ electrons form and are pinned as the filling increases. In fact, the feature of an isotropic metallic phase around $`\nu _2=1/6`$ can be identified in Fig.1 of Ref.. In addition, the $`\nu _2=1/5`$ state in our calculations is also slightly lower in energy than the corresponding CDW state while no FQHE was observed in this region. This may be explained as the ground state energies of compressible states (CDW and CF liquid) being lowered by disorder while the incompressible FQHE state hardly being affected. For $`n3`$, the RIQH effects are even weaker to be identified in the experiment . Our calculation also shows that it is very delicate. The CF cohesive energy is slightly higher than the CDW at $`\nu _3=1/6`$ (Fig. 3). However, the disorder should lower the CF ground state energy a little bit more than that of the CDW because of its gapless ground state.
In this work, we have presented a systematic computations of the cohesive energies of the CF liquid, the Laughlin liquid and the CDW in the lowest as well as in higher LLs. We find that the CF liquid and Laughlin liquid generally dominate at the lowest and second LLs while the CDW dominates at the higher LLs. The energy of the CF liquid at $`\nu =4+1/6`$ anomalously lower than that of the CDW may lead to the remarkable phenomenon of RIQH. Since the difference of energies in this region is very slight and the treatment of the screenings from the lower LLs is based on the Hartree-Fock approximation in the large $`n`$ limit, the results may not be completely decisive for $`n=2,3`$. Further numerical experiments are required to justify our conclusion.
One of the authors (SJY) thanks J. K. Jain for help. This work is supported in part by the NSF of China.

FIGURES

Fig.1 The pair-distribution functions of Laughlin states at $`n=0`$ (the solid line), $`n=1`$ (the dash-dotted line) and $`n=2`$ (the dotted line). (a) $`\nu _n=1/3`$. (b) $`\nu _n=1/5`$.
Fig.2 The pair-distribution functions of CF liquid states in the lowest LL at $`\nu =1/2`$ (the solid line), $`\nu =1/4`$ (the dash-dotted line) and $`\nu =1/6`$ (the dotted).
Fig.3 The cohesive energy per electron versus $`\nu _n`$ at the $`n`$-th LLs. The bubble state and the stripe state are represented by the solid line and dashed line, respectively. The triangle and square denote the Laughlin state and the CF liquid state, respectively. The dotted line is a guide to the eye. (a) $`n=0`$. (b) $`n=1`$. (c) $`n=2`$. (d) $`n=3`$.

TABLES

| $`\nu _n`$ | $`E_{coh}^B`$ | $`E_{coh}^{UCDW}`$ | $`E_{coh}^L/E_{coh}^{CF}`$ |
| --- | --- | --- | --- |
| 1/2 | -0.1320 | -0.1120 | -0.1551 |
| 1/3 | -0.1796 | -0.1469 | -0.2005 |
| 1/4 | -0.1945 | -0.1574 | -0.2040 |
| 1/5 | -0.1967 | -0.1584 | -0.2018 |
| 1/6 | -0.1941 | -0.1547 | -0.1971 |
| 1/8 | -0.1850 | -0.1466 | -0.1857 |

Table 1. The cohesive energy per electron for $`n=0`$ LL. The last column are energies for Laughlin states (odd-denominator fillings) or CF liquids (even-denominator fillings), respectively. The energy is in units of $`e^2/\kappa l_0`$.

| $`\nu _n`$ | $`E_{coh}^B`$ | $`E_{coh}^{UCDW}`$ | $`E_{coh}^L/E_{coh}^{CF}`$ |
| --- | --- | --- | --- |
| 1/2 | -0.0721 | -0.0687 | -0.0768 |
| 1/3 | -0.0996 | -0.0918 | -0.1062 |
| 1/4 | -0.1145 | -0.1014 | -0.1169 |
| 1/5 | -0.1229 | -0.1055 | -0.1256 |
| 1/6 | -0.1272 | -0.1069 | -0.1283 |
| 1/8 | -0.1295 | -0.1063 | -0.1288 |

Table 2. The same as Table 1 for $`n=1`$ LL.

| $`\nu _n`$ | $`E_{coh}^B`$ | $`E_{coh}^{UCDW}`$ | $`E_{coh}^L/E_{coh}^{CF}`$ |
| --- | --- | --- | --- |
| 1/2 | -0.0448 | -0.0456 | -0.0387 |
| 1/3 | -0.0635 | -0.0598 | -0.0603 |
| 1/4 | -0.0713 | -0.0659 | -0.0681 |
| 1/5 | -0.0762 | -0.0690 | -0.0764 |
| 1/6 | -0.0798 | -0.0707 | -0.0802 |
| 1/8 | -0.0834 | -0.0719 | -0.0829 |

Table 3. The same as Table 1 for $`n=2`$ LL.

| $`\nu _n`$ | $`E_{coh}^B`$ | $`E_{coh}^{UCDW}`$ | $`E_{coh}^L/E_{coh}^{CF}`$ |
| --- | --- | --- | --- |
| 1/2 | -0.0318 | -0.0329 | -0.0290 |
| 1/3 | -0.0441 | -0.0435 | -0.0420 |
| 1/4 | -0.0509 | -0.0480 | -0.0466 |
| 1/5 | -0.0531 | -0.0502 | -0.0490 |
| 1/6 | -0.0560 | -0.0514 | -0.0530 |
| 1/8 | -0.0593 | -0.0523 | -0.0587 |

Table 4. The same as Table 1 for $`n=3`$ LL.


text:
# Ring dark solitary waves: experiment versus theory

## I Introduction

Mathematically, optical dark spatial solitons (DSSs) are exact solutions of the one-dimensional nonlinear Schrödinger equation (NLSE) for negative nonlinearity and nonvanishing boundary conditions . Physically, they form on background beams of finite width as self-supported intensity dips due to the counterbalance between beam self-defocusing and diffraction. The required negative nonlinearity of the medium causes an inevitable reduction of the beam intensity. Losses, saturation and high transverse dimensionality result in nonintegrable model equations. Despite certain adiabatic relaxation characteristics, the solitary solutions of these equations have a large number of characteristics  in common with the soliton solution of the one-dimensional NLSE and are widely denoted by the term “dark soliton”. DSSs are generated as dark stripes , whereas the only known truly two-dimensional optical DSS is the optical vortex soliton (OVS) . Characteristic of the phase portraits of these beams are a 1D $`\pi `$-phase jump and an on-axis $`2\pi `$ helical phase ramp, also denoted as edge-phase and screw-phase dislocations, respectively.

Optical ring dark solitary waves (RDSWs) were first introduced by Kivshar and Yang . Their evolution was studied in the frame of the adiabatic approximation of perturbation theory. The quasi-one-dimensional treatment allowed the authors to obtain an expression for the RDSW’s transverse velocity $`dR/dz`$ as a function of the ring radius and the contrast. In addition, a linear stability analysis was performed. The authors showed that in the small-amplitude limit these ring dark waves are described by the cylindrical Korteweg- de Vries equation, which has ring soliton solutions . Recently, Frantzeskakis and Malomed derived new equations of evolution for small-amplitude solitary waves on a finite background  by means of a multiscale expansion method applied to the generalized cylindrical NLSE. Their long-wave solitary solutions propagate “on top” of the continuous-wave solution (i.e. on the background) as a dark perturbation in Kerr defocusing media. Depending on the background beam intensity, they can appear both dark and anti-dark for Kerr-like saturable defocusing nonlinearity.

The RDSWs may turn out to be of practical interest because of their ability to induce waveguides in which multiple signal beams could be guided parallel . This requires detailed analysis of the dynamics of these solitary structures. The first experimental generation of optical RDSWs  was conducted with pure amplitude modulation in front of the nonlinear medium (reflective micrometer-sized dots illuminated by the background beam). Under such conditions, pairs of phases oppositely shifted by about $`\pi /2`$ were measured in each diametrical slice of the RDSWs . In a subsequent experiment , RDSWs were generated from odd initial conditions by using binary computer-generated holograms (CGHs) . Inside a circle the phase was shifted by $`\pi `$ vs. the outlying area by shifting the photolithographically produced interference lines by half a grating period. In both experiments the phase profiles were measured by using the multiple-frame interferometric technique for optical wavefront reconstruction . It was found that along the NLM the transverse velocity of the RDSWs born from even initial conditions is higher than that of the odd ones. Numerical simulations have shown that, besides the ring radius , the transverse velocity of the RDSW can be effectively controlled by initial phase modulation of the background beam (inside, outside or on both sides of the phase dislocation) or by nonlinear interaction with a second coaxial dark formation (RDSW or OVS) . In saturable media, the RDSW’s transverse velocity is lower than that in Kerr NLM. At high saturation levels, in formal contradiction to but in actual agreement with the one-dimensional case, decay of the dark rings into a loop of optical vortices of alternate helicities should be expected .

We have measured experimentally the dynamics of ring dark solitary waves generated by different methods and compared those results with the predictions of the analytical method and numerical data. In particular, we analyzed the non-monotonic transverse velocity of RDSWs vs. background beam intensity in saturable NLM. The RDSWs are generated in the $`+1`$ and $`1`$ diffraction order of the background beams. Comparative numerical simulations that account for the estimated contrast of the RDSWs and medium saturation are carried out. An intentional reduction of the magnitude of the phase shift encoded in the CGHs was found to have a substantial influence on the dark ring radii and the contrast at the entrance of the NLM. The differences in the nonlinear evolution of the RDSWs enable us to verify qualitatively the analytical result of Kivshar and Yang (see Eqs. 6 and 7 in Ref. ).

## II RDSW transverse velocity vs. intensity

The experimental arrangement is shown in Fig. 1a. In order to obtain a dark ring of a desired radius and a phase jump of a certain magnitude, a single-line $`Ar^+`$-laser ($`\lambda =488nm`$) is used to reconstruct the respective CGH. The latter are produced photolithographically with a grating period of $`18\mu m`$ and ensure a diffraction efficiency of nearly $`10\%`$ in first order. The nonlinear medium is ethylene glycol dyed with DODCI (Lambdachrome). The background beam, diffracted in first order with the dark ring nested in, is transmitted through a slit placed about 15 cm behind the CGH and is gently focused on the entrance of the NLM. After passing the desired nonlinear propagation path length ($`0.5cm`$ to $`8.5cm`$) the dark beam is partially reflected by a prism immersed in the liquid and is projected directly on a co-moving CCD array with a resolution of $`13\mu m`$. Alternatively, the dark beam can be recorded at the exit of the NLM. For further evaluation the images are stored by an image- capturing card.

The data presented in this section refer to an absorption of $`\alpha =0.107cm^1`$ at $`\lambda =488nm`$. In a calibration measurement we generated a 1D DSS by using a CGH of the respective type. The soliton constant $`Ia^2`$ (i.e. the product of the background-beam intensity $`I`$ and the square of the dark beam width $`a`$ measured at the $`1/e^2`$-level) is found to reach its asymptotically constant value for input powers $`P_{\mathrm{sol}}^{1D}90mW`$ (Fig. 1b, dashed curve). Since the saturation of the (thermal) nonlinearity is expected to have important consequences on the RDSW’s nonlinear dynamics , we realized a background-beam self-bending scheme similar to that used in . The required asymmetry is introduced by tilting the prism immersed in the NLM. This results in different propagation path lengths for different parts of the beam. The strength of the self-bending effect (Fig. 1b, solid curve) measured in the near field gives an estimated saturation power $`P_{\mathrm{sat}}35mW`$. The choice of a saturation model for absorptive nonlocal nonlinearity is not trivial . Here we use a sigmoidal fit of the type $`f(I)I/(1+I/I_{sat})^3`$.

When a RDSW is generated by pure amplitude modulation of the background beam at the entrance of the NLM, its transverse velocity increases monotonically with the intensity (up to $`2.5P_{\mathrm{sol}}^{1D}`$; see Fig. 7 in Ref. ). By using the CGH of a black ring (with a phase jump of $`\pi `$), we observe a non-monotonic change of the ring radius vs. background beam power for a nonlinear propagation length of $`z=8.5cm`$ (Fig. 2a). The corresponding data for the contrast of the RDSW $`A^2=(I_0I_{min})/I_0`$ show reversed non-monotonic behaviour (Fig. 2b, solid curve). Generally, this is consistent with physical intuition and the theory  that the smaller the ring, the darker the solitary wave and the higher the total phase shift. As will be discussed later (see Fig. 4, the region on the left-hand side of the vertical dashed line), the dark rings initially generated with a $`\pi `$-phase shift (“black”) become wider by diffraction and thus arrive as “dark-grey” at the entrance of the NLM. In our case the contrast is lowered to $`A^2=0.88`$ (see Fig. 2b). We measured the change of the RDSW’s radius at $`z=8.5cm`$ as a function of the laser power. An evident minimum is observed for a power slightly below $`P_{sol}^{1D}`$. For this value the ring contrast is highest. Further insight can be gained by evaluating the analytical results for the phase of the solitary wave in the models of competing cubic-quintic and threshold nonlinearities (see Eqs. 6-9 in Ref. ). With the data for $`A^2`$, both models show well-pronounced minima in the RDSW transverse velocity for powers slightly lower than $`P_{\mathrm{sol}}^{1D}`$ (Fig. 2, dashed curve). It should be noted that the experimental data in Fig. 2 were obtained by compensating losses during the propagation with an appropriate change of the input power. However, even for fixed input power the RDSW changes its transverse velocity along the medium. In that sense Fig. 2a presents the effect of the enhancement of the phase shift and the contrast of the RDSW by saturating the nonlinearity in an integral form. This result, however, does not explain why in our previous experiments only a monotonic increase of the RDSW radii was observed .

In order to explain this discrepancy we look at the evolution of the RDSW before it enters the NLM. Numerical simulations are carried out for the evolution of “grey” solitary waves with initial radii $`R_0=R(z=0)`$ twice as large as the ring width $`r_0=r(z=0)`$. The numerical procedure used to solve the 2D NLSE is based on the beam propagation method over a $`1024\times 1024`$ grid. In agreement with our measurement, initially “black” rings generated by CGHs become “grey” after passing four Rayleigh diffraction lengths $`L_{Diff}=kr_0^2`$ and reach the estimated contrast $`A^2=0.88`$. The evolution of the RDSW’s radius for different intensities is calculated for up to $`10`$ nonlinear lengths ($`L_{NL}=(k|n_2|I_0)^1`$, $`k`$ is the wavenumber, and $`|n_2|I_0`$ the nonlinear refractive index correction). The behaviour of the RDSW radius in the NLM can be summarized as follows:
i) For a fixed input intensity in the range of $`(0.22.2)I_{\mathrm{sol}}^{1D}`$ all RDSWs have their minimal ring radius for propagation distances between $`1.5L_{NL}`$ and $`3L_{NL}`$.
ii) For nonlinear propagation distances shorter than $`2L_{NL}`$ the grey rings monotonically reduce their radii for increasing intensity (Fig. 3, lowest box).
iii) For long propagation distances ($`>6L_{NL}`$) this tendency is reversed (Fig. 3, upper box).
iv) In between, the radius of the RDSW as a function of the intensity shows a characteristic minimum around $`4.5L_{NL}`$ (Fig. 3, middle box). This propagation distance is in reasonable agreement with the NLM length of $`8.5cm`$ corresponding to $`4L_{NL}`$, for which Fig. 2 is obtained.

Generally, we found qualitative confirmation of the existence of a minimal radius $`R_{min}`$ of the RDSW, which depends on the initial contrast $`A^2`$ . The evolution of the contrast itself is strongly influenced by the saturation of the nonlinearity. In view of the above it is interesting to recall the observation which we made when modeling the decay of multiple-charged OVSs under radial perturbations : The small-amplitude even-phase initial conditions for perturbation and the saturation of the nonlinearity resulted in a negative transverse velocity of the rings. Their contrast was found to increase with decreasing radius and the radial phase variation reached a maximum of $`0.75\pi `$ at $`5L_{NL}`$. The initially decaying and still strongly overlapping optical vortices partially recovered their multiple-charged state (Fig. 7 in Ref. ). Periodic sign changes in the transverse velocity of the dark rings in the propagation direction were observed in the simulations up to $`10L_{NL}`$ with a period of approximately $`2L_{NL}`$.

## III RDSW transverse dynamics vs. initial phase shift

It should be pointed out that initially “black” RDSWs with high transverse dynamics are not likely to be present in a real experiment, at least when generated by holograms. The reason lies in the necessity to select one of the diffracted beams, which requires a certain free-space propagation. After some diffraction is “accumulated”, the dark rings with small radii spread out substantially and their contrast decreases. Dark rings with large radii are less affected and they evolve more slowly inside the NLM. Figure 4 is intended to visualize this before and inside the NLM. The coordinate $`z=0`$ is fixed on the nonlinear interface. Solid and dashed lines correspond to Kerr and saturable Kerr nonlinearity ($`s=I/I_{sat}=0.4`$). The simulation presented is carried out for $`I=I_{\mathrm{sol}}^{1D}`$ and for an initial phase jump of $`\pi `$, which will be flattened in the course of propagation. Generally, the saturation of the nonlinearity leads to a reduction of the dynamics of the RDSW. For our further considerations it is important to note that there is no difference whether the phase inside the ring $`\mathrm{\Phi }_{in}`$ is bigger or lower than the phase $`\mathrm{\Phi }_{out}`$ outside the ring, provided $`\mathrm{\Delta }\mathrm{\Phi }=|\mathrm{\Phi }_{in}\mathrm{\Phi }_{out}|=\pi `$. Equivalently, there should be no difference in the transverse dynamics of the RDSWs generated by the $`+1`$st and $`1`$st diffraction order beam when a phase jump of $`\pi `$ is encoded in the CGHs. This was carefully proved by numerical simulations and experimentally.

### A Linear evolution for $`|\mathrm{\Delta }\mathrm{\Phi }|<\pi `$

The situation changes when we produce dark ring waves encoded on CGHs with phase shifts smaller than $`\pi `$. Unfortunately, in binary CGHs only $`\pi /N`$ steps are possible, where one grating period consists of $`2N`$ elementary stripes. In our case $`N=3`$ and the grey rings are designed for $`|\mathrm{\Delta }\mathrm{\Phi }|=2\pi /3`$. The beam of a single-mode He-Ne laser is expanded 3.5 times and reproduces the respective CGH. By means of a lens ($`f=4cm`$) the beam transmitted directly and the $`\pm 1`$st order diffracted beams are imaged on a screen located about $`100cm`$ behind. By varying the lens position we are able to follow the dark ring evolution. Figure 5a shows the images of the still overlapping first-order diffracted beams, which interfere with the zeroth order one. It can be seen clearly that the rings with $`|\mathrm{\Delta }\mathrm{\Phi }|=2\pi /3`$ have smaller (larger) radius than each of the dark rings with an initial phase shift $`|\mathrm{\Delta }\mathrm{\Phi }|=\pi `$.

The ring radii observed for lens positions from $`4cm`$ to $`16cm`$ are shown in the inset of Fig. 5b. In order to get better spatial resolution, the measurement was repeated with a CCD array illuminated directly after the CGH (Fig. 5b). Blank circles denote the dark ring radii for $`|\mathrm{\Delta }\mathrm{\Phi }|=\pi `$, the solid circles refer to the ring with encoded phase jump $`|\mathrm{\Delta }\mathrm{\Phi }|=2\pi /3`$. Note that in the two graphs of Fig. 5b the ring radius $`R`$ is denoted in CCD-camera pixels, but these units are not directly comparable because of the different focusing conditions. On the CGHs the ring radius is encoded in 20 pix. (i.e. $`R_0=60\mu m`$). It is easy to understand the observed differences: The diffraction tends to broaden the dark rings and flatten the phase jumps. As a result, for $`\mathrm{\Phi }_{in}<\mathrm{\Phi }_{out}`$ an effective concave phasefront is formed which “focuses” the dark ring. In the other diffraction order $`\mathrm{\Phi }_{out}`$ appears to be bigger than $`\mathrm{\Phi }_{in}`$. The phase change induced by the diffraction results in an effective convex phasefront that increases the dark ring radius. For $`|\mathrm{\Delta }\mathrm{\Phi }|=\pi `$, the diffraction influences both first-order beams in the same way, but more weakly.

### B Nonlinear evolution of initially grey RDSWs

Although the nonlinear regime was analytically analyzed by Kivshar and Yang , to the best of our knowledge it has not been investigated experimentally. The basic results (Eqs. 5-7 in Ref ) can be summarized as follows:
i) For each RDSW there is a minimum ring radius $`R_{min}`$ for which the RDSW has its highest contrast.
ii) Depending on the initial value of the phase jump, the RDSW either collapses along the NLM to reach $`R_{min}`$ and diverge monotonically thereafter, or it diverges monotonically from the beginning.

Figure 6 is intended to clarify these conclusions. We solved the NLSE numerically by keeping the following parameters close to the experimental values: RDSW radii ($`(R/r)|_{z=0}=2,4,`$ and $`6`$), background-beam intensity ($`I=I_{\mathrm{sol}}^{1D}`$), initial contrast $`A^2=0.88`$, phase jump $`|\mathrm{\Delta }\mathrm{\Phi }|=0.837\pi `$, and saturation of the nonlinearity $`s=0.4`$. In Fig. 6a the corresponding results for the case of a higher phase inside the rings ($`\mathrm{\Phi }_{in}>\mathrm{\Phi }_{out}`$) are presented. In this case, the RDSWs start diverging immediately after entering the NLM. The data shown in Fig. 6b refer to the case of a lower phase inside the dark ring and higher phase outside ($`\mathrm{\Phi }_{in}<\mathrm{\Phi }_{out}`$). The effective concave wavefront forces the RDSWs to collapse to a minimum ring diameter. Thereafter they start to broaden faster. The main tendency, that the broader the RDSW initially is, the longer the nonlinear propagation needed to reach $`R_{min}`$ will be, is well expressed. The other general tendency, that the smaller the RDSW diameter is, the higher its transverse dynamics is, can also be clearly seen.

Both tendencies and the influence of the saturation seem to be understood theoretically, but still need experimental confirmation. We performed additional measurements with the experimental setup shown in Fig. 1, in which the telescope is replaced by a single lens ($`f=12cm`$) in order to focus the beams near the entrance of the NLM. Moderate and strong saturation $`s`$ can be achieved by raising the absorption $`\alpha `$ to $`0.2cm^1`$ ($`P_{\mathrm{sol}}^{1D}=25mW;P_{\mathrm{sat}}=40mW;s=0.6`$) and $`0.4cm^1`$ ($`P_{\mathrm{sol}}^{1D}=20mW;P_{\mathrm{sat}}=14(\pm 5)mW;s1.4`$). In order to keep the conditions for both first-order beams as similar as possible, the beam transmitted directly was blocked at the entrance of the NLM and both $`\pm 1`$st order beams were recorded simultaneously. Unfortunately, due to the interaction of the two background beams, the interaction-free nonlinear propagation path length is limited. The results of the measurements are shown in Fig. 7a for moderate and high saturation (host graph and inset, respectively). Blank circles denote “black” RDSWs radii, solid ones RDSW radii for $`\mathrm{\Delta }\mathrm{\Phi }=\pm 2\pi /3`$. The best results for moderate saturation are achieved with CGHs containing phase dislocations with radii of $`R=30\mu m`$. For high saturation, rings with dislocation radii of $`60\mu m`$ are used. As seen, after some $`5mm`$ of nonlinear propagation the “grey” RDSWs reach $`R_{min}`$ and diverge monotonically afterwards. The beam of the other diffraction order and the RDSW generated for $`|\mathrm{\Delta }\mathrm{\Phi }|=\pi `$ start diverging from the beginning. The interpolated curves in the figure are intended to guide the eye. The interpolated curves for $`z<2mm`$ and $`z>17mm`$ are not shown for an obvious physical reason$`:`$ As mentioned, the grey RDSWs do not enter the NLM with equal radii and equal (absolute) phase shifts. For larger propagation distances, the saturation-like behaviour is caused mainly by interaction with the second co-propagating beam. The reduced transverse velocities of the RDSWs for higher saturation are also strongly pronounced. The grey-scale images in Fig. 7b show the RDSWs at $`z=2mm`$ and $`8mm`$ inside the NLM.

By evaluating the data corresponding to the largest RDSWs in our experiment ($`R=90\mu m`$) we observed internal ring splitting for moderate saturation (Fig. 8). Despite their large radii, these RDSWs have non-zero transverse velocity. They broaden and flatten due to the saturation (see Fig. 8, frame at $`z=2mm`$) while evolving into grey rings. It seems that such a grey ring can serve as a background on which two coaxial rings can form, subsequently narrow and thus get a high modulation depth for increased propagation distance (Fig. 8, frames at $`z=5,8,11,14,`$ and $`17mm`$). The mechanism for their formation includes effects from increasing a phase shift and decreasing a transverse velocity due to the saturation of the nonlinearity , as well as repulsion between opposite quasi-2D phase dislocations. The relation of these coaxial rings to the small-amplitude dark solitary waves described by Frantzeskakis and Malomed  should be clarified.

## IV Conclusion

It is worth noting that it is difficult (if possible at all) to ensure perfect odd initial conditions (radial phase jumps equal to $`\pi `$ and contrast equal to unity) for the generation of RDSWs in cubic nonlinear media. When such a phase jump is encoded in a CGH, there is no difference between the evolution of the two beams in first diffraction order. This changes gradually for smaller phase jumps. In qualitative agreement with the theory  we presented evidence of the monotonic increase of the dark ring radius for $`\mathrm{\Phi }_{in}>\mathrm{\Phi }_{out}`$ and the corresponding non-monotonicity in the opposite case. For grey initial conditions and saturation of the nonlinearity we found characteristic regions for the dependence of the RDSW radius on the background-beam intensity. The radius can increase, decrease monotonically or go through a minimum. The results presented are in qualitative agreement with the analytical ones regarding the phase of the solitary wave in the presence of saturation .

###### Acknowledgements.

The authors are pleased to acknowledge highly useful discussions with Prof. Yu. Kivshar and his suggestions for improving the manuscript. A. D. would like to thank the Alexander von Humboldt Foundation for the award of a fellowship and the opportunity to work in the stimulating atmosphere of Max-Planck-Institut für Quantenoptik (Garching, Germany). This work was also supported by the National Science Foundation of Bulgaria and by the Science Fund of Sofia University.


text:
# References

Gauge choice and geodetic deflection in conformal gravity

A. Edery1, A. A. Méthot2 and M. B. Paranjape3

1Center for High Energy Physics, Department of Physics,
McGill University, 3600 University Street,
Montréal, Québec, Canada, H3A 2T8

2,3Groupe de Physique des Particules, Département de Physique,
Université de Montréal, C.P. 6128, succ. centre-ville,
Montréal, Québec, Canada, H3C 3J7

3National Center for Theoretical Sciences, Physics Division,
National Tsing Hua University
Hsinchu, Taiwan 300, ROC

11footnotetext: edery@hep.physics.mcgill.ca22footnotetext: methot@lps.umontreal.ca33footnotetext: paranj@lps.umontreal.ca

Abstract

Conformal gravity has been proposed as an alternative theory of gravity which can account for flat galactic rotation curves without recourse to copious quantities of dark matter. However it was shown that for the usual choice of the metric, the result is catastrophic for null or highly relativistic geodesics, the effect is exactly the opposite yielding an effective repulsion and less deflection in this case. It is the point of this paper, that any result for massive geodesics depends on the choice of conformal gauge, in contradistinction to the case of null geodesics. We show how it is possible to choose the gauge so that the theory is attractive for all geodesics.

One of the most important problems facing astrophysics is that of missing matter. All observations indicate that particles, massive or massless, are attracted more strongly by cosmological sources than would be expected on the basis of the observable (luminous) mass of the source, together with the use of Newton/Einstein gravity. There are two evident potential solutions to this problem. We can postulate the existence of (copious) quantities of non-luminous matter (dark matter), in just the right distribution, to account for the excess attraction, or we can reject Einstein gravity and propose an alternative theory of gravity chosen to better reflect the cosmological phenomenology. Either solution implies a radical departure from our present understanding of cosmology. Conformal gravity, in recent years, has been proposed as such an alternative theory of gravitation.

Kazanas and Mannheim have had reasonable success in fitting galactic rotation curves. They consider the spherically symmetric vacuum solution in conformal gravity

$$d\tau ^2=B(r)dt^2\frac{dr^2}{B(r)}r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2\right)$$
(1)
with

$$B(r)=1\frac{2\beta }{r}+\gamma rkr^2$$
(2)
where $`\beta ,\gamma `$ and $`k`$ are constants. $`kr^2`$ corresponds to a cosmological solution (conformal to a Robertson-Walker background) so that $`k`$ is chosen small enough not to have any effect at galactic scales while $`\beta `$ and $`\gamma `$ are chosen such that they do influence the gravitation exactly at these scales. Then for weak fields, at galactic scales, the effective Newtonian potential for non-relativistic matter is given by

$$\varphi (r)=\frac{2\beta }{r}+\gamma r.$$
(3)
Clearly for $`\gamma >0`$ this gives an additional attractive linear potential. Fixing $`\gamma `$ phenomenologically yields reasonable fits to galactic rotation curves. However, it was recently pointed out that a positive sign for $`\gamma `$ yields replusion for null and highly relativistic geodesics in blatant contradiction with observation. The deflection is given by $`\mathrm{\Delta }\phi =\frac{4\beta }{r_0}\gamma r_0`$ where $`r_0`$ is the radius of closest approach. This seems to arrest the development of the Kazanas-Mannheim program.

Our point in this paper is that it is possible to circumvent this difficulty by an appropriate choice of the conformal gauge. The key observation is that while massless geodesics are conformally invariant, massive geodesics are not. A massless geodesic satisfies

$$d\tau ^2=0=B(r)dt^2\frac{dr^2}{B(r)}r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2\right)$$
(4)
hence conformal rescaling $`d\tau ^2\mathrm{\Omega }d\tau ^2`$, evidently, has no effect on the geodesic equation. Massive geodesics are, however, sensitive to the conformal factor $`\mathrm{\Omega }`$.

Take the sign of $`\gamma `$ in (2) to be negative, $`\gamma <0`$, so that the conformal gravity implies additional attraction for the null geodesics ($`\mathrm{\Delta }\phi =\frac{4\beta }{r_0}+|\gamma |r_0`$). The metric

$$d\tau ^2=\mathrm{\Omega }(r)\left(B(r)dt^2\frac{dr^2}{B(r)}r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2\right)\right)$$
(5)
will continue to be a vacuum solution of conformal gravity with the same massless geodesics. We will show how to find $`\mathrm{\Omega }(r)`$ so that non-relativistic massive geodesics will also exhibit additional attraction. We will work to first order in the weak field approximation, however, it is clear that our idea is not restricted to this domain. Massless geodesics are insensitive to the conformal gauge whatever the gravitational field, while massive geodesics are always sensitive to the conformal factor, since the coupling of mass to the theory necessarily breaks the conformal invariance.

To identify the effective Newtonian potential that non-relativistic matter will feel, we must make a change of coordinates bringing the metric back to the standard spherically symmetric form. We make the transformation

$$r^{}=r\sqrt{\mathrm{\Omega }(r)}$$
(6)
which is imposed to yield

$$d\tau ^2=B^{}(r^{})dt^2\frac{dr_{}^{}{}_{}{}^{2}}{B^{}(r^{})}r_{}^{}{}_{}{}^{2}\left(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2\right)$$
(7)
with

$$B^{}(r^{})=1\frac{2\beta }{r^{}}\gamma r^{}kr_{}^{}{}_{}{}^{2},$$
(8)
since $`\gamma `$ is now negative. This yields the set of equations

$`\mathrm{\Omega }(r)B(r)`$ $`=`$ $`B^{}(r^{})`$ (9)
$`{\displaystyle \frac{\mathrm{\Omega }(r)dr^2}{B(r)}}`$ $`=`$ $`{\displaystyle \frac{dr_{}^{}{}_{}{}^{2}}{B^{}(r^{})}}.`$ (10)
These equations can be easily solved in first order, weak field perturbation theory. The weak field limit is phenomenologically, abundantly justified. We take $`\mathrm{\Omega }(r)=1+\psi (r)`$ and assume that all $`r`$ dependent functions are small compared to 1. Then expanding (9) to first order immediately yields the solution

$$\psi (r)=2\gamma r.$$
(11)
Using

$$\frac{dr^{}}{dr}=\frac{d}{dr}\sqrt{12\gamma r}r\frac{d}{dr}(1\gamma r)r=12\gamma r$$
(12)
and using (6) and noting that $`r`$ may be replaced with $`r^{}`$ in terms that are already of first order, yields,

$`{\displaystyle \frac{\mathrm{\Omega }(r)dr^2}{B(r)}}`$ $``$ $`{\displaystyle \frac{(12\gamma r)dr^2}{1\frac{2\beta }{r}+\gamma rkr^2}}`$ (13)
$``$ $`{\displaystyle \frac{dr_{}^{}{}_{}{}^{2}}{1\frac{2\beta }{r^{}}\gamma r^{}kr_{}^{}{}_{}{}^{2}}}={\displaystyle \frac{dr_{}^{}{}_{}{}^{2}}{B^{}(r^{})}}`$ (14)
verifying (10).

Thus we see that the sign of $`\gamma `$ can be reversed with a coordinate change and a concomitant change of the conformal gauge in first order weak field perturbation. The null geodesics are insensitive to these transformations, indeed, if $`r=r(\theta )`$ is the original orbit equation for the null geodesic, then it is just transformed to $`r^{}(\theta )=r^{}(r(\theta ))`$. The range of $`\theta `$ in the orbit equation, which determines the scattering angle, is unchanged, hence neither is the scattering angle or the deflection. The sign of $`\gamma `$ being negative, this choice corresponds to the null geodesics behaving as if there is additional attraction in comparison to the usual Newton/Einstein case. There are of course local changes in the trajectory, but these do not affect the global behaviour of the null orbits.

The massive geodesics on the other hand are sensitive to conformal transformations. Indeed the non-relativistic massive geodesics for the metric (7) are also deflected as if there is additional attraction since the sign of $`\gamma `$ is negative, as observed in  and which consequently served as the cornerstone of the conformal gravity program.

In conclusion, it is possible to choose a metric and conformal gauge for which conformal gravity will yield additional attraction above the Newton/Einstein result for both massless and non-relativistic, massive particles. Hence the possibility that conformal gravity could give a solution of the missing matter problem is not closed. Much work should be done, however, before embracing this theory as a viable phenomenological alternative. The exact scenario of the spontaneous or explicit breaking of the conformal invariance must be examined in detail. Unless we understand the mechanism, the theory will remain of only nominal theoretical interest.

Acknowledgements

We thank NSERC of Canada for partial financial support and the Department of Physics, National Central University, Chungli, Taiwan, where this paper was written.


text:
# Supersymmetry solution for finitely extensible dumbbell model

## I Introduction

Simple mechanical models are helpful for understanding the complex dynamical behaviour of polymer molecules. The so-called elastic dumbbell constitutes the simplest model that captures the effects of stretching and orientation of the polymer . In this model the polymer is represented by two beads which are connected by a spring. The classical Hookean dumbbell model  is characterised by a linear spring force. Its relaxation times decrease like the inverse of the mode number and the relaxation modulus is of a single exponential form. As noted by many authors, the finite extensibility of the molecules has to be included in order to avoid unphysical behaviour such as infinite viscosities predicted by the Hookean model. It is even argued “that taking into account the nonlinear stretching of the polymer molecules is the most important correction that should be made to the Hookean model in order to describe real systems” (Ref. , Chap. 13, p. 81). Molecular arguments suggest the finitely extensible spring force law to be given by the inverse Langevin function . Very often, rational  or Padé approximations  to the inverse Langevin function are used since they are analytically more tractable. All these approximations show the desired physical characteristics: a linear regime for small extensions of the spring and a divergence at a finite extension. These forces therefore prohibit the molecule to stretch beyond a certain value which is related to the total length of the chain.

In general, the finite extensibility of the molecules prohibits analytical expressions for the spectrum of relaxation times and the relaxation modulus. The shear relaxation modulus for the Warner spring force (so-called FENE model, ) has been determined numerically . The results are very close to those obtained for the Peterlin approximation to the Warner force (FENE-P model, ). The effect of the finite extensibility on the spectrum of relaxation times remains less studied. One might speculate that the finite extensibility decreases the slowest relaxation time as described by the Peterlin approximation. However, knowledge of the whole spectrum of relaxation times would be very desirable for understanding the dynamics of dumbbell models. It would allow, for example, to identify the number of relevant slow modes that contribute to the polymer contribution to the stress. Work in this direction was performed in  where the universal constitutive equation of dilute polymeric solutions was found. There, it was assumed that in the limit of low Deborah and Weissenberg number the modes corresponding to the two lowest eigenvalues dominate.

Here, we compute the relaxation times, eigenfunctions and the relaxation modulus of a finitely extensible dumbbell model exactly. The model is a one-dimensional version of the finitely extensible dumbbell models that uses trigonometric functions to approximate the inverse Langevin force. It is demonstrated that this force approximates the inverse Langevin function as good or even better than the Warner force frequently used in Brownian dynamics simulations. The present study therefore clarifies the effect of the finitely extensible spring force on the spectrum of relaxation times, eigenfunctions and the relaxation modulus.

## II Model description

Let $`f(𝒒;t)`$ denote the distribution function of the connector vector $`𝒒`$ of the dumbbell at time $`t`$. The time evolution of $`f`$ is given by the well-known Fokker-Planck equation

$$_tf=D\left\{[U]+\right\}fL_{\mathrm{FP}}f$$
(1)
where we have assumed that the spring force $`𝑭`$ can be derived from the dimensionless potential $`U`$ via $`𝑭(𝒒)=k_\mathrm{B}TU(𝒒)`$. The diffusion constant $`D=2k_\mathrm{B}T/\zeta `$ is assumed to be independent of $`𝒒`$ which means we neglect hydrodynamic interactions. Note that the equilibrium distribution $`f_{\mathrm{eq}}`$ is given by $`f_{\mathrm{eq}}=N_{\mathrm{eq}}e^U`$, $`N_{\mathrm{eq}}`$ being a normalization constant. As is well-known the Fokker-Planck equation (1) can be transformed into an equivalent imaginary time Schrödinger equation . If $`\psi `$ is a solution to

$$_t\psi =DH\psi $$
(2)
with the operator $`H=(1/D)e^{U/2}L_{\mathrm{FP}}e^{U/2}`$, then the function $`f(𝒒;t)=\psi (𝒒;t)f_{\mathrm{eq}}^{1/2}(𝒒)`$ solves the original Fokker-Planck equation (1). In particular, solving the eigenvalue problem $`H\psi _n=ϵ_n\psi _n`$ leads immediately to the solution of the eigenvalue problem of the original Fokker-Planck operator $`L_{\mathrm{FP}}f_n=\lambda _n^1f_n`$ with $`\lambda _n^1=Dϵ_n`$ and $`f_n=\psi _nf_{\mathrm{eq}}^{1/2}`$, where factors $`D`$ are introduced for later convenience. Inserting (1) for the Fokker-Planck operator $`H`$ yields the form of the Schrödinger operator

$$H=^2+V(𝒒),$$
(3)
where the Schrödinger potential $`V`$ is given by $`V=\mathrm{\Phi }^2\mathrm{\Phi }`$, with $`\mathrm{\Phi }=U/2`$. Therefore, each solvable potential for the Schrödinger equation serves as a solvable model for the Fokker-Planck equation. In the following, we benefit from the fact that simple Schrödinger potentials $`V`$ can have quite complicated counterparts $`U`$ for the Fokker-Planck equation. Note also that operator $`H`$ is Hermitian which is in general not the case for $`L_{\mathrm{FP}}`$.

The simplest exactly solvable potential for the Schrödinger equation related to finite extensibility is the three dimensional spherical symmetric infinite well potential. Up to a constant, the corresponding potential $`U`$ of the Fokker-Planck equation is given by $`U=2\mathrm{ln}\psi _{\mathrm{eq}}`$, where $`\psi _{\mathrm{eq}}`$ is the ground state of the Schrödinger problem. In the present case, the ground state is known to be $`\psi _{\mathrm{eq}}\mathrm{sin}(\pi q/q_0)/q`$ with $`q=|𝒒|`$ leading to the force law

$$𝑭(𝒒)=2k_\mathrm{B}T\frac{1(\pi q/q_0)\mathrm{cot}(\pi q/q_0)}{q}\frac{𝒒}{q}.$$
(4)
Since the force (4) is one-parametric, the “spring constant” is no adjustable parameter but determined by the maximum extension $`q_0`$. If we define the spring constant $`h`$ in the linear regime of (4), the dimensionless finite extensibility parameter $`b=hq_0^2/(k_\mathrm{B}T)`$ is fixed: $`b=b_3=2\pi ^2/3`$. Parameter $`b`$ is known to denote roughly the number of monomer units represented by the dumbbell and should therefore be a large number. Thus, (4) is of limited use. However, if $`b=b_3`$ is acceptable the force (4) compares very well to the inverse Langevin force and Cohen’s Padé approximation, as shown in the inset of Fig. 1. In this case, exact eigenfunctions of the Schrödinger operator are known to be the spherical Bessel functions and the exact eigenvalues are given by the zeros of this functions.

Let us now consider the one-dimensional case. One-dimensional nonlinear dumbbell models serve as toy models for various approximations that have been proposed to obtain closed constitutive equations for polymer solutions . In addition, the one-dimensional models are believed to describe the elongational behaviour of the polymer molecule reasonably . The one-dimensional version of the operators $`L_{\mathrm{FP}}`$ and $`H`$ is obtained by replacing $`𝒒`$ with $`q`$ and $``$ with $`d/dq`$.

To account for the finite extensibility we propose the following two parameter family of force laws

$$F(q)=hq_1\mathrm{tan}(q/q_1),\text{for}q_0<q<q_0$$
(5)
where $`q_1=2q_0/\pi `$ is determined by the maximum extension $`q_0`$ and $`h`$ denotes the “spring constant” in the sense $`F(q)hq`$ for $`q0`$. In Fig. 1 we plot the force $`F/h`$ as a function of the reduced extension $`q/q_0`$. Fig. 1 also shows the Warner force  $`F(q)=hq/(1q^2/q_0^2)`$, the inverse Langevin force $`(h/3)L^1(q/q_0)`$, with $`L(x)=\mathrm{coth}(x)1/x`$ and Cohen’s Padé approximation $`F(q)=(hq/3)(3q^2/q_0^2)/(1q^2/q_0^2)`$ . It is easily verified that these approximations give the correct limiting behaviour not only for small but also near the maximum extension. In addition, as can be seen from Fig. 1, the force (5) serves as an even better approximation to the inverse Langevin force as does the Warner force. Therefore we accept (5) as a reasonable force law to substitute for the inverse Langevin force in one dimension.

The potential $`U`$ from which the force (5) can be derived is $`U(q)=b_1\mathrm{ln}\mathrm{cos}(q/q_1)`$, where we introduced the dimensionless parameter $`b_1=hq_1^2/(k_\mathrm{B}T)=(2/\pi )^2b`$. The equilibrium distribution function of the Fokker-Planck equation is $`f_{\mathrm{eq}}(q)=N_{\mathrm{eq}}\mathrm{cos}^{b_1}(q/q_1)`$, with $`N_{\mathrm{eq}}^1=q_12^{b_1}B[(b_1+1)/2,(b_1+1)/2]`$, where $`B[x,y]`$ is the Beta function. The equilibrium distribution is very close to the one corresponding to the Warner force $`f_{\mathrm{eq}}(q)=N_{\mathrm{eq}}(1q^2/q_0^2)^{b/2}`$, with $`N_{\mathrm{eq}}^1=q_0B[1/2,(b+2)/2]`$.

## III Supersymmetry solution

To obtain exact eigenvalues and eigenfunctions of the Fokker-Planck equation (1) we exploit supersymmetry methods for the corresponding Schrödinger problem. Note, that the model force (5) belongs to the class of so-called Poeschl-Teller potentials which may also be solved exactly by Schrödinger’s factorization method . Consider the original potential $`U`$ together with the inverted potential $`U`$ and denote by $`f_{}`$ the solution for the potential $`U`$ and $`f_+`$ for $`U`$. As done before, defining $`f_\pm =\sqrt{N_{\mathrm{eq}}}e^{\pm U/2}\psi _\pm `$ we arrive at $`(1/D)_t\psi _\pm =H_\pm \psi _\pm `$ with operators $`H_\pm =^2+V_\pm (q)`$. The “partner potentials” are given by $`V_\pm (q)=\mathrm{\Phi }^2(q)\pm \mathrm{\Phi }^{}(q)`$ with $`\mathrm{\Phi }=U^{}/2`$, where $`\mathrm{\Phi }^{}`$ and $`U^{}`$ denote the derivative of $`\mathrm{\Phi }`$ resp. $`U`$ with respect to $`q`$. Schrödinger operators $`H_\pm `$ are very well studied in the literature since they occur in the so-called “Witten model” , which is the simplest model that shows all typical features of supersymmetric quantum mechanics.

Rewriting the force (5) as $`F=2k_\mathrm{B}T\mathrm{\Phi }`$ we obtain $`\mathrm{\Phi }_\alpha (q)=\alpha \mathrm{tan}(q/q_1)`$, for $`q_0<q<q_0`$. We denote explicitly the dependence on the parameter $`\alpha =b_1/(2q_1)`$ which is proportional to the spring constant $`h`$. To proceed further we take advantage of the concept of “shape invariance” introduced by Gendenshteîn . The partner potentials $`V_\pm (\alpha ,q)=\mathrm{\Phi }_\alpha ^2(q)\pm \mathrm{\Phi }_\alpha ^{}(q)`$ are called shape invariant if they are related by $`V_+(\alpha _k,q)=V_{}(\alpha _{k+1},q)+R(\alpha _{k+1})`$ where the new parameter $`\alpha _{k+1}`$ is a function of $`\alpha _k`$. In our case, $`\mathrm{\Phi }_\alpha =\alpha \mathrm{tan}(q/q_1)`$, shape invariance is easily verified for $`\alpha _{k+1}=\alpha _k+1/q_1`$, $`\alpha _0=\alpha `$ and $`R(\alpha _{k+1})=\alpha _{k+1}^2\alpha _k^2`$. For shape invariant potentials it is possible to define a series of operators that are isospectral except for the lowest eigenvalue which is $`_{k=1}^nR(\alpha _k)`$. Going back to the original Schrödinger operator with lowest eigenvalue zero the complete spectrum is found to be $`ϵ_n=_{k=1}^nR(\alpha _k)`$ . Inserting the special form of $`R`$ and multiplying by $`D`$, the exact (inverse) eigenvalues of the Fokker-Planck operator (1) are given by

$$\lambda _n=2\lambda _\mathrm{H}\left[n+n^2/b_1\right]^1,n=1,2,\mathrm{}$$
(6)
where $`\lambda _\mathrm{H}=\zeta /(4h)`$ is the time constant of the Hookean dumbbell. Note, that the eigenvalue zero corresponds to the equilibrium distribution. Corresponding eigenfunctions of operator $`H`$ are

$$\psi _n(q)=C_na_0^{}a_1^{}\mathrm{}a_{n1}^{}\mathrm{cos}^{n+b_1/2}(q/q_1)$$
(7)
where $`a_k^{}=d/dq+\mathrm{\Phi }_{\alpha _k}`$ are generalized creation operators and $`C_n`$ denote normalization constants . Remember that eigenfunctions of the original Fokker-Planck operator are obtained by $`f_n=\psi _nf_{\mathrm{eq}}^{1/2}`$. The first eigenfunctions read: $`f_0=f_{\mathrm{eq}}`$, $`f_1=N_1\mathrm{sin}(q/q_1)f_{\mathrm{eq}}`$, $`f_2=N_2[(b_1+1)\mathrm{sin}^2(q/q_1)\mathrm{cos}^2(q/q_1)]f_{\mathrm{eq}}`$ with $`N_1^2=b_1+2`$ and $`N_2^2=(b_1+4)/[2(b_1+1)]`$. Note, that the eigenfunctions are orthonormal, $`f_n,f_m=\delta _{nm}`$, in the scalar product

$$g,h=_{q_0}^{q_0}f_{\mathrm{eq}}^1(q)g(q)h(q)𝑑q.$$
(8)
The linear viscoelastic behaviour can be obtained from linear response theory. For a given flow kinematics the perturbation of the Fokker-Planck operator is given by $`L_{\mathrm{ext}}(t)=_q\kappa (t)q`$ where $`\kappa (t)`$ denotes the velocity gradient.

The polymer contribution to the stress, $`\tau /(k_\mathrm{B}T)=\sigma ,f1`$, $`\sigma =f_{\mathrm{eq}}qU^{}`$, is found to be

$$\tau /(k_\mathrm{B}T)=_{\mathrm{}}^tG(tt^{})\kappa (t^{})𝑑t^{},$$
(9)
where the relaxation modulus $`G(t)`$ is given by the equilibrium correlation function $`G(t)=\sigma ,e^{L_{\mathrm{FP}}t}\sigma 1`$. Inserting the completeness relation of the eigenfunctions, $`\delta (qq^{})=_n\psi _n(q)\psi _n(q^{})=f_{\mathrm{eq}}^1(q^{})_nf_n(q)f_n(q^{})`$, we obtain

$$G(t)=\underset{n=1}{\overset{\mathrm{}}{}}\left|\sigma ,f_{2n}\right|^2e^{t/\lambda _{2n}},$$
(10)
which is of the form of the “generalized Maxwell model”. Note, that only even terms are included in (10) since they correspond to symmetric eigenfunctions whereas odd terms, corresponding to antisymmetric eigenfunctions, vanish by symmetry. Analytical expressions for $`G(t)`$ resulting from dumbbell models are known for linear springs (Hookean dumbbell) where $`G(t)=e^{t/\lambda _\mathrm{H}}`$. For the Warner spring force no such expression is available. The Peterlin approximation  leads also to a single-exponential form with relaxation time $`\lambda _\mathrm{H}[1+1/b]^2`$. The relaxation modulus (10) in our model shows a spectrum of relaxation times $`\lambda _n`$, the lowest one being $`\lambda _2=\lambda _\mathrm{H}[1+2/b_1]^1`$. The finite extensibility of the dumbbell therefore reduces the longest relaxation time when compared to the Hookean dumbbell. Only in the limit $`b\mathrm{}`$ the result for the Hookean dumbbell is recovered. As can be seen from Fig. 2, $`G(t)`$ is dominated by the lowest relaxation time since the relative weight of the higher modes decreases rapidly. For $`b=20`$, the relative weight of the first two modes is $`99\%`$, for $`b=50`$ the first mode carries $`99\%`$ and for $`b\mathrm{}`$ all weights except for the first vanish, reflecting the single-exponential form of the Hookean dumbbell.

The relaxation modulus determines the zero-elongational viscosity $`\overline{\eta }_0=_0^{\mathrm{}}G(t)𝑑t`$. As shown in , $`\overline{\eta }_0`$ can be expressed in terms of second moments of the equilibrium distribution function for arbitrary dumbbell models. In our model, the resulting integral cannot be done analytically.

Knowledge of the exact relaxation times and eigenfunctions (6) and (7) can further be used to integrate the kinetic equation (1) in the presence of a given flow field. Expanding the distribution function into eigenfunctions, $`f(q;t)=_nc_n(t)f_n(q)`$, leads to

$$\dot{c}_n=\lambda _n^1c_n\kappa (t)\underset{k}{}A_{nk}c_k,$$
(11)
with $`A_{nk}=f_n,_qqf_k`$ a constant matrix independent of the flow. Coefficients $`c_n(t)`$ are defined by $`c_n(t)=f_n,f(q;t)`$. Eqs. (11) are an equivalent formulation of the kinetic equation (1). Conservation of total probability is guaranteed by $`\dot{c}_0=0`$, which is easily verified from (11). The stress $`\tau `$ is determined by the coefficients $`c_n`$ via $`\tau /(k_\mathrm{B}T)=_{n=1}^{\mathrm{}}\sigma ,f_{2n}c_{2n}`$. Because of symmetry, only even terms contribute to $`\tau `$. Moreover, Eqs. (11) for even coefficients decouple from those for odd coefficients since $`A_{nk}`$ is nonzero only if $`n`$ and $`k`$ are both even or both odd. Therefore, if we are interested only in the polymer contribution to the stress, only even coefficients have to be included in (11). Given $`A_{nk}`$, Eqs. (11) can be integrated for a finite set of eigenfunctions quite efficiently using standard integrators for ordinary differential equations. How many eigenfunctions have to be included in order to obtain a desired accuracy is a delicate question. A simple estimation is the following. To a first approximation, coefficients $`c_n`$ are still of exponential form but with a different time constant $`\lambda _n^1+\kappa A_{nn}`$. Modes that decrease quadratically with number $`n`$ will then be negligible. This leads to a minimal number of modes $`n_c`$ that satisfies $`n_c^2/b_1=n_c+2\lambda _\mathrm{H}|\kappa (t)A_{n_cn_c}|`$. If this number is not too large, method (11) is much more efficient than Brownian dynamics simulation.

## IV Conclusions

We have obtained the spectrum of relaxation times, the eigenfunctions and the relaxation modulus for a one-dimensional finitely extensible dumbbell model exactly. It is observed that the relaxation times are reduced when compared to the Hookean dumbbell model. Moreover, the decrease of the relaxation times with mode number $`n`$ shows a crossover from $`n^1`$ to $`n^2`$ at the value of the finite extensibility parameter $`b`$. The relaxation modulus is obtained in the form of the “generalized Maxwell model”. For rather large values of the finite extensibility parameter $`b`$, which are commonly believed to be physically meaningful (see, e.g. ), the linear viscoelastic behaviour is well described by a single relaxation time. This is because the relative weight of the higher modes is decreasing rapidly with increasing $`b`$. Therefore, it is not surprising that the linear viscoelastic regime is well described by the Peterlin approximation as found in . To calculate the non-linear dynamical response to a given flow kinematics we propose a numerical integration scheme using the expansion of the distribution function into the exact eigenfunctions in the absence of flow.

All these results are limited to the one-dimensional case. In three dimensions the class of solvable potentials is much more restricted. In particular we did not find any relevant two parameter potential that can be solved exactly. However, there is some evidence, that “a simple one-dimensional version of the FENE theory captures qualitatively and quantitatively the elongational behaviour of the actual three-dimensional theory” (Ref. , Introduction).

## Acknowledgment

The authors gratefully acknowledge valuable discussions with H. C. Öttinger.