text:
# Beyond Fermi pseudopotential: a modified GP equation

## I Introduction

The Gross-Pitaevskii equation (GPE) has in recent years played a central role in our understanding of gaseous Bose-Einstein condensates (BEC) dal99 ; leg01 . Its astonishing success has to do with the fact that in most experiments to date, both the gaseous parameter, $`na^3`$ (here $`n`$ is the number density and $`a`$ is the scattering length), and the parameter $`a/L_0`$, which characterizes the strength of confinement (with $`L_0`$ being a measure of the size of the trap), have been very small. Recently, a new breed of experiments feshbe  have emerged that take advantage of Feshbach resonances in cold atomic collision tie92  to make the scattering length $`a`$ tunable through a magnetic field. Such experiments promise to take us into parameter regimes in which the validity of GPE and its underlying assumptions, such as the shape-independent pseudopotential approximation dal99 ; leg01 , have to be carefully examined.

Most of existing theories going beyond the GPE bra97 ; fab99 ; blu01  have focused on quantum fluctuation that is responsible for the lowest order correction in $`na^3`$. But it has become clear that eventually, a better description of atom-atom interaction, beyond the pseudopotential approximation that is the basis of the GP theory, will be required. For homogeneous systems, this has been demonstrated in a recent work by Cowell et al. cow02 , which shows that different potentials having the same scattering length $`a`$ can lead to vastly different ground state energy for $`na^3`$ of the order of 0.05 or greater. For inhomogeneous systems, recent works on two atoms in a trap have shown that the shape-independent approximation becomes less valid under strong confinement tie00 ; bol02 ; blu02 .

Neither of these results are surprising. The shape-independent approximation assumes the two-body scattering amplitude to be a constant over the momentum (or energy) range of interest. For larger $`na^3`$ or stronger confinement, a greater range of momentum states become involved and this assumption becomes less applicable. This point has already been illustrated for two-atom systems under strong confinement, where it has been shown that inclusion of the energy-dependence of scattering amplitude leads to much improved results for the energy levels of the system tie00 ; bol02 ; blu02 . The goal of this article is to show how a similar idea can be generalized to a N-body system.

In light of the success of the GP theory dal99 ; leg01 , we focus here on a perturbative approach around the GPE. An effective potential is proposed that incorporates the energy dependence of the two-body scattering amplitude through an effective-range expansion. A corresponding modified GP equation is derived. By taking advantage of the relationship between effective range and scattering length as implied by the quantum-defect theory (QDT) gao98b ; gao01 , no new inputs from scattering calculations are required other than the standard scattering length and the $`C_6`$ coefficient. For the ground state energy of an inhomogeneous system, our numerical results show that the modified equation leads to better agreements with diffusion Monte Carlo (DMC) results blu01  than theories that consider only the quantum fluctuation correction. More importantly, the theory introduces the concept of shape-dependent confinement correction, an origin of shape-dependence on the interaction potential that is absent in a homogeneous system.

## II Modified GP Equation

In a second-quantization formulation, a system of N bosons in an external trapping potential $`V_{\text{ext}}(𝐫)`$ is described by a Hamiltonian

$$\widehat{H}=𝑑𝐫\widehat{\mathrm{\Psi }}^{}(𝐫)\left[\frac{\mathrm{}^2}{2m}^2+V_{\text{ext}}(𝐫)\right]\widehat{\mathrm{\Psi }}(𝐫)+\frac{1}{2}𝑑𝐫_1𝑑𝐫_2\widehat{\mathrm{\Psi }}^{}(𝐫_1)\widehat{\mathrm{\Psi }}^{}(𝐫_2)V(𝐫_1𝐫_2)\widehat{\mathrm{\Psi }}(𝐫_2)\widehat{\mathrm{\Psi }}(𝐫_1).$$
(1)
Here $`V`$ is the interaction between particles. $`\widehat{\mathrm{\Psi }}(𝐫)`$ and $`\widehat{\mathrm{\Psi }}^{}(𝐫)`$ are the bosonic annihilation and creation operators that satisfy the commutation relations

$`[\widehat{\mathrm{\Psi }}^{}(𝐫_1),\widehat{\mathrm{\Psi }}^{}(𝐫_2)]`$ $`=[\widehat{\mathrm{\Psi }}(𝐫_1),\widehat{\mathrm{\Psi }}(𝐫_2)]=0,`$ (2)
$`[\widehat{\mathrm{\Psi }}(𝐫_1),\widehat{\mathrm{\Psi }}^{}(𝐫_2)]`$ $`=\delta (𝐫_1𝐫_2).`$ (3)
Taking the inter-particle potential $`V`$ to be that of a mean field

$$V_{\text{MF}}=g\delta (𝐫_1𝐫_2),$$
(4)
where

$$g=\frac{4\pi \mathrm{}^2a}{m},$$
(5)
and ignoring quantum fluctuation, one arrives at the usual GP equation dal99 ; leg01

$$i\mathrm{}\frac{}{t}\mathrm{\Phi }(𝐫,t)=\left[\frac{\mathrm{}^2}{2m}^2+V_{\text{ext}}(𝐫)+g\left|\mathrm{\Phi }\right|^2\right]\mathrm{\Phi }(𝐫,t).$$
(6)
Here $`\mathrm{\Phi }=\widehat{\mathrm{\Psi }}`$ and is nomalized by $`𝑑𝐫|\mathrm{\Phi }|^2=N`$, with $`N`$ being the total number of particles.

From a two-body scattering point of view, the mean-field potential, given by Eq. (4), is such that it gives, in the first Born approximation, the correct foward scattering amplitude at zero energy. It is a shape-independent approximation that ignores completely the energy-dependence of the scattering amplitude. While this may be a good approximation for small values of $`na^3`$ and for weak confinement, it becomes less valid for strong confinement and/or greater $`na^3`$. This is true even at zero temperature, as confinement and/or quantum depletion invariably involve nonzero momentum states.

For a better treatment of atomic interaction, while preserving much of the structure of the GP theory, an effective interaction is proposed here that incorporates the energy dependence of the scattering amplitude through an effective-range expansion. Specifically, we look for a $`\widehat{V}_{\text{eff}}`$ that gives, in the first Born approximation, the correct real part of the forward scattering amplitude,

$$\frac{m}{4\pi \mathrm{}^2}𝑑𝐫e^{i𝐤𝐫}\widehat{V}_{\text{eff}}e^{+i𝐤𝐫}=\text{Re}[f(k,\theta =0)],$$
(7)
to the order of $`k^2`$ (the first order in energy). This requirement basically ensures a better optical potential for an atom moving in the medium of others rod67 .

For the range of energies of interest here, only the $`s`$ wave scattering is important. In this case, the scattering amplitude can be written in terms of the $`s`$ wave phase shift $`\delta _0(k)`$ as

$$f(k,\theta )=\frac{1}{k}\frac{1}{\mathrm{cot}\delta _0(k)i}.$$
(8)
From the standard effective-range expansion

$$k\mathrm{cot}\delta _0(k)=\frac{1}{a}+\frac{1}{2}r_ek^2+\mathrm{},$$
(9)
where $`r_e`$ is the effective range, we have, to the order of $`k^2`$,

$$\text{Re}[f(k,\theta )]=a+a^2(a\frac{1}{2}r_e)k^2.$$
(10)
It is straightforward to verify that the $`\widehat{V}_{\text{eff}}`$ satisfying Eq. (7) can be written as

$$\widehat{V}_{\text{eff}}=V_{\text{MF}}(𝐫_1𝐫_2)+\widehat{V}_{\text{mod}}\text{ }(𝐫_1𝐫_2),$$
(11)
where

$$\widehat{V}_{\text{mod}}=\frac{g_2}{2}\left[\delta (𝐫_1𝐫_2)_{𝐫_1𝐫_2}^2+_{𝐫_1𝐫_2}^2\delta (𝐫_1𝐫_2)\right],$$
(12)
and

$$g_2=\frac{4\pi \mathrm{}^2a^2(a\frac{1}{2}r_e)}{m}.$$
(13)
This potential is consistent with the pseudopotential expansion of Huang and Yang hua57 , except that it is Hermitian and is represented in the coordinate space.

With this effective interaction, the many-body Hamiltonian becomes

$`\widehat{H}`$ $`\widehat{H}_{\text{MF}}+\widehat{H}_{\text{mod}}`$
$`\left\{{\displaystyle 𝑑𝐫\widehat{\mathrm{\Psi }}^{}(𝐫)\left[\frac{\mathrm{}^2}{2m}^2+V_{\text{ext}}(𝐫)\right]\widehat{\mathrm{\Psi }}(𝐫)}+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫_1𝑑𝐫_2\widehat{\mathrm{\Psi }}^{}(𝐫_1)\widehat{\mathrm{\Psi }}^{}(𝐫_2)V_{\text{MF}}(𝐫_1𝐫_2)\widehat{\mathrm{\Psi }}(𝐫_2)\widehat{\mathrm{\Psi }}(𝐫_1)}\right\}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫_1𝑑𝐫_2\widehat{\mathrm{\Psi }}^{}(𝐫_1)\widehat{\mathrm{\Psi }}^{}(𝐫_2)\widehat{V}_{\text{mod}}(𝐫_1𝐫_2)\widehat{\mathrm{\Psi }}(𝐫_2)\widehat{\mathrm{\Psi }}(𝐫_1)}.`$ (14)
The term associated with $`\widehat{V}_{\text{mod}}`$ can be simplified in the two-body centre-of-mass frame \[$`𝐫=𝐫_1𝐫_2`$, $`𝐑=(𝐫_1+𝐫_2)/2`$\], in which $`\widehat{H}_{\text{mod}}`$ takes the form:

