text:
# 1 Abstract

## 1 Abstract

In these talks we review some of the recent results on open strings and noncommutative gauge theories, starting from the early calculations of open strings in a constant electromagnetic background. We discuss both the neutral string and the charged string. In the latter case, the scaling limit that leads to noncommutative abelian gauge theory can be generalized to a scaling limit in which multiple noncommutativity parameters enter. Our approach corresponds to expanding a theory with $`U(N)`$ Chan-Paton factors around a background $`U(1)^N`$ gauge field with different magnetic fields in each $`U(1)`$. This scaling limit can be interpreted in terms of a matrix model. We also describe an open string model with a time-dependent noncommutativity parameter. This model is the open string version of a WZW model based on a non-semi-simple group. It has a time-dependent background, and a spacetime metric of the plane wave type supported by a Neveu-Schwarz two-form potential.

## 2 Introduction

There are many ways how string theory seems to challenge our current understanding of spacetime. Probably, string dualities are the most striking example of this statement, not to mention supersymmetry and extra dimensions. More recently, a new fact has emerged, the noncommutativity of spacetime variables at the string scale, a fact that obviously might have important implications for the structure of spacetime. Noncommutativity originally emerged in the context of open strings, starting from the treatment of open string field theory in . More recently, it reappeared in the context of matrix theory compactified on a torus . Finally, it showed up in the low energy description of strings in an electromagnetic background . In these notes we will concentrate in the last approach. Our interest is to investigate the emergence of noncommutativity in various open string models, and compute the noncommutativity parameters there. The question behind our investigation is to see what novel noncommutative space geometries can emerge from string theory. We will go from the single noncommutativity parameter of the neutral string to the multiple noncommutative parameters of the open string in a $`U(N)`$ background, and finally to the time-dependent noncommutativity parameters in a model with a time-dependent background.

A very interesting question of course is what is the generalization of the Moyal-Weyl product corresponding to more general noncommutativity situations. A related question is what is the corresponding Born-Infeld action. In the $`U(1)`$ case the noncommutativity can be seen to emerge directly from a rewriting/reinterpretation of the Born-Infeld action. How does that work in the non-abelian case? Is there an equivalent in the time-dependent case? These are still open questions, and we will not have much to say about them.

## 3 Neutral String

The propagator for the open neutral string in a slowly varying background $`U(1)`$ gauge field was computed in  by solving the equations of motion

$$\text{ }\text{ }\text{ }\text{ }<x^i(z,\overline{z})x^j(\zeta ,\overline{\zeta })>=2\pi \alpha ^{}\delta ^2(z\zeta )G^{ij}$$
(3.1)
with the boundary conditions

$$(_z_{\overline{z}})<x^i(z,\overline{z})x^j(\zeta ,\overline{\zeta })>+2\pi \alpha ^{}g^{ik}B_k\mathrm{}(_z+_{\overline{z}})<x^{\mathrm{}}(z,\overline{z})x^j(\zeta ,\overline{\zeta })>|_{z=\overline{z}}=0$$
(3.2)
where $`\text{ }\text{ }\text{ }\text{ }4_z_{\overline{z}}`$. Here the open string worldsheet $`\mathrm{\Sigma }`$ will denote the disc with Euclidean metric $`\gamma ^{\alpha \beta }=\delta ^{\alpha \beta }`$, and the complex coordinates $`z,\overline{z}`$ are related to the original strip coordinates $`\sigma ,\tau `$ with $`\tau `$ rotated to $`ti\tau `$ and $`\mathrm{}t\mathrm{},0\sigma \pi ,`$ by $`ze^{t+i\sigma },\overline{z}e^{ti\sigma }`$ with $`\mathrm{Im}z0`$.

The propagator was found to be

$`<x^i(z,\overline{z})x^j(\zeta ,\overline{\zeta })>=\alpha ^{}[\frac{1}{2}g^{ij}\mathrm{ln}(z\zeta )+\frac{1}{2}g^{ij}\mathrm{ln}(\overline{z}\overline{\zeta })`$
$`+(\frac{1}{2}g^{ij}+G^{ij}+{\displaystyle \frac{\theta ^{ij}}{2\pi \alpha ^{}}})\mathrm{ln}(z\zeta )`$
$`+(\frac{1}{2}g^{ij}+G^{ij}{\displaystyle \frac{\theta ^{ij}}{2\pi \alpha ^{}}})\mathrm{ln}(\overline{z}\overline{\zeta }){\displaystyle \frac{i}{2\alpha ^{}}}\theta ^{ij}]`$ (3.3)
where

$`G^{ij}`$ $`=`$ $`[(g+2\pi \alpha B)^1g(g2\pi \alpha B)^1]^{ij}`$
$`G_{ij}`$ $``$ $`g_{ij}(2\pi \alpha ^{})^2(Bg^1B)_{ij}=[(g2\pi \alpha B)g^1(g+2\pi \alpha B)]_{ij}`$
$`\theta ^{ij}`$ $``$ $`(2\pi \alpha ^{})^2[(g+2\pi \alpha B)^1B(g2\pi \alpha B)^1]^{ij}.`$ (3.4)
In  the above propagator was used to compute the “equal time” commutator of the string operators on the boundary via a short distance expansion procedure , and to define a noncommutativity parameter

$$[X^i,X^j]=i\theta ^{ij},$$
(3.5)
where $`\theta ^{ij}`$ is given in (3). This shows that the $`X^i`$ are coordinates in a noncommutative space with noncommutativity parameter $`\theta `$. Indeed in  it was also shown that, in the scaling limit $`\alpha ^{}=\sqrt{ϵ}0`$ and $`g_{ij}=ϵ0`$, a noncommutative gauge theory emerged with field strength

$$\widehat{F}_{ij}=_i\widehat{A}_j_j\widehat{A}_ii\widehat{A}_i\widehat{A}_j+i\widehat{A}_j\widehat{A}_i$$
(3.6)
where $`\widehat{A}_i`$ is related to an ordinary gauge potential $`A_i`$ via the relation

$$\widehat{F}=\frac{1}{1+F\theta }F,$$
(3.7)
known as the Seiberg-Witten map . In (3.6) the $``$ product is the Moyal-Weyl product

$$(fg)(x)=e^{\frac{i}{2}\theta ^{ij}_i_{}^{}{}_{j}{}^{}}f(x)g(x^{})|_{x=x^{}}$$
(3.8)
The relation between this noncommutative field strength and the commutative one we started from is encapsulated in the Born-Infeld action $`\sqrt{det(g+2\pi \alpha ^{}(B+F))},`$ the effective target space action for an open string with slowly varying background fields . It turns out that

$$\sqrt{det(g+2\pi \alpha ^{}(B+F))}\sqrt{det(G+2\pi \alpha ^{}\widehat{F})}$$
(3.9)
where $`g_{ij}`$ is the ‘closed string metric’ and $`G_{ij}`$ is the ‘open string metric’ which appear in (3).

## 4 Charged Open String

We will now find the noncommutativity parameter in the charged open string case . Rather than studying the scaling limit of an open string theory in a fully non-abelian $`U(N)`$ background (which is still difficult), we will restrict ourselves to backgrounds that reside in the $`U(1)^N`$ Cartan subgroup, but with different background $`U(1)`$ fields on each brane. So we will discuss the charged string, i.e. a string with different $`U(1)`$ backgrounds at each end. We consider magnetic backgrounds with constant field strength. In this case we have the exact mode expansion to start with, originally derived in  and , and we use it to compute the propagator on the disk. Starting from this propagator, we repeat the argument in  to compute two noncommutativity parameters, one at each end of the string. This emphasizes the interpretation that the noncommutativity of the $`D`$-brane worldvolume in the presence of a background $`B`$-field along the brane is really a property of the endpoint of the string, rather than a feature of the worldvolume itself. For $`N`$ $`D`$-branes, there are $`N`$ different noncommutativity parameters, even when the $`D`$-branes are coincident. Hence, the short distance behavior of the operator products of tachyon vertex operators which are inserted on the boundaries defines star products with different noncommutativity parameters. These enter into the computation of the scattering amplitudes of the scaling limit theory.

The worldsheet action for the “charged” string with different magnetic fields at each end is

$`S={\displaystyle \frac{1}{4\pi \alpha ^{}}}{\displaystyle _\mathrm{\Sigma }}g_{\mu \nu }_\alpha X^\mu ^\alpha X^\nu {\displaystyle \frac{i}{2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t(B_{ij}^{(1)}X^i_tX^j|_{\sigma =0}+B_{ij}^{(2)}X^i_tX^j|_{\sigma =\pi }).`$ (4.10)
Here $`0\mu ,\nu 25`$ and the open strings end on $`Dp`$-branes in the $`(0,i)`$ directions for $`1ip`$. Variation of (4.10) gives the equations of motion for the worldsheet field

$$(_\sigma ^2+_t^2)X^\mu (\sigma ,t)=0$$
(4.11)
together with the boundary conditions at each end of the string

$`g_{ij}_\sigma X^j+2\pi i\alpha ^{}B_{ij}^{(1)}_tX^j|_{\sigma =0}=0`$
$`g_{ij}_\sigma X^j2\pi i\alpha ^{}B_{ij}^{(2)}_tX^j|_{\sigma =\pi }=0`$ (4.12)
and

$`_\sigma X^0|_{\sigma =0,\sigma =\pi }=0`$
$`X^I|_{\sigma =0}=0,p+1I25.`$ (4.13)
Note that the worldsheet action for the neutral string, whose propagator in the directions along the $`Dp`$-brane is given in (LABEL:anprop), is a special case of (4.10) with $`B_{ij}^{(1)}=B_{ij}^{(2)}`$. For simplicity, we now specialize to the case where the ends of the string live on $`D2`$-branes, i.e. $`i=1,2`$. In this case the magnetic fields have only one component and we relabel them as $`B_{12}^{(1)}=q_1B_{12}`$ and $`B_{12}^{(2)}=q_2B_{12}`$, where $`q_1+q_20`$.

In a basis given by $`X^\pm (\sigma ,t)X^1(\sigma ,t)\pm iX^2(\sigma ,t)`$, the charged string normal mode expansion  can be written as

$$X^+(z,\overline{z})=x^++\frac{i}{2}\sqrt{2\alpha ^{}}\underset{r𝒵+A}{}\frac{a_r}{r}(z^r+\overline{z}^r)\frac{1}{2}\sqrt{2\alpha ^{}}B\underset{r𝒵+A}{}\frac{a_r}{r}(z^r\overline{z}^r)$$
(4.14)
$$X^{}(z,\overline{z})=x^{}+\frac{i}{2}\sqrt{2\alpha ^{}}\underset{s𝒵A}{}\frac{\stackrel{~}{a}_s}{s}(z^s+\overline{z}^s)+\frac{1}{2}\sqrt{2\alpha ^{}}B\underset{s𝒵A}{}\frac{\stackrel{~}{a}_s}{s}(z^s\overline{z}^s).$$
(4.15)
with commutation relations for the zero modes given by

$`[a_r,\stackrel{~}{a}_s]=2Gr\delta _{r,s};[a_r,a_r^{}]=0=[\stackrel{~}{a}_s,\stackrel{~}{a}_s^{}];`$
$`[x^+,x^{}]={\displaystyle \frac{2}{(q_1+q_2)B_{12}}};[a_r,x^\pm ]=0=[\stackrel{~}{a}_s,x^\pm ].`$ (4.16)
Here $`g_{ij}=g^1\delta _{ij}`$, $`G_{ij}=G^1\delta _{ij}`$, where now

$$G=\frac{g}{1+B^2},B=gq_12\pi \alpha ^{}B_{12},A=\frac{1}{\pi }(\mathrm{arctan}B+\mathrm{arctan}\frac{q_2}{q_1}B).$$
(4.17)
Unlike the neutral string case, oscillators $`a_r,\stackrel{~}{a}_s`$ are non-integrally moded, since $`r=n+A`$, $`s=nA`$, for $`n𝒵`$. We use the mode expansions computed above to compute the charged propagator. The details of this calculation can be found in . For $`|z|>|\zeta |`$, we find it is given by

$`<X^+(z,\overline{z})X^{}(\zeta ,\overline{\zeta })>\beta |X^+(z,\overline{z})X^{}(\zeta ,\overline{\zeta })|\alpha `$
$`={\displaystyle \frac{2\alpha ^{}\pi g}{B+\frac{q_2}{q_1}B}}2\alpha ^{}G\frac{1}{A}(\zeta ^A+\overline{\zeta }^A1)+2i\alpha ^{}GB\frac{1}{A}(\zeta ^A\overline{\zeta }^A)`$
$`+\alpha ^{}G[f(\frac{\zeta }{z})+f(\frac{\overline{\zeta }}{\overline{z}})+f(\frac{\zeta }{\overline{z}})+f(\frac{\overline{\zeta }}{z})]`$
$`+\alpha ^{}GB^2[f(\frac{\zeta }{z})+f(\frac{\overline{\zeta }}{\overline{z}})f(\frac{\zeta }{\overline{z}})f(\frac{\overline{\zeta }}{z})]`$
$`+2i\alpha ^{}GB[f(\frac{\zeta }{\overline{z}})+f(\frac{\overline{\zeta }}{z})]`$ (4.18)
where

