101 lines
3.4 KiB
C++
101 lines
3.4 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: ln_Density_ME.cc
|
|
|
|
Purpose: Computes the density operator \rho(x = 0) matrix element
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
/**
|
|
Given a index \f$j\f$ together with a bra and ket respectively parametrized
|
|
by sets of rapidities \f$\{ \mu \}\f$ and \f$\{ \lambda \}\f$,
|
|
this function returns the product (see e.g.
|
|
[J.-S. Caux et al, JSTAT P01008 (2007)](https://doi.org/10.1088/1742-5468/2007/01/P01008))
|
|
\f[
|
|
V_j^+ \equiv \prod_l
|
|
\frac{\tilde{\mu}_l - \tilde{\lambda}_j + i}{\tilde{\lambda}_l - \tilde{\lambda}_j + i}.
|
|
\f]
|
|
*/
|
|
complex<DP> Fn_V (int j, LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate)
|
|
{
|
|
complex<DP> result = 1.0;
|
|
|
|
for (int m = 0; m < lstate.N; ++m) {
|
|
result *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II)
|
|
/(rstate.lambdaoc[m] - rstate.lambdaoc[j] + II);
|
|
}
|
|
|
|
return(result);
|
|
}
|
|
|
|
/**
|
|
Computes the log of the density operator \f$\hat{\rho}(x = 0)\f$ matrix element
|
|
between lstate (bra) and rstate (ket). If lstate == rstate, density matrix element = N/L;
|
|
if momentum difference is zero but states are different, then matrix element is zero
|
|
(but this function then returns the log artificially set to -200.0).
|
|
*/
|
|
complex<DP> ln_Density_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate)
|
|
{
|
|
if (lstate.Ix2 == rstate.Ix2) return(log(lstate.N/lstate.L));
|
|
|
|
else if (lstate.iK == rstate.iK) return(-200.0); // so exp(.) is zero
|
|
|
|
SQMat_DP one_plus_U (0.0, lstate.N);
|
|
|
|
Vect_CX Vplus (0.0, lstate.N); // contains V^+_j; V^-_j is the conjugate
|
|
// prod_{m\neq j} (\mu_m - \lambdaoc_j)/(\lambdaoc_m - \lambdaoc_j):
|
|
Vect_CX ln_Fn_Prod (0.0, lstate.N);
|
|
Vect_DP rKern (0.0, lstate.N); // K(lambdaoc_j - lambdaoc_(p == arbitrary))
|
|
|
|
// "Phantom" rapidity in ME expression. Choice doesn't matter, see 1990_Slavnov_TMP_82
|
|
// after (3.8). Choose rapidity around the middle.
|
|
//int p = 0;
|
|
int p = rstate.N/2-1;
|
|
|
|
DP Kout = lstate.K - rstate.K;
|
|
|
|
for (int a = 0; a < lstate.N; ++a) {
|
|
Vplus[a] = Fn_V (a, lstate, rstate);
|
|
ln_Fn_Prod[a] = 0.0;;
|
|
for (int m = 0; m < lstate.N; ++m)
|
|
if (m != a) ln_Fn_Prod[a] += log(complex<DP>(lstate.lambdaoc[m] - rstate.lambdaoc[a])
|
|
/(rstate.lambdaoc[m] - rstate.lambdaoc[a]));
|
|
rKern[a] = rstate.Kernel (a, p);
|
|
}
|
|
|
|
for (int a = 0; a < lstate.N; ++a)
|
|
for (int b = 0; b < lstate.N; ++b)
|
|
one_plus_U[a][b] = (a == b ? 1.0 : 0.0)
|
|
+ 0.5 * ((lstate.lambdaoc[a] - rstate.lambdaoc[a])/imag(Vplus[a]))
|
|
* real(exp(ln_Fn_Prod[a])) * (rstate.Kernel(a,b) - rKern[b]);
|
|
|
|
complex<DP> ln_ddalpha_sigma = lndet_LU_dstry(one_plus_U);
|
|
|
|
complex<DP> ln_prod_V = 0.0;
|
|
for (int a = 0; a < lstate.N; ++a) ln_prod_V += log(2.0 * II * imag(Vplus[a]));
|
|
|
|
complex<DP> ln_prod_2 = 0.0;
|
|
for (int a = 0; a < lstate.N; ++a)
|
|
for (int b = 0; b < lstate.N; ++b)
|
|
ln_prod_2 += log((rstate.lambdaoc[a] - rstate.lambdaoc[b] + II)/(lstate.lambdaoc[a] - rstate.lambdaoc[b]));
|
|
|
|
ln_ddalpha_sigma += ln_prod_V + ln_prod_2 - log(2.0 * II * imag(Vplus[p]));
|
|
|
|
return (log(-II * Kout) + ln_ddalpha_sigma - 0.5 * (lstate.lnnorm + rstate.lnnorm));
|
|
}
|
|
|
|
}
|