$$\widehat{H}_{\text{mod}}=\frac{1}{4}g_2𝑑𝐑𝑑𝐫\widehat{\mathrm{\Psi }}^{}(𝐑+\frac{𝐫}{2})\widehat{\mathrm{\Psi }}^{}(𝐑\frac{𝐫}{2})\left[\delta (𝐫)_𝐫^2+_𝐫^2\delta (𝐫)\right]\widehat{\mathrm{\Psi }}(𝐑\frac{𝐫}{2})\widehat{\mathrm{\Psi }}(𝐑+\frac{𝐫}{2}).$$
(15)
After integrating over $`𝐫`$ and some integrations by parts, $`\widehat{H}_{\text{mod}}`$ can be written as

$$\widehat{H}_{\text{mod}}=\frac{1}{4}g_2𝑑𝐑\widehat{\mathrm{\Psi }}^{}(𝐑)\left[^2\left(\widehat{\mathrm{\Psi }}^{}(𝐑)\widehat{\mathrm{\Psi }}(𝐑)\right)\right]\widehat{\mathrm{\Psi }}(𝐑).$$
(16)
Note that the Laplacian operator operates only on the number-density operator $`\widehat{\mathrm{\Psi }}^{}(𝐑)\widehat{\mathrm{\Psi }}(𝐑)`$. Simplification of $`\widehat{H}_{\text{MF}}`$ follows the usual procedure, and leads finally to

$`\widehat{H}`$ $`={\displaystyle 𝑑𝐑\widehat{\mathrm{\Psi }}^{}(𝐑)\left[\frac{\mathrm{}^2}{2m}^2+V_{\text{ext}}(𝐑)\right]\widehat{\mathrm{\Psi }}(𝐑)}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle }d𝐑\widehat{\mathrm{\Psi }}^{}(𝐑)[g\left(\widehat{\mathrm{\Psi }}^{}(𝐑)\widehat{\mathrm{\Psi }}(𝐑)\right)`$
$`+{\displaystyle \frac{1}{2}}g_2^2\left(\widehat{\mathrm{\Psi }}^{}(𝐑)\widehat{\mathrm{\Psi }}(𝐑)\right)]\widehat{\mathrm{\Psi }}(𝐑).`$ (17)
From this Hamiltonian, both the static and the dynamic properties of a many-atom Bose system can be studied. We focus here on the ground state energy of the system.

Ignoring both quantum fluctuation and $`\widehat{H}_{\text{mod}}`$, Eq. (17) leads to the standard GP energy functional

$$E_{\text{GP}}\left[\mathrm{\Phi }\right]=𝑑𝐫\left[\frac{\mathrm{}^2}{2m}\left|\mathrm{\Phi }\right|^2+V_{\text{ext}}(𝐫)\left|\mathrm{\Phi }\right|^2+\frac{1}{2}g\left|\mathrm{\Phi }\right|^4\right],$$
(18)
which corresponds to a time-independent GP equation for the ground-state wave function with a chemical potential $`\mu `$

$$\left[\frac{\mathrm{}^2}{2m}^2+V_{\text{ext}}(𝐫)+g\left|\mathrm{\Phi }\right|^2\right]\mathrm{\Phi }(𝐫,t)=\mu \mathrm{\Phi }(𝐫,t).$$
(19)
Inclusion of the lowest-order term due to quantum fluctuation leads to an energy functional lhy57 ; fab99

$`E_{\text{MGPI}}\left[\mathrm{\Phi }\right]`$ $`={\displaystyle }d𝐫\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\right|\mathrm{\Phi }|^2+V_{\text{ext}}(𝐫)|\mathrm{\Phi }|^2`$
$`+{\displaystyle \frac{1}{2}}g\left|\mathrm{\Phi }|^4(1+{\displaystyle \frac{128}{15}}{\displaystyle \frac{1}{\sqrt{\pi }}}a^{3/2}|\mathrm{\Phi }\left|\right)\right],`$ (20)
where the extra term due to quantum fluctuation was first derived by Lee, Huang and Yang lhy57 , and will be called the LHY term. This functional corresponds to a modified GP equation that we call MGPI fab99 ; blu01

$$\left[\frac{\mathrm{}^2}{2m}^2+V_{\text{ext}}(𝐫)+g\left|\mathrm{\Phi }\right|^2+g_1\left|\mathrm{\Phi }\right|^3\right]\mathrm{\Phi }(𝐫,t)=\mu \mathrm{\Phi }(𝐫,t),$$
(21)
where

$$g_1=g\frac{32}{3\sqrt{\pi }}a^{3/2}.$$
(22)
Inclusion of both quantum fluctuation and $`\widehat{H}_{\text{mod}}`$ to the lowest order leads to an energy functional

$`E_{\text{MGPII}}\left[\mathrm{\Phi }\right]`$ $`={\displaystyle }d𝐫\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\right|\mathrm{\Phi }|^2+V_{\text{ext}}(𝐫)|\mathrm{\Phi }|^2`$
$`+{\displaystyle \frac{1}{2}}g|\mathrm{\Phi }|^4+{\displaystyle \frac{2}{5}}g_1|\mathrm{\Phi }|^5+{\displaystyle \frac{1}{4}}g_2\left|\mathrm{\Phi }|^2^2\left(\right|\mathrm{\Phi }|^2)\right].`$ (23)
The corresponding modified GP equation, which we call MGPII, is

$`[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_{\text{ext}}(𝐫)+g|\mathrm{\Phi }|^2+g_1|\mathrm{\Phi }|^3`$
$`+{\displaystyle \frac{g_2}{2}}^2\left(\right|\mathrm{\Phi }|^2)]\mathrm{\Phi }(𝐫,t)=\mu \mathrm{\Phi }(𝐫,t).`$ (24)
MGPII incorporates the shape dependence on the interaction potential through the effective range $`r_e`$. This additional parameter does not, however, add much extra complexity to the theory, as $`r_e`$ and $`a`$ are generally related. For a hard sphere potential, we have the well-know relation

$$r_e=(2/3)a.$$
(25)
For atoms with van der Waals interaction, $`r_e`$ can be determined from $`a`$ by gao98b

$$r_e/\beta _6=\left(\frac{2}{3x_e}\right)\frac{1}{(a/\beta _6)^2}\left\{1+\left[1x_e(a/\beta _6)\right]^2\right\},$$
(26)
where $`\beta _6=\left(mC_6/\mathrm{}^2\right)^{1/4}`$ is a length scale associated with the van der Waals interaction, and $`x_e\left[\mathrm{\Gamma }(1/4)\right]^2/\left(2\pi \right)`$. Similar relations for other long-range potentials should also exist, an assertion that can be deduced from a more general QDT consideration gao01 .

The physical meaning of the correction due to $`\widehat{H}_{\text{mod}}`$ can be better understood by estimating the order-of-magnitude of each term in Eq. (23) or (24). Consider a system of $`N`$ bosons confined to a length scale $`L_0`$. The order-of-magnitude for the number density is $`nN/(L_0)^3`$. Let $`ϵ_{\text{MF}}na^3(\mathrm{}^2/2m)(1/a)^2`$ to represent the energy per particle due to the mean field. It is easy to show that the LHY quantum fluctuation term is of the order of $`ϵ_{\text{LHY}}ϵ_{\text{MF}}(na^3)^{1/2}`$, while the $`\widehat{H}_{\text{mod}}`$ term is the order of $`ϵ_{\text{mod}}ϵ_{\text{MF}}(a/L_0)^2`$. It is thus clear that the two corrections are of different origin. The $`\widehat{H}_{\text{mod}}`$ correction has an order-of-magnitude that is determined primarily by the strength of confinement. It depends also on the shape of the interaction potential through $`g_2`$ and $`r_e`$. We will call this correction the shape-dependent confinement correction. It is a new source of shape dependence that is absent in a homogeneous system ($`L_0\mathrm{}`$ with $`n`$ being fixed). In this regard, it differs substantially from the shape dependence that can come at higher densities from two-body correlation at shorter length scales cow02 ; gao03 .

A further comparison of the magnitudes of the LHY and the shape-dependent confinement corrections reveals other interesting features. From the estimates above, we have $`ϵ_{\text{mod}}/ϵ_{\text{LHY}}(a/L_0)^2/(na^3)^{1/2}(a/L_0)^{1/2}/N^{1/2}`$. This means, for instance, that relative to the LHY correction, the shape-dependent confinement correction is more important for small $`N`$ than for large $`N`$, a point that will be further demonstrated numerically.