$$f(\rho )\underset{r=n+A;n0}{}\frac{\rho ^r}{r};\underset{A0}{lim}f(\rho )=\mathrm{ln}(1\rho )+\underset{A0}{lim}\frac{\rho ^A}{A}.$$
(4.19)
As a first step toward computing the commutator on the boundary, we need to compute the propagator on the boundary. We distinguish the two boundary regions of the open string disk as follows. On the boundary $`\sigma =0`$, we have $`z=|z|=\tau `$ and $`\zeta =|\zeta |=\tau ^{}`$ so $`\tau ,\tau ^{}>0`$; while on $`\sigma =\pi `$, then $`z=|z|e^{i\pi }=\tau `$ and $`\zeta =|\zeta |e^{i\pi }=\tau ^{}`$ so here $`\tau ,\tau ^{}<0`$.

For $`|z|>|\zeta |`$, and on the boundary $`\sigma =0`$, the propagator is

$`<X^+(z,\overline{z})X^{}(\zeta ,\overline{\zeta })>|_{\sigma =0}`$ $`={\displaystyle \frac{2\alpha ^{}\pi g}{B+\frac{q_2}{q_1}B}}+4\alpha ^{}G{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n+A}}(\frac{\zeta }{z})^{n+A}`$ (4.20)
$`+{\displaystyle \frac{2\alpha ^{}G}{A}}{\displaystyle \frac{4\alpha ^{}G}{A}}\zeta ^A,`$
$`<X^{}(z,\overline{z})X^+(\zeta ,\overline{\zeta })>|_{\sigma =0}`$ $`={\displaystyle \frac{2\alpha ^{}\pi g}{B+\frac{q_2}{q_1}B}}+4\alpha ^{}G{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{nA}}(\frac{\zeta }{z})^{nA}`$ (4.21)
$`+{\displaystyle \frac{2\alpha ^{}G}{A}}{\displaystyle \frac{4\alpha ^{}G}{A}}z^A.`$
By using the above equations, we compute the commutator at the $`\sigma =0`$ end of the string.

$`[X^+(\tau ),X^{}(\tau )]=T(X^+(\tau )X^{}(\tau ^{})X^+(\tau )X^{}(\tau ^+))`$
$`\underset{ϵ0}{lim}(<X^+(\tau )X^{}(\tau ϵ)><X^{}(\tau +ϵ)X^+(\tau )>),(\mathrm{for}ϵ>0)`$
$`=\underset{ϵ0}{lim}({\displaystyle \frac{4\alpha ^{}\pi g}{B+\frac{q_2}{q_1}B}}+4\alpha ^{}G[{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n+A}}(\frac{\tau ϵ}{\tau })^{n+A}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{nA}}(\frac{\tau }{\tau +ϵ})^{nA}]`$
$`{\displaystyle \frac{4\alpha ^{}G}{A}}(\tau ϵ)^A+{\displaystyle \frac{4\alpha ^{}G}{A}}(\tau +ϵ)^A)`$
$`={\displaystyle \frac{4\alpha ^{}\pi g}{B+\frac{q_2}{q_1}B}}+4\alpha ^{}G\pi \mathrm{cot}\pi A`$
$`=4\alpha ^{}\pi (g)^2{\displaystyle \frac{q_12\pi \alpha ^{}B_{12}}{1+(g)^2(q_12\pi \alpha ^{}B_{12})^2}}`$
$`=2\mathrm{\Theta }^{12}.`$ (4.22)
Notice that $`\mathrm{\Theta }^{12}`$ is the same expression that appears in the neutral string. The analogous calculation at $`\sigma =\pi `$ provides a different commutator.

$`[X^+(\tau ),X^{}(\tau )]=T(X^+(\tau )X^{}(\tau ^{})X^+(\tau )X^{}(\tau ^+))`$
$`\underset{ϵ0}{lim}(<X^+(\tau )X^{}(\tau +ϵ)><X^{}(\tau ϵ)X^+(\tau )>),(\mathrm{for}ϵ>0)`$
$`=4\alpha ^{}\pi (g)^2{\displaystyle \frac{q_2\pi \alpha ^{}B_{12}}{1+(g)^2(q_22\pi \alpha ^{}B_{12})^2}}`$
$`=2\stackrel{~}{\mathrm{\Theta }}^{12}.`$ (4.23)
In the limit $`q_1q_2`$, then $`\stackrel{~}{\mathrm{\Theta }}^{12}\mathrm{\Theta }^{12}`$. Indeed in the neutral string case, where both ends of the string are on the same $`D`$-brane, the noncommutativity parameter at one end of the string is equal to minus that of the other end. For $`U(N)`$ Chan-Paton factors, the background magnetic fields can take on N possible values, giving rise to $`N`$ noncommutativity parameters.

As in , one can show that in the scaling limit ($`\alpha ^{}0`$, keeping $`G`$ and $`\mathrm{\Theta }^{12}`$ fixed), the operator product expansion reduces to the star product

$`e^{ipX(\tau )}e^{iqX(\tau ^{})}`$ $``$ $`e^{\mathrm{\Theta }^{12}(p_{}q_+p_+q_{})}:e^{i(p+q)X(\tau ^{})}:`$ (4.24)
$``$ $`e^{ipX(\tau ^{})}e^{iqX(\tau ^{})}.`$
For $`\sigma =\pi `$ the same equations will hold with $`\mathrm{\Theta }^{12}`$ replaced by $`\stackrel{~}{\mathrm{\Theta }}^{12}`$.

The same results about noncommutativity parameters can be derived by starting with the charged string propagator on the annulus, as calculated in  from a charged string annulus propagator. Since the result is a short distance effect, it is independent of the topology of the worldsheet.

## 5 The Spectrum in the Scaling Limit and a Matrix Model

Our explicit knowledge of the mode expansion allows us to compute the spectrum of the gauge theory obtained in the scaling limit, and compare it to one predicted by the matrix model. In addition to the $`U(1)^N`$ massless gauge bosons, we find charged vector states that survive for each Landau level. The states of the limiting non-abelian noncommutative gauge theory are no longer massless, but rather tachyonic or massive.

As in , we consider the scaling limit $`g^1ϵ`$ and $`\alpha ^{}\sqrt{ϵ}`$, for $`ϵ0`$. Actually this limit means letting the dimensionless quantity $`\alpha ^{}B_{12}\sqrt{ϵ}`$, while keeping $`B_{12}`$ fixed. Then we have that the noncommutativity parameters are finite in the scaling limit and are given by $`\mathrm{\Theta }^{12}(q_1B_{12})^1`$, $`\stackrel{~}{\mathrm{\Theta }}^{12}(q_2B_{12})^1`$.

One can see from (4.17) that $`\mathrm{tan}\pi A=\frac{B+\frac{q_2}{q_1}B}{1\frac{q_2}{q_1}B^2}`$. In the scaling limit,

$$\mathrm{tan}\pi A\frac{(q_1+q_2)\sqrt{ϵ}}{2\pi q_1q_2}A\frac{(q_1+q_2)\sqrt{ϵ}}{2\pi ^2q_1q_2}.$$
(5.25)
Using these formulae, one finds  that the mass formulae for the two polarization states of the charged vectors are finite in the scaling limit and behave as

$`\mathrm{mass}^2={\displaystyle \frac{1}{2\alpha ^{}}}A(1+A){\displaystyle \frac{1}{2}}{\displaystyle \frac{(q_1+q_2)B_{12}}{2\pi ^2q_1q_2}}={\displaystyle \frac{(q_1+q_2)G}{q_2\mathrm{\Theta }^{12}}},`$
$`\mathrm{mass}^2={\displaystyle \frac{1}{2\alpha ^{}}}A(3A)){\displaystyle \frac{3}{2}}{\displaystyle \frac{(q_1+q_2)B_{12}}{2\pi ^2q_1^2q_2\mathrm{\Theta }^{12}}}=3{\displaystyle \frac{(q_1+q_2)G}{q_2\mathrm{\Theta }^{12}}}.`$ (5.26)
Also in the spectrum are the charged vectors at different Landau levels. For each charged boson, the two polarizations have different masses. They differ from those of (LABEL:mtach) by integer multiples of $`\frac{2(q_1+q_2)G}{q_2\mathrm{\Theta }^{12}}`$. (Notice that each Landau level has infinite degeneracy). All other states in the charged string spectrum become infinitely heavy and decouple in the scaling limit. So the complete spectrum of our $`U(1)^N`$ noncommutative field theory, which is derived from $`N`$ neutral and $`N^2N`$ charged string sectors, is described by $`N`$ massless neutral gluons and the charged vectors above.

We have seen that our theory, which is a $`U(2)`$ gauge theory expanded around a $`U(1)\times U(1)`$ background, has a scaling limit. The next question is to find the corresponding action.

It turns out that the noncommutative gauge theory action for the $`U(1)\times U(1)`$ case can be described by a matrix model. A matrix theory action in terms of two infinite-dimensional matrices $`X^1,X^2`$ is

$$S=\frac{1}{2g_{YM}^2}Tr_{}[X^1,X^2]^2$$

where $`X^1,X^2`$ are some operators acting in a Hilbert space $``$. Its equations of motion are

$$\delta _k\mathrm{}[X^k,[X^{\mathrm{}},X^j]]=0.$$

A solution  is $`X_{\mathrm{sol}}^i`$ given by a $`2\times 2`$ block matrix of the form

$$X_{\mathrm{sol}}^i=\left(\begin{array}{cc}y^i& 0\\ 0& z^i\end{array}\right)$$

where $`y,z`$ satisfy the Heisenberg algebras

$$[y^i,y^j]=i\theta _1^{ij}\mathrm{𝟏}$$

$$[z^i,z^j]=i\theta _2^{ij}\mathrm{𝟏},$$

i.e. $`X_{\mathrm{sol}}^1+iX_{\mathrm{sol}}^2=\left(\begin{array}{cc}\sqrt{2\theta _1}a_1& 0\\ 0& \sqrt{2\theta _2}a_2\end{array}\right).`$ Above, $`a_1,a_2`$ are the creation and annihilation operators in the Hilbert space $``$. Then it can be shown  by expanding the equations of motions around $`X^i=X_{\mathrm{sol}}^i+y^i`$ and computing the coefficients of the terms quadratic in $`y^i`$, that the spectrum of fluctuations about this classical solution coincides with the $`\alpha ^{}0`$ limit of the string spectrum described above. One needs to identify

$$\theta _1\theta _2=\frac{(q_1+q_2)G}{q_2\mathrm{\Theta }^{12}},$$

where the RHS follows the notation from (LABEL:mtach) .

## 6 Time-dependent Background

In most of the examples currently known, the noncommutativity parameter is constant. An obvious task is to look for time-dependent noncommutativity parameters. In  we studied an open string model, whose target space has a plane wave metric supported by a time-dependent Neveu-Schwarz two-form potential. This background was studied in  for closed strings, while in  we looked at the open string version. We quantized the sigma-model in light-cone gauge, computed the worldsheet propagator, and used it to derive a time-dependent noncommutativity parameter. Indeed, for large values of the time parameter, this model resembles a neutral string in a constant background $`B`$ field, hence it is a good candidate for spacetime noncommutativity.

We do not solve the exact open mode expansion in closed form, but we compute it as a power series in a suitable parameter $`\mu `$. This expansion is adequate to show noncommutativity. To do that, we follow the same strategy adopted for the neutral and charged string. We compute a mode expansion, derive from it a worldsheet propagator on the disk, and evaluate the commutator on the boundaries to find a time-dependent noncommutativity parameter.

The worldsheet action coupling a string to a general metric and background Neveu-Schwarz field is

$$S=_\mathrm{\Sigma }𝑑\tau 𝑑\sigma [\sqrt{\gamma }\gamma ^{\alpha \beta }G_{MN}_\alpha X^M_\beta X^M+B_{MN}ϵ^{\alpha \beta }_\alpha X^M_\beta X^N]$$
(6.27)
where we choose the string worldsheet $`\mathrm{\Sigma }`$ with Lorentz signature, and have rescaled the scalar worldsheet fields by $`(2\sqrt{\pi \alpha ^{}})^1`$ so that the $`X^M`$ are dimensionless. We consider the time-dependent background provided by the model based on a non-semi-simple group discussed in , and adopt the same notation, with $`X^M=(a_1,a_2,u,v)`$, and $`u`$ being identified with the time in the target space.

The background field $`G_{MN}`$ and $`B_{MN}`$ are given by

$$G_{MN}=\left(\begin{array}{cccc}1& 0& \frac{a_2}{2}& 0\\ 0& 1& \frac{a_1}{2}& 0\\ \frac{a_2}{2}& \frac{a_1}{2}& b& 1\\ 0& 0& 1& 0\end{array}\right),B_{MN}=\left(\begin{array}{cccc}0& u& 0& 0\\ u& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$
(6.28)
The Lorentz signature target space metric $`G_{MN}`$ can be recognized as a plane wave metric . The background is time-dependent because of the the $`u`$-dependence of $`B_{12}`$. In  it was shown that this model is exactly conformally invariant (i.e. to all orders in $`\alpha ^{}`$) by checking the one-loop $`\beta `$ function equations for the closed string backgrounds, and then proving that there were no higher loop contributions.

Here, since we are interested in noncommutativity, we consider open string boundary conditions. This case is not conformally invariant. <sup>1</sup><sup>1</sup>1A conformally invariant version of (6.28) is studied in , but its noncommutativity parameter although non-constant, is not time-dependent. The background (6.28) satisfies the the Born-Infeld field equations for $`Nu`$,

$$(D_MF_{NL})(1F^2)^{1LM}=0$$
(6.29)
where $`(1F^2)^{1LM}=(1+F)^{1LP}G_{PN}(1F)^{1NM}`$ and $`(1F)_{MN}G_{MN}2\pi \alpha ^{}F_{MN}.`$ In our case $`F_{MN}=B_{MN}`$. Using (6.28), one can check that the nonvanishing components of the Ricci tensor and affine connections are $`R_{uu}=\frac{1}{2}`$, $`\mathrm{\Gamma }_{uj}^i=\frac{1}{2}ϵ_j^i`$, $`\mathrm{\Gamma }_{ui}^v=\frac{a^i}{4}.`$ It follows that $`(D_MF_{NL})(1F^2)^{1LM}=ϵ_{ij}(1F^2)^{1ju}=0,`$ for $`Nu`$. But

$$(D_MF_{uL})(1F^2)^{1LM}=\frac{u}{1+(2\pi \alpha ^{}u)^2}.$$
(6.30)
In terms of the $`a_i,u,v`$ variables, the sigma model action is

$$S=_\mathrm{\Sigma }𝑑\tau 𝑑\sigma [\sqrt{\gamma }\gamma ^{\alpha \beta }(_\alpha a^i_\beta a^i+2_\alpha u_\beta v+b_\alpha u_\beta u+ϵ_{ij}_\alpha u_\beta a^ia^j)+ϵ^{\alpha \beta }ϵ_{ij}u_\alpha a^i_\beta a^j].$$
(6.31)
Although our background is not conformally invariant, we will consider a light-cone version of the sigma model in order to study open string propagators in a B-field with linear time dependence. We identify the target space time $`u`$ with the worldsheet time $`\tau `$ via $`u=\mu \tau `$, where $`\mu `$ is a dimensionless parameter. The equations of motion and boundary conditions for the transverse fields $`a^i`$ written in terms of $`Xa^1+ia^2`$ and $`\stackrel{~}{X}a^1ia^2`$ become:

$`\text{ }\text{ }\text{ }Xi\mu (_\sigma X_\tau X)=0,\text{ }\text{ }\text{ }\stackrel{~}{X}+i\mu (_\sigma \stackrel{~}{X}_\tau \stackrel{~}{X})=0,`$
$`[_\sigma X+i\mu \tau _\tau X]|_{\sigma =0,\pi }=0,[_\sigma \stackrel{~}{X}i\mu \tau _\tau \stackrel{~}{X}]|_{\sigma =0,\pi }=0,`$ (6.32)
where $`\text{ }\text{ }\text{ }\text{ }_\tau ^2+_\sigma ^2=4z\overline{z}_z_{\overline{z}}.`$

For large $`\tau `$ (so that $`\tau `$ can be considered constant), notice the similarity of the boundary condition in (LABEL:eomlg) with the boundary condition for an open string in a background $`B`$ field. Since in the latter case the noncommutativity parameter is proportional to the background, this suggests we should expect here a noncommutativity parameter which depends on time.

The solution of (LABEL:eomlg) is given by the normal mode expansion for the transverse coordinates $`X`$ and $`\stackrel{~}{X}`$, to first order in $`\mu `$:

$`X(\sigma ,\tau )=x_0+a_0[\tau +\mu (i\tau \sigma +{\displaystyle \frac{i}{2}}\tau ^2)]`$
$`+{\displaystyle \underset{n0}{}}a_ne^{in\tau }[{\displaystyle \frac{i}{n}}\mathrm{cos}n\sigma +\mu (({\displaystyle \frac{1}{2n^2}}i{\displaystyle \frac{\tau }{n}})\mathrm{sin}n\sigma +({\displaystyle \frac{i}{2n^2}}+{\displaystyle \frac{(\sigma \tau )}{2n}})\mathrm{cos}n\sigma )]`$ (6.33)
$`\stackrel{~}{X}(\sigma ,\tau )=\stackrel{~}{x}_0+\stackrel{~}{a}_0[\tau \mu (i\tau \sigma +{\displaystyle \frac{i}{2}}\tau ^2)]`$
$`+{\displaystyle \underset{n0}{}}\stackrel{~}{a}_ne^{in\tau }[{\displaystyle \frac{i}{n}}\mathrm{cos}n\sigma \mu (({\displaystyle \frac{1}{2n^2}}i{\displaystyle \frac{\tau }{n}})\mathrm{sin}n\sigma +({\displaystyle \frac{i}{2n^2}}+{\displaystyle \frac{(\sigma \tau )}{2n}})\mathrm{cos}n\sigma )].`$ (6.34)
## 7 Time-dependent Noncommutativity

Having found a mode expansion, the propagator on the disk can be computed along the lines of . We will use the notations $`z=e^{i(\tau +\sigma )}`$, $`\overline{z}=e^{i(\tau \sigma )}`$, $`\zeta =e^{i(\tau ^{}+\sigma ^{})}`$ and $`\overline{\zeta }=e^{i(\tau ^{}\sigma ^{})}`$ Then for $`|z|>|\zeta |`$, the propagator to order $`\mu `$ is

$`<X(z,\overline{z})\stackrel{~}{X}(\zeta ,\overline{\zeta })>`$ $`=`$ $`i4\alpha ^{}(\tau +\mu (i\tau \sigma +{\displaystyle \frac{i}{2}}\tau ^2))`$ (7.35)
$`+`$ $`4\alpha ^{}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}e^{in(\tau \tau ^{})}[{\displaystyle \frac{1}{n}}\mathrm{cos}n\sigma \mathrm{cos}n\sigma ^{}`$
$`+`$ $`i\mu \mathrm{cos}n\sigma ^{}(({\displaystyle \frac{1}{2n^2}}+{\displaystyle \frac{i\tau }{n}})\mathrm{sin}n\sigma ({\displaystyle \frac{i}{2n^2}}+{\displaystyle \frac{(\sigma \tau )}{2n}})\mathrm{cos}n\sigma )`$
$``$ $`i\mu \mathrm{cos}n\sigma (({\displaystyle \frac{1}{2n^2}}{\displaystyle \frac{i\tau ^{}}{n}})\mathrm{sin}n\sigma ^{}+({\displaystyle \frac{i}{2n^2}}{\displaystyle \frac{(\sigma ^{}\tau ^{})}{2n}})\mathrm{cos}n\sigma ^{})`$
$``$ $`{\displaystyle \frac{\mu }{n^2}}\mathrm{cos}n\sigma \mathrm{cos}n\sigma ^{}]`$
$`+`$ $`\mu (c_1\tau +c_0).`$
We are free to add the function $`\mu (c_1\tau +c_0)`$ to the expression since it does not affect the equation of motion or the boundary condition for the propagator to first order in $`\mu `$. For $`|z|>|\zeta |`$, the expression for $`<\stackrel{~}{X}(z,\overline{z})X(\zeta ,\overline{\zeta })>`$ is given by letting $`\mu \mu `$ in the above propagator. In the $`\mu 0`$ limit, these propagators reduce to the open bosonic string propagator $`lim_{\mu 0}<X(z,\overline{z})\stackrel{~}{X}(\zeta ,\overline{\zeta })>=2\alpha ^{}(\mathrm{ln}|z\zeta |+\mathrm{ln}|z\overline{\zeta }|)`$.

We evaluate the propagator on the worldsheet boundary at $`\sigma =0`$ and $`\sigma =\pi `$. We denote the points on the boundary at $`\sigma =0`$ with $`z=e^{i\tau }𝒯`$, and $`\zeta =e^{i\tau ^{}}=𝒯^{}`$. We get on the $`\sigma =0`$ boundary

$`<X(z,\overline{z})\stackrel{~}{X}(\zeta ,\overline{\zeta })>|_{\sigma =0}`$ $`=`$ $`i4\alpha ^{}(\tau +\mu {\displaystyle \frac{i}{2}}\tau ^2)+\mu (c_1\tau +c_0)`$
$``$ $`4\alpha ^{}\mathrm{ln}(1e^{i(\tau \tau ^{})})2\alpha ^{}\mu i(\tau \tau ^{})\mathrm{ln}(1e^{i(\tau \tau ^{})})`$
$`=`$ $`4\alpha ^{}\mathrm{ln}(𝒯𝒯^{})`$
$`+`$ $`\mu (2\alpha ^{}\mathrm{ln}^2𝒯2\alpha ^{}\mathrm{ln}({\displaystyle \frac{𝒯}{𝒯^{}}})\mathrm{ln}(1{\displaystyle \frac{𝒯^{}}{𝒯}})+(c_1i\mathrm{ln}𝒯+c_0))`$
$`<\stackrel{~}{X}(z,\overline{z})X(\zeta ,\overline{\zeta })>|_{\sigma =0}`$ $`=`$ $`i4\alpha ^{}(\tau \mu {\displaystyle \frac{i}{2}}\tau ^2)\mu (c_1\tau +c_0)`$ (7.36)
$``$ $`4\alpha ^{}\mathrm{ln}(1e^{i(\tau \tau ^{})})+2\alpha ^{}\mu i(\tau \tau ^{})\mathrm{ln}(1e^{i(\tau \tau ^{})})`$
Then at $`\sigma =0`$ the commutator is

$`[X(𝒯),\stackrel{~}{X}(𝒯)]`$ $`=`$ $`T(X(𝒯)\stackrel{~}{X}(𝒯^{})X(𝒯)\stackrel{~}{X}(𝒯^+))`$ (7.37)
$``$ $`\underset{ϵ0}{lim}(<X(𝒯)\stackrel{~}{X}(𝒯ϵ)><\stackrel{~}{X}(𝒯+ϵ)X(𝒯)>),(\mathrm{for}ϵ>0)`$
$`=`$ $`\mu (4i\alpha ^{})(\pi \mathrm{ln}𝒯i\mathrm{ln}^2𝒯)`$
$`=`$ $`\mu \mathrm{\hspace{0.17em}4}\alpha ^{}(\pi \tau +\tau ^2)\mathrm{\Theta },`$
where we chose $`c_1=2\pi \alpha ^{}`$, $`c_0=0`$, and use $`lim_{ϵ0}(\mathrm{ln}(1+ϵ)\mathrm{ln}ϵ)=0`$. The noncommutativity parameter $`\mathrm{\Theta }`$ is time-dependent.

Similarly, at $`\sigma =\pi `$ the commutator is

$`[X(𝒯),\stackrel{~}{X}(𝒯)]`$ $`=`$ $`T(X(𝒯)\stackrel{~}{X}(𝒯^{})X(𝒯)\stackrel{~}{X}(𝒯^+))`$ (7.38)
$``$ $`\underset{ϵ0}{lim}(<X(𝒯)\stackrel{~}{X}(𝒯+ϵ)><\stackrel{~}{X}(𝒯ϵ)X(𝒯)>),(\mathrm{for}ϵ>0)`$
$`=`$ $`(i4\alpha ^{})\mu [\pi \mathrm{ln}𝒯i\mathrm{ln}^2𝒯]`$
$`=`$ $`\mu \mathrm{\hspace{0.17em}\hspace{0.17em}4}\alpha ^{}(\pi \tau +\tau ^2)`$
Thus for small $`\mu `$, we have

$`\mathrm{\Theta }`$ $`=`$ $`\mu \mathrm{\hspace{0.17em}\hspace{0.17em}4}\alpha ^{}(\pi \tau +\tau ^2)\mathrm{at}\sigma =0,`$
$`\mathrm{\Theta }`$ $`=`$ $`\mu \mathrm{\hspace{0.17em}\hspace{0.17em}4}\alpha ^{}(\pi \tau +\tau ^2)\mathrm{at}\sigma =\pi .`$ (7.39)
For small $`\tau `$, the theta parameter at the $`\sigma =0`$ end of the string is minus that at the $`\sigma =\pi `$ end. This is the case for the neutral string in a constant background $`B`$ field as well. In fact, although we have worked only to lowest order in $`\mu `$, we can see directly from the equations of motion and boundary conditions (in $`z`$, $`\overline{z}`$) variables in the limit of large $`z`$, i.e. large $`i\tau `$, a limit for which $`z^10`$, that the system reduces to the neutral string with the identification $`\mu \tau =B`$, a constant. (In the large $`\tau `$ limit, we note that $`\mathrm{ln}|z|`$ is approximately constant, in the sense that it is changing slowly, i.e. its derivative $`|z|^1`$ is small. Therefore, for large $`\tau `$ the noncommutativity parameter becomes constant, and our model is similar to the neutral string.) For large $`\tau `$, using the neutral string expressions, we find the noncommutativity parameter be time-dependent

$`\mathrm{\Theta }`$ $`=`$ $`4\alpha ^{}\pi B=4\alpha ^{}\mu \pi \tau \mathrm{at}\sigma =0,`$
$`\mathrm{\Theta }`$ $`=`$ $`4\alpha ^{}\pi B=4\alpha ^{}\mu \pi \tau \mathrm{at}\sigma =\pi .`$ (7.40)
We have shown that our model exhibits noncommutativity for both small and large $`\tau `$. The expectation is that the model will remain noncommutative with a time-dependent noncommutativity parameter for all times.

## 8 Conclusions

The noncommutativity properties of string theories are interesting at many levels. In string theory, noncommutativity has lead to new insights and new techniques. From a mathematical point of view, noncommutativity in string theory is interesting because of its connection to noncommutative geometries and noncommutative algebras. Many open questions remain. One is what is the generalization of the Moyal-Weyl product in the case of multiple non commutativity parameters, as in the case of non-abelian gauge theories. What is the generalization of the map in (3.7)? Partial answers have been attempted in , but much more needs to be understood. In fact, we might be just at the beginning of our understanding of the relation between strings and noncommutativity. We suspect that the inclusion of additional terms in the worldsheet action, such as mass terms or Ramond backgrounds , will lead to novel types of noncommutativities. Probably there is a sort of “noncommutative democracy” among the various pieces of the worldsheet action. This set of noncommutativity parameters, one for each of the additional terms in the worldsheet action, could be reflected in the interpretation of the open string field theory star product as a continuous tensor product of Moyal products.

Finally, the most interesting aspect is that spacetime noncommutativity discussed above offers an insight on the structure of spacetime, as it implies a spacetime uncertainty principle

$$\mathrm{\Delta }X^\mu \mathrm{\Delta }X^\nu \alpha ^{}.$$

This relation suggests the existence of a spacetime cutoff at the Planck scale, and hence possible deviations from the smoothness of spacetime at small distances. It might well be that these aspects of Planck scale physics lead to distinctive signatures observable in cosmology. In particular, they could leave their mark on the spectrum of density fluctuations in the early universe, and modify the inflationary perturbation spectrum.


format: text

text:
# Incompressible Couette Flow Thanks for Grzegorz Juraszek (for English languague checking).

## 1 Introduction

Main problem is shown in figure (1). There is viscous flow between two parallel plates. Upper plate is moving in x direction with constans velocity $`(U=U_e)`$. Lower one is not moving $`(U=0)`$. We are looking for a solution to describe velocity vector field in the model (between two plates).

## 2 Fundamental Equations

Most of incompressible fluid mechanics (dynamics) problems are described by simple Navier-Stokes equation for incompressible fluid velocity, which can be written with a form:

$$\frac{\stackrel{}{u}}{t}=(\stackrel{}{u}\stackrel{}{})\stackrel{}{u}\stackrel{}{}\phi +\upsilon \stackrel{}{}^2\stackrel{}{u}+\stackrel{}{g},$$
(1)
where $`\phi `$ is defined is defined as the relation of pressure to density:

$$\phi =\frac{p}{\rho }$$
(2)
and $`\upsilon `$ is a kinematics viscosity of the fluid.

We will also use a continuity equation, which can be written as follows 111I assume constans density of the fluid.:

$$D=\stackrel{}{}\stackrel{}{v}=0$$
(3)
Of course, in a case of couette incompressible flow we will use several simplifications of (1).

## 3 Mathematical Formulation of the Couette Problem

Incompressible couette problem is not needed to solve full Navier-Stokes equations. There is no external force, so first simplification of (1) will be:

$$\frac{\stackrel{}{u}}{t}=(\stackrel{}{u}\stackrel{}{})\stackrel{}{u}\stackrel{}{}\phi +\upsilon \stackrel{}{}^2\stackrel{}{u},$$
(4)
In  there can be found easy proof that in couette problem there are no pressure gradients, which means that:

$$\stackrel{}{}\phi =0,$$
(5)
We will ignore a convection effects so, equation (4) can be written with a form:

$$\frac{\stackrel{}{u}}{t}=\upsilon \stackrel{}{}^2\stackrel{}{u}$$
(6)
Now we have simple differential equation for velocity vector field. That equation is a vector type and can be simplified even more. Let us write continuity equation (3) in differential form. Let $`\stackrel{}{u}=(u,v)`$, then continuity equation can be expanded as follows222Only with assumption of non compressible fluid.:

$$\frac{u}{x}+\frac{v}{y}=0$$
(7)
We know that there is no $`u`$ velocity component gradient along x axis (symmetry of the problem), so:

$$\frac{u}{x}=0$$
(8)
Evaluation of Taylor series333Those simple expressions can be found in , chapter 9.2. at points $`y=0`$ and $`y=D`$ gives us a proof that only one possible and physically correct value for y component of velocity $`\stackrel{}{u}`$ is:

$$v=0$$
(9)
Because of (9) equation (6) can be written as follows:

$$\frac{u}{t}=\upsilon \frac{^2u}{y^2}$$
(10)
Our problem is now simplified to mathematical problem of solving equations like (10). That is now a governing problem of incompressible couette flow analysis.

### 3.1 Analytical soulution

Analytical solution for velocity profile of steady flow, without time-changes (steady state) can be found very easily in an equation:

$$\upsilon \frac{^2u}{y^2}=0$$
(11)
And without any changes of viscosity $`\upsilon `$ it can be written in form:

$$\frac{^2u}{y^2}=0$$
(12)
After simple two times integration of equation (12) we have analytical solution function of (12):

$$u=c_1y+c_2$$
(13)
where $`c_1`$ and $`c_2`$ are integration constans.

### 3.2 Boundary Conditions for The Analytical Solution

Simple boundary conditions are provided in that problem. We know that:

$$u=\{\begin{array}{cc}0\hfill & \text{y=0}\hfill \\ u_e\hfill & \text{y=D}\hfill \end{array}$$
(14)
Simple applying it to our solution (13) gives a more specified one, where $`c_1=\frac{u_e}{D}`$ and $`c_2=0`$:

$$u=\frac{u_e}{D}y$$
(15)
It means, that a relationship between $`u`$ and $`y`$ is linear. A Better idea to write that with mathematical expression is:

$$\frac{u}{y}=\frac{u_e}{D}$$
(16)
Where $`\frac{u_e}{D}`$ is a constans for the problem (initial $`u`$ velocity vs size of the model).

## 4 Numerical Solution

### 4.1 Non-dimensional Form

Let us define some new non-dimensional variables444Exacly the same, like in .:

$$u=\{\begin{array}{ccc}u^{^{}}\hfill & =\hfill & \frac{u}{u_e}\hfill \\ y^{^{}}\hfill & =\hfill & \frac{y}{D}\hfill \\ t^{^{}}\hfill & =\hfill & \frac{t}{D/u_e}\hfill \end{array}$$
(17)
Now let us place these variables into equation (10) and we have now a non-dimensional equation written as follows:

$$\rho \frac{u/u_e}{(tu_e)/D}\left(\frac{u_e^2}{D}\right)=\upsilon \frac{^2(u/u_e)}{(y/D)^2}\left(\frac{u_e}{D^2}\right)$$
(18)
Now we replace all the variables to nondimensional, like defined in (17):

$$\rho \frac{u^{^{}}}{t^{^{}}}\left(\frac{u_e^2}{D}\right)=\upsilon \frac{^2u^{^{}}}{y^{}_{}{}^{}2}\left(\frac{u_e}{D^2}\right)$$
(19)
Now we will remove all $`^{^{}}`$ chars from that equation (only for simplification of notation), and the equation becomes555Constans simplification also implemented to:

$$\frac{u}{t}=\left(\frac{\upsilon }{D\rho u_e}\right)\frac{^2u}{y^2}$$
(20)
In that equation Reynold’s number $`Re`$ appears, and is defined as:

$$Re=\frac{D\rho u_e}{\upsilon }$$
(21)
Where $`Re`$ is Reynold’s number that depends on $`D`$ height of couette model. Finally, the last form of the equation for the couette problem can be written as follows:

$$\frac{u}{t}=\left(\frac{1}{Re}\right)\frac{^2u}{y^2}$$
(22)
We will try to formulate numerical solution of the equation (22).

### 4.2 Finite - Difference Representation

In our solution of equation (22) we will use Crank-Nicolson technique, so discrete representation of that equation can be written 666That representation is based on central discrete differential operators. as:

$$u_j^{n+1}=u_j^n+\frac{\mathrm{\Delta }t}{2(\mathrm{\Delta }y)^2Re}(u_{j+1}^{n+1}+u_{j+1}^n2u_j^{n+1}+u_{j1}^{n+1}+u_{j1}^n)$$
(23)
Simple grouping of all terms which are placed in time step (n+1) on the left side and rest of them - on right side, gives us an equation which can be written as:

$$Au_{j1}^{n+1}+Bu_j^{n+1}+Au_{j+1}^{n+1}=K_j$$
(24)
Where $`K_j`$ is known and depends only on values $`u`$ at $`n`$ time step:

$$K_j=\left(1\frac{\mathrm{\Delta }t}{(\mathrm{\Delta }y)^2Re}\right)u_j^n+\frac{\mathrm{\Delta }t}{2(\mathrm{\Delta }y)^2Re}(u_{j+1}^n+u_{j1}^n)$$
(25)
Constans $`A`$ and $`B`$ are defined as follows777Directly from equation (23).:

$$A=\frac{\mathrm{\Delta }t}{2(\mathrm{\Delta }y)^2Re}$$
(26)
$$B=1+\frac{\mathrm{\Delta }t}{(\mathrm{\Delta }y)^2Re}$$
(27)
### 4.3 Crank - Nicolson Implicit scheme

For numerical solution we will use one-dimensional grid points $`(1,\mathrm{},N+1)`$ where we will keep calculated $`u`$ velocities. That means $`u`$ has values from the range $`(u_1,u_2,\mathrm{},u_{N+1})`$. We know (from fixed boundary conditions), that: $`u_1=0`$ and $`u_{N+1}=0`$. Simple analysis of the equation (24) gives us a system of equations, which can be described by matrix equation:

$$𝐀\stackrel{}{X}=\stackrel{}{Y}$$
(28)
Where $`𝐀`$ is tridiagonal $`[N1][N1]`$ matrix of constant $`A`$ and $`B`$ values:

$$𝐀=\left[\begin{array}{cccccc}B\hfill & A\hfill & & & & \\ A\hfill & B\hfill & A\hfill & & & \\ & & \mathrm{}\hfill & & & \\ & & & \mathrm{}\hfill & & \\ & & & A\hfill & B\hfill & A\hfill \\ & & & & A\hfill & B\hfill \end{array}\right]$$
(29)
$`\stackrel{}{X}`$ vector is a vector of $`u`$ values:

$$\stackrel{}{X}=[u_1,u_2,\mathrm{},u_{N+1}]$$
(30)
$`\stackrel{}{Y}`$ vector is a vector of constans $`K_j`$ values:

$$\stackrel{}{Y}=[K_1,\mathrm{},K_{N+1}]$$
(31)
## 5 Solving The System of Linear Equations

Now the problem comes to solving the matrix - vector equation (28). There are a lot of numerical methods for that888Especially for tridiagonal matrices like $`𝐀`$, and we will try to choose two of them: Thomas and Gauss method. Both are very similar, and I will start with a description of my implementation with the simple Gauss method.

### 5.1 Gauss Method

Choice of the Gauss method for solving system of linear equations is the easiest way. This simple algorithm is well known, and we can do it very easily by hand on the paper. However, for big matrices (big $`N`$ value) a computer program will provide us with a fast and precise solution. A very important thing is that time spent on writing (or implementing, if Gauss procedure was written before) is very short, because of its simplicity.

I used a Gauss procedure with partial choice of a/the general element. That is a well known technique for taking the first element from a column of a matrix for better numerical accuracy.

The whole Gauss procedure of solving a system of equations contains three steps. First, we are look- ing for the general element.

After that, when a general element is in the first row (we make an exchange of rows999We made it for $`𝐀`$ matrix, and for $`\stackrel{}{X},\stackrel{}{Y}`$ too.) we make some simple calculations (for every value in every row and column of the matrix) for the simplified matrix to be diagonal (instead of a tridiagonal one which we have at the beginning). That is all, because after diagonalization I implement a simple procedure (from the end row to the start row of the matrix) which calculates the whole vector $`\stackrel{}{X}`$. There are my values of $`u_i`$ velocity in all the model101010More detailed description of Gauss method can be found in a lot of books on numerical methods, like ..

### 5.2 Thomas Method

Thomas’ method, described in  is simplified version of Gauss method, created especially for tridiagonal matrices. There is one disadvantage of Gauss method which disappears when Thomas’ method is implemented. Gauss method is rather slow, and lot of computational time is lost, because of special type of matrix. Tridiagonal matrices contain a lot of free (zero) values. In the Gauss method these values joins the calculation, what is useless.

Thomas’ simplification for tridiagonal matrices is to get only values from non-zero tridiagonal part of matrix. Simply applying a Thomas’ equations for our governing matrix equation (28) gives us:

$$d_i^{^{}}=B\frac{A^2}{d_{i1}^{^{}}}$$
(32)
$$u_i^{^{}}=u_i\frac{u_i^{^{}}A}{d_{i1}^{^{}}}$$
(33)
We know that exact value of $`u_M`$ is defined as follows:

$$u_M=\frac{u_M^{^{}}}{B}$$
(34)
Now solution of the system of equations will be rather easy. We will use recursion like that:

$$u_{i1}=\frac{u_{i1}^{^{}}Au_i}{B}$$
(35)
That easy recursion provides us a solution for the linear system of equations.

## 6 Results

Main results are provided as plots of the function:

$$\frac{y}{D}=f\left(\frac{u}{u_e}\right)$$
(36)
In figure (2) there is drawn an analytical solution to the problem of couette flow. That is linear function, and we expect that after a several time steps of numerical procedure we will have the same configuration of velocity field.

### 6.1 Different Time Steps

In figure (3) there are results of velocity $`u`$ calculation for several different time steps. Analytical solution is also drawn there.

As we can see in the figure (3) - the solution is going to be same as analytical one. Beginning state (known from boundary conditions) is changing and relaxing.

### 6.2 Results for Different Reynolds Numbers

In the figure (4) there is plot of numerical calculations for different Reynold’s numbers. For example Reynold’s number depends on i.e. viscosity of the fluid, size of couette model. As it is shown on the plot there is strong relationship between the speed of the velocity field changes and Reynold’s number. In a couple of words: when Reynolds number increases - frequency of changes also increases.

### 6.3 Results for Different Grid Density

In figure (5) there is an example of calculations of velocity field for different grid density ($`N`$ number). We see that there is also strong correlation between grid density, and speed of changes on the grid. Also, very interesting case $`N=10`$ shows, that for low density of the grid changes are very fast, and not accurate.

### 6.4 Conclusion

Solving of Incompressible Couette Problem can be good way to check numerical method, because of existing analytical solution. In that report there were presented two methods of solving system of equations: Gauss and Thomas’ method. System of equations was taken from Crank-Nicolson Implicit scheme. Well known linear relationships were observed.

## 7 APPENDIX A

```

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define N (40)
#define NN (N+1)

void Zamien(double *a, double *b) {
    double c;
    c=*a;   *a=*b;  *b=c;
}

void WypiszMacierz(double A[NN][NN], int n) {
    int i,j;

    for(j=1;j<n;j++)
    {
        for(i=1;i<n;i++)            // show matrix
        {
            printf("%2.4f   ",A[i][j]);
        }
        printf("\n");
    }
}

void Gauss(double A[NN][NN], double *b, double *x, int n) {
    int i,j,k;
    double m;

// Gauss Elimination

    for(i=0;i<n;i++)
    {
        // Step #1: Change governing element

        m=fabs(A[i][i]);
        k=i;

        for(j=i+1;j<n;j++)
        if(fabs(A[i][j])>m)
        {
            m=fabs(A[i][j]);
            k=j;
        }
        if(k!=i)
        for(j=0;j<n;j++)
        {
            Zamien(&A[j][i],&A[j][k]);
            Zamien(&b[i+1],&b[k+1]);
        }

        // Step #2: make it triangle
        for(j=i+1;j<n;j++)
        {
            m = A[i][j]/A[i][i];

            for(k=i;k<n;k++)
                A[k][j] = A[k][j] - m*A[k][i];

            b[j+1] = b[j+1] - m*b[i+1];
        }
    }
// Step#3: Solve now
    for(i=n-1;i>=1;i--)
    {
        for(j=i+1;j<n;j++)
            b[i+1] = b[i+1]-A[j][i]*x[j+1];
        x[i+1] = b[i+1]/A[i][i];
    }
}

int main(void)
{
    double U[N*2+2]={0},A[N*2+2]={0},B[N*2+2]={0},C[N*2+2]={0},D[N*2+2]={0},Y[N*2+2]={0};
// initialization
    double OneOverN = 1.0/(double)N;
    double Re=5000;                     // Reynolds number
    double EE=1.0;                      // dt parameter
    double t=0;
    double dt=EE*Re*OneOverN*2;         // delta time
    double AA=-0.5*EE;
    double BB=1.0+EE;
    int KKEND=1122;
    int KKMOD=1;
    int KK;                             // for a loop
    int i,j,k;                          // for loops too
    int M;                              // temporary needed variable
    double GMatrix[NN][NN]={0};             // for Gauss Elimination
    double test;
    Y[1]=0; // init

// apply boundary conditions for Couette Problem

    U[1]=0.0;
    U[NN]=1.0;

// initial conditions (zero as values of vertical velocity inside of the couette model)

    for(j=2;j<=N;j++)
        U[j]=0.0;

    A[1]=B[1]=C[1]=D[1]=1.0;

    for(KK=1;KK<=KKEND;KK++)
    {
        for(j=2;j<=N;j++)
        {
            Y[j]=Y[j-1]+OneOverN;
            A[j]=AA;
            if(j==N)
                A[j]=0.0;
            D[j]=BB;
            B[j]=AA;
            if(j==2)
                B[j]=0.0;
            C[j]=(1.0-EE)*U[j]+0.5*EE*(U[j+1]+U[j-1]);
            if(j==N)
                C[j]=C[j]-AA*U[NN];
        }

// Gauss
// C[]       -   free
// A[]B[]D[] -   for matrix calculation
// U[]       -   X

// calculate matrix for Gauss Elimination

        GMatrix[0][0]=D[1];
        GMatrix[1][0]=A[1];

        for(i=1;i<N-1;i++)
        {
            GMatrix[i-1][i]=B[i+1];     // GMatrix[1][2]=B[2]
            GMatrix[i][i]=D[i+1];       // GMatrix[2][2]=D[2]
            GMatrix[i+1][i]=A[i+1];     // GMatrix[3][2]=A[2]
        }

        GMatrix[N-2][N-1]=B[N];
        GMatrix[N-1][N-1]=D[N];

        Gauss(GMatrix,C,U,N);           // Gauss solving function

        Y[1]=0.0;
        Y[NN]=Y[N]+OneOverN;

        t=t+dt;                         // time increment
        test=KK % KKMOD;

       if(test < 0.01)                  // print the results
        {
            printf("KK,TIME\n");        // info 1
            printf("%d,%f\n",KK,t);

            printf(",J,Y[J],U[j]\n");   // info 2
            for(j=1;j<=NN;j++)
                printf("%d , %f, %f\n",j,U[j],Y[j]);
                printf("\n \n \n \n");      // for nice view of several datas
        }
     }

    return (1);
}

```

## 8 APPENDIX B

```

#include <stdlib.h>
#include <stdio.h>

#define N (50)
#define NN (N+1)

int main(void)
{
    double U[N*2+2],A[N*2+2],B[N*2+2],C[N*2+2],D[N*2+2],Y[N*2+2];
// initialization

    double OneOverN = 1.0/(double)N;
    double Re=7000;                     // Reynolds number
    double EE=1.0;                      // dt parameter
    double t=0;
    double dt=EE*Re*OneOverN*2;         // delta time
    double AA=-0.5*EE;
    double BB=1.0+EE;
    int KKEND=122;
    int KKMOD=1;
    int KK;                             // for a loop
    int j,k;                            // for loops too
    int M;                              // temporary needed variable

    double test;
    Y[1]=0; // init

// apply boundary conditions for Couette Problem
    U[1]=0.0;
    U[NN]=1.0;

// initial conditions (zero as values of vertical velocity inside of the couette model)
    for(j=2;j<=N;j++)
        U[j]=0.0;
    A[1]=B[1]=C[1]=D[1]=1.0;

    printf("dt=%f,    Re=%f,   N=%d \n",dt,Re, N);

    for(KK=1;KK<=KKEND;KK++)
    {
        for(j=2;j<=N;j++)
        {
            Y[j]=Y[j-1]+OneOverN;
            A[j]=AA;
            if(j==N)
                A[j]=0.0;

            D[j]=BB;
            B[j]=AA;

            if(j==2)
                B[j]=0.0;
            C[j]=(1.0-EE)*U[j]+0.5*EE*(U[j+1]+U[j-1]);

            if(j==N)
                C[j]=C[j]-AA*U[NN];
        }

        // upper bidiagonal form
        for(j=3;j<=N;j++)
        {
            D[j]=D[j]-B[j]*A[j-1]/D[j-1];
            C[j]=C[j]-C[j-1]*B[j]/D[j-1];
        }

        // calculation of U[j]
        for(k=2;k<N;k++)
        {
            M=N-(k-2);
            U[M]=(C[M]-A[M]*U[M+1])/D[M];   // Appendix A
        }

        Y[1]=0.0;
        Y[NN]=Y[N]+OneOverN;
        t=t+dt;                             // time increment
        test=KK % KKMOD;

        if(test < 0.01)                     // print the results
        {
            printf("KK,TIME\n");            // info 1
            printf("%d,%f\n",KK,t);

            printf(",J,Y[J],U[j]\n");       // info 2
            for(j=1;j<=NN;j++)
                printf("%d , %f, %f\n",j,U[j],Y[j]);
                printf("\n \n \n \n");      // for nice view of several datas
        }
    }
    return (1);
}

```


text:
# Optimum Monte Carlo Simulations: Some Exact Results

## 1 Introduction

Since the original Metropolis algorithm appeared five decades ago, countless studies have employed the technique to evaluate the thermodynamic properties of model systems\[Allen and Tildesley, 1987, Binder, 1997, Frenkel and Smit, 2002\]. The essence of the method is to generate a sequence of configurations that represent a given thermodynamic ensemble, often the canonical ensemble. Properties of interest are then obtained as averages over the configurations. At each step of the simulation a trial configuration is obtained from the current one by making a random displacement in the configuration space. This might correspond to, for example, displacing a randomly selected particle. The trial configuration is either accepted or rejected with a probability given by the appropriate Boltzmann factor for the ensemble. In case of rejection, the current configuration is retained for use in evaluating the properties of interest.

Since many applications of MC are computationally intensive, a much addressed issue has been the optimization of the simulation with respect to one or more control parameters so that the configuration space is sampled in the most efficient way. Bouzida, Kumar and Swendsen (BKS)\[Bouzida et al., 1992, Swendsen, 2002\], as part of a program aimed at improving the efficiency of MC simulations of biomolecules, performed numerical studies of the simple harmonic oscillator (SHO) where the convergence of the simulation depends on the maximum displacement. There is no unique measure of efficiency, but two simple choices are the mean squared and mean absolute displacements. As the maximum displacement tends to zero or infinity it is clear that the average value of both of these quantities tends to zero since, in the first case the particle does not move, while in the second all attempted moves are rejected. Thus both quantities have a maximum for some intermediate value of the maximum displacement.

It is also useful to consider the dynamical process associated with an MC simulation, even though it does not correspond to the actual dynamics of the system. One can, for example, calculate various time correlation functions that can be used to develop alternative efficiency criteria\[Kolafa, 1988, Mountain and Thirumalai, 1994\].

BKS\[Bouzida et al., 1992\] performed numerical studies of the SHO in one, two and three dimensions and examined the acceptance ratio $`P_{acc}`$ (i.e., the fraction of accepted trial configurations), mean squared, $`<(\mathrm{\Delta }x)^2)>`$, and mean absolute, $`<|\mathrm{\Delta }x|>`$, displacements as a function of the maximum displacement, $`\delta `$. They found that the acceptance ratio decreases approximately exponentially for small to intermediate values of $`\delta `$ and then inversely for larger values. In one dimension they found that the maxima in $`<(\mathrm{\Delta }x)^2>`$ and $`<|\mathrm{\Delta }x|>`$ occur at $`P_{acc}=0.42`$ and $`P_{acc}=0.56`$, respectively. In higher dimensions the results depend on how the jump is made. BKS considered two cases: in one the jumps are performed uniformly to any point in a spherical volume of radius $`\delta `$ centered on the current position, a choice that favors larger radial displacements for dimensions $`D>1`$. In a second method, the jumps were sampled uniformly in the radius (and randomly in the orientation) so that all radial displacements are equally probable. In the former case BKS observed that, for a given $`\delta `$, $`P_{acc}`$ decreases as a function of $`D`$, while for uniform radius sampling the numerical results suggested that $`P_{acc}`$ is independent of $`D`$.

In addition to these static properties, BKS also examined the correlation time, $`\tau `$, of the energy-energy correlation function. They observed a minimum correlation time for an acceptance ratio of approximately $`50\%`$.

Here we present an analytical study of the SHO in arbitrary dimension. We obtain exact expressions for the acceptance ratio and the mean squared and mean absolute displacements as functions of the maximum displacement $`\delta `$. We show that when the trial jump is selected uniformly in the radius, the results are independent of the dimension. We also present an analysis of the dynamics of the process.

## 2 ONE DIMENSION

We first investigate the case of a SHO in one dimension, whose potential energy is given by $`V(x)=kx^2/2`$ where $`x`$ is the position and $`k`$ is the stiffness constant.

In a standard Metropolis Monte Carlo Simulation one makes a trial move with a uniform random displacement selected between $`\delta `$ and $`\delta `$. The dynamical process generated by the successive trial moves of the Monte Carlo simulation can be written

$`{\displaystyle \frac{dP(x,t)}{dt}}`$ $`={\displaystyle \frac{1}{2\delta }}{\displaystyle _\delta ^\delta }𝑑hW(xx+h)P(x,t)+{\displaystyle \frac{1}{2\delta }}{\displaystyle _\delta ^\delta }𝑑hW(x+hx)P(x+h,t)`$ (1)
where $`W(xx+h)`$ denotes the transition rate from the state $`x`$ to the state $`x+h`$ and $`P(x,t)`$ is the probability to find the oscillator at position $`x`$ at time $`t`$. To ensure the convergence towards equilibrium, a sufficient condition is given by detailed balance, which is expressed as

$$\frac{W(xx+h)}{W(x+hx)}=\frac{P_{eq}(x+h)}{P_{eq}(x)}$$
(2)
where

$$P_{eq}(x)=c\mathrm{exp}(\beta V(x))$$
(3)
and $`c=\sqrt{\beta k/2\pi }`$ is a normalization constant ensuring that $`_{\mathrm{}}^{\mathrm{}}𝑑xP_{eq}(x)=1`$. One solution of Eq. (2) is the Metropolis rule,

$$W(xx+h)=\mathrm{min}(1,\mathrm{exp}(\beta (V(x+h)V(x))).$$
(4)
Most of the properties of the $`1D`$ SHO can be obtained analytically. For instance, the acceptance ratio, which is the number of accepted trials over the total number of trials can be expressed as

$$P_{acc}(\delta )=c_{\mathrm{}}^+\mathrm{}𝑑xe^{\beta V(x)}\frac{1}{2\delta }_\delta ^\delta 𝑑hW(xx+h).$$
(5)
In this equation $`\mathrm{exp}(\beta V(x))dx`$ is the probability that the oscillator is between $`x`$ and $`x+dx`$, $`dh/2\delta `$ is the probability of selecting a random displacement between $`h`$ and $`h+dh`$ and $`W(xx+h)`$ is the probability of accepting the trial displacement (given by Eq. (4)). Integration over the allowed values of $`x`$ and $`h`$ then gives the average acceptance probability. Since displacements to the left and right are symmetric, we need consider only one direction. For displacements to the right, Eq. (5), can be written as

$`{\displaystyle \frac{d(\delta P_{acc}(\delta ))}{d\delta }}`$ $`=c{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xe^{\beta V(x)}W(xx+\delta ).`$ (6)
For $`x<\delta /2`$, $`W(xx+\delta )=1`$ and for $`x>\delta /2`$, $`W(xx+\delta )=e^{\beta k((x+\delta )^2x^2)/2}`$. One thus obtains

$`{\displaystyle \frac{d(\xi P_{acc}(\xi ))}{d\xi }}`$ $`=1erf\left({\displaystyle \frac{\xi }{2}}\right)`$ (7)
where $`erf(x)`$ is the error function and $`\xi =\sqrt{\frac{\beta k}{2}}\delta `$. Using the initial condition, i.e., $`P_{acc}(0)=1`$, the solution of the differential equation (7) is

$$P_{acc}(\xi )=1erf\left(\frac{\xi }{2}\right)+\frac{2}{\sqrt{\pi }\xi }\left(1e^{\xi ^2/4}\right).$$
(8)
The function is plotted in Fig. 1.

The mean squared displacement, $`<(\mathrm{\Delta }x)^2>`$, and the mean absolute displacement, $`<|\mathrm{\Delta }x|>`$, defined as

$`<(\mathrm{\Delta }x)^2>`$ $`=c{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x\mathrm{exp}(\beta V(x)){\displaystyle \frac{1}{2\delta }}{\displaystyle _\delta ^\delta }𝑑hh^2W(xx+h),`$ (9)
$`<|\mathrm{\Delta }x|>`$ $`=c{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x\mathrm{exp}(\beta V(x)){\displaystyle \frac{1}{2\delta }}{\displaystyle _\delta ^\delta }𝑑h|h|W(xx+h).`$ (10)
can be quite simply obtained from the generating function

$$Z(\lambda )=c_{\mathrm{}}^+\mathrm{}𝑑x\mathrm{exp}(\beta V(x))\frac{1}{2\delta }_\delta ^\delta 𝑑h\mathrm{exp}(\lambda |h|)W(xx+h),$$
(11)
by the derivatives with respect to $`\lambda `$ and set $`\lambda =0`$, i.e., $`<|\mathrm{\Delta }x|>=(Z(\lambda )/\lambda )_{\lambda =0}`$, $`<(\mathrm{\Delta }x)^2>=(^2Z(\lambda )/\lambda ^2)_{\lambda =0}`$ (and, of course, $`P_{acc}=Z(\lambda =0)`$).

By multiplying both sides of Eq. (11) by $`\xi `$ and next differentiating with respect to $`\xi `$, one obtains

$$\frac{(\xi Z(\lambda ,\xi ))}{\xi }=\mathrm{exp}\left(\sqrt{\frac{2}{\beta k}}\lambda \xi \right)\left(1erf\left(\frac{\xi }{2}\right)\right)$$
(12)
which, after defining $`\stackrel{~}{\lambda }=\sqrt{\frac{2}{\beta k}}\lambda `$, leads to

$$Z(\lambda ,\xi )=\frac{1}{\stackrel{~}{\lambda }\xi }\left[1\mathrm{exp}(\stackrel{~}{\lambda }\xi )\left(1erf\left(\frac{\xi }{2}\right)\right)\mathrm{exp}(\stackrel{~}{\lambda }^2)\left(erf\left(\frac{\xi }{2}+\stackrel{~}{\lambda }\right)erf(\stackrel{~}{\lambda })\right)\right]$$
(13)
By taking the derivatives of the above formula with respect to $`\lambda `$ and evaluating the resulting expressions at $`\lambda =0`$, it is easy to show that the mean squared and mean absolute displacements are given by

$$\sqrt{\frac{\beta k}{2}}<|\mathrm{\Delta }x|>=\left[\frac{\xi }{2}\left(1erf\left(\frac{\xi }{2}\right)\right)\frac{1}{\sqrt{\pi }}\mathrm{exp}\left(\frac{\xi ^2}{4}\right)+\frac{erf\left(\frac{\xi }{2}\right)}{\xi }\right],$$
(14)
$$\frac{\beta k}{2}<(\mathrm{\Delta }x)^2>=\frac{1}{3}\left[\xi ^2\left(1erf\left(\frac{\xi }{2}\right)\right)+\frac{8}{\sqrt{\pi }\xi }\left(1\left(1+\frac{\xi ^2}{4}\right)\mathrm{exp}\left(\frac{\xi ^2}{4}\right)\right)\right].$$
(15)
The maximum in the mean squared and mean absolute displacements occur for $`\xi =2.61648`$ and $`\xi =1.76332`$, which corresponds to acceptance ratio values of $`P_{acc}=0.41767`$ and $`P_{acc}=0.558239`$ in agreement with the numerical results of BKS\[Bouzida et al., 1992\]: see Figs. 2 and 3.

## 3 D dimensions

We show here that the acceptance ratio, mean squared displacement and other quantities of interest can be obtained exactly in any dimension. Note that the term “volume” should be interpreted as the hypervolume in $`D`$ dimensions, e.g., area in 2D and volume in 3D. For simplicity, we derive exact expressions in odd dimensions, but similar results can be obtained in even dimensions.

### 3.1 Uniform volume sampling

In $`D`$ dimensions, the acceptance ratio is expressed as

$$P_{acc}(\delta )=\frac{c_D}{\delta ^DV_D}d^D𝐫\mathrm{exp}(\beta kr^2/2)_{|h|\delta }d^D𝐡\mathrm{min}(1,\mathrm{exp}(\beta k(|𝐫+𝐡|^2r^2)/2))$$
(16)
where $`V_D=\pi ^{D/2}/\mathrm{\Gamma }(D/2+1)`$ is volume of the sphere of unit radius in $`D`$ dimensions and

$`c_D`$ $`=\left(DV_D{\displaystyle _0^{\mathrm{}}}𝑑rr^{D1}e^{\beta kr^2/2}\right)^1=\left({\displaystyle \frac{\beta k}{2\pi }}\right)^{D/2}.`$ (17)
In odd dimensions, the derivative of Eq. (16) with respect to $`\delta `$ can be written explicitly by using generalized spherical coordinates

$`{\displaystyle \frac{d(\delta ^DP_{acc}(\delta ))}{d\delta }}`$ $`={\displaystyle \frac{c_D}{V_D}}{\displaystyle d^D𝐫\mathrm{exp}(\beta kr^2/2)}`$ (18)
$`\delta ^{D1}{\displaystyle 𝑑\mathrm{\Omega }\mathrm{min}(1,\mathrm{exp}(\beta k(r\delta \mathrm{cos}\varphi _1+\delta ^2/2)))},`$
where $`d\mathrm{\Omega }=(_{j=1}^{D2}(\mathrm{sin}(\varphi _j))^{D1j}d\varphi _j)d\varphi _{D1}`$ such that $`𝑑\mathrm{\Omega }=DV_D`$. The first $`D2`$ variables $`\varphi _j`$ are integrated from $`0`$ to $`\pi `$, whereas $`\varphi _{D1}`$ is integrated from $`0`$ to $`2\pi `$. If we denote $`u=\mathrm{cos}(\varphi _1)`$, perform integration over $`\varphi _2\mathrm{}\varphi _{D1}`$, and introduce the variable $`v=r\sqrt{\frac{\beta k}{2}}`$ Eq. (18) can be rewritten as

$`{\displaystyle \frac{d(\xi ^DP_{acc}(\xi ))}{d\xi }}`$ $`=D\xi ^{D1}{\displaystyle \frac{2}{\mathrm{\Gamma }\left(\frac{D}{2}\right)_1^1𝑑u(1u^2)^{(D3)/2}}}\times `$
$`\times `$ $`{\displaystyle _0^+\mathrm{}}v^{D1}𝑑ve^{v^2}{\displaystyle _1^1}𝑑u(1u^2)^{(D3)/2}\mathrm{min}(1,\mathrm{exp}((2\xi vu+\xi ^2)))`$ (19)
where $`\mathrm{min}(1,\mathrm{exp}((2\xi vu+\xi ^2)))=\mathrm{exp}(((2\xi vu+\xi ^2))`$ for $`v<\xi /2`$ with $`1<u<1`$ and for $`v>\xi /2`$ with $`\xi /(2v)<u<1`$ and $`\mathrm{min}(1,\mathrm{exp}((2\xi vu+\xi ^2)))=1`$ for $`r>\xi /2`$ with $`1<u<\xi /(2v)`$. Using that $`_1^1𝑑u(1u^2)^{(D3)/2}=\frac{\mathrm{\Gamma }\left(\frac{D1}{2}\right)}{\mathrm{\Gamma }\left(\frac{D}{2}\right)}\sqrt{\pi }`$, Eq. (3.1) then becomes

$`{\displaystyle \frac{d(\xi ^DP_{acc}(\xi ))}{d\xi }}=`$ $`{\displaystyle \frac{2D\xi ^{D1}}{\sqrt{\pi }\mathrm{\Gamma }\left(\frac{(D1)}{2}\right)}}({\displaystyle _0^{\xi /2}}dvv^{D1}e^{v^2}{\displaystyle _1^1}du(1u^2)^{(D3)/2}\mathrm{exp}((\xi ^2+2\xi vu))`$
$`+`$ $`{\displaystyle _{\xi /2}^+\mathrm{}}v^{D1}dve^{v^2}[{\displaystyle _1^{\xi /(2v)}}du(1u^2)^{(D3)/2}+`$
$`+`$ $`{\displaystyle _{\xi /(2v)}^1}du(1u^2)^{1/2}\mathrm{exp}(((\xi ^2+2\xi vu)))]).`$ (20)
After some calculation (see Appendix A) one obtains,

$`{\displaystyle \frac{d(\xi ^DP_{acc}(\xi ))}{d\xi }}`$ $`=D\xi ^{D1}\left(1erf\left({\displaystyle \frac{\xi }{2}}\right)\right)`$ (21)
which gives, for instance, in three dimensions

$$P_{acc}(\xi )=1erf\left(\frac{\xi }{2}\right)+\frac{8}{\sqrt{\pi }\xi ^3}\left(1\left(1+\frac{\xi ^2}{4}\right)\mathrm{exp}\left(\frac{\xi ^2}{4}\right)\right)$$
(22)
Figure 1 shows the acceptance ratio $`P_{acc}`$ versus $`\delta `$ in one, two and three dimensions.

A similar calculation for generating function $`Z_D(\lambda ,\xi )`$ leads to

$$\frac{(\xi ^DZ_D(\stackrel{~}{\lambda },\xi ))}{\xi }=D\xi ^{D1}e^{\stackrel{~}{\lambda }\xi }\left(1erf\left(\frac{\xi }{2}\right)\right),$$
(23)
and to

$$Z_D(\stackrel{~}{\lambda },\xi )=D(\xi )^{1D}\frac{^{D1}}{\stackrel{~}{\lambda }^{D1}}Z_{D=1}(\stackrel{~}{\lambda },\xi ),$$
(24)
where the expression for $`Z_{D=1}(\stackrel{~}{\lambda },\xi )`$ is given in Eq. (13). Since $`P_{acc}(\xi )=Z_D(\stackrel{~}{\lambda }=0,\xi )`$, it follows from Eq. (24) that the acceptance ratio in $`D`$ dimensions is equal, up a factor $`D(\xi )^{1D}\sqrt{\frac{\beta k}{2}}`$, to the mean squared displacement $`<(\mathrm{\Delta }x)^2>`$ in one dimension: compare Eq. (22) to Eq. (15).

Although straightforward, the algebra rapidly becomes tedious, and we only illustrate the results by giving the expression of the mean squared displacement $`<(\mathrm{\Delta }r)^2>`$ in three dimensions

$`{\displaystyle \frac{\beta k}{2}}<(\mathrm{\Delta }r)^2>`$ $`={\displaystyle \frac{^2}{\stackrel{~}{\lambda }^2}}Z_{D=3}(\stackrel{~}{\lambda },\xi )|_{\stackrel{~}{\lambda }=0}`$ (25)
$`={\displaystyle \frac{3}{\xi ^2}}\left({\displaystyle \frac{\beta k}{2}}\right)^2<(\mathrm{\Delta }x)^4>_{D=1}`$
$`={\displaystyle \frac{12}{5}}\left[{\displaystyle \frac{\xi ^2}{4}}\left(1erf\left({\displaystyle \frac{\xi }{2}}\right)\right)+{\displaystyle \frac{16}{\sqrt{\pi }\xi ^2}}\left(1\left(1+{\displaystyle \frac{\xi ^2}{4}}+{\displaystyle \frac{\xi ^4}{32}}\right)e^{\left(\frac{\xi ^2}{4}\right)}\right)\right].`$
The mean squared and mean absolute displacement are plotted versus the acceptance ratio $`P_{acc}`$ in one, two and three dimensions in Figs. 2 and 3. Note that the maximum is shifted to the left, i.e., to the smallest values of the acceptance ratio, when the space dimension increases.

### 3.2 Uniform radius sampling

The acceptance ratio $`P_{acc,w}`$ in $`D`$ dimensions can be expressed as

$`P_{acc,w}(\delta )`$ $`=c_D{\displaystyle _D}d^D𝐫\mathrm{exp}(\beta k𝐫^2/2)`$ (26)
$`{\displaystyle _{|h|\delta }}d^D𝐡P_w(h)\mathrm{min}(1,\mathrm{exp}(\beta k(|𝐫+𝐡|^2r^2)/2))`$
where $`P_w(h)`$ is the weighted probability. For a uniform distribution in radius, $`h^{D1}P_w(h)=(DV_D\delta )^1`$. Using the method developed in the above section, it is straightforward to obtain that

$$\frac{d(\xi P_{acc,w}(\delta ))}{d\xi }=\left(1erf\left(\frac{\xi }{2}\right)\right)$$
(27)
which shows that the acceptance ratio is the same whatever the dimension and explains the data collapse observed in \[Bouzida et al., 1992\]. Similarly, the generating function $`Z(\lambda ,\xi )`$ can be shown to obey to the differential equation

$$\frac{Z(\lambda ,\xi )}{\xi }=\mathrm{exp}(\sqrt{\frac{2}{\beta k}}\lambda \xi )\left(1erf\left(\frac{\xi }{2}\right)\right),$$
(28)
independently of the dimension $`D`$. This proves that $`Z(\lambda ,\xi )`$ and all moments such as $`<|\mathrm{\Delta }r|>`$ and $`<(\mathrm{\Delta }r)^2>`$ are independent of dimension, as numerically found by BKS\[Bouzida et al., 1992\].

## 4 Dynamic behavior

In addition to the exact results for the static properties presented above, we have investigated the dynamic behavior of the SHO by using numerical and analytical approaches. For simplicity, we discuss only the unidimensional case, but the approach can be generalized to $`D`$ dimensions.

The master equation describing the dynamical evolution of the system during the Monte Carlo simulation, Eq. (1), can be formally written

$$\frac{d𝐏(t)}{dt}=L𝐏(t)$$
(29)
where $`L`$ is a linear operator acting on $`P`$ and the Metropolis rule, Eq (4), is used for the transition rate. We consider the spectrum of eigenvalues $`\lambda `$ of $`L`$. Denoting $`P_\lambda (x)`$ the eigenfunction associated with $`\lambda `$ and introducing $`f_\lambda (x)`$ via $`P_\lambda (x)=P_{eq}(x)f_\lambda (x)`$, where $`P_{eq}(x)`$ is given in Eq. (3), one can express the eigenvalue equation

$$\lambda P_{eq}(x)f_\lambda (x)=L(P_{eq}(x)f_\lambda (x))$$
(30)
as

$$\lambda f_\lambda (x)=\frac{1}{2\delta }_\delta ^\delta 𝑑hMin(1,e^{\beta (V(x+h)V(x))})(f_\lambda (x)f_\lambda (x+h)).$$
(31)
Multiplying both sides by $`\mathrm{exp}(\beta V(x))f_\lambda ^{}(x)`$, where the star denotes a complex conjugate, and integrating over $`x`$ then gives

$$\lambda =\frac{1}{4\delta }\frac{_{\mathrm{}}^+\mathrm{}𝑑x_\delta ^\delta 𝑑hMin(e^{\beta V(x)},e^{\beta V(x+h)})|f_\lambda (x+h)f_\lambda (x)|^2}{_{\mathrm{}}^+\mathrm{}𝑑xe^{\beta V(x)}|f_\lambda (x)|^2}.$$
(32)
As anticipated for a Markov process satisfying detailed balance, one deduces from the above formula that all eigenvalues are real and positive; the smallest eigenvalue is $`\lambda _0=0`$ and it is associated with $`f_0(x)=constant0`$. One need consider only real eigenfunctions. Moreover, the eigenvalues can be sorted according to the symmetry of the associated eigenfunctions: it is easy to check that the eigenfunctions are either even or odd functions of $`x`$, due to the fact that the potential $`V(x)`$ is an even function of $`x`$.

Any solution of the master equation can be expanded as

$$P(x,t)=P_{eq}(x)(1+\underset{\lambda >0}{}c_\lambda f_\lambda (x)e^{\lambda t}),$$
(33)
and a similar expansion applies to the conditional probability $`P(x,t|x_0,0)`$ from which one can compute any time-dependent correlation function. The long-time kinetics governing the approach to equilibrium in $`P(x,t)`$ and in any correlation function is characterized by the smallest non-zero eigenvalue for which the amplitude, i.e., the projection of $`P(x,0)`$, or of the dynamic observable, onto the relevant eigenfunction, does not vanish.

Since, according to Eq. (32), the eigenvalues are expressed as the ratio of two positive quadratic functionals (that in the denominator being also definite), one can use the Rayleigh-Ritz procedure to find a variational upper-bound for the eigenvalues\[Dettman, 1962, Arfken, 1985\]. Consider first the smallest non-zero eigenvalue $`\lambda _1`$. For any real function $`\varphi (x)`$ which is both normalized and orthogonal to $`f_0(x)`$, i.e., satisfies

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}P_{eq}(x)\varphi (x)^2𝑑x`$ $`=1`$ (34)
$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}P_{eq}(x)\varphi (x)𝑑x`$ $`=0,`$ (35)
one has the inequality

$$\lambda _1\lambda _1[\varphi ]=\frac{1}{4\delta }_{\mathrm{}}^+\mathrm{}𝑑x_\delta ^\delta 𝑑hMin(e^{\beta V(x)},e^{\beta V(x+h)})|\varphi (x+h)\varphi (x)|^2.$$
(36)
A convenient choice of trial functions is provided by linear combinations of Hermite polynomials $`H_n(\xi )`$ (where, as in the previous sections $`\xi =\sqrt{\frac{\beta k}{2}}\delta `$), with $`n>0`$, since these form, up to a trivial multiplicative factors, an orthonormal basis with respect to the weight function $`\mathrm{exp}(\beta V(x))`$ with $`V(x)=(1/2)kx^2`$. For $`\lambda _1`$, which is associated with an odd eigenfunction, one need consider only the odd polynomials $`H_{2n+1}(\xi )`$, $`n0`$.

The simplest estimate of $`\lambda _1`$ is provided by taking

$$\varphi (\xi )=\frac{H_1(\xi )}{\sqrt{2\sqrt{\pi }}}=\sqrt{\frac{2}{\sqrt{\pi }}}\xi ,$$
(37)
which gives

$$\lambda _1[\varphi ]=\frac{4}{3}\left(\frac{\xi ^2}{4}(1erf\left(\frac{\xi }{2}\right))\right)+\frac{2}{\sqrt{\pi }\xi }(1(1+\frac{\xi ^2}{4}\mathrm{exp}\left(\frac{\xi ^2}{4}\right))$$
(38)
With this choice of $`\varphi (\xi )`$, $`\lambda _1[\varphi ]`$ simply reduces to the mean squared displacement $`<(\mathrm{\Delta }x)^2>`$ multiplied by $`\left(\frac{\beta k}{2}\right)`$ (see Eq. (13)). One then derives from section 2 that $`\lambda _1[\varphi ]`$ versus $`\xi `$ passes though a maximum for $`\xi 2.611648`$, which corresponds to an acceptance ratio of $`P_{acc}=0.41767`$.

A better estimate of $`\lambda _1`$ can be obtained by using a linear combination of $`H_1(\xi )`$ and $`H_3(\xi )`$:

$$\varphi (\xi ;\theta )=\left[\frac{\mathrm{cos}(\theta )}{\sqrt{2\sqrt{\pi }}}H_1(\xi )+\frac{\mathrm{sin}(\theta )}{\sqrt{48\sqrt{\pi }}}H_3(\xi )\right]$$
(39)
where only one independent parameter $`\theta `$ appears due to the normalization condition. The best bound is determined by minimizing the expression $`\lambda _1[\varphi ]`$ with respect to $`\theta `$:

$$\frac{\lambda _1[\varphi ]}{\theta }=0.$$
(40)
The result is a lengthly algebraic formula that is plotted in Figs. 4 and 5, together with the expression in Eq. (38).

An improved estimate of $`\lambda _1`$ can be derived by noting that at large $`\xi `$, $`\lambda _1`$ is inversely proportional to $`\xi `$. Actually, one can show that this is true for all eigenvalues except $`\lambda _0=0`$. By considering Eq. (31) in the limit where $`x`$ goes to zero, one arrives at the result (see Appendix B)

$$\lambda (\xi )\frac{\sqrt{\pi }}{2\xi }+0(e^{\xi ^2}),$$
(41)
valid for large $`\xi `$. Note that the correction terms are very small as soon as $`\xi 3`$. One can build an estimate of $`\lambda _1(\xi )`$ by using the piecewise function that is equal to $`\frac{\sqrt{\pi }}{2\xi }`$ for $`\xi \xi ^{}`$ and is equal to $`\lambda _1[\varphi ]`$ obtained for a linear combination of $`H_1`$ and $`H_3`$ (see above) for $`\xi \xi ^{}`$, where $`\xi ^{}`$ is the value at which $`\lambda _1[\varphi ]=\frac{\sqrt{\pi }}{2\xi }`$. $`\lambda (\xi )`$ is then maximum for $`\xi =\xi ^{}2.35`$; the corresponding value of the acceptance ratio is $`P_{acc}0.56`$. The estimate is shown in Figs. 4 and 5.

In order to compare our results to the BKS paper \[Bouzida et al., 1992\], it is necessary to calculate the second eigenvalue $`\lambda _2`$. The trial function is then chosen in a subspace orthogonal not only to $`f_0(\xi )=constant`$ but also to the eigenfunction associated with $`\lambda _1`$. A convenient choice is provided by (normalized) linear combinations of the even Hermite polynomials $`H_{2n}(\xi )`$ with $`n1`$. An estimate of $`\lambda _2`$ can be obtained by using a linear combination of $`H_2(\xi )`$ and $`H_4(\xi )`$:

$$\varphi (\xi ;\theta )=\left[\frac{\mathrm{cos}(\theta )}{2\sqrt{2\sqrt{\pi }}}H_2(\xi )+\frac{\mathrm{sin}(\theta )}{8\sqrt{6\sqrt{\pi }}}H_4(\xi )\right]$$
(42)
and the result is shown in Figs. 6 and 7, together with the zeroth-order approximation obtained with only $`H_2(\xi )`$ and the improved estimate taking into account the large-$`\xi `$ behavior.

In addition to the above analytical estimates, we have also performed a numerical study of the spectrum of eigenvalues of the master equation, Eq. (1). The latter has been discretized in $`x`$-space by taking a constant step size $`\mathrm{\Delta }`$, which leads to a matrix form,

$`{\displaystyle \frac{d𝐏(t)}{dt}}=\mathrm{𝐖𝐏}(t)`$ (43)
where $`𝐏`$ is a vector with components $`P_i(t)=P(x_i,t)`$ and the elements $`W_{ij}`$ of the matrix $`𝐖`$ are such that

$`W_{ii}`$ $`=\mathrm{\Delta }{\displaystyle \underset{j=iNh,ji}{\overset{i+Nh}{}}}W(ij)`$
$`W_{ij}`$ $`=\mathrm{\Delta }W(ji)`$ (44)
and $`W(ij)=Min(1,e^{\beta (V(x_j)V(x_i))})`$. The eigenvalues $`\lambda `$ can then be obtained via an exact numerical diagonalization of the matrix $`W`$. In practice, convergence is obtained for a unidimensional lattice of $`400`$ sites, where the $`x`$-range is $`[10,10]`$, and $`\xi `$ goes from $`0`$ to $`5`$. One checks that the lowest eigenvalue is equal to zero and corresponds to the equilibrium state, and that all other eigenvalues are real and strictly positive and behave as $`\sqrt{\pi }/(2\xi )`$ for large enough $`\xi `$. The results for the first non-zero eigenvalue $`\lambda _1`$ and $`\lambda _2`$ are displayed in Figs. 4 and 6, respectively. One can see that the best analytical estimate described above (linear combination of two Hermite polynomials plus exact asymptotic behavior at large $`\xi `$), is in excellent agreement with the numerical value in both cases.

The above analysis allows us to derive by analytical means the simulation result obtained by BKS\[Bouzida et al., 1992\] for the dependence of the characteristic time $`\tau `$ of the energy-energy correlation function on the acceptance ratio: approximating $`\tau `$ by $`1/\lambda _2`$ (since the energy $`V(x)`$ is an even function of $`x`$, its projection on the first eigenfunction associated with $`\lambda _1`$ vanishes), using for $`\lambda _2`$ our best analytical estimate, and combining this with the exact result for the acceptance ratio $`P_{acc}`$ in section 2 lead to the full curve plotted in Fig. 8; the time $`\tau `$ is minimum for the acceptance ratio close to $`0.47`$, as found by BKS\[Bouzida et al., 1992\] ($`0.50`$).

In Fig. 8, we have also plotted $`1/\lambda _1`$ versus the acceptance ratio: it is minimum for $`P_{acc}=0.45`$. Note that using the results of this and the preceding sections, one can rigorously show that the correlation time $`\tau `$, no matter how it is precisely defined, diverges as $`1/P_{acc}`$ when $`P_{acc}0`$ and as $`1/(1P_{acc})`$ when $`P_{acc}1`$.

## 5 Conclusion

We have obtained exact results concerning Metropolis Algorithms for the displacement of a particle in the simple harmonic potential. Our analysis provides a theoretical explanation of the numerical results obtained by BKS\[Bouzida et al., 1992\]. In particular, we show that the results become independent of the space dimension when the successive trial moves are sampled according to a Metropolis algorithm with a uniform distribution in radius (instead of volume). This rationalizes the search for efficient Monte Carlo methods for the simulation of systems with intrinsic inhomogeneity and anisotropy such as biological molecules\[Bouzida et al., 1992, Swendsen, 2002\]

## Appendix A Acceptance probability

The derivation of Eq. (21) can be done from Eq. (3.1) after some manipulations whose details are given here. Let us denote $`I_D`$ and $`J_D`$ as

$$I_D=_0^{\xi /2}𝑑vv^{D1}e^{v^2}_1^1𝑑u(1u^2)^{(D3)/2}\mathrm{exp}((\xi ^2+2\xi vu))$$
(45)
$`J_D`$ $`={\displaystyle _{\xi /2}^+\mathrm{}}dvv^{D1}e^{v^2}[{\displaystyle _1^{\xi /(2v)}}du(1u^2)^{(D3)/2}+`$ (46)
$`+{\displaystyle _{\xi /(2v)}^1}du(1u^2)^{(D3)/2}\mathrm{exp}(((\xi ^2+2\xi vu)))].`$
Changing the variable $`u`$ to $`y=vu`$ and $`v`$ to $`t=v^2y^2`$ leads to the following relations

$$I_D=_0^{\xi /2}𝑑y\frac{\mathrm{exp}((\xi +y)^2)+\mathrm{exp}((y\xi )^2)}{2}_0^{\xi ^2/4y^2}𝑑tt^{(D3)/2}e^t$$
(47)
$`J_D`$ $`={\displaystyle _{\xi /2}^+\mathrm{}}𝑑y{\displaystyle \frac{\mathrm{exp}(y^2)}{2}}{\displaystyle _0^+\mathrm{}}𝑑tt^{(D3)/2}e^t`$ (48)
$`+{\displaystyle _0^{\xi /2}}𝑑y{\displaystyle \frac{\mathrm{exp}((y\xi )^2)}{2}}{\displaystyle _{\xi ^2/4y^2}^+\mathrm{}}𝑑tt^{(D3)/2}e^t`$
$`+{\displaystyle _0^{\xi /2}}𝑑y{\displaystyle \frac{\mathrm{exp}((y+\xi )^2)}{2}}{\displaystyle _{\xi ^2/4y^2}^+\mathrm{}}𝑑tt^{(D3)/2}e^t`$
$`+{\displaystyle _{\xi /2}^+\mathrm{}}𝑑y{\displaystyle \frac{\mathrm{exp}((y+\xi )^2)}{2}}{\displaystyle _0^+\mathrm{}}𝑑tt^{(D3)/2}e^t`$
Using that $`_0^+\mathrm{}𝑑tt^{(D3)/2}e^t=\mathrm{\Gamma }((D1)/2)`$, one obtains that

$$I_D+J_D=\frac{\sqrt{\pi }}{2}\mathrm{\Gamma }\left(\frac{D1}{2}\right)\left(1erf\left(\frac{\xi }{2}\right)\right)$$
(49)
Inserting Eqs. (49) in Eq. (3.1) leads to Eq. (21).

## Appendix B Asymptotic behavior of the eigenvalues for large $`\delta `$

Consider Eq. (31) in the limit $`x0`$ ; one obtains after rearranging the various terms:

$`\lambda f_\lambda (x)`$ $`={\displaystyle \frac{1}{2\delta }}{\displaystyle _\delta ^\delta }𝑑h\mathrm{exp}({\displaystyle \frac{\beta k}{2}}((x+h)^2x^2))(f_\lambda (x)f_\lambda (x+h))`$ (50)
$`+{\displaystyle \frac{1}{2\delta }}{\displaystyle _{2x}^0}𝑑h(1\mathrm{exp}({\displaystyle \frac{\beta k}{2}}((x+h)^2x^2)))(f_\lambda (x)f_\lambda (x+h)).`$
The second term of the right-hand side of Eq. (50) is at most of order $`x^3|f_\lambda (x)|`$ and is always negligible so that one can rewrite Eq. (50) as

$$\left(\lambda \delta _0^\delta 𝑑he^{\frac{\beta k}{2}h^2}+O(x^2)\right)f_\lambda (x)\frac{e^{\frac{\beta k}{2}x^2}}{2}_\delta ^\delta 𝑑h\mathrm{exp}(\frac{\beta k}{2}(x+h)^2)f_\lambda (x+h).$$
(51)
Shifting the variable from $`h`$ to $`x+h`$ in the integral of the r.h.s of Eq. (51) and using the orthogonality of $`f_\lambda (x)`$ to $`f_0(x)=constant`$ for all non-zero eigenvalues $`\lambda `$, i.e., $`_{\mathrm{}}^+\mathrm{}𝑑x\mathrm{exp}(\frac{\beta k}{2}x^2)f_\lambda (x)=0`$, leads to the following expression,

$$\left(\lambda \delta \sqrt{\frac{\pi }{2\beta k}}+O(x^2)\right)f_\lambda (x)\frac{1+O(x^2)}{2}\left[_{x+\delta }^+\mathrm{}𝑑he^{\frac{\beta kh^2}{2}}f_\lambda (h)+_{\mathrm{}}^{x\delta }𝑑he^{\frac{\beta kh^2}{2}}f_\lambda (h)\right]$$
(52)
for any non-zero $`\lambda `$. The r.h.s. of Eq. (52) can be Taylor expanded, which gives

$`\left(\lambda \delta \sqrt{{\displaystyle \frac{\pi }{2\beta k}}}+O(x^2)\right)f_\lambda (x)`$ $`{\displaystyle \frac{1}{2}}[{\displaystyle _\delta ^+\mathrm{}}dhe^{\frac{\beta kh^2}{2}}(f_\lambda (h)+f_\lambda (h))`$ (53)
$`xe^{\frac{\beta kh^2}{2}}(f_\lambda (h)f_\lambda (h))+0(x^2)].`$
If the eigenfunctions $`f_\lambda `$ is an even function of $`x`$, one then derives that $`f_\lambda (x)=f_\lambda (0)+O(x^2)`$ with $`f_\lambda (0)0`$ and

$$\left(\lambda \delta \sqrt{\frac{\pi }{2\beta k}}\right)=_\delta ^+\mathrm{}𝑑h\mathrm{exp}(\frac{\beta kh^2}{2})\frac{f_\lambda (h)}{f_\lambda (0)},$$
(54)
which after introducing $`\xi =\sqrt{\frac{\beta k}{2}}\delta `$ can be rewritten as

$$\lambda =\frac{\sqrt{\pi }}{2\xi }+\sqrt{\frac{2}{\beta k}}_\xi ^+\mathrm{}𝑑h\mathrm{exp}(\frac{\beta kh^2}{2})\frac{f_\lambda (\sqrt{\frac{2}{\beta k}}h)}{f_\lambda (0)}.$$
(55)
If the eigenfunction is an odd function of $`x`$, one has that $`f_\lambda (x)=f_\lambda ^{}(0)x(1+O(x^2)`$ with $`f_\lambda ^{}(0)0`$ and

$$\lambda =\frac{\sqrt{\pi }}{2\xi }\frac{f_\lambda (\sqrt{\frac{2}{\beta k}}\xi )}{f_\lambda (0)}\mathrm{exp}(\xi ^2).$$
(56)
From Eqs. (55) and (56), one immediately obtains that all non-zero eigenvalues behave as

$$\lambda \frac{\sqrt{\pi }}{2\xi }+O(\mathrm{exp}(\xi ^2))$$
(57)
when $`\xi +\mathrm{}`$, since $`f_\lambda (x)`$ diverges more slowly than $`e^{x^2}`$ when $`x+\mathrm{}`$.

## References