text:
# Introduction

## Introduction

TeV scale Supersymmetry (SUSY) is theoretically the most attractive way of stabilizing the mass of the Higgs boson against radiative corrections, around the electroweak symmetry (EW) breaking scale. Hence, SUSY searches, along with those for the Higgs boson, form the focus of physics studies at the current and future colliders. In the simplest version of SUSY models one assumes conservation of a quantum number $`R_p(1)^{3B+L+2S}`$, where $`S`$ stands for the spin of the particle. The conservation of this quantum number is guaranteed only if the superpotential is invariant under the discrete $`R_p`$ symmetry. However, there really is no deep theoretical reason for this symmetry. As a matter of fact, it is possible to have $`R_p`$-violating terms in the superpotential while respecting supersymmetry as well as invariance under the SM gauge transformation. These are given by

$$W_\mathit{}_p=\frac{1}{2}\lambda _{ijk}L_iL_jE_k^c+\lambda _{ijk}^{}L_iQ_jD_k^c+\frac{1}{2}\lambda _{ijk}^{\prime \prime }U_i^cD_j^cD_k^c+\kappa _iL_iH_2.$$
(1)
Here $`L_i`$,$`Q_i`$ are the doublet lepton, quark superfields, and $`E_i,U_i,D_i`$ the singlet lepton, quark superfields, $`i`$ being the generation index. The above necessarily mean non-conservation of the baryon number $`B`$ and/or the lepton number $`L`$. The $`\mathit{}`$ couplings allow the generation of $`\nu `$ masses in an economical way without introducing any new fields. The contribution from the trilinear $`\lambda ,\lambda ^{}`$ comes at the one- or two-loop level and the bilinear $`\kappa _i`$ generate it at tree level . The Kamioka and SNO experiments now give an unambiguious evidence for $`\nu `$ masses. In SUSY models with $`R/_p`$, there exists enough freedom to generate observed mass patterns with sparticle masses in the TeV range, and hence with clean predictions for collider searches. This interplay between indirect constraints, coming from $`\nu `$ physics and the collider signals of $`R/_p`$ SUSY is a very interesting aspect of the $`R/_p`$ SUSY studies. Of course the simultaneous existence of the $`\lambda ^{}`$ and $`\lambda ^{\prime \prime }`$ couplings and TeV scale masses for the sfermions will mean a very rapid proton decay. This can be cured by adopting $`B`$ conservation, i.e. $`\lambda ^{\prime \prime }=0`$. This choice is actually preferred if $`R/_p`$ terms are not to wash out baryogenesis. Unified string theories actually prefer models with $`B`$ conservation and $`R/_p`$ . These models treat the lepton and the quark fields differently and have two discrete symmetries. $`B`$ conservation and $`R/_p`$ eliminate not just the dimension-$`4`$ operators for proton decay but also the dimension-$`5`$ operators. Needless to say that all this makes $`R/_p`$ theoretically very interesting.

Of course the large number, 48 as per Eq. (1), of these essentially Yukawa-type couplings, with no theoretical indications about their sizes, is not a very satisfactory feature of $`R/_p`$ SUSY. Luckily many of these unknown couplings are constrained by low energy processes such as $`p/\mu `$ decay, $`\nu `$ masses or cosmological arguments like generating the right amount of baryon asymmetry in the early Universe. There exist a host of constraints on $`R/_p`$ couplings . Most of the constraints come from virtual effects caused by sparticles in loops and some from direct collider searches. Many of these get less severe with increasing number of the third-generation indices of the involved couplings. The constraints on the $`\mathit{}`$ couplings coming from $`\nu `$ masses can be particularly severe but may also depend on model assumptions. Hence a study of the same couplings in a collider environment, can certainly help to clarify model building.

## $`R_p`$ violation at the colliders

Effects of $`R/_p`$ at colliders fall into qualitatively different classes, depending on the size of the $`R/_p`$ couplings. The one effect that is always present is an unstable lightest supersymmetric particle (LSP), which need not even be the $`\stackrel{~}{\chi }_1^0`$. If at least one of the $`R/_p`$, $`\lambda ,\lambda ^{}`$ couplings is $`>10^6`$, the LSP will decay in the detector. In our work we consider the case of these $`\mathit{},R/_p`$ couplings and take the $`\stackrel{~}{\chi }_1^0`$ to be the LSP. The decays of the LSP $`\stackrel{~}{\chi }_1^0f\overline{f}_1f_2`$ can give rise to strikingly different final states  as compared to the SUSY signal with $`R_p`$ conservation. For masses of $`\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_1^\pm `$ of interest at the LHC and NLC, even a $`t`$ in the final state will be allowed. For larger $`R/_p`$ couplings, decays of particles and sparticles other than the LSP via $`R/_p`$ interactions  can take place. For example, $`\stackrel{~}{t}bl,tb\stackrel{~}{l}`$. These then may compete with the $`R_p`$-conserving decays. The $`R/_p`$ couplings can also give contributions to tree level processes due to virtual sparticle exchange, such as $`pp\mu ^+\mu ^{}`$  or $`pp(p\overline{p})t\overline{t}X`$. In the last case the interference terms give rise to non-vanishing $`t`$ polarization , which can be probed by the decay lepton energy/angular distribution . In addition, for large enough $`R/_p`$ couplings, one also has resonant or non-resonant production of a single sparticle via the $`R/_p`$ couplings . For example, much as $`H^{}`$ produced via $`t\overline{b}H^{}`$ coupling, one can also have $`ppt\overline{b}\stackrel{~}{\tau }`$ via $`\lambda _{333}^{}`$ as shown in

Figure 1. A few points are worth noting. The production cross-section is appreciable even, for $`m_{\stackrel{~}{t}}>m_t`$ and for $`\lambda _{333}^{}`$ as small as $`0.05`$. $`\stackrel{~}{\tau }`$ so produced will have both $`R/_p`$ and $`R_p`$-conserving decays. The latter will produce $`\stackrel{~}{\chi }_1^0`$, which will further decay through some $`R/_p`$ coupling. Some of these decays can fake the charged Higgs signal. This can be seen in a little more detail as follows. For large $`\lambda _{333}^{},\stackrel{~}{\tau }^{}`$s can have the $`R_p`$ conserving decays, $`\stackrel{~}{\tau }\tau \stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{\tau }\nu _\tau \stackrel{~}{\chi }_1^{}`$ for $`m_t>m_{\stackrel{~}{\tau }}`$, whereas for $`m_t<m_{\stackrel{~}{\tau }}`$, it can also have the $`R/_p`$ decay $`\stackrel{~}{\tau }b\overline{t}`$. The net final state produced will then be decided by relative branching ratios of $`\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_1^+`$ into different channels. Production of $`\stackrel{~}{\tau }`$ through $`R/_p`$ couplings and its decay via the same will give rise to $`ppt\overline{b}\stackrel{~}{\tau }Xt\overline{b}t\overline{b}X`$. Note that this is the same final state as the $`H^\pm `$. A comprehensive study of the $`H^+`$ signal at the LHC, in this case, will thus require full analysis of the three-body decays of the $`\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_1^\pm `$. The $`R_p`$-conserving decays of the $`\stackrel{~}{\tau }`$ followed by the $`R/_p`$ decays of the $`\stackrel{~}{\chi }_1^0`$ from the $`\stackrel{~}{\tau }`$ decay can also produce

$`ppt\overline{b}\stackrel{~}{\tau }X`$ $`(2t)(2b)(2\tau )X`$ (2)
$`tb(2\overline{b})\tau \nu _\tau X,`$
etc. These will then give rise to the characteristic $`\mathit{}`$ signal of like-sign fermion pairs.

As mentioned already, of particular importance are the possibilities of testing the low energy constraints, in ‘direct’ signals at the colliders. Since the third generation sfermions are expected to be lighter, they may give rise to larger virtual effects. For $`\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_1^+`$ with masses of interest at the LHC, NLC, final states with third-generation fermions including $`t`$ are possible. One needs a study of the $`R/_p`$ decays of $`\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_1^\pm `$ retaining effects of the mass of the third-generation fermions, for $`\mathit{}`$ coupling.

We therefore calculated three-body $`R/_p`$ decays of the $`\stackrel{~}{\chi }_1^0`$ with dominant $`\lambda ,\lambda ^{}`$ couplings involving at least two third-generation indices. We kept the mass effects of the third-generation fermions and the sfermion $`L`$$`R`$ mixing terms complex and analysed numerically for cases with unified/ununified gaugino masses, including the effect of subdominant $`\lambda ,\lambda ^{}`$.

## Three-body decays of the $`\stackrel{~}{\chi }_1^0`$

In this section we present the results for the above-mentioned case of dominant $`\lambda _{333}^{}`$ coupling. The calculation and the results can be trivially generalized  to cases of other couplings, dominant or subdominant.

Figure 2 shows the diagrams contributing to the three-body decay $`\stackrel{~}{\chi }_1^0b\overline{b}\nu _\tau `$. If $`m_{\stackrel{~}{\chi }_1^0}>m_t`$, then the decay into a massive final state containing $`t`$, viz. $`\stackrel{~}{\chi }_1^0\overline{b}t\nu _\tau `$ \+ C.C., is also possible. Our expressions for the decay widths  are the most general in that they include the effect of the finite fermion masses and are applicable to the case of complex $`L`$$`R`$ mass term. In the limit of real $`L`$$`R`$ sfermion mass terms, these agree with the earlier calculations  with finite fermion mass and, in the limit of zero fermion mass, they agree with those of Refs.  up to the corrections pointed out in Ref. .

The decay widths clearly depend on the $`L`$$`R`$ mixing in the sfermion sector as well as on the composition of the $`\stackrel{~}{\chi }_1^0`$. For example, for a Wino-like $`\stackrel{~}{\chi }_1^0`$, in the absence of $`L`$$`R`$ mixing, the second diagram with a $`\stackrel{~}{b}_R`$ exchange will not contribute. A substantial gaugino/higgsino mixing causes the width into the massive mode to be large at low $`\mathrm{tan}\beta `$, whereas for large $`\mathrm{tan}\beta `$ it is the massless mode that gets enhanced, compared to the corresponding widths for unmixed cases.

As an illustration, the results for the $`\stackrel{~}{\chi }_1^0`$ widths into massive and massless modes, for $`\lambda _{333}^{}=1`$ and $`M_{\stackrel{~}{\chi _1^0}}=600`$ GeV, are presented in