We have chosen, up to this point, a formulation that most closely parallels the standard GP formulation dal99 ; leg01 . For a small particle number $`N`$, a number-conserving Hartree-Fock formulation is in fact more appropriate esr97 ; leg01 . The corresponding results are easily derived. For example, at the level of MGPII, we have

$`E_{\text{MGPII}}\left[\varphi \right]`$ $`=N{\displaystyle }d𝐫\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\right|\varphi |^2+V_{\text{ext}}(𝐫)|\varphi |^2`$
$`+{\displaystyle \frac{1}{2}}g(N1)\left|\varphi \right|^4+{\displaystyle \frac{2}{5}}g_1(N1)^{3/2}\left|\varphi \right|^5`$
$`+{\displaystyle \frac{1}{4}}g_2(N1)\left|\varphi |^2^2\left(\right|\varphi |^2)\right].`$ (27)
and

$`[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_{\text{ext}}(𝐫)+g(N1)|\varphi |^2+g_1(N1)^{3/2}|\varphi |^3`$
$`+{\displaystyle \frac{g_2}{2}}(N1)^2\left(\right|\varphi |^2)]\varphi (𝐫,t)=\mu \varphi (𝐫,t),`$ (28)
where $`\varphi `$ is normalized according to $`𝑑𝐫|\varphi |^2=1`$. It is the Hartree-Fock formulation that we actually implement in our numerical calculations, but the difference is small for large $`N`$ (between $`N`$ and $`N1`$). Setting $`g_2`$, or both $`g_2`$ and $`g_1`$, to zero in Eqs. (27) and (28) leads to Hartree-Fock results in MGPI and GP, respectively.

## III Numerical Results

We present here our numerical results for the ground state energy of a BEC in different approximations. They are obtained using an imaginary time evolution method with ADI algorithm pre92 . The code has been tested by reproducing the GPE results of Ref. dal96 .

### III.1 Results for hard spheres

A many-body system of hard spheres is an ideal test-ground for various theoretical models, as accurate results for its ground-state energy are available from diffusion Monte Carlo (DMC) calculations for both homogeneous gio99  and inhomogeneous cases blu01 .

Figure 1 shows the relative differences of the energy functional between GP and DMC for $`N`$ bosons in a spherical harmonic trap characterized by a length scale $`L_0=\sqrt{\mathrm{}/2m\omega }`$. It is clear from these results that the ground state energy of the system can differ substantially from those predicated by the GP theory, for either strong confinement (large $`a/L_0`$) and/or high density (large $`N`$).

As shown in Ref. blu01 , most of these differences, at least for the range of parameters considered, can be accounted for by the LHY quantum fluctuation correction lhy57 , at a level corresponding to MGPI. This is illustrated in Figure 2(a).

In the same figure, it is also shown that a consistently better agreement with DMC is achieved by MGPII, which includes the shape-dependent confinement correction. This means that of the remaining difference between MGPI and DMC, at least part of it is shape-dependent and requires a better description of atom-atom interaction. Note that in the case of hard spheres, $`g_2`$ is always positive. The sign of the shape-dependent confinement correction is then determined by the expectation value of $`^2\left(\left|\varphi \right|^2\right)`$, which is negative for the ground state.

More results of MGPII are presented in Figure 2(b). It shows that at the level of MGPII, the relative difference from DMC has been reduced from up to 20% for the GP theory (see Figure 1) to less than 2% for the range of parameters considered.

### III.2 Results for atoms with van der Waals interaction

With the introduction of $`\widehat{H}_{\text{mod}}`$, the ground state energy becomes dependent upon the shape of the long-range potential through the relationship between $`r_e`$ and $`a`$. For atoms with van der Waals interaction, $`r_e`$ and $`a`$ are related by Eq. (26). Figure 3 is a graphic illustration of this relation. Note that $`r_e`$ diverges for $`a=0`$, corresponding to the fact that in this particular case, the energy dependence of the scattering amplitude around zero energy cannot be described by an effective-range expansion gao00 .

Also note that unlike the case of hard spheres, the $`g_2`$ for a van der Waals potential may become negative for sufficiently small $`a/\beta _6`$. Consequently, the shape-dependent confinement correction to the ground state energy may in principle be positive. \[Both points suggest that something interesting happens around $`a=0`$ that may deserve a separate investigation.\]

Figure 4(a) gives a comparison of ground state energies predicted by MGPII and by the GP theory for a realistic experimental configuration. Specifically, it is for a 85Rb condensate in a cylindrical trap with an aspect ratio of $`\epsilon =6.8/17.5`$ don01 . A $`C_6`$ value of 4698 a.u. is used van02 . It shows that deviations from the GP theory can become substantial for either high density and/or strong confinement. With the tunability of the scattering length via a Feshbach resonance don01 , such many-body effects beyond the GP theory may soon become observable.

Figure 4(b) shows more explicitly the shape-dependent confinement correction. Its magnitude increases with $`a/L_0`$, as expected. For a fixed $`a/L_0`$, its relative contribution is more significant for small $`N`$ than for large $`N`$, consistent with our earlier discussion. The corrections are mostly negative (except for very small $`a`$) since $`g_2`$ is greater than zero for most of the data points shown. \[Because $`L_0=34850`$ a.u. is much greater than $`\beta _6=164.2`$ a.u., most of the data points correspond to $`a/\beta _61`$. The $`r_e`$ is then approximately a constant $`r_e(3\pi )^1[\mathrm{\Gamma }(1/4)]^2\beta _6`$, and the corresponding $`g_2`$’s are positive.\]

## IV Conclusion

In conclusion, an effective interaction and the corresponding modified GP equation have been proposed that take into account the energy dependence of the two-body scattering amplitude through an effective-range expansion. The resulting theory, called MGPII, leads to better agreements with diffusion Monte Carlo calculations blu01  than either the GP theory, or MGPI that considers only the quantum fluctuation correction. The theory expands, considerably, the parameter space (specified by $`na^3`$ and $`a/L_0`$), in which a GP type of formulation can be applied. It introduces the concept of shape-dependent confinement correction and shows that for a fixed confinement, its relative contribution is more significant for small $`N`$ than for large $`N`$.

Finally, we point out that both MGPI and MGPII are perturbative expansions around the GP theory. They will eventually fail for sufficiently large $`na^3`$ and/or sufficiently strong confinement. For high density, shape dependence of a different origin will eventually emerge from two-body correlations at short length scales cow02 ; gao03 . For strong confinement, the local density approximation, implicit in both MGPI and MGPII, will eventually fail. Exploration into those regimes will require theories that differ much more substantially from the GP formulation.

###### Acknowledgements.

The authors are deeply indebted to Dr. Doerte Blume for providing her results of DMC calculations. This work was supported by the National Natural Science Foundation of China No. 19834060 and the Key Project of Knowledge Innovation Program of Chinese Academy of Science No. KJCX2-W7. BG is also supported in part by the US National Science Foundation.


format: text

text:
# Reverse Integration for Computing Stationary Points of Unstable Stiff Systems

## 1 Introduction

We consider the problem of computing a stationary point, $`y_0:f(y_0)=0`$, of the differential equation

$$\frac{dy}{dt}=f(y)$$

when $`y_0`$ is unstable. Here $`f`$ is $`\mathrm{}^n\mathrm{}^n`$.

If the differential equation is exponentially stable to a stationary point in reverse time then an integration backwards in time will approach that stationary point. In a neighborhood of such a point the eigenvalues of

$$J=\frac{f}{y}$$

will be in the positive half plane. In that case we can simply numerically integrate in the reverse direction (negative t), and, if we start sufficiently close to the stationary point and use a small enough step size, we will asymptotically approach it.

We are interested in various situations in which backward integration is not feasible. If all we have is a “black box” time-stepper, appropriate only for integration in the forward direction (this could be a legacy code, or a stochastic kinetic Monte Carlo simulator which only works in forward time) we cannot integrate backwards directly. If there are eigenvalues in both half planes, then the differential equation is unstable in either direction and an accurate integration will have the same properties. We are particularly interested in the case that the instability in the forward direction is “mild” - that is, the eigenvalues in the positive half plane are fairly close to the origin, while the instability in the reverse direction is “severe” - that is, the remaining eigenvalues have quite negative real parts and the corresponding solution components are strongly stable in the forward direction. This type of situation arises in many physical problems where strongly damped components are coupled with mildly growing components - that is, in stiff system where the slow components are unstable. This type of problem could arise, for example, in a singular perturbation context, or when a Differential-Algebraic equation is regularized by converting the algebraic terms into stiffly stable differential components.