Figure 3 as a function of $`m_{}`$, the squark mass scale. We use two different choices of $`L`$$`R`$ mixing in the squark sector and different compositions of the $`\stackrel{~}{\chi }_1^0`$. The massive decay mode has a large width for low $`\mathrm{tan}\beta `$ and large higgsino–gaugino mixing. Otherwise the massless mode is larger, although the massive mode is non-negligible. From the left panel we see that, for a $`\stackrel{~}{\chi _1^0}\stackrel{~}{W}`$, the widths are enhanced by an order of magnitude with respect to the case of a $`\stackrel{~}{\chi }_1^0\stackrel{~}{B}`$. The right panel shows that with moderate $`L`$$`R`$ squark mixing, even for a $`\stackrel{~}{\chi }_1^0\stackrel{~}{W},`$ the $`\stackrel{~}{b}`$-mediated diagram contributes to the massive mode. Further the figures also show that smaller $`\stackrel{~}{t_1},\stackrel{~}{b_1}`$, masses enhance the width for the massive mode for a $`\stackrel{~}{B}`$-like $`\stackrel{~}{\chi }_1^0`$. For the same reasons the widths into the massive and massless final state also depend on the value of $`\mathrm{tan}\beta `$.

Figure 4 shows the $`b\overline{b}\nu _\tau ,t\overline{b}\nu `$ decay widths, as a function of the squark mass scale $`m_{}`$. The massive mode dominates at lower $`\mathrm{tan}\beta `$ and the massless one for the higher values. Thus, to summarize, we find that with $`\lambda _{333}^{}`$ dominant, the massive decay has a large width for low $`\mathrm{tan}\beta `$ and large higgsino–gaugino mixing. Otherwise the width for the massless mode is larger, though that for massive mode is also non-negligible. All the plots have been shown for $`\lambda _{333}^{}=1`$. In the resonant sfermion region one finds very little dependence of the width on $`\lambda _{333}^{}`$ and otherwise the widths for other values of the couplings simply scale down like $`\lambda _{333}^{}_{}{}^{}2`$.

The analysis can be easily generalized to the case where there are more than one large $`R/_p`$ couplings, and the results are shown in Figure 5. We notice that the widths into different final states are not too dissimilar, and more importantly the branching ratio into a final state with $`t`$ is not too small.

## Conclusions

It is important to have a collider probe of the $`R/_p\lambda `$ and $`\lambda ^{}`$ couplings with more than one third-generation index. Further, single charged slepton production via the $`R/_p`$ coupling $`\lambda _{i33}^{}`$ can produce final states, which can be confused with the $`H^\pm `$ signal. Even for smaller but dominant $`R/_p\lambda `$ and $`\lambda ^{}`$ couplings with more than one third-generation index, the three-body $`R/_p`$ decays of $`\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_1^+`$ can have important phenomenological consequences for new particle searches at the future colliders. If the $`\stackrel{~}{\chi }_1^0`$ is ‘Wino’-like, its decay widths increase by over an order of magnitude. If the $`\stackrel{~}{\chi }_1^0`$ is lighter than the $`t`$, its $`R/_p`$ decays will be dominantly into massless fermions: $`b\overline{b}\nu _\tau `$ for dominant $`\lambda _{333}^{}`$ coupling, $`c\overline{b}\tau ,s\overline{b}\nu _\tau `$ and the charge conjugate modes $`\overline{c}b\overline{\tau },\overline{s}b\overline{\nu }_\tau `$ for dominant $`\lambda _{323}^{}`$, and for $`\lambda _{233}`$ into $`\overline{\mu }\tau \nu _\tau `$ and $`\mu \overline{\tau }\overline{\nu }_\tau `$. For a $`\stackrel{~}{\chi }_1^0`$ heavier than the $`t`$, the massive decay modes can become competitive for large $`L`$$`R`$ sfermion mass term and/or for substantial mixing in the higgsino/gaugino sector, at not too large $`\mathrm{tan}\beta `$. All the $`R/_p`$ decays of the $`\stackrel{~}{\chi }_1^0`$ produce final states with like-sign dileptons as the tell-tale signature of the $`\mathit{}`$, as long as the final state $`t`$ decays hadronically.

Acknowledgements
I wish to thank my collaborators F. Borzumati, J.L. Kneur and F. Takayama for a nice collaboration, which resulted in the publication on which this talk is based.


format: text

text:
# hep-th/0210199 Representation of Small Conformal Algebra in 𝜅-basis

## 1 Introduction

The basic ingredients in the construction of the covariant string field theory are Witten’s star product  and BRST operator. Recent progress in the diagonalization  of the Neumann matrices $`M^{rs}`$ defining star product allows one to identify this product with the continuous Moyal product . We have

$$\underset{n=1}{\overset{\mathrm{}}{}}M_{mn}^{rs}v_n^{(\kappa )}=\mu ^{rs}(\kappa )v_m^{(\kappa )}$$

where $`\mathrm{}<\kappa <\mathrm{}`$, the eigenvalues $`\mu ^{rs}(\kappa )`$ are

$$\mu ^{rs}(\kappa )=\frac{1}{1+2\mathrm{cosh}\frac{\pi \kappa }{2}}\left[12\delta _{r,s}+e^{\frac{\pi \kappa }{2}}\delta _{r+1,s}+e^{\frac{\pi \kappa }{2}}\delta _{r,s+1}\right].$$

The eigenvectors $`v_n^{(\kappa )}`$ are given by their generating function

$$f^{(\kappa )}(z)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{v_n^{(\kappa )}}{\sqrt{n}}z^n=\frac{1}{\kappa \sqrt{𝒩(\kappa )}}(1e^{\kappa \mathrm{tan}^1z})$$
(1.1)
where

$$𝒩(\kappa )=\frac{2}{\kappa }\mathrm{sinh}\left(\frac{\pi \kappa }{2}\right).$$

In spite of the fact that we understand the nature of the star product, it is still unclear how BRST operator acts in the basis in which star product is simple. Moreover we do not know a representation of Virasoro generators in this basis. The authors of  tried to construct generator $`L_0`$ us a kernel operator in the $`\kappa `$-basis. They show that this kernel is not defined as a distribution. In this paper I show that small conformal algebra $`\{L_1,L_0,L_1\}`$ is represented by a certain difference operators in the $`\kappa `$-basis.

Let me briefly remind a construction of the Fock space at hand . In the discrete basis we have creation and annihilation operators $`a_n^{}`$ and $`a_n`$, $`n=1,2,\mathrm{}`$. Introduce the corresponding operators in the continuous basis via

$$a_\kappa ^{}=\underset{n=1}{\overset{\mathrm{}}{}}v_n^{(\kappa )}a_n^{}\text{and}a_\kappa =\underset{n=1}{\overset{\mathrm{}}{}}v_n^{(\kappa )}a_n.$$
(1.2)
Let $`f_1(\kappa ),\mathrm{},f_m(\kappa )`$ be functions from Schwartz space, then the states

$$a^{}(f_1)\mathrm{}a^{}(f_m)|0,\text{where}a^{}(f)=_{\mathrm{}}^{\mathrm{}}𝑑\kappa f(\kappa )a_\kappa ^{}$$
(1.3)
form a dense subset in this Fock space. Since the operator

$$L_0=\underset{n=1}{\overset{\mathrm{}}{}}na_n^{}a_n$$
(1.4)
contains only one creation and one annihilation operator it is completely determined by specifying its action on one particle states. Therefore we can translate the action of operator $`L_0`$ on states to its action on functions from the Schwartz space via

$$[L_0,a^{}(f)]=a^{}(L_0[f]).$$
(1.5)
In the next section I show that

$$L_0[f](\kappa )=\frac{1}{4}\left[\sqrt{\kappa (\kappa +i2^{})}f(\kappa +i2^{})+\sqrt{\kappa (\kappa i2^{})}f(\kappa i2^{})\right],$$
(1.6)
where $`2^{}=20`$. The operator (1.6) is defined on a class of holomorphic functions in the strip $`2<\mathrm{}\kappa <2`$. These functions may have poles on the boundary, and the term $`i2^{}`$ shows how one approaches these poles. I.e. one should approach the pole on the line $`\kappa +2i`$ from the bottom and the pole on the line $`\kappa 2i`$ from the top. If it happens that the function $`f`$ is holomorphic in strip $`4<\mathrm{}\kappa <4`$ one can apply operator (1.6) once again, etc.

Notice that the operators $`L_1`$ and $`L_1`$ also contain only one creation and one annihilation operator, and therefore they are completely determined by specifying their action on one particle subspace. We can write

$$[L_{\pm 1},a^{}(f)]=a^{}(L_{\pm 1}[f]),$$
(1.7)
where

$`L_1[f]`$ $`={\displaystyle \frac{\kappa }{2}}f(\kappa )+{\displaystyle \frac{i}{4}}\left[\sqrt{\kappa (\kappa +i2^{})}f(\kappa +i2^{})\sqrt{\kappa (\kappa i2^{})}f(\kappa i2^{})\right]`$ (1.8a)
$`L_1[g]`$ $`={\displaystyle \frac{\kappa }{2}}f(\kappa ){\displaystyle \frac{i}{4}}\left[\sqrt{\kappa (\kappa +i2^{})}f(\kappa +i2^{})\sqrt{\kappa (\kappa i2^{})}f(\kappa i2^{})\right]`$ (1.8b)
The paper is organized as follows. In Sections 2 and 3 we give a derivation of the equations (1.6) and (1.8) correspondingly. In Section 4 we check that the generators (1.6) and (1.8) defining in a class of difference operators have indeed correct commutation relations. The Appendix contains the technical information.

## 2 Operator $`L_0`$

Using the fact that

$$[L_0,a_n^{}]=na_n^{}.$$
(2.1)
one can easily obtain that $`L_0`$ acts on the generating function (1.1) (it should be considered as function of $`\kappa `$; $`z`$ is an external parameter) in the following way

$$L_0[\widehat{f}^{(\kappa )}]=z\frac{d}{dz}\widehat{f}^{(\kappa )}(z).$$
(2.2)
From this one can easily obtain a formal expression for the kernel of the operator $`L_0`$

$$L_0(\kappa ,\kappa ^{})\underset{n=1}{\overset{\mathrm{}}{}}nv_n^{(\kappa )}v_n^{(\kappa ^{})}.$$
(2.3)
It was argued in  that this kernel is not defined as a distribution. This means that the operator $`L_0`$ is not a kernel operator in the $`\kappa `$-basis.

It was explained in the Introduction that operator $`L_0`$ is completely determined by its restriction onto the one-particle Fock space. To actually find it we will use the following strategy. First, find the inverse of the operator $`L_0`$ on one particle Fock subspace. In other words we are looking for an operator $`G_0`$ such that

$$G_0L_0|\psi =L_0G_0|\psi =|\psi $$

for any state $`|\psi =a^{}(f)|0`$ specifying by function $`f`$ from Schwartz space. Using relation (1.5) one can formulate this statement in the Schwartz space as

$$G_0[L_0[f]](\kappa )=L_0[G_0[f]](\kappa )=f(\kappa )$$
(2.4)
for any Schwartz function $`f(\kappa )`$. From (2.3) one can easily obtain a formal expression for the kernel of operator $`G_0`$

$$G_0(\kappa ,\kappa ^{})=\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n}v_n^{(\kappa )}v_n^{(\kappa ^{})}.$$
(2.5)
Straight forward calculations of this kernel, which I present in Appendix A, show that operator $`G_0`$ is indeed a kernel operator and

$$G_0(\kappa ,\kappa ^{})=\left[\frac{\theta (\kappa )}{\kappa }\right]^{1/2}\left[\frac{\theta (\kappa ^{})}{\kappa ^{}}\right]^{1/2}\frac{1}{4\mathrm{cosh}\left[\frac{\pi }{4}(\kappa \kappa ^{})\right]}.$$
(2.6)
Here $`\theta (\kappa )=2\mathrm{tanh}\frac{\pi \kappa }{4}`$ is a non-commutativity parameter specifying the continuous Moyal algebra .

Second, we find the inverse of operator $`G_0`$ on the one-particle Fock space. In other words we are going to solve equation (2.4). The resulting operator will coincide with the operator $`L_0`$ on the one-particle Fock space. But $`L_0`$ is completely determined by its action on this subspace and therefore we will actually find operator $`L_0`$ on the whole Fock space.

So we need to solve the equation $`G_0[L_0[f]]=f`$:

$$_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}\left[\frac{\theta (\kappa )}{\kappa }\right]^{1/2}\frac{1}{4\mathrm{cosh}\left[\frac{\pi }{4}(\kappa \kappa ^{})\right]}\left[\frac{\theta (\kappa ^{})}{\kappa ^{}}\right]^{1/2}g(\kappa ^{})=f(\kappa ),$$
(2.7)
where $`g=L_0[f]`$. It is convenient to introduce new functions $`G`$ and $`F`$ via

$$G(\kappa ^{})=\left[\frac{\theta (\kappa ^{})}{\kappa ^{}}\right]^{1/2}g(\kappa ^{})\text{and}F(\kappa )=\left[\frac{\theta (\kappa )}{\kappa }\right]^{1/2}f(\kappa ).$$
(2.8)
Then equation (2.7) takes a from

$$_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}\frac{1}{4\mathrm{cosh}\left[\frac{\pi }{4}(\kappa \kappa ^{})\right]}G(\kappa ^{})=F(\kappa ).$$
(2.9)
This equation can be easily solved by using Fourier transformation method:

$$F(\kappa ^{})=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega e^{i\omega \kappa ^{}}\stackrel{~}{F}(\omega )$$

Using that

$$_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}\frac{e^{i\omega \kappa ^{}}}{4\mathrm{cosh}\left[\frac{\pi }{4}(\kappa \kappa ^{})\right]}=\frac{e^{i\omega \kappa }}{\mathrm{cosh}2\omega }$$

one obtains

$$\stackrel{~}{G}(\omega )=\mathrm{cosh}2\omega \stackrel{~}{F}(\omega ).$$

Finally the solution to (2.9) is of the form

$$G(\kappa )=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega _{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}e^{i\omega (\kappa \kappa ^{})}\mathrm{cosh}2\omega F(\kappa ^{}).$$
(2.10)
The integral over $`\omega `$ does not exist (this once again confirms that $`L_0`$ is not a kernel operator) and therefore we can not simply switch integrals over $`\omega `$ and over $`\kappa ^{}`$. One can easily see that the integral operator in the right hand side is well defined at least on functions on $`\kappa ^{}`$ with Gaussian fall off at infinity. One can give a precise mathematical meaning to this operator if instead of Fourier transform one will use a Laplace transform and define $`0<\kappa <\mathrm{}`$. In this case we can rigorously prove that222Of course, expression (2.11) can be obtained by a formal integration of (2.10).

$$G(\kappa )=\frac{1}{2}\left[F(\kappa +2ii0)+F(\kappa 2i+i0)\right].$$
(2.11)
The function $`F`$ here is supposed to be holomorphic function in the strip $`2<\mathrm{}\kappa <2`$. It is allowed to have poles on the boundary of the strip. The $`\pm i0`$ in the formula above tells us how one should approach this poles. Namely one has to approach the pole on the line $`\kappa +2i`$ from the bottom and the one on the line $`\kappa 2i`$ from the top. Substitution of (2.8) into (2.11) yields expression (1.6) for the operator $`L_0`$.

Let us now check that (2.11) is indeed a solution to equation (2.9):

$$\begin{array}{c}_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}\frac{1}{4\mathrm{cosh}\left[\frac{\pi }{4}(\kappa \kappa ^{})\right]}\frac{1}{2}\left[F(\kappa +2ii0)+F(\kappa 2i+i0)\right]\hfill \\ \hfill =\frac{1}{8}_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}F(\kappa ^{})\left[\frac{i}{\mathrm{sinh}\frac{\pi }{4}(\kappa \kappa ^{})+i\epsilon \mathrm{cosh}\frac{\pi }{4}(\kappa \kappa ^{})}+\frac{i}{\mathrm{sinh}\frac{\pi }{4}(\kappa \kappa ^{})i\epsilon \mathrm{cosh}\frac{\pi }{4}(\kappa \kappa ^{})}\right]\\ \hfill =F(\kappa ).\end{array}$$

Here in the middle line we have made a shift of integration variable, which is allowed due to the fact that $`F(\kappa )`$ is holomorphic function on the strip $`2<\mathrm{}\kappa <2`$. This completes the derivation of formula (1.6) defining the operator $`L_0`$ in the class of difference operators.

Let me show now that the operator $`L_0`$ defining through (1.6) acts onto generating function $`f^{(\kappa )}`$ in the way it requires by (2.2). Indeed,

$$\begin{array}{c}L_0[f^{(\kappa )}(z)]=\frac{1}{4}[\frac{\sqrt{\kappa (\kappa +2i)}}{(\kappa +2i)\sqrt{\frac{2}{\kappa +2i}\mathrm{sinh}(\frac{\pi \kappa }{2})e^{i\pi }}}(1e^{(\kappa +2i)\mathrm{arctan}z})\hfill \\ \hfill +\frac{\sqrt{\kappa (\kappa 2i)}}{(\kappa 2i)\sqrt{\frac{2}{\kappa 2i}\mathrm{sinh}(\frac{\pi \kappa }{2})e^{i\pi }}}(1e^{(\kappa 2i)\mathrm{arctan}z})]\\ \hfill =\frac{i}{4}e^{\kappa \mathrm{arctan}z}\left(e^{2i\mathrm{arctan}z}e^{2i\mathrm{arctan}z}\right)=\frac{z}{1+z^2}e^{\kappa \mathrm{arctan}z}=z\frac{d}{dz}f^{(\kappa )}(z).\end{array}$$
(2.12)
## 3 Operators $`L_1`$ and $`L_1`$

The operators $`L_1`$ and $`L_1`$ are related by hermitian conjugation. Therefore it is enough to find the operator $`L_1`$. As it is done in the previous section we restrict our considerations to the zero momentum sector. The action of $`L_1`$ on the generating function $`f^{(\kappa )}(z)`$ is given by the formula

$$L_1[f^{(\kappa )}(z)]=\frac{1}{𝒩(\kappa )}+\frac{d}{dz}f^{(\kappa )}(z).$$
(3.1)
Therefore formally its kernel can be written as

$$L_1(\kappa ,\kappa ^{})\underset{m=1}{\overset{\mathrm{}}{}}v_{m+1}^{(\kappa )}\sqrt{(m+1)m}v_m^{(\kappa ^{})}$$
(3.2)
As it is explained in the Introduction the operator $`L_1`$ is completely determined by its action onto the one-particle Fock space. Therefore in order to find it we will use the same method as the one used in the previous section. There will be only one modification of the method related to the fact that the operator $`L_1`$ (conjugated to the operator $`L_1`$) has a one-dimensional kernel. First, we find an operator $`G_1`$ such that it “inverts” $`L_1`$ on the one-particle Fock space. More precisely $`G_1`$ satisfies the equation

$$G_1[L_1[f]](\kappa )=f(\kappa )\text{and}L_1[G_1[f]](\kappa )=P_1[f](\kappa ),$$
(3.3a)
where $`P_1`$ is a projector on the subspace on which $`L_1`$ is a non-degenerate operator:
$$P_1[f](\kappa )=f(\kappa )v_1^{(\kappa )}_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}v_1^{(\kappa ^{})}f(\kappa ^{}).$$
(3.3b)
From equation (3.2) one can easily obtain the formal expression for the kernel of operator $`G_1`$:

$$G_1(\kappa ,\kappa ^{})=\underset{m=1}{\overset{\mathrm{}}{}}v_m^{(\kappa )}\frac{1}{\sqrt{m(m+1)}}v_{m+1}^{(\kappa ^{})}.$$
(3.4)
Straightforward calculations (very similar to the ones presented in Appendix A) show that this expression does indeed define a distribution:

$$G_1(\kappa ,\kappa ^{})=\frac{1}{\kappa \sqrt{𝒩(\kappa )𝒩(\kappa ^{})}}\frac{1}{2\kappa }\sqrt{\frac{𝒩(\kappa )}{𝒩(\kappa ^{})}}\frac{\kappa \kappa ^{}}{\mathrm{sinh}\frac{\pi }{2}(\kappa \kappa ^{})}$$
(3.5)
and therefore $`G_1`$ is a kernel operator.

Second, we find the “inverse” of operator $`G_1`$ on the one particle Fock space. In other words we are going to solve the first equation in (3.3a). The resulting operator will coincide with the operator $`L_1`$ on the one-particle Fock space. Since $`L_1`$ is completely determined by its action on this subspace we will actually find operator $`L_1`$ on the whole Fock space.

So we need to solve the equation

$$_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}\frac{1}{\kappa \sqrt{𝒩(\kappa )𝒩(\kappa ^{})}}g(\kappa ^{})\frac{1}{2\kappa }_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}\sqrt{\frac{𝒩(\kappa )}{𝒩(\kappa ^{})}}\frac{\kappa \kappa ^{}}{\mathrm{sinh}\frac{\pi }{2}(\kappa \kappa ^{})}g(\kappa ^{})=f(\kappa ),$$
(3.6)
where $`g(\kappa ^{})=L_1[f](\kappa ^{})`$. Notice now that by the construction $`g`$ is a function that belongs to the image of $`L_1`$, therefore it satisfies the equation $`P_1[g]=g`$. This equation just reflects the fact that there is no state $`a_1^{}|0`$ in the image of operator $`L_1`$ restricted to the one-particle subspace. Eventually the first term in the equation is identically zero333 We use the fact that $`v_1^{(\kappa ^{})}=\frac{1}{\sqrt{𝒩(\kappa ^{})}}`$. It is useful to introduce new functions $`F`$ and $`G`$

$$G(\kappa ^{})=\sqrt{𝒩(\kappa ^{})}g(\kappa ^{})\text{and}F(\kappa )=\kappa \sqrt{𝒩(\kappa )}f(\kappa )$$
(3.7)
for which the equation takes extremely simple form

$$\frac{1}{2}_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}\frac{\kappa ^{}\kappa }{\mathrm{sinh}\frac{\pi }{2}(\kappa ^{}\kappa )}G(\kappa ^{})=F(\kappa ).$$

This equation can be solved using method of Fourier transform. Using that

$$_{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}e^{i\omega \kappa ^{}}\frac{\kappa ^{}}{\mathrm{sinh}\frac{\pi }{2}\kappa ^{}}=\frac{2}{\mathrm{cosh}^2\omega }$$

one can easily obtain the solution

$$G(\kappa )=\frac{1}{2}F(\kappa )\frac{1}{4\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega _{\mathrm{}}^{\mathrm{}}𝑑\kappa ^{}e^{i\omega (\kappa \kappa ^{})}\mathrm{cosh}2\omega F(\kappa ^{}).$$

Using the arguments we presented equation (2.10) one obtains

$$G(\kappa )=\frac{1}{2}F(\kappa )\frac{1}{4}\left[F(\kappa +i2^{})+F(\kappa i2^{})\right].$$

Here the function $`F(\kappa )`$ is supposed to be holomorphic on the strip $`2<\mathrm{}\kappa <2`$. Now using the relations (3.7) one obtains expression (1.8) defining the operator $`L_1`$ in a class of difference operators. An expression for the operator $`L_1`$ is easily obtained by hermitian conjugation.

Using the same technique as we used in deriving (2.12) we can show that operators $`L_1`$ and $`L_1`$ defined via (1.8) act on the generating function (1.1) as follows

$`L_1[f^{(\kappa )}(z)]`$ $`=z^2{\displaystyle \frac{d}{dz}}f^{(\kappa )}(z)`$ (3.8a)
$`L_1[f^{(\kappa )}(z)]`$ $`={\displaystyle \frac{1}{𝒩(\kappa )}}+{\displaystyle \frac{d}{dz}}f^{(\kappa )}(z)`$ (3.8b)
## 4 Commutation relations

In this section we are going to show that operators $`L_0`$, $`L_{\pm 1}`$ represented by the difference operators (1.6) and (1.8) satisfy the standard commutation relations.

First, let me demonstrate how to calculate the commutator of $`L_1`$ and $`L_1`$. From (1.8) we obtain the identities

$`L_1[L_1[g]](\kappa )`$ $`={\displaystyle \frac{\kappa }{2}}L_1[g]+{\displaystyle \frac{i}{4}}\left[\sqrt{\kappa (\kappa +2i)}L_1[g](\kappa +2i)\sqrt{\kappa (\kappa 2i)}L_1[g](\kappa 2i)\right];`$ (4.1a)
$`L_1[L_1[g]](\kappa )`$ $`={\displaystyle \frac{\kappa }{2}}L_1[g]{\displaystyle \frac{i}{4}}\left[\sqrt{\kappa (\kappa +2i)}L_1[g](\kappa +2i)\sqrt{\kappa (\kappa 2i)}L_1[g](\kappa 2i)\right].`$ (4.1b)
By subtracting the last equation from the first one one obtains

$$\begin{array}{c}[L_1,L_1][g](\kappa )=\frac{\kappa }{2}\left(L_1[g]L_1[g]\right)(\kappa )\hfill \\ \hfill +\frac{i}{4}\left[\sqrt{\kappa (\kappa +2i)}\left(L_1[g]+L_1[g]\right)(\kappa +2i)\sqrt{\kappa (\kappa 2i)}\left(L_1[g]+L_1[g]\right)(\kappa 2i)\right]\end{array}$$

Substitution of (1.8) and simplification yield

$$[L_1,L_1][g]=2L_0[g].$$
(4.2)
Second , we check commutation relations between $`L_0`$ and $`L_1`$.

$$\begin{array}{c}[L_0,L_1][g](\kappa )=\frac{\kappa }{2}L_0[g](\kappa )\hfill \\ \hfill +\frac{1}{4}\left[\sqrt{\kappa (\kappa +2i)}\left(L_1[g]iL_0[g]\right)(\kappa +2i)+\sqrt{\kappa (\kappa 2i)}\left(L_1[g]+iL_0[g]\right)(\kappa 2i)\right]\end{array}$$
(4.3)
Substitution of (1.6) and (1.8) and simplification yield

$$[L_0,L_1][g]=L_1[g].$$
(4.4)
## Acknowledgments

I would like to thank A. Konechny, H. Liu, S. Lukyanov and G. Moore for useful discussions. I would like to thank H. Liu and G. Moore for showing me the draft of paper . I would like to thank A. Konechny and M. Kroyter for the comments on this paper. The work of was supported in part by RFBR grant 02-01-00695.

Note added: While this paper was nearing completion, the paper  appeared, which contains in Section 5 similar expressions for $`L_0`$ and $`L_{\pm 1}`$. The reason for publishing this paper is mainly to give another point of view on the Virasoro operators in the $`\kappa `$-basis.

## Appendix

## Appendix A Calculations of $`G_0`$

The kernel for operator $`G_0`$ can be written as the following contour integral

$$G_0(\kappa ,\kappa ^{})=\frac{1}{\kappa \kappa ^{}\sqrt{𝒩(\kappa )𝒩(\kappa ^{})}}\frac{1}{2\pi i}_C\frac{dz}{z}(1e^{\kappa \mathrm{arctan}\frac{1}{z}})(1e^{\kappa ^{}\mathrm{arctan}z}).$$
(A.1)
Notice that this integral has two logarithmic singularities at the points $`\pm i`$. Therefore the complex plane has two cuts starting at those points. We choose them to lie on the imaginary axis (for more details on the contour see Appendix in ). To calculate the integral we deform the contour $`C`$ in such way that it will go along the cuts. We introduce $`z=ixϵ`$, $`z=ix+ϵ`$, $`z=ix+ϵ`$ and $`z=ixϵ`$ on contours $`C_+`$, $`C_{}`$, $`C_+^{}`$ and $`C_{}^{}`$ respectively

$`\kappa \kappa ^{}\sqrt{𝒩(\kappa )𝒩(\kappa ^{})}G_0(\kappa ,\kappa ^{})={\displaystyle \frac{1}{2\pi i}}{\displaystyle _1^{\mathrm{}}}`$ $`[{\displaystyle \frac{idx}{ixϵ}}(1e^{\kappa \mathrm{arctan}\frac{1}{ixϵ}})(1e^{\kappa ^{}\mathrm{arctan}(ixϵ)})`$
$`+{\displaystyle \frac{idx}{ix+ϵ}}\left(1e^{\kappa \mathrm{arctan}\frac{1}{ix+ϵ}}\right)\left(1e^{\kappa ^{}\mathrm{arctan}(ix+ϵ)}\right)`$
$`+{\displaystyle \frac{idx}{ix+ϵ}}\left(1e^{\kappa \mathrm{arctan}\frac{1}{ix+ϵ}}\right)\left(1e^{\kappa ^{}\mathrm{arctan}(ix+ϵ)}\right)`$
$`+{\displaystyle \frac{idx}{ixϵ}}(1e^{\kappa \mathrm{arctan}\frac{1}{ixϵ}})(1e^{\kappa ^{}\mathrm{arctan}(ixϵ)})]`$
Using that $`arctanh(ix\pm ϵ)=\pm \frac{\pi }{2}+i\mathrm{coth}^1x`$ and $`arctanh(ix\pm ϵ)^1=i\mathrm{coth}^1x`$ one obtains

$`={\displaystyle \frac{1}{2\pi i}}{\displaystyle _1^{\mathrm{}}}{\displaystyle \frac{dx}{x}}`$ $`[(1e^{\kappa (i\mathrm{coth}^1x)})(1e^{\kappa ^{}(\frac{\pi }{2}+i\mathrm{coth}^1x)})`$
$`\left(1e^{\kappa (i\mathrm{coth}^1x)}\right)\left(1e^{\kappa ^{}(\frac{\pi }{2}+i\mathrm{coth}^1x)}\right)`$
$`+\left(1e^{\kappa i\mathrm{coth}^1x}\right)\left(1e^{\kappa ^{}(\frac{\pi }{2}i\mathrm{coth}^1x)}\right)`$
$`(1e^{\kappa i\mathrm{coth}^1x})(1e^{\kappa ^{}(\frac{\pi }{2}i\mathrm{coth}^1x)})]`$
Now the change of variable

$$x=\mathrm{coth}u\text{and}dx=\frac{du}{\mathrm{sinh}^2u}$$

yields

$`={\displaystyle \frac{1}{\pi i}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{\mathrm{sinh}2u}}`$ $`[(1e^{i\kappa u})(1e^{\frac{\pi \kappa ^{}}{2}i\kappa ^{}u})`$
$`\left(1e^{i\kappa u}\right)\left(1e^{\frac{\pi \kappa ^{}}{2}i\kappa ^{}u}\right)`$
$`+\left(1e^{i\kappa u}\right)\left(1e^{\frac{\pi \kappa ^{}}{2}+i\kappa ^{}u}\right)`$
$`(1e^{i\kappa u})(1e^{\frac{\pi \kappa ^{}}{2}+i\kappa ^{}u})]`$
Simplification yields

$$\begin{array}{c}=\frac{2}{\pi i}\mathrm{sinh}\frac{\pi \kappa ^{}}{2}_0^{\mathrm{}}\frac{du}{\mathrm{sinh}2u}\left[e^{i\kappa ^{}u}e^{i(\kappa ^{}\kappa )u}e^{i\kappa ^{}u}+e^{i(\kappa ^{}\kappa )u}\right]\hfill \\ \hfill =\frac{2}{\pi i}\mathrm{sinh}\frac{\pi \kappa ^{}}{2}_{\mathrm{}}^{\mathrm{}}𝑑u𝒫\frac{1}{\mathrm{sinh}2u}\left[e^{i\kappa ^{}u}+e^{i(\kappa ^{}\kappa )u}\right]\end{array}$$
(A.2)
Using the fact that

$$\frac{1}{\pi i}_{\mathrm{}}^{\mathrm{}}𝑑u𝒫\frac{1}{\mathrm{sinh}2u}e^{i\beta u}=\frac{1}{2}\mathrm{tanh}\frac{\pi \beta }{4}$$
(A.3)
one obtains

$$=\mathrm{sinh}\frac{\pi \kappa ^{}}{2}\left[\mathrm{tanh}\frac{\pi \kappa ^{}}{4}+\mathrm{tanh}\frac{\pi (\kappa \kappa ^{})}{4}\right]=2\frac{\mathrm{sinh}\frac{\pi \kappa }{4}\mathrm{sinh}\frac{\pi \kappa ^{}}{4}}{\mathrm{cosh}\frac{\pi (\kappa \kappa ^{})}{4}}.$$
(A.4)
Finally one gets

$$G_0(\kappa ,\kappa ^{})=\left[\frac{\theta (\kappa )}{\kappa }\right]^{1/2}\left[\frac{\theta (\kappa ^{})}{\kappa ^{}}\right]^{1/2}\frac{1}{4\mathrm{cosh}\left[\frac{\pi }{4}(\kappa \kappa ^{})\right]}.$$
(A.5)


text:
# Asymptotic Behavior of the Emptiness Formation Probability in the Critical Phase of XXZ Spin Chain

## Abstract

We study the Emptiness Formation Probability (EFP) for the spin 1/2 $`XXZ`$ spin chain. EFP $`P(n)`$ detects a formation of ferromagnetic string of the length $`n`$ in the ground state. It is expected that EFP decays in a Gaussian way for large strings $`P(n)n^\gamma C^{n^2}`$. Here, we propose the explicit expressions for the rate of Gaussian decay $`C`$ as well as for the exponent $`\gamma `$. In order to confirm the validity of our formulas, we employed an ab initio simulation technique of the density-matrix renormalization group to simulate $`XXZ`$ spin chain of sufficient length. Furthermore, we performed Monte-Carlo integration of the Jimbo-Miwa multiple integral for $`P(n)`$. Those numerical results for $`P(n)`$ support our formulas fairly definitely.

We consider interacting spin 1/2 on the one-dimensional infinite lattice. The model is described by the Hamiltonian

$$=\underset{j=\mathrm{}}{\overset{\mathrm{}}{}}\left\{S_j^xS_{j+1}^x+S_j^yS_{j+1}^y+\mathrm{\Delta }\left(S_j^zS_{j+1}^z\frac{1}{4}\right)\right\},$$
(1)
where $`S=\sigma /2`$ and $`\sigma `$ are Pauli matrices. It exhibits diverse physics with varying the anisotropy parameter $`\mathrm{\Delta }`$. For $`\mathrm{\Delta }1`$ the ground state has a long-range ferromagnetic order, and a finite excitation gap opens in the low-lying spectrum. For $`\mathrm{\Delta }>1`$ the ground state develops antiferromagnetic order, the gap also exists. For moderate anisotropy $`1<\mathrm{\Delta }1`$, the antiferromagnetic order is dissolved by quantum fluctuations. The gap closes, and the long-range asymptotic of correlation function $`S_iS_j`$ decays as a power law at zero temperature.

In this article we study the Emptiness Formation Probability (EFP),

$$P(n)=\underset{j=1}{\overset{n}{}}\left(S_j^z+\frac{1}{2}\right),$$
(2)
in the domain of critical phase. The importance of EFP was emphasized in the book . It reflects the nature of the ground-state and it is a good indicator of the phase-separation.

In spite of its seemingly complicated definition, EFP has emerged quite naturally in the course of studies of the quantum integrability of the $`XXZ`$ spin chain. The first progress in calculating of EFP was achieved in Ref. , where $`P(3)`$ was obtained for the isotropic antiferromagnetic chain $`(\mathrm{\Delta }=1)`$. Subsequently, Jimbo and Miwa derived multiple integral representation for all correlation functions by means of vertex operators approach . Their integral formula for EFP is:

$`P(n)`$ $`=`$ $`\left(\nu \right)^{\frac{n(n+1)}{2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}\lambda _1}{2\pi }}\mathrm{}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}\lambda _n}{2\pi }}{\displaystyle \underset{a>b}{}}{\displaystyle \frac{\mathrm{sinh}(\lambda _a\lambda _b)}{\mathrm{sinh}\left(\left(\lambda _a\lambda _b\mathrm{i}\pi \right)\nu \right)}}`$ (3)
$`\times {\displaystyle \underset{k=1}{\overset{n}{}}}{\displaystyle \frac{\mathrm{sinh}^{nk}\left(\left(\lambda _k+\mathrm{i}\pi /2\right)\nu \right)\mathrm{sinh}^{k1}\left(\left(\lambda _k\mathrm{i}\pi /2\right)\nu \right)}{\mathrm{cosh}^n(\lambda _k)}},`$
where the parameter $`\nu `$ is related to the anisotropy $`\mathrm{\Delta }`$ as

$$\nu =\frac{1}{\pi }\mathrm{cos}^1(\mathrm{\Delta }).$$
(4)
The integral representation was reproduced in the framework of algebraic Bethe ansatz in . A good deal of new developments have been reported for EFP. For example, in the case of $`\mathrm{\Delta }=1/2`$ the simple formula for $`P(n)`$ has been conjectured in  and then proved in . For $`\mathrm{\Delta }=0`$ the asymptotic form of $`P(n)`$ as $`n1`$ was found in . In either case the large-$`n`$ asymptotic behavior is given by $`P(n)n^\gamma C^{n^2}`$ with

$`C`$ $`=`$ $`\sqrt{2},\gamma ={\displaystyle \frac{1}{4}},(\nu =1/2),`$
$`C`$ $`=`$ $`{\displaystyle \frac{8}{3\sqrt{3}}},\gamma ={\displaystyle \frac{5}{36}},(\nu =1/3).`$ (5)
Moreover, $`P(n)`$ was calculated recently for $`\mathrm{\Delta }=1`$ for strings of the length $`n6`$, see papers . Their result suggests that $`P(n)`$ at $`\mathrm{\Delta }=1`$ decays in a Gaussian way as well. Meanwhile, the asymptotic behavior of EFP for the whole critical regime was analyzed in field-theoretical framework, by Abanov and Korepin . The Gaussian decay naturally appears in this approach.

Based on those developments, it would appear reasonable that the following asymptotic form of EFP

$$P(n)An^\gamma C^{n^2},$$
(6)
holds for all over the critical regime $`(1<\mathrm{\Delta }1)`$. We propose the explicit expressions for the rate of Gaussian decay $`C`$ and the power-law exponent $`\gamma `$:

$`C`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }^2(1/4)}{\pi \sqrt{2\pi }}}\mathrm{exp}\left\{{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}t}{t}}{\displaystyle \frac{\mathrm{sinh}^2(t\nu )\mathrm{e}^t}{\mathrm{cosh}(2t\nu )\mathrm{sinh}(t)}}\right\},`$ (7)
$`\gamma `$ $`=`$ $`{\displaystyle \frac{1}{12}}+{\displaystyle \frac{\nu ^2}{3(1\nu )}}.`$ (8)
One can confirm that above formulas reproduce the exact results for $`\mathrm{\Delta }=0`$ and $`\mathrm{\Delta }=1/2`$ (5).

In order to confirm the validity of the formulas for general anisotropies $`0\nu <1`$, we have performed extensive numerical calculations of two kinds: One is the first-principle simulation method, namely, the density-matrix renormalization group (DMRG) , and the other is the Monte-Carlo numerical integration of the multiple integration formula (3) of Jimbo and Miwa. As noted afterwards, those methods are compensative, and we were able to perform reliable simulations for over all the critical regime. As a result, we could confirm fairly definitely that the above general formulas are indeed correct.

Let us turn to addressing the numerical-simulation results. In Tables 1-4, we have listed the DMRG results of EFP for $`\nu =0.2,\mathrm{},0.8`$. (The methodological details will be explained afterwards.)

In addition, in each table, we have presented the logarithm of the ratio of two adjacent EFPs, which should behave in the form

$$\mathrm{ln}\left(\frac{P(n)}{P(n+1)}\right)\gamma \mathrm{ln}(1+\frac{1}{n})+(2n+1)\mathrm{ln}C,$$
(9)
according to the long-distance asymptotic formula (6). Note that, after taking the ratio, we are able to kill the contribution of the constant factor $`A`$ of (6), for which, at present, we have no analytical prediction. (We, however, could estimate $`A`$ numerically from our DMRG data as is shown in Table 5.) In this way, the resultant processed data can be directly comparable with the analytical conjecture.

We see that our processed DMRG data $`\mathrm{ln}(P(n)/P(n+1))`$ are extremely close to the analytical prediction (9). Actually we find that they coincide up to two or three digits for general values of $`n`$ and $`\nu `$. Therefore we conclude that our formulas for the asymptotic form are valid indeed for all over the regime $`0\nu <1`$. Especially see the Table 6, where we compare the asymptotic formula with the known exact values for $`\mathrm{\Delta }=1`$ . Note that, in this case, we have $`C=\mathrm{\Gamma }^2(1/4)/(\pi \sqrt{2\pi })=1.66925\mathrm{},\gamma =1/12`$, and $`A=0.841.`$

We, however, remark that rather large discrepancies are seen for such cases either very large $`n`$ for small $`\nu `$ or small $`n`$ for large $`\nu `$. The former deviations are merely due to the numerical round off errors. Since in computers, real numbers are stored in 8-byte space, and the precision is of the order of $`10^{15}`$ at best. Hence it is in principle difficult to calculate the correlations less than the magnitude $`<10^{12}`$ reliably. The latter deviations are not so surprising, because our general formulas should be justified for long distances of EFP. This short-range deviation will be further exploited in the succeeding Monte-Carlo numerical integration analyses.

So far, we had analyzed processed EFP data. In Fig. 1, we present the raw (unprocessed) EFP data obtained with the DMRG method. The dotted lines are the analytical conjecture with the factors given in Table 5. We see that our first-principle data are well-fitted by our general formulas (7) and (8) for various values of $`\nu `$. However, as noted above, for some cases, there appear deviations.

Let us turn to addressing details of the DMRG method . DMRG is a sort of computer-aided real-space decimation, where the number of states for a renormalized block is retained within a bound tractable in computers. In this way, through the successive applications of DMRG, we are able to access very large system sizes. Hence, this method is suitable for surveying the long-distance asymptotic behavior of EFP. As a matter of fact, efficiency of the method was demonstrated in the preceding works : This success was rather unexpected, because in general, DMRG is not very efficient for such systems exhibiting criticality. In particular, long-distance asymptotic behavior of two-point correlation function is deteriorated rather severely. Those pathologies would be due to the truncation of bases through numerical renormalizations. However, as far as EFP is concerned, DMRG is proven to be free from such difficulties. Below, we will overview technical points of our numerical simulation. Our algorithm is standard , and we refer readers to a recent proceeding  for full account of methodological details. We have employed the so-called infinite algorithm, which is adequate to investigate bulk properties at the ground state. (For the purpose of studying finite-size scaling behavior, the finite algorithm would be more suitable.) We have remained, at most, $`m=300`$ bases for a renormalized block. The density matrix eigenvalue $`\{w_\alpha \}`$ of remaining bases indicates the statistical weight of each remained state. We found $`w_\alpha >510^{12}`$: That is, we have remained almost all relevant states with appreciable statistical weight $`w_\alpha >510^{12}`$, which may indicate error of the present simulation. We have repeated 300 renormalizations, and hence, total system size extends to $`L=600`$. The DMRG data alternates in turn through renormalizations: Note that the number of spins consisting a renormalized block increases by one after another through renormalizations. The problem is that the Hilbert-space structures are incompatible with respect to those cases whether the block contains even or odd spins. (Haldane system ($`S=1`$) is not affected by this difficulty.) Therefore, we have taken arithmetic mean over those two cases.

The DMRG method is not efficient in the close vicinity of the ferromagnetic isotropic point $`\nu =1`$. Because in the vicinity of that point, the spin-wave velocity tends to vanish, and there appear numerous nearly-degenerated low-lying levels. Those nearly-degenerated levels are very hard to resolve in the process of numerical diagonalization. The diagonalization is a significant part in the DMRG procedure, and hence, DMRG becomes hardly applicable. In order to compensate this drawback, we have employed another numerical method, that is, the Monte-Carlo integration. We will explain it in the following.

In Fig. 2, we have presented numerical data with Monte-Carlo integration for the multiple integral formula (3) of Jimbo and Miwa. The technical details will be explained afterwards. From the plot, we see that the decay of EFP is gradually modified as for $`\nu 1`$, and surprisingly enough, the decay approaches the simple exponential formula of $`(1/2)^n`$ eventually. This feature is precisely in accord with the aforementioned finding that DMRG data deviate from the Gaussian formula as for $`\nu 1`$. Therefore, we see that at least for short-distance asymptotic form of EFP is governed by the simple exponential decay $`0.5^n`$ as $`\nu 1`$. As a matter of fact, right at $`\nu =1`$, one may easily verify the pure exponential decay $`P(n)`$ from the Jimbo-Miwa formula.

The crossover from Gaussian to exponential decay suggests that the quantum fluctuations are suppressed as $`\nu 1`$. As as matter of fact, field theoretical consideration  reveled that EFP measures the probability of formations of ferromagnetic islands with size $`n`$ surrounded by antiferromagnetic background from the viewpoint of Euclidean space-time. Therefore, our result suggests that this Euclidean space-time picture is deteriorated for $`\nu 1`$, because the spin-wave velocity tends to vanish, and thus, the size of the islands grows along the imaginary-time direction abruptly.

In the following, we explain technical details of the Monte-Carlo integration. We have used the subroutine package described in the textbook . As a random number generator, we have employed “Mersenne Twistor”, MT19937 in Ref. . For each plot, we have performed eight-billion Monte-Carlo steps. This main Monte-Carlo procedure is preceded by preliminary five million Monte-Carlo steps which are aimed to improve the efficiency of Monte-Carlo sampling by surveying the integrand in the multi-dimensional space. Irrespective of $`n`$ and $`\nu `$, we found that the statistical errors are of the order of $`10^8`$.

To summarize, we have advocated compact explicit formulas (7) and (8) for the long-distance asymptotic form behavior of EFP. The formulas reproduce the presently-available exact results obtained at special solvable points of $`\mathrm{\Delta }=0`$ and $`1/2`$, and therefore, we expect that our formulas are valid for over all the critical regime. We have performed extensive simulations of DMRG and Monte-Carlo integration. As a result, we found that our general formulas are indeed correct over the critical regime; see Tables and Fig. 1. In other words, EFP decays obeying the Gaussian formula (6) in the domain of criticality, and the decay rate becomes slower in the ferromagnetic side. This tendency is understood by a physical picture that EFP detects a formation of condensed particles (ferromagnetic string). Apparently, for ferromagnetic side, the condensation becomes promoted, and EFP gets enhanced, albeit the asymptotic form is still maintained to be Gaussian. However, as for the extreme limit to the ferromagnetic isotropic point, the Gaussian decay is smeared out by an exponential form at least for short distances. This fact reflects the suppression of quantum fluctuations in the ferromagnetic side. In other words, this crossover may be regarded as the precursor of the onset of true (long-range order) phase-separation for $`\mathrm{\Delta }<1`$, where EFP should not decay, but may saturate to a certain finite resident value. Finally we like to remark that the rate of Gaussian decay $`C`$, Eq.(7), can be evaluated analytically for some other rational values of $`\nu `$. For example, we have found

$`C|_{\nu =1/4}`$ $`=`$ $`2^{\frac{1}{4}}\mathrm{e}^{\frac{\mathrm{G}}{\pi }},(\mathrm{G}:\mathrm{Catalan}\mathrm{number}),`$
$`C|_{\nu =1/6}`$ $`=`$ $`2\sqrt{{\displaystyle \frac{2}{3}}},C|_{\nu =2/3}={\displaystyle \frac{16(\sqrt{2}1)}{3\sqrt{3}}},`$ (10)
by use of the integral representations for the logarithm of (multiple) Gamma functions. These facts may indicate a possibility to prove our analytic formulas (6) rigorously for these values of $`\nu `$ just in the similar way as $`\nu =1/2`$ and $`1/3`$.

Acknowledgment We thank A. Abanov, H.Boos, F.Essler, F.Smirnov, M.Takahashi and A. Zamolodchikov for useful discussions. V. Korepin was supported by by NSF Grant PHY-9988566, S. Lukyanov was supported in part by DOE Grant DE-FG02-96ER40959, Y. Nishiyama was supported by Grant-in-Aid for Young Scientists No. 13740240, M. Shiroishi was supported by Grant-in-Aid for Young Scientists No. 14740228.


text:
# Exploring Defective Eigenvalue Problems with the Method of Lifting

## 1 Introduction

In this paper, we present a method for improving the accuracy of computed eigenvalues and eigenvectors of defective matrices . This algorithm, which we refer to as lifting, involves construction of a nondefective matrix in a higher dimensional space. Standard linear algebra procedures are shown to yield more accurate numerical results for the lifted matrix. Our scheme formally resembles the construction of an augmented system in pseudo-arclength continuation . Its formulation was motivated by the structure of the linear stability analysis of a particular chemical reaction-diffusion system, known as the chlorine-dioxide–iodine–malonic-acid system  (henceforth referred to as CDIMA), which exhibits the Turing instability.

First, we present the Lifting Theorem. We demonstrate its numerical implementation on “small” and “large” defective eigenproblems. We describe the connection with the CDIMA linear stability problem, and discuss the general applicability of this scheme.

## 2 Lifting

An eigenvalue is defective if its algebraic multiplicity is greater than its geometric multiplicity. (The geometric multiplicity of an eigenvalue equals the dimension of the space spanned by the eigenvectors corresponding to that eigenvalue.) Consider an eigenvalue with geometric multiplicity equal to one: It is nondefective if the corresponding left and right eigenvectors are nonorthogonal. For a defective eigenvalue, lifting corresponds to a geometric procedure of casting the corresponding eigenvectors of the matrix and the adjoint matrix into a higher dimensional space in which these vectors are no longer orthogonal. In this section, we show how to construct such a well-conditioned matrix in the higher dimensional space from a defective matrix.

For simplicity, we define the notation:

algebraic multiplicity of eigenvalue $`\mu `$ of matrix $`𝑴`$ $``$ $`AM(𝑴;\mu )`$
geometric multiplicity of eigenvalue $`\mu `$ of matrix $`𝑴`$ $``$ $`GM(𝑴;\mu ).`$
Although in the following we consider zero eigenvalues and corresponding nullvectors, the discussion holds generally for nonzero eigenvalues, with $`𝑨=𝑴\mu 𝑰_N`$, where $`𝑰_N`$ is the appropriate $`N\times N`$ identity matrix.

Theorem: Consider an $`N\times N`$ matrix $`𝑨`$, such that zero is a defective eigenvalue of $`𝑨`$ with $`AM(𝑨;0)=n`$ and $`GM(𝑨;0)=1`$, $`n>1`$. The right and left unit nullvectors of $`𝑨`$ are given by $`\mathit{\varphi }`$ and $`𝝍^{\mathbf{}}`$, respectively. An $`(N+1)\times (N+1)`$ matrix $`𝓛`$ constructed according to

$$𝓛=\left(\begin{array}{cc}𝑨& \mathrm{𝟎}\\ \mathrm{𝟎}^{}& 0\end{array}\right)+𝒗𝒘^{\mathbf{}}.$$
(1)
where

$$𝒗=\left(\begin{array}{c}𝐯\\ \eta \end{array}\right),𝒘=\left(\begin{array}{c}𝐰\\ \omega \end{array}\right),$$
(2)
will have a simple zero eigenvalue if

$`\text{(i)}𝐰^{}\mathit{\varphi }`$ $``$ $`0,`$ (3)
$`\text{(ii)}𝝍^{}𝐯`$ $``$ $`0,`$ (4)
$`\text{(iii)}\eta \omega `$ $``$ $`0.`$ (5)
$`𝐯`$ and $`𝐰`$ are $`N`$-dimensional vectors, and $`\eta `$ and $`\omega `$ are scalars.

Proof: First, we show there exists a nontrivial vector $`𝚽`$, given by

$$𝚽=\left(\begin{array}{c}𝐱\\ \xi \end{array}\right),$$
(6)
such that

$$𝓛𝚽=\mathrm{𝟎}.$$
(7)
For $`𝚽`$ to be a right nullvector of $`𝓛`$, the $`N`$ components of $`\left(𝓛𝚽\right)`$ in the original space, as well as the $`(N+1)^{th}`$ component in the inflated space must be zero:

$`𝑨𝐱+\alpha 𝐯`$ $`=`$ $`0`$ (8)
$`\alpha \eta `$ $`=`$ $`0,`$ (9)
where

$$\alpha =𝐰^{}𝐱+\omega \xi .$$
(10)
Since $`\eta 0`$, for Eq. 9 to be satisfied, we must have $`\alpha =0`$. Then, Eq. 8 implies $`𝐱=r\mathit{\varphi }`$, where $`r`$ is a constant, and $`\xi `$ must be given by

$$\xi =r\frac{𝐰^{}\mathit{\varphi }}{\omega }.$$
(11)
Hence, the pointwise projection of the right nullvector of $`𝓛`$ onto the original $`N`$-dimensional space is the right nullvector of $`𝑨`$, up to a multiplicative constant.111If $`𝚽`$ is normalized to give the unit right nullvector, then $`r^2=\left(1+𝐰^{}\mathit{\varphi }/\omega \right)^1`$. Additionally, its $`(N+1)^{th}`$ component, $`\xi `$, is nonzero.

Likewise, we can show that $`𝓛`$ has a nontrivial left nullvector $`𝚿^{\mathbf{}}`$ given by

$$𝚿=\left(\begin{array}{c}𝐲\\ \zeta \end{array}\right),$$
(12)
where

$$𝐲=s\text{ }𝝍,\zeta =s\frac{𝝍^{\mathbf{}}𝐯}{\eta }.$$
(13)
Again, $`s`$ is a constant, and we have $`\zeta 0`$.

Next, we must show that zero is a simple eigenvalue of $`𝓛`$. This will be done in two steps:

$`(\mathrm{i})`$ $`GM(𝓛;0)=1,`$ (14)
$`(\mathrm{ii})`$ $`𝚿^{\mathbf{}}𝚽\mathbf{=}𝝍^{\mathbf{}}\mathit{\varphi }+\zeta \xi 0.`$
First, we prove $`GM(𝓛;0)=1`$ by contradiction. Assume $`GM(𝓛;0)>1`$. Then, there exists $`𝚽^{\mathbf{}}\mathbf{}𝚽`$ such that

$$𝓛𝚽^{\mathbf{}}=\mathrm{𝟎},𝚽^{\mathbf{}}=\left(\begin{array}{c}𝐱^{}\\ \xi ^{}\end{array}\right).$$
(15)
Eqs. 8 and 9 imply $`𝑨𝐱^{}=\mathrm{𝟎}`$. We distinguish two cases:

* If $`𝐱^{}r\mathit{\varphi }`$, then $`𝐱^{}`$ is a right nullvector of $`𝑨`$, different from $`r\mathit{\varphi }`$. This violates the assumption $`GM(𝑨;0)`$=1.
* If $`𝐱^{}=r\mathit{\varphi }`$, we must have $`\xi ^{}\xi `$ for $`𝚽^{\mathbf{}}\mathbf{}𝚽`$ . However, from Eqs. 9 and 10, $`\xi ^{}=r𝐰^{}\mathit{\varphi }/\omega =\xi `$. Therefore, $`𝚽^{\mathbf{}}`$ is not different from $`𝚽`$.

Hence, we have shown that $`GM(𝓛;0)=GM(𝑨;0)=1`$. (This proof is easily generalized to $`GM(𝓛;0)=GM(𝑨;0)=\nu `$.) To show (ii), we note that zero being a defective eigenvalue of $`𝑨`$ with $`GM(𝑨;0)=1`$ implies $`𝝍^{\mathbf{}}\mathit{\varphi }=0`$. However, the inner product given by Eq. 14 will still be nonzero, since $`\zeta \xi 0`$. This completes the proof that zero is a simple eigenvalue of $`𝓛`$.

It is important to consider the generic case in numerical applications, where $`𝑨`$ is not exactly defective but almost-defective, $`𝝍^{\mathbf{}}\mathit{\varphi }0`$. For zero to be a nondefective eigenvalue of $`𝓛`$, the inner product given by Eq. 14 must be nonzero:

$$𝝍^{\mathbf{}}\mathit{\varphi }+\left(\frac{𝝍^{\mathbf{}}𝐯}{\eta }\right)\left(\frac{𝐰^{}\mathit{\varphi }}{\omega }\right)0.$$
(16)
For example, constructing $`𝓛`$ with $`𝒗`$ and $`𝒘`$ such that $`𝐯=\mathit{\varphi }`$, $`\eta =1`$, and $`\omega =𝐰^{}\mathit{\varphi }`$, would violate the condition above. Therefore, when $`𝑨`$ is almost-defective, Eq. 16 must be additionally satisfied to guarantee that $`𝓛`$ will not be defective. Below, we give a general prescription for choosing the lifting vectors, which unlike the pathological special case just presented, should avoid construction of a defective lifted matrix.

Consider choosing (for simplicity) $`\eta =\omega =1`$, and $`(𝐯,𝐰)=(\beta 𝐯_𝐫,\gamma 𝐰_𝐫)`$:

$$\mathrm{v}_{\mathrm{r}}^{}{}_{i}{}^{}=\mathrm{𝚛𝚊𝚗𝚍}(i)/N_\mathrm{v},\mathrm{w}_{\mathrm{r}}^{}{}_{j}{}^{}=\mathrm{𝚛𝚊𝚗𝚍}(j)/N_\mathrm{w},$$
(17)
where $`\mathrm{𝚛𝚊𝚗𝚍}(i)`$ are random numbers on $`[1,1]`$, $`N_\mathrm{v}`$ and $`N_\mathrm{w}`$ are normalizing factors leading to unit random vectors, and $`\gamma `$ and $`\beta `$ are constants. This construction should guarantee that when $`𝑨`$ is singular and exactly defective (with $`AM(𝑨;0)=n`$, $`GM(𝑨;0)=1`$), or close to defective, then zero will be a simple eigenvalue of $`𝓛`$. We will refer to the constants $`\gamma `$ and $`\beta `$ as lifting parameters, since the $`(N+1)^{th}`$ components of the left/right nullvectors of $`𝓛`$ are proportional to these constants, respectively, and measure the projections of these vectors onto the new, inflated direction.

To extend the above theorem and proof to the case where $`𝑨`$ is defective and $`AM(𝑨;0)>GM(𝑨;0)>1`$, we must show that $`𝓛`$ can be constructed such that zero is a nondefective, multiple eigenvalue with $`AM(𝓛;0)=GM(𝓛;0)=GM(𝑨;0)`$. With multiple zero eigenvalues, $`𝓛`$ will be defective if there exists a vector $`𝐠`$ such that

$$𝓛𝐠=𝚽_𝒊,$$
(18)
for some $`𝚽_𝒊𝒩(𝓛)`$, where $`i=1\mathrm{}\nu `$, and $`𝒩(𝓛)`$ is the nullspace of $`𝓛`$. (If such a vector $`𝐠`$ existed, then it would be a generalized eigenvector of $`𝓛`$, associated with an $`m\times m`$ Jordan block, where $`m>1`$.) Taking the inner product of Eq. 18 with $`𝚿_𝒋^{\mathbf{}}`$, where $`𝚿_𝒋^{\mathbf{}}𝒩(𝓛^{})`$, $`j=1\mathrm{}\nu `$, gives

$$𝚿_𝒋^{\mathbf{}}𝚽_𝒊=0.$$
(19)
Therefore, for $`𝓛`$ to be nondefective, it must be constructed such that no vector in the right nullspace of $`𝓛`$ is orthogonal to its entire left nullspace.

## 3 Numerical Results

### 3.1 $`2\times 2`$ Defective Eigenproblem

Consider the following matrix, $`𝑴`$:

$$𝑴=\left(\begin{array}{cc}\pi & 1\\ \pi ^2/4& ϵ\end{array}\right).$$
(20)
The eigenvalues of $`𝑴`$ are given by

$$\mu _\pm =\left[\pi +ϵ\pm \sqrt{\left(ϵ^22\pi ϵ\right)}\right]/2.$$
(21)
The two components of an eigenvector $`\mathit{\varphi }^{\mathbf{}}`$$`=(\varphi _1,\varphi _2)`$ satisfy

$$\varphi _2/\varphi _1=\mu \pi .$$
(22)
where $`\mu =\mu _\pm `$. For $`ϵ=0`$, $`𝑴`$ is defective. We use the lifting algorithm to explore the error in computing the eigenvector associated with $`\mu _+`$ when $`𝑴`$ is close to being defective, $`ϵ1`$. The lifted matrix $`𝓛`$ is constructed with $`𝑨=𝑴\mu _+𝑰_2`$ according to

$$𝒗=\left(\begin{array}{c}\beta 𝐯_𝐫\\ 1\end{array}\right),𝒘=\left(\begin{array}{c}\beta 𝐰_𝐫\\ 1\end{array}\right),$$
(23)
where $`(𝐯_𝐫,𝐰_𝐫)`$ are two-dimensional random unit vectors, and we set $`\gamma =\beta `$ without loss of generality.

The right unit nullvector $`𝚽`$ of $`𝓛`$ is computed using MATLAB’s eig() function . We explore the lifting error, constructed according to

$$(ϵ,\beta )=\left|\mathrm{\Phi }_2/\mathrm{\Phi }_1(\mu _+\pi )\right|.$$
(24)
as a function of parameters $`ϵ`$ and $`\beta `$.

Figure 1 shows $`(ϵ,\beta )`$ as a function of the lifting parameter $`\beta `$ for different values of $`ϵ`$. The basic trend is that of decreasing error as a function of increasing $`\beta `$ for small to moderate values of $`\beta `$. In this figure, each point corresponds to the mean error computed using $`1000`$ unit random lifting vectors, $`(𝐯_𝐫,𝐰_𝐫)`$. The root-mean-square (rms) of the distribution of errors is of the order of the mean error for small values of lifting parameter. For $`\beta \stackrel{>}{}𝒪(1)`$, the rms is at most an order of magnitude greater than the mean error for values of $`ϵ`$ close to machine precision; for larger values of $`ϵ`$, the rms is equal to a few times the mean error. We have verified that the larger values of rms are a result of a few outliers. In Figure 2, for lifting parameter $`\beta =1.0`$, we compare the mean error with lifting (open circles) and without lifting (open triangles) as a function of $`ϵ`$. In this figure, we also show the error using lifting vectors given by (v, w)$`=(𝝍,\mathit{\varphi })`$, such that the nonorthogonality conditions Eqs. 3 and 4 are satisfied by definition (open squares). We note that the lifting error is within an order of magnitude of machine precision for all values of $`ϵ`$.

The condition number of the simple singular value of the lifted matrix is given by the reciprocal of

$$s(0)=\left|𝚿^{}𝚽\right|,$$
(25)
where $`𝚿`$ and $`𝚽`$ are the left and right unit nullvectors of $`𝓛`$. In Figure 3, we show the condition number as a function of the lifting parameter $`\beta `$ for different values of $`ϵ`$, using a single pair of unit random lifting vectors: For $`ϵ1`$, the condition number clearly improves with increasing lifting parameter.

### 3.2 $`N\times N`$ Defective Eigenproblem

For applications in which the matrix is a discrete approximation to a continuous operator, $`𝑨`$ will have dimension $`N1`$. To explore the method of lifting in computing the defective eigenvector of larger matrices, we constructed such an example from the $`2\times 2`$ matrix above: with $`𝑴`$ in the top left block, along with the $`(N2)\times (N2)`$ block tridiagonal matrix from Poisson’s equation in the bottom right, we applied a similarity transformation to the resulting $`N\times N`$ matrix to obtain a general matrix with no special properties (such as symmetry or bandedness). The eigenvalue, $`\mu _+`$, of this matrix is exactly known, and the lifted matrix $`𝓛`$ is constructed as in the $`2\times 2`$ case.

Once the right nullvector of $`𝓛`$, given by $`𝚽^{\mathbf{}}=(𝐱,\xi )`$ is numerically computed, the lifting error is constructed as before:

$$(ϵ,\beta )=\left|\right(𝑺𝐱\left){}_{2}{}^{}/\right(𝑺𝐱\left){}_{1}{}^{}(\mu _+\pi )\right|,$$
(26)
(where $`𝑺`$ is the similarity transformation). In Figure 4, we show the mean error and its rms computed using $`50`$ unit random lifting vectors, $`(𝐯_𝐫,𝐰_𝐫)`$. We find that for small values of $`\beta `$, the lifting error is dominated by the error in computing the “zero” eigenvalue, $`\lambda _0`$, of $`𝓛`$. For $`\beta \stackrel{>}{}𝒪(1)`$, the magnitude of this eigenvalue remains at machine precision, while the lifting error is larger. For values of $`\beta `$ such that the magnitude of the “zero” eigenvalue is at machine precision, the lifting error is dominated by the error in computing the associated eigenvector. We note that there is an “optimal” value of $`\beta `$ giving the smallest lifting error, which is of order unity. For values of $`\beta `$ much larger than unity, the error increases as a function of increasing $`\beta `$.

## 4 CDIMA Eigenproblem

The symmetry-breaking instability of a system from a homogeneous steady state to a patterned state, predicted in 1952 by Alan Turing , was observed for the first time in the chlorite-iodide-malonic acid ClO$`{}_{2}{}^{}{}_{}{}^{}`$-I--MA (CIMA) reaction-diffusion system . In practice, due to boundary chemical feeds, the steady state is not homogeneous, but rather depends on the single space coordinate along the feed gradient. A realistic model of the simpler chlorine dioxide-iodine-malonic acid ClO2-I2-MA (CDIMA) reaction, which is similar to CIMA in terms of its stationary pattern-forming and dynamical behavior, was put forth by Lengyel, Rabai, and Epstein (henceforth referred to as LRE) . In particular, they demonstrated that the fortuitous choice of starch as the color indicator for visualizing the patterns in this system provides the necessary rescaling of the diffusion coefficient of the activator (iodide) relative to the inhibitor (chlorite) species required for this instability to occur. Starch is a large molecule that binds to the gel matrix of the reactor and is thus effectively immobile. The reversible binding of iodide with starch results in the purple starch-triiodide complex and an effective slowing down of the diffusion rate of iodide.

We begin with the general formulation of an $`(N+1)`$-component reaction-diffusion system with $`N`$ mobile and a single immobile species (as in the CDIMA system), following :

$$\frac{𝓒}{t}=𝓓^2𝓒+𝓕(𝑼)+𝜻g.$$
(27)
In the evolution equation above,

$$𝓒=\left(\begin{array}{c}𝑼\\ W\end{array}\right),𝓓=\left(\begin{array}{cc}𝑫& \mathrm{𝟎}\\ \mathrm{𝟎}^{}& 0\end{array}\right),𝓕=\left(\begin{array}{c}𝒇\\ 0\end{array}\right),𝜻=\left(\begin{array}{c}𝑹\\ 1\end{array}\right),$$
(28)
$`𝑼`$ represents the concentrations of mobile species and $`W`$ that of the immobile species, $`𝒇=𝒇(𝑼;𝜶)`$ describes the reactions among the mobile species with control parameters $`𝜶`$, and $`g=g(𝑼,W;𝜷)`$ describes the reaction between the mobile species and the immobile species with control parameters $`𝜷`$. $`𝑫`$ is the diffusion matrix for the mobile subsystem, and $`𝑹`$ is a constant stiochiometry vector. We refer the interested reader to previous works for the description of the LRE model of the CDIMA system.

The steady state of this system in the presence of externally imposed boundary feed gradients satisfies:

$`𝒇(𝑼_s;\alpha )+𝑫{\displaystyle \frac{^2𝑼_s}{z^2}}`$ $`=`$ $`0,`$ (29)
$`g(𝑼_s,W_s;\beta )`$ $`=`$ $`0`$ (30)
with boundary conditions $`𝓒_s(0)=𝓒_0`$ and $`𝓒_s(\mathrm{})=𝓒_{\mathrm{}}`$. We will focus on the physically relevant special case of scalar diffusion, $`𝑫=D\text{ }𝑰\text{ }`$, where  $`𝑰`$  is the identity matrix, as the diffusion constants of ions in solution are comparable. It is clear from the formulation given by Eq. 27, that the steady state of the mobile subsystem is independent of the immobile reactions. However, as illustrated in previous works , the immobile subsystem affects the stability of the mobile system. The stability of the steady state to a mode with transverse wavenumber $`k`$ ($`\widehat{k}\widehat{z}`$) leads to the following linear operator:

$$𝓛_k=\left(\begin{array}{cc}𝒅𝒇& \mathrm{𝟎}\\ \mathrm{𝟎}^{}& 0\end{array}\right)D\left(k^2+_z^2\right)\left(\begin{array}{cc}𝐈& \mathrm{𝟎}\\ \mathrm{𝟎}^{}& 0\end{array}\right)+𝜻𝒅𝒈^{},$$
(31)
where $`𝒅𝒇=𝒇_𝐔`$ is the Jacobian derivative of $`𝒇`$ with respect to $`𝑼`$, and $`𝒅𝒈^{}=(g_𝐔,g_W)`$ is the Jacobian derivative of $`g`$ with respect to $`(𝑼,W)`$, evaluated at the steady state. The eigenvalue problem determining the stability of $`𝓒_s(z)`$ to the mode with wavenumber $`k`$ is given by

$$𝓛_k𝚿_k=\lambda _k𝚿_k$$
(32)
where $`𝚿=(𝚿^M,\mathrm{\Psi }^I)`$. The space $`S`$ on which the linear operator acts can be written as $`S=S_MS_I`$, where $`S_M`$ and $`S_I`$ are mobile and immobile subspaces, respectively. In gradient systems, the basis set spanning $`S_M`$ is given by the Cartesian product of $`N`$ orthogonal unit vectors $`\{\widehat{e}_i\}`$, $`i=1,\mathrm{},N`$, corresponding to the N mobile species, and $`N_z`$ basis functions, defined on the interval in the $`z`$-direction, for example given by $`\{\mathrm{sin}m\pi z/L_z\}`$ for Dirichlet boundary conditions, where $`m=1,\mathrm{},N_z`$. Similarly, $`S_I`$ would be spanned by $`\widehat{e}_{N+1}\{\mathrm{sin}m\pi z/L_z\}`$. We make the physical assumption that no instability occurs in the immobile subsystem, $`𝚿^M0`$.

A Turing bifurcation occurs when a single $`\lambda _k`$ goes through zero and becomes positive as the parameter $`𝜶`$ goes through the critical value $`𝜶_𝒄`$. These conditions

$$\lambda _{k_c}=\frac{d\lambda _k}{dk}|_{k=k_c}=0$$
(33)
are satisfied at the onset of the instability, and define the critical wavenumber $`k_c`$. Note that in the case of scalar diffusion, no such instability can occur in the absence of reaction with the immobile species. From the above,

$$𝓛_{k_c}𝚿_{k_c}=0,$$
(34)
and

$$\frac{𝓛}{k}|_{k_c}𝚿_{k_c}+𝓛_{k_c}\frac{𝚿}{k}|_{k_c}=0.$$
(35)
Taking the inner product of Eq. 35 on the left by $`𝚿_{k_c}^{}`$, the nullvector of $`𝓛_{k_c}^{}`$, we find

$$𝚿_{k_c}^{M}{}_{}{}^{}𝚿_{k_c}^M=0,$$
(36)
for $`k_c0`$. In the absence of the immobile reaction, this implies that zero is a defective eigenvalue of $`𝓛_{k_c}`$ at the Turing bifurcation. However, it follows that

$$𝚿_{k_c}^{}𝚿_{k_c}=\mathrm{\Psi }_{k_c}^{I}{}_{}{}^{}\mathrm{\Psi }_{k_c}^I0$$
(37)
since there must exist a nontrivial solution $`\mathrm{\Psi }^I`$ to

$$g_𝐔𝚿_{k_c}^M+g_W\mathrm{\Psi }_{k_c}^I=0.$$
(38)
(Otherwise a nontrivial null vector $`𝚿_{k_c}`$ would not exist.) Hence, zero is not a defective eigenvalue of $`𝓛_{k_c}`$. Furthermore, it can be easily shown that the singularity of $`𝓛_k`$ as a function of parameters $`𝜶`$ and wavenumber $`k`$ is independent of the immobile reaction .

In this way, the structure of the linear stability of the $`(N+1)`$-component reaction-diffusion system at the Turing bifurcation motivated the formulation of a general scheme for appropriately lifting a defective eigenvalue problem to a higher dimensional space in which it is no longer defective. Despite this connection, our numerical examples of the lifting algorithm in the previous section did not include linear stability analysis of the CDIMA system. This is because the improvement in accurate calculation of the eigenvector associated with the defective eigenvalue is best illustrated when the eigenvalue is exactly known. For the gradient CDIMA system, the defective eigenvalue $`\mu _cDk_c^2`$ is not known apriori.

## 5 Concluding Remarks

The stability analysis of the experimental CDIMA reaction-diffusion system has led to the formulation of a new scheme in computational linear algebra for more accurate numerical solution of defective eigenvectors. We have demonstrated improvement in numerical accuracy by several orders of magnitude using illustrative toy examples. This scheme should be of utility in physical applications where accurate eigenvectors of defective matrices are sought, as well as of general use in improving accuracy of matrix computations involving such matrices.