In an earlier paper ( we considered projective methods for stiff problems with gaps in their spectra. In the projective method, a numerical solution is computed at a sequence of relatively closely spaced points in time using a conventional integrator with small time steps, and then a “giant” step is taken using polynomial extrapolation from the last few of the points computed by the small steps. This giant step was called the projective step. It was taken forward in time. The small steps stabilized the fast (“strongly stable”) components; the large, projective step had a stability region associated with explicit, large step methods that are stable for slowly damped components. The combined method had a stability region for linear problems that had two components, one that caused damping of the fast components and one that caused damping of the slow.

In the discussion below we will continue to talk about the “rapidly damped” components, referring to those that are damped in the forward time direction, and to unstable components as those that grow in the forward time direction. We will do this even when we discuss integrating in the reverse direction, so as to avoid always having to qualify the terms “stable components” and “unstable components.”

In this note we consider using projective methods in which the projective step is taken backwards in time while the small steps using the conventional integrator remain forward in time. Thus the overall reverse projective integration step consists of small regular forward steps followed by a “giant” reverse step, giving a net integration backwards in time. We will show that the resulting stability regions consist, once again, of two components. The first is due to the conventional small step forward integrator and as before leads to stability of the rapidly damped components. The second, due to the projective step, now leads to stability of the slowly growing, unstable components of the differential equation.

## 2 Analysis

We will analyze the simplest of these methods following the technique used in . The reverse projective integration step we will discuss here consists of $`k+1`$ inner integration steps of size $`h`$ forward from $`t_n`$ to $`t_{n+k+1}`$, followed by one projective step of size $`Mh`$ to arrive at $`t_{k+1M}`$. The projective step takes the form

$$y_{k+1M}^p=My_k(M1)y_{k+1}$$

As in  we assume that the error amplification of an inner integrator step is $`\rho (h\lambda )`$. Assuming that the inner integrator is first order accurate we have

$$\rho (h\lambda )=1+h\lambda +\mathrm{O}(h\lambda )^2$$

and if the inner integrator is the forward Euler method, the last term can be dropped. Then we see that the reverse projective method has an error amplification $`\sigma `$ given by

$$\sigma =\rho ^k(M(M1)\rho ).$$

As before, we define the stability region in the $`\rho `$-plane as the set of $`\rho `$ for which $`\sigma `$ is not outside the unit disk, and plot its boundary by finding the set of $`\rho `$ such that $`|\sigma |=1`$. Figure 1 shows the plots for $`k=2`$ and four different values of $`M`$. Note that, because the process goes forward $`k+1=3`$ steps before going backwards $`M`$ steps, the values actually correspond to net reverse steps of $`h,2h,4h`$, and $`8h`$ respectively. The method is stable inside the regions shown. As $`M`$ gets large, these regions asymptotically tend to a pair of disks. One, centered at $`1+1/M`$ and of radius $`1/M`$, corresponds to the stability region of the forward Euler method because the reverse projective step is the equivalent of a forward Euler method in the reverse direction. The second is centered at the origin and has radius $`M^{1/k}`$. It represents the region where the damping of $`\rho ^k`$ is sufficient to overcome the growth proportional to $`M`$. In other words, the regions are essentially similar to those for forward projective methods, except that the stability region due to the projective step corresponds to $`h\lambda `$ in the positive half plane since that step is taken in the reverse direction.

## 3 Stationary Methods

In the above we have taken $`M>k+1`$ so that the net progress is negative in time. If we choose $`M=k+1`$ the total step length of the reverse projective integration method will be zero; we will call this a stationary projective method. It has the curious stability region shown in Figure 2. It is stable along a section of the real axis that includes $`\rho =1`$ (which is $`h\lambda =0`$). The boundary crosses itself at $`\rho =1`$. The intersection is at right angles because in the neighborhood of $`\rho =1`$ simple algebra shows that $`\sigma `$ is given by

$$\sigma =1k(k+1)/2(\rho 1)^2k^2(k1)/2(\rho 1)^3$$
(1)
so the stability in the neighborhood of $`\rho =1`$ is equivalent to requiring that

$$\mathrm{Real}((\rho 1)^2)0$$

which happens when $`\mathrm{Arg}(\rho )[\pi /4,\pi /4]`$ or $`[3\pi /4,5\pi /4]`$. The form of eq. (1) arises because the first-order accuracy of the method for a net step size of length zero means that

$$\sigma (0\lambda )=1+0\lambda +\mathrm{O}(h\lambda )^2=1+\mathrm{O}(\rho 1)^2$$

It is possible that this procedure will be of value in the context of finding consistent initial conditions for differential algebraic equations (see ). If the algebraic components of a DAE were converted to stiff components (as has been advocated by some), the initialization problem is to find an initial condition that is “on the slow manifold” (i.e., the solution manifold of the original DAE.)

## 4 Comments

The procedure we just outlined will converge to saddle-type points whose linearization is characterized by a gap between strong stable modes and weak unstable ones. An existing forward simulation code will not in general converge to such points (the numerical trajectories will move away forward in time, and “explode away” backward in time). We can therefore think of our procedure as a computational superstructure that transforms a forward simulation code into a contraction mapping capable of converging to such saddle unstable points. The zero-net backward movement process described at the end can be related to consistent initialization algorithms for differential-algebraic equations.

Beyond the computation of saddle-type fixed points, however, if the forward in time dynamics possess such a separation of time scales globally, and an attracting, forward-invariant, low-dimensional slow manifold exists, the procedure becomes an “on manifold” reverse time integration scheme. That is, the regularizing action of the forward integration allows us to follow the on-manifold well-behaved trajectories backward in time. This may then provide a meaningful -and very simple to implement- way of regularizing the reverse, “on the slow manifold”, dynamics of stiff sets of ODEs and even discretizations of disspative PDEs. For example, in contexts where a low-dimensional inertial manifold exists for a dissipative PDE , our superstructure enables a direct simulation of an accurate discretization of the PDE to follow backward trajectories on the inertial manifold without ever having to explicitly derive an inertial form (or approximate inertial form).

In the spirit of the last observation, it is interesting to consider the implications of such a process for a meaningful reverse integration of systems described by microscopic/stochastic simulators. In many practically relevant cases, the coarse-grained behavior of such simulators can be described by the evolution in time of a few “master” moments of microscopically evolving distributions. The remaining, higher moments, become quickly slaved by forward simulation to the slow, master moments. We have recently proposed coarse projective integration schemes that use short bursts of appropriately initialized microscopic/stochastic simulation to estimate the time-derivative of the unavailable coarse equations for the master modes, and “project” these modes forward in time  If the coarse projection is performed backward in time, the procedure will allow us to follow the regularized reverse time behavior of the coarse variables. This is done using the microscopic/stochastic forward-in-time simulator directly, circumventing the necessity of deriving an explicit macroscopic closure. It becomes therefore possible to use a forward-in-time molecular dynamics simulator to extract regularized reverse-time information of coarse system variables. We have already demonstrated the feasibility of this technique in the case of coarse molecular dynamics simulations of a dipeptide folding kinetics in water . Coarse reverse integration allows us to use microscopic simulators to quickly escape free-energy minima, to converge on certain transition states (saddles on the free-energy surface) and, more generally, may enhance our ability to explore the structure of free energy surfaces.

In summary, the technique holds promise towards (a) the computation of unstable, saddle-type fixed points using existing simulators; (b) the regularized, “on manifold” backward-in-time integration of certain dissipative PDEs possessing low-dimensional, exponentially attracting slow manifolds; and (c) the use of microscopic/stochastic simulators to track coarse-grained behavior backward in time, enhancing the ability to escape free energy minima and to locate saddle-type coarse-grained “transition states”.

Acknolwedgements. This work was partially supported by AFOSR (Dynamics and Control, Dr. B. King) and an NSF ITR grant. Discussions with Profs. P. G. Kevrekidis (UMass), Ju Li (OSU) and Dr. G. Hummer (NIH) are gratefully acknowledged.


text:
# Integration of twisted Poisson structures

## 1. Introduction

Poisson manifolds may be regarded as the infinitesimal form of symplectic groupoids , i.e., Lie groupoids endowed with a multiplicative symplectic form. Up to singularities, Poisson manifolds may be integrated to symplectic groupoids as described in  (conditions under which integration with no singularities is possible are given in ). In this paper we generalize this result to the case when the two structures (of symplectic groupoid and of Poisson manifold) are twisted by a closed $`3`$-form.

Let $`M`$ be a smooth manifold. A pair $`(\pi ,\varphi )`$, where $`\pi `$ is a bivector field and $`\varphi `$ is a closed $`3`$-form, is called a twisted Poisson structure if it satisfies the equation

(1.1)
$$[\pi ,\pi ]=\frac{1}{2}^3\pi ^\mathrm{\#}\varphi ,$$

where $`[,]`$ denotes the Schouten–Nijenhuis bracket and $`\pi ^\mathrm{\#}`$ is the vector bundle homomorphism $`T^{}MTM`$ induced by $`\pi `$ (viz., $`\pi ^\mathrm{\#}(x)(\sigma ):=\pi (x)(\sigma ,)`$, with $`xM`$, $`\sigma T_x^{}M`$). According to , one also says that $`\pi `$ is a $`\varphi `$-Poisson tensor. In the case $`\varphi =0`$ one recovers the usual notions of Poisson tensor and Poisson manifold. Twisted Poisson structures have been extensively studied in the physics literature, e.g., .

As explained in , a twisted Poisson structure induces a Lie algebroid structure on $`T^{}M`$ with anchor map $`\pi ^\mathrm{\#}`$ and Lie bracket of sections $`\sigma `$ and $`\tau `$ defined by

(1.2)
$$[\sigma ,\tau ]:=\mathrm{L}_{\pi ^\mathrm{\#}\sigma }\tau \mathrm{L}_{\pi ^\mathrm{\#}\tau }\sigma \mathrm{d}\pi (\sigma ,\tau )+\varphi (\pi ^\mathrm{\#}\sigma ,\pi ^\mathrm{\#}\tau ,).$$

In particular, $`f,gC^{\mathrm{}}(M)`$ we have:

(1.3) $`[\mathrm{d}f,\mathrm{d}g]`$ $`=\mathrm{d}\{f,g\}+\varphi (X_f,X_g,)`$
and
(1.4) $`[X_f,X_g]`$ $`=X_{\{f,g\}}+\pi ^\mathrm{\#}(\varphi (X_f,X_g,)),`$
with $`\{f,g\}=\pi (\mathrm{d}f,\mathrm{d}g)`$ and $`X_f=\pi ^\mathrm{\#}\mathrm{d}f`$.

We will denote this Lie algebroid by $`T^{}M_{(\pi ,\varphi )}`$. Sections of its exterior algebra are ordinary differential forms. One may define a derivation $`\delta `$ deforming the de Rham differential by $`\varphi `$; viz., we define a graded derivation $`\delta :\mathrm{\Omega }^{}(M)\mathrm{\Omega }^{+1}(M)`$ by setting $`\delta f=\mathrm{d}f`$ if $`fC^{\mathrm{}}(M)`$ and

$$\delta \sigma =\mathrm{d}\sigma \iota _{\pi ^\mathrm{\#}\sigma }\varphi ,$$

if $`\sigma \mathrm{\Omega }^1(M)`$. It turns out that

$$\delta [\sigma ,\tau ]=[\delta \sigma ,\tau ]+[\sigma ,\delta \tau ],\sigma ,\tau \mathrm{\Omega }^1(M),$$

and that $`\delta ^2=[\varphi ,]`$ (where we have extended the Lie bracket to the whole of $`\mathrm{\Omega }^{}(M)`$ as a biderivation). So $`(T^{}M_{(\pi ,\varphi )},\delta )`$ constitutes an example of a quasi Lie bialgebroid , a generalization of Drinfeld’s quasi Lie bialgebras .

If $`T^{}M_{(\pi ,\varphi )}`$ may be integrated to a Lie groupoid $`(GM,\alpha ,\beta )`$ (i.e., if it exists a Lie groupoid $`G`$ whose Lie algebroid is $`T^{}M_{(\pi ,\varphi )}`$), the differential $`\delta `$ induces extra structure on $`G`$. Namely, denoting by $`\alpha `$ and $`\beta `$ the source and target maps of $`G`$, then $`G`$ may be endowed with a non-degenerate, multiplicative $`2`$-form $`\omega `$ that satisfies

$$\mathrm{d}\omega =\alpha ^{}\varphi \beta ^{}\varphi .$$

In other words, $`(\omega ,\varphi )`$ is a $`3`$-cocycle for the double complex $`\mathrm{\Omega }^{}(G^{()})`$, where $`G^{(0)}=M`$, $`G^{(1)}=G`$ and elements of $`G^{(k)}`$ are $`k`$-tuples of elements of $`G`$ that may be multiplied (in the given order). One differential is de Rham and the other is the groupoid-complex differential. Observe that in the true Poisson case (i.e., $`\varphi =0`$), $`\omega `$ is closed, so $`G`$ is an ordinary symplectic groupoid. In the general case, $`G`$ is called a non-degenerate twisted symplectic groupoid, and the non-degenerate $`2`$-form $`\omega `$ is said to be relatively $`\varphi `$-closed. The main theorem of the paper (conjectured in ) is

###### Theorem.

There is a bijection between integrable twisted Poisson structures and source-simply connected non-degenerate twisted symplectic groupoids.

Here “integrable twisted Poisson structure” means that the associated Lie algebroid is integrable.

In Sect. 2 we give an introduction to non-degenerate twisted symplectic groupoids and prove that they induce twisted Poisson structures on the base manifolds (Theorem 2.6 on page 2.6).

In Sect. 5 we prove the Theorem, though in a more general setting. In fact, as shown in the generaliztion  (see also ) of the construction given in ), to any Lie algebroid $`A`$ one can associate a topological source-simply-connected groupoid $`G(A)`$, which is the Lie groupoid integrating $`A`$ whenever $`A`$ is integrable. The topological groupoid $`G(A)`$ is defined as the leaf space of a smooth foliation, as we recall in Sect. 3; so it makes sense to define on it a notion of smooth functions and forms. In the case when $`A`$ is $`T^{}M_{(\pi ,\varphi )}`$, we prove that $`G(A)`$ may always be endowed with a non-degenerate, multiplicative, relatively $`\varphi `$-closed $`2`$-form $`\omega `$. The construction is a modification, described in Sect. 4, of the method developed in , where the true Poisson case (i.e., $`\varphi =0`$) was dealt with.

As a final remark, we mention that general multiplicative $`2`$-forms, their infinitesimal counterparts and their integrations are being treated in .

###### Acknowledgment.

We thank Jim Stasheff and Alan Weinstein for useful discussions and comments. We acknowledge Zurich University (P. X.) and Penn State University (A. S. C.) for their kind hospitality during the preparation of the work.

## 2. Non-degenerate twisted symplectic groupoids

###### Definition 2.1.

A non-degenerate twisted symplectic groupoid is a Lie groupoid $`(GM,\alpha ,\beta )`$ equipped with a non-degenerate $`2`$-form $`\omega \mathrm{\Omega }^2(G)`$ and a $`3`$-form $`\varphi \mathrm{\Omega }^3(M)`$ such that

1. $`\mathrm{d}\varphi =0`$;
2. $`\mathrm{d}\omega =\alpha ^{}\varphi \beta ^{}\varphi `$;
3. $`\omega `$ is multiplicative, i.e., the $`2`$-form $`(\omega ,\omega ,\omega )`$ vanishes when being restricted to the graph of the groupoid multiplication $`\mathrm{\Lambda }G\times G\times G`$.

Let $`\pi _G`$ denote the bivector field on $`G`$ corresponding to $`\omega `$. Then $`(\pi _G,\mathrm{\Omega })`$, where $`\mathrm{\Omega }=\alpha ^{}\varphi \beta ^{}\varphi `$, defines a twisted Poisson structure on $`G`$ in the sense of .

For any $`\xi \mathrm{\Gamma }(A)`$, by $`\stackrel{}{\xi }`$ and $`\stackrel{}{\xi }`$ we denote its corresponding right and left invariant vector fields on the groupoid $`G`$, respectively. The following properties can be easily verified.

###### Proposition 2.2.

1. $`ϵ^{}\omega =0`$, where $`ϵ:MG`$ is the natural embedding;
2. $`i^{}\omega =\omega `$, where $`i:GG`$ is the groupoid inversion;
3. for any $`\xi ,\eta \mathrm{\Gamma }(A)`$, $`\omega (\stackrel{}{\xi },\stackrel{}{\eta })`$ is a right invariant function on $`G`$, and $`\omega (\stackrel{}{\xi },\stackrel{}{\eta })`$ is a left invariant function on $`G`$;
4. $`\omega (\stackrel{}{\xi },\stackrel{}{\eta })=0`$;
5. $`\omega (\stackrel{}{\xi },\stackrel{}{\eta })(x)=\omega (\stackrel{}{\xi },\stackrel{}{\eta })(x^1)`$.

###### Proof.

The proof is standard, and essentially follows from the multiplicativity of $`\omega `$.

1. For any $`\delta _m^{},\delta _m^{\prime \prime }T_mM`$, since $`(\delta _m^{},\delta _m^{},\delta _m^{}),(\delta _m^{\prime \prime },\delta _m^{\prime \prime },\delta _m^{\prime \prime })T\mathrm{\Lambda }`$, it follows that $`\omega (\delta _m^{},\delta _m^{\prime \prime })=0`$.
2. $`xG`$ and $`\delta _x^{},\delta _x^{\prime \prime }T_xG`$, it is clear that $`(\delta _x^{},i_{}\delta _x^{},\alpha _{}\delta _x^{}),(\delta _x^{\prime \prime },i_{}\delta _x^{\prime \prime },\alpha _{}\delta _x^{\prime \prime })T\mathrm{\Lambda }`$. Thus, by (1), we have

$$\omega (\delta _x^{},\delta _x^{\prime \prime })+\omega (i_{}\delta _x^{},i_{}\delta _x^{\prime \prime })=0,$$

and (2) follows.
3. For any $`\xi ,\eta \mathrm{\Gamma }(A)`$, $`(\stackrel{}{\xi }(x),0_y,\stackrel{}{\xi }(xy)),(\stackrel{}{\eta }(x),0_y,\stackrel{}{\eta }(xy))T\mathrm{\Lambda }`$. Thus

$$\omega (\stackrel{}{\xi }(x),\stackrel{}{\eta }(x))\omega (\stackrel{}{\xi }(xy),\stackrel{}{\eta }(xy))=0.$$

Hence $`\omega (\stackrel{}{\xi },\stackrel{}{\eta })`$ is a right invariant function on $`G`$. Similarly, $`\omega (\stackrel{}{\xi },\stackrel{}{\eta })`$ is a left invariant function on $`G`$.
4. By considering the vectors $`(\stackrel{}{\xi }(x),0_{\beta (x)},\stackrel{}{\xi }(x))`$ and $`(0_x,\stackrel{}{\eta }(\beta (x)),\stackrel{}{\eta }(x))T\mathrm{\Lambda }`$, we obtain $`\omega (\stackrel{}{\xi }(x),\stackrel{}{\eta }(x))=0`$.
5. follows from (2) and the fact that $`i_{}\stackrel{}{\xi }=\stackrel{}{\xi }`$.

Define a section $`\gamma \mathrm{\Gamma }(^2A^{})`$ and a bundle map: $`\lambda :AT^{}M`$ by

(2.1)
$$\omega (\stackrel{}{\xi },\stackrel{}{\eta })=\alpha ^{}\gamma (\xi ,\eta ),\xi ,\eta \mathrm{\Gamma }(A),$$

and

(2.2)
$$\lambda (\xi ),v=\omega (\stackrel{}{\xi }(m),v),\xi A|_m,vT_mM$$

###### Lemma 2.3.

1. $`\omega (\stackrel{}{\xi },\stackrel{}{\eta })=\beta ^{}\gamma (\xi ,\eta ),\xi ,\eta \mathrm{\Gamma }(A)`$;
2. for all $`\xi ,\eta \mathrm{\Gamma }(A)`$,

(2.3)
$$\gamma (\xi ,\eta )=\rho (\xi ),\lambda (\eta );$$
3. $`\lambda :AT^{}M`$ is a vector bundle isomorphism.

###### Proof.

1. follows from Proposition 2.2 (5).
2. We have

$$\omega (\stackrel{}{\xi },\stackrel{}{\eta })=\omega (\stackrel{}{\xi }\stackrel{}{\xi },\stackrel{}{\eta })=\omega (\stackrel{}{\eta },\rho (\xi ))=\rho (\xi ),\lambda (\eta )$$
3. Assume that $`\lambda (\xi )=0`$. That is, $`\omega (\stackrel{}{\xi }(m),v)=0,vT_mM`$, which implies that $`\stackrel{}{\xi }(m)\text{ }\text{ }\omega =0`$ by Proposition 2.2 (4). Hence $`\xi =0`$ since $`\omega `$ is non-degenerate. This means that $`\lambda `$ is injective. On the other hand, assume that $`v(\lambda (A|_m))^{}`$. Then $`\omega (\stackrel{}{\xi }(m),v)=0,\xi A|_m`$. Thus $`v\text{ }\text{ }\omega =0`$ using Proposition 2.2 (1), which implies that $`v=0`$. Therefore $`\lambda `$ is surjective.

###### Lemma 2.4.

For any $`fC^{\mathrm{}}(M)`$,

(2.4)
$$\stackrel{}{\lambda ^1(\mathrm{d}f)}=X_{\alpha ^{}f};\stackrel{}{\lambda ^1(\mathrm{d}f)}=X_{\beta ^{}f}.$$

###### Proof.

First, one shows that $`X_{\alpha ^{}f}`$ is a right invariant vector field on $`G`$ and $`X_{\beta ^{}f}`$ is a left invariant vector field. This can be shown using the same argument as in the case of symplectic groupoids . Namely the multiplicativity of $`\omega `$ together with dimension counting implies that the graph $`\mathrm{\Lambda }`$ is coisotropic with respect to $`(\pi _G,\pi _G,\pi _G)`$. The later implies that $`X_{\alpha ^{}f}`$ is a right invariant vector field on $`G`$ and $`X_{\beta ^{}f}`$ is a left invariant vector field.

Second, for any $`vT_mM`$, we have

$$\omega (X_{\alpha ^{}f}(m),v)=\alpha ^{}\mathrm{d}f(m),v=\mathrm{d}f(m),\alpha _{}v=\mathrm{d}f(m),v$$

It thus follows that $`\lambda (X_{\alpha ^{}f})=\mathrm{d}f`$, or $`\stackrel{}{\lambda ^1(\mathrm{d}f)}=X_{\alpha ^{}f}`$. The other equation can be proved similarly. ∎

By pulling back the $`2`$-form $`\gamma \mathrm{\Gamma }(^2A^{})`$ via $`\lambda ^1`$, one obtains a bivector field $`\pi \mathrm{\Gamma }(^2TM)`$. We introduce a bracket and Hamiltonian vector fields by the usual definitions; i.e., $`\{f,g\}=\pi (\mathrm{d}f,\mathrm{d}g)`$ and $`X_f=\pi ^\mathrm{\#}(\mathrm{d}f)`$.

###### Corollary 2.5.

(2.5)
$$\alpha _{}\pi _G=\pi ;\beta _{}\pi _G=\pi ;$$

or equivalently

(2.6)
$$\alpha _{}X_{\alpha ^{}f}=X_f;\beta _{}X_{\beta ^{}f}=X_f,fC^{\mathrm{}}(M).$$

###### Proof.

For any $`f,gC^{\mathrm{}}(M)`$,

$$\begin{array}{c}\{\alpha ^{}f,\alpha ^{}g\}=\omega (X_{\alpha ^{}f},X_{\alpha ^{}g})=\omega (\stackrel{}{\lambda ^1(\mathrm{d}f)},\stackrel{}{\lambda ^1(\mathrm{d}g)})=\hfill \\ \hfill =\alpha ^{}(\pi (\mathrm{d}f,\mathrm{d}g))=\alpha ^{}\{f,g\}.\end{array}$$

Similarly, we have $`\{\beta ^{}f,\beta ^{}g\}=\beta ^{}\{f,g\}`$. ∎

We are now ready to prove the main result of the section.

###### Theorem 2.6.

1. $`\pi `$ is a $`\varphi `$-Poisson tensor in the sense of —i.e., it satisfies (1.1).
2. The bundle map $`\lambda :AT^{}M`$ establishes a Lie algebroid isomorphism, where the Lie algebroid on $`T^{}M`$ is induced by the twisted Poisson tensor $`\pi `$ as given by Eq. (1.2).

###### Proof.

Let $`\mathrm{\Omega }=\alpha ^{}\varphi \beta ^{}\varphi `$. Thus $`f,gC^{\mathrm{}}(M)`$

$`(X_{\alpha ^{}f}X_{\alpha ^{}g})\text{ }\text{ }\mathrm{\Omega }`$ $`=`$ $`(X_{\alpha ^{}f}X_{\alpha ^{}g})\text{ }\text{ }\alpha ^{}\varphi `$
$`=`$ $`\alpha ^{}[(\alpha _{}X_{\alpha ^{}f}\alpha _{}X_{\alpha ^{}g})\text{ }\text{ }\varphi ]`$
$`=`$ $`\alpha ^{}[X_fX_g\text{ }\text{ }\varphi ].`$
Thus by Eq. (1.4)

$`[X_{\alpha ^{}f},X_{\alpha ^{}g}]X_{\{\alpha ^{}f,\alpha ^{}g\}}=\pi _G^\mathrm{\#}(\mathrm{\Omega }(X_{\alpha ^{}f},X_{\alpha ^{}g},)`$
$`=\pi _G^\mathrm{\#}(\alpha ^{}\varphi (X_f,X_g,))`$
Thus it follows that

$$\lambda [X_{\alpha ^{}f},X_{\alpha ^{}g}]=\mathrm{d}\{f,g\}+\varphi (X_f,X_g,).$$

Note that $`\lambda `$ intertwines the anchors: $`\pi ^\mathrm{\#}\lambda =\rho `$, according to Eq. (2.3). Therefore, using Lie algebroid properties, one shows that the push forward Lie algebroid on $`T^{}M`$ via $`\lambda `$ is given by Eq. (1.2). This forces, by the Jacobi identity, $`\pi `$ to be $`\varphi `$-Poisson, and $`\lambda `$ is a Lie algebroid isomorphism between $`A`$ and $`(T^{}M)_{\pi ,\varphi }`$. ∎

## 3. Integration of $`T^{}M_{(\pi ,\varphi )}`$

We briefly describe the integration procedure for Lie algebroids of , adapted to the case of $`T^{}M_{(\pi ,\varphi )}`$. First one defines the manifold $`P(T^{}M_{(\pi ,\varphi )})`$ of $`C^1`$-Lie algebroid morphisms $`TIT^{}M_{(\pi ,\varphi )}`$, where $`I`$ is the interval $`[0,1]`$ and $`TI`$ is given its canonical Lie algebroid structure. An element of $`PT^{}M_{(\pi ,\varphi )}`$ consists of a $`C^2`$-path $`X:IM`$ together with a section $`\eta `$ of $`T^{}IX^{}T^{}M`$ satisfying

$$\mathrm{d}X=\pi ^\mathrm{\#}(X)\eta .$$

On this manifold one may consider as equivalent two elements which are related by a Lie algebroid morphism $`T(I\times I)T^{}M_{(\pi ,\varphi )}`$ that fixes the endpoints. The quotient space $`G(T^{}M_{(\pi ,\varphi )})`$ may be given a groupoid structure. For our purposes it is however better to use a different description of $`G(T^{}M_{(\pi ,\varphi )})`$, i.e., as the leaf space of a foliation. Namely, let $`P_0\mathrm{\Gamma }(T^{}M_{(\pi ,\varphi )})`$ be the space of $`C^2`$-paths in the Lie algebra of sections of $`T^{}M_{(\pi ,\varphi )}`$ with endpoints at zero. We give this space the structure of a Lie algebra by the pointwise Lie bracket. One may then define an infinitesimal action of this Lie algebra on $`P(T^{}M_{(\pi ,\varphi )})`$. To describe it, we prefer to introduce local coordinates $`\{x^i\}`$ on $`M`$ (alternatively, one may use a torsion-free connection). Since $`\{\mathrm{d}x^i\}`$ is a local basis of sections of $`T^{}M_{(\pi ,\varphi )}`$, we may define structure functions $`f`$ by

$$[\mathrm{d}x^i,\mathrm{d}x^j]=f_k^{ij}\mathrm{d}x^k,$$

where a sum over repeated indices is understood. If we write locally $`\pi =\pi ^{ij}_i_j`$ and $`\varphi =\varphi _{ijk}\mathrm{d}x^i\mathrm{d}x^j\mathrm{d}x^k`$, we may compute:

$$f_k^{ij}=_k\pi ^{ij}+\pi ^{mi}\pi ^{nj}\varphi _{mnk}.$$

The action is then as follows. To $`BP_0\mathrm{\Gamma }(T^{}M_{(\pi ,\varphi )})`$ we associate a vector field $`\xi _B`$ on $`P(T^{}M_{(\pi ,\varphi )})`$. We can always write $`\xi _B=\xi _B^h+\xi _B^v`$ with $`\xi _B^h(X,\eta )\mathrm{\Gamma }(I,X^{}TM)`$ and $`\xi _B^v(X,\eta )\mathrm{\Gamma }(I,T^{}IX^{}T^{}M)`$. We set then

(3.1a) $`(\xi _B^h(X,\eta ))^i`$ $`=\pi ^{ij}(X)(B_X)_j,`$
(3.1b) $`(\xi _B^v(X,\eta ))_i`$ $`=\mathrm{d}(B_X)_if_i^{rs}(X)\eta _r(B_X)_s,`$
where $`B_X`$ is the section of $`X^{}T^{}M`$ defined by $`B_X(t)=B(t)(X(t))`$.

Thus, the infinitesimal action of $`P_0\mathrm{\Gamma }(T^{}M_{(\pi ,\varphi )})`$ defines a foliation on $`P(T^{}M_{(\pi ,\varphi )})`$ and $`G(T^{}M_{(\pi ,\varphi )})`$ is its quotient space. Let us briefly recall its groupoid structure. The target map $`\alpha `$ associates to a class of morphisms $`(X,\eta )`$ the value of $`X`$ at $`0`$, while the source map $`\beta `$ associates to it the value of $`X`$ at $`1`$ (observe that the infinitesimal action preserves the endpoints of $`X`$). The identity section associates to a point $`m`$ in $`M`$ the class $`ϵ(m)`$ of the constant path at $`m`$ with $`\eta =0`$. The product is obtained by joining the base paths and restricting the fiber maps consequently (the product is more precisely defined on smooth representatives such that $`\eta `$ vanishes with its derivatives at the endpoints).

## 4. Quasi-symplectic reduction

In this section we describe how to obtain $`G(T^{}M_{(\pi ,\varphi )})`$ by some sort of symplectic reduction, though our replacement for a symplectic form will be a non-degenerate but not necessarily closed $`2`$-form.

Let $`T^{}PM`$ denote the manifold of $`C^1`$-bundle maps $`TIT^{}M`$ (over $`C^2`$-maps). This space is morally a cotangent bundle and as such it has a canonical symplectic structure $`\mathrm{\Omega }_0`$. Explicitly, a point in $`T^{}PM`$ is a pair $`(X,\eta )`$ where $`X`$ is a $`C^2`$-path $`IM`$ and $`\eta `$ is a $`C^1`$-section of $`T^{}IX^{}T^{}M`$. The tangent space at $`(X,\eta )`$ is the direct sum of $`T_{(X,\eta )}^hT^{}PM=\mathrm{\Gamma }(I,X^{}TM)`$ and $`T_{(X,\eta )}^vT^{}PM=\mathrm{\Gamma }(I,T^{}IX^{}T^{}M)`$. Using this splitting, we write

(4.1)
$$\mathrm{\Omega }_0(X,\eta )(\xi _1e_1,\xi _2e_2)=_Ie_1,\xi _2e_2,\xi _1,$$

where $`,`$ denotes the canonical pairing between tangent and cotangent fibers of $`M`$.

Using the $`3`$-form $`\varphi `$ on $`M`$ we may also define a second $`2`$-form on $`T^{}PM`$:

(4.2)
$$\mathrm{\Omega }_1(X,\eta )(\xi _1e_1,\xi _2e_2)=\frac{1}{2}_I\varphi (X)(\pi ^\mathrm{\#}(X)\eta ,\xi _1,\xi _2).$$

The $`2`$-form $`\mathrm{\Omega }=\mathrm{\Omega }_0+\mathrm{\Omega }_1`$ is still non-degenerate but no longer closed.

The manifold $`P(T^{}M_{(\pi ,\varphi )})`$ introduced in the previous Section may be regarded as a submanifold of $`T^{}PM`$. If we introduce “momentum maps” $`H:T^{}PMP_0\mathrm{\Gamma }(T^{}M_{(\pi ,\varphi )})^{}`$ by

$$H_B(X,\eta )=_IB_X,\mathrm{d}X\pi ^\mathrm{\#}(X)\eta ,$$

then $`P(T^{}M_{(\pi ,\varphi )})`$ is $`H^1(0)`$. One may check that $`\mathrm{d}H_B`$ lies in the image of $`\mathrm{\Omega }`$ for any $`BP_0\mathrm{\Gamma }(T^{}M_{(\pi ,\varphi )})`$; so, since $`\mathrm{\Omega }`$ is non-degenerate, one may define a map $`B\widehat{\xi }_B`$ that associates a vector field $`\widehat{\xi }_B`$ on $`T^{}PM`$ to $`B`$ by

(4.3)
$$\iota _{\widehat{\xi }_B}\mathrm{\Omega }=\mathrm{d}H_B.$$

One may easily check that the restriction of $`\widehat{\xi }_B`$ to $`P(T^{}M_{(\pi ,\varphi )})`$ is tangent to it. More to the point, one may check that the vector field on $`P(T^{}M_{(\pi ,\varphi )})`$ so obtained is precisely the $`\xi _B`$ of (3.1) which defines the infinitesimal action of $`P_0\mathrm{\Gamma }(T^{}M_{(\pi ,\varphi )})`$ on $`P(T^{}M_{(\pi ,\varphi )})`$.

## 5. Proof of the Theorem

In the setting of the previous Section, we want to prove that the restriction $`\underset{¯}{\mathrm{\Omega }}`$ of $`\mathrm{\Omega }`$ to $`P(T^{}M_{(\pi ,\varphi )})`$ is basic w.r.t. to the projection $`p:P(T^{}M_{(\pi ,\varphi )})G(T^{}M_{(\pi ,\varphi )})`$, viz., $`\underset{¯}{\mathrm{\Omega }}=p^{}\omega `$; moreover, we want to prove that $`\omega `$ satisfies all the required conditions.

Observe that $`\underset{¯}{\mathrm{\Omega }}`$ is automatically horizontal by (4.3). On the other hand, unlike the usual symplectic case, it is not clear that $`\underset{¯}{\mathrm{\Omega }}`$ is also invariant; in fact, at first, we may only see that $`\mathrm{L}_{\xi _B}\underset{¯}{\mathrm{\Omega }}=\iota _{\xi _B}\mathrm{d}\underset{¯}{\mathrm{\Omega }}=\iota _{\xi _B}\mathrm{d}\underset{¯}{\mathrm{\Omega }}_1`$, where $`\underset{¯}{\mathrm{\Omega }}_1`$ denotes the restriction of $`\mathrm{\Omega }_1`$ to $`P(T^{}M_{(\pi ,\varphi )})`$. To proceed, we must understand $`\underset{¯}{\mathrm{\Omega }}_1`$ better.

Let $`PM`$ be the manifold of $`C^2`$-paths in $`M`$. Let $`\mathrm{ev}:I\times PMM`$ be the evaluation map and $`\mathrm{pr}:I\times PMPM`$ the projection to the second factor. Define $`\mathrm{\Phi }=\mathrm{pr}_{}\mathrm{ev}^{}\varphi \mathrm{\Omega }^2(PM)`$, where $`\mathrm{pr}_{}`$ denotes integration along the fiber. If we finally denote by $`q:P(T^{}M_{(\pi ,\varphi )})PM`$ the map that retains only the base map of the Lie algebroid morphism, we realize immediately that

$$\underset{¯}{\mathrm{\Omega }}_1=q^{}\mathrm{\Phi }.$$

By the generalized Stokes’ Theorem and the fact that $`\varphi `$ is closed, we obtain $`\mathrm{d}\mathrm{\Phi }=\underset{¯}{\alpha }^{}\varphi \underset{¯}{\beta }^{}\varphi `$, where $`\underset{¯}{\alpha }`$ and $`\underset{¯}{\beta }`$ are the maps $`PMM`$ that assign to a path its values at $`0`$ and at $`1`$, respectively. Thus,

$$\mathrm{d}\underset{¯}{\mathrm{\Omega }}=q^{}(\underset{¯}{\alpha }^{}\varphi \underset{¯}{\beta }^{}\varphi ).$$

Since the vector field $`\xi _B`$ does not move the endpoints, we conclude that $`\iota _{\xi _B}\mathrm{d}\underset{¯}{\mathrm{\Omega }}=0`$, viz., that $`\underset{¯}{\mathrm{\Omega }}`$ is invariant as well. We write then $`\underset{¯}{\mathrm{\Omega }}=p^{}\omega `$ as at the beginning of the Section. The $`2`$-form $`\omega `$ on $`G(T^{}M_{(\pi ,\varphi )})`$ is clearly multiplicative since the product is defined by joining the paths and $`\mathrm{\Omega }`$ is defined as an integral. Moreover, recalling the definition of the source and target map $`\beta `$ and $`\alpha `$, we observe that $`\underset{¯}{\alpha }q=\alpha p`$ and $`\underset{¯}{\beta }q=\beta p`$. So we may write the equation above as

$$\mathrm{d}\underset{¯}{\mathrm{\Omega }}=p^{}(\alpha ^{}\varphi \beta ^{}\varphi ).$$

Since $`\mathrm{d}\underset{¯}{\mathrm{\Omega }}=p^{}\mathrm{d}\omega `$ and $`p`$ is a surjection, this shows that $`\omega `$ is relatively $`\varphi `$-closed.

Finally, we need to prove that the $`2`$-form $`\omega `$ is non-degenerate. It is clear from the construction that $`\omega `$ is non-degenerate along the identity $`M`$. The claim thus follows from the following:

###### Lemma 5.1.

A multiplicative $`2`$-form $`\omega \mathrm{\Omega }^2(G)`$ on a Lie groupoid $`GM`$ is non-degenerate if and only if it is non-degenerate along the identity $`M`$.

###### Proof.

First of all, note that for any $`\delta _xT_xG`$, and $`\xi \mathrm{\Gamma }(A)`$, we have

(5.1) $`\omega (\stackrel{}{\xi }(x),\delta _x)=\omega (\stackrel{}{\xi }(v),\beta _{}\delta _x)`$
(5.2) $`\omega (\stackrel{}{\xi }(x),\delta _x)=\omega (\stackrel{}{\xi }(u),\alpha _{}\delta _x),`$
where $`u=\alpha (x)`$ and $`v=\beta (x)`$. Eq. (5.1), for instance, follows from the fact that both $`(\delta _x,\delta _x,\beta _{}\delta _x)`$, and $`(0,\stackrel{}{\xi }(x),\stackrel{}{\xi }(v))`$ are tangent to the graph of the groupoid multiplication $`\mathrm{\Lambda }G\times G\times G`$. Eq. (5.2) can be proved similarly. Now assume that $`\delta _x\mathrm{ker}\omega _x`$. It follows from Eq. (5.1) that $`\beta _{}\delta _x\mathrm{ker}\omega _v`$ since $`M`$ is isotropic with respect to $`\omega `$. Therefore $`\beta _{}\delta _x=0`$ by assumption. Hence $`\delta _x=\stackrel{}{\eta }(x)`$. On the other hand, according to Eq. (5.2), one has $`\omega (\stackrel{}{\eta }(u),T_uM)=0`$ since $`\alpha `$ is a submersion. This implies that $`\stackrel{}{\eta }(u)\mathrm{ker}\omega _u`$. Therefore $`\stackrel{}{\eta }(u)=0`$ by assumption. This implies that $`\delta _x=\stackrel{}{\eta }(x)=0`$. This concludes the proof. ∎

We need now to prove that the correspondence between $`\varphi `$-twisted Poisson structures and twisted symplectic groupoids is a bijection. The proof is divided into two steps.

### Step 1

By construction (see ) the Lie algebroid of $`G(T^{}M_{(\pi ,\varphi )})`$ is $`T^{}M_{(\pi ,\varphi )}`$. As discussed in Sect. 2, the relatively $`\varphi `$-closed, multiplicative, non-degenerate $`2`$-form $`\omega `$ determines an automorphism $`\lambda `$ of $`T^{}M`$ and a bivector field $`\gamma `$ on $`M`$ as in Eqs. (2.1) and (2.2). We have to show that $`\lambda `$ is the identity and that $`\gamma =\pi `$. First of all we observe that it is enough to consider (2.1) at the unit element $`ϵ(m)G(T^{}M_{(\pi ,\varphi )})`$ corresponding to $`mM`$:

$$\omega (ϵ(m))(\stackrel{}{\xi _1}(ϵ(m)),\stackrel{}{\xi _2}(ϵ(m)))=\gamma (m)(\xi _1,\xi _2),\xi _1,\xi _2A|_m.$$

By construction $`ϵ(m)`$ is the equivalence class of the path $`X(t)=m`$, $`\eta (t)=0`$, $`tI=[0,1]`$. The vector field $`\stackrel{}{\xi _i}`$, $`i=1,2`$, evaluated at $`ϵ(m)`$ is the projection to $`T_{ϵ(m)}G(T^{}M_{(\pi ,\varphi )})`$ of the vector $`\widehat{\xi }_iT_{(m,0)}P(T^{}M_{(\pi ,\varphi )})`$ defined by $`\widehat{\xi }_i(t)=(\pi ^\mathrm{\#}(m)\xi _it,\xi _i\mathrm{d}t)`$. Observing then that for $`\eta =0`$ the $`2`$-form $`\mathrm{\Omega }_1`$ of Eq. (4.2) vanishes, we get, also using (4.1),

$$\begin{array}{c}\omega (ϵ(m))(\stackrel{}{\xi _1}(ϵ(m)),\stackrel{}{\xi _2}(ϵ(m)))=\mathrm{\Omega }_0(m,0)(\widehat{\xi }_1,\widehat{\xi }_2)=\hfill \\ \hfill =2_0^1\pi (m)(\xi _1,\xi _2)tdt=\pi (m)(\xi _1,\xi _2),\end{array}$$

which shows $`\gamma =\pi `$. As for (2.2), observe that $`\omega (ϵ(m))(\stackrel{}{\xi _1}(ϵ(m)),v)`$ is just $`\mathrm{\Omega }_0(m,0)(\widehat{\xi }_1,\widehat{v})`$ with $`\widehat{v}(t)=(v,0)`$. As a consequence,

$$\omega (ϵ(m))(\stackrel{}{\xi _1}(ϵ(m)),v)=_0^1\xi _1,vdt=\xi _1,v,$$

which shows that $`\lambda `$ is the identity.

### Step 2

Assume that $`(GM,\omega +\varphi )`$ is an $`\alpha `$-simply connected non-degenerate twisted symplectic groupoid. Let $`\pi `$ be its induced $`\varphi `$-twisted Poisson structure on $`M`$. Then the above integration process integrates the Lie algebroid $`T^{}M_{(\pi ,\varphi )}`$ into a Lie groupoid, which is known to be isomorphic to $`GM`$, and a multiplicative $`2`$-form $`\omega ^{}`$ on that groupoid. By identifying this groupoid with $`GM`$, therefore one may think $`\omega ^{}`$ as a multiplicative $`2`$-form on $`G`$. One needs to show that $`\omega ^{}=\omega `$. By Step 1, we conclude that $`\omega ^{}`$ and $`\omega `$ must coincide along the identity space $`M`$. Let $`\stackrel{~}{\omega }=\omega \omega ^{}`$. Then $`\stackrel{~}{\omega }`$ is a multiplicative closed $`2`$-form on $`G`$ and $`\stackrel{~}{\omega }|_M=0`$. Given any $`\xi \mathrm{\Gamma }(A)`$, it is easy to see that $`(\stackrel{}{\xi }(\alpha (x)),0,\stackrel{}{\xi }(x))`$ is tangent to the graph $`\mathrm{\Lambda }`$ of groupoid multiplication. On the other hand, for any $`\delta _xT_xG`$, it is also clear that $`(\alpha _{}\delta _x,\delta _x,\delta _x)T\mathrm{\Lambda }`$. It thus follows that

$$\stackrel{~}{\omega }(\stackrel{}{\xi }(\alpha (x)),\alpha _{}\delta _x)\stackrel{~}{\omega }(\stackrel{}{\xi }(x),\delta _x)=0.$$

Therefore we have $`\stackrel{}{\xi }\text{ }\text{ }\stackrel{~}{\omega }=0`$. Thus

$$L_\stackrel{}{\xi }\stackrel{~}{\omega }=(di_\stackrel{}{\xi }+i_\stackrel{}{\xi }d)\stackrel{~}{\omega }=0,$$

which implies that $`\stackrel{~}{\omega }=0`$ since any point in $`G`$ can be reached by a product of (local) bisections generated by $`\stackrel{}{\xi }`$. This concludes the proof.