You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

111 lines
5.5KB

  1. /*************************************************************************
  2. * Copyright (C) 2015 by Andrius Štikonas <andrius@stikonas.eu> *
  3. * *
  4. * This program is free software; you can redistribute it and/or modify *
  5. * it under the terms of the GNU General Public License as published by *
  6. * the Free Software Foundation; either version 3 of the License, or *
  7. * (at your option) any later version. *
  8. * *
  9. * This program is distributed in the hope that it will be useful, *
  10. * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  11. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  12. * GNU General Public License for more details. *
  13. * *
  14. * You should have received a copy of the GNU General Public License *
  15. * along with this program. If not, see <http://www.gnu.org/licenses/>.*
  16. *************************************************************************/
  17. // Calculation of mutual information in the setup of http://arxiv.org/abs/1503.08161
  18. #include <cmath>
  19. #include <complex>
  20. #include <iostream>
  21. static double alpha, beta;
  22. std::complex<double> crossRatio (std::complex<double>, std::complex<double>, std::complex<double>, std::complex<double>);
  23. double Fitzpatrick(std::complex<double>, std::complex<double>, double);
  24. int main()
  25. {
  26. // Parameters that can be changed:
  27. double tPlus = 0; // time on the left boundary
  28. double tMinus = 0; // time on the right boundary
  29. alpha = 0.4; // encodes the conformal dimention of local operator h_Psi
  30. double y = 3; // endpoint of interval A
  31. double L = 15; // length of interval A
  32. double epsilon = 0.01; // smearing parameter
  33. beta = 10; // inverse temperature
  34. // double tOmega = /*5.4418*/; // thermal state is perturbed by operator inserted at time -tOmega
  35. double c = 600; // central charge. Must be large in our approximation
  36. // End of parameters
  37. for(unsigned int i = 0; i < 1000; ++i)
  38. {
  39. double tOmega = i*2*L/1000;
  40. // Operator insertion points: Left boundary
  41. std::complex<double> x1 (0, -epsilon), x4 (0, epsilon), x1bar, x4bar;
  42. x1bar = conj(x1);
  43. x4bar = conj(x4);
  44. double x2 = y - tOmega - tMinus;
  45. double x2bar = y + tOmega + tMinus;
  46. double x3 = L + x2;
  47. double x3bar = L + x2bar;
  48. // Operator insertion points: Right boundary
  49. std::complex<double> x6 (y - tPlus - tOmega, beta/2);
  50. std::complex<double> x6bar (y + tPlus + tOmega, -beta/2);
  51. std::complex<double> x5 = L + x6;
  52. std::complex<double> x5bar = L + x6bar;
  53. // Cross-ratios for S_A
  54. std::complex<double> zA = crossRatio(x1, x2, x3, x4);
  55. std::complex<double> zAbar = crossRatio(x1bar, x2bar, x3bar, x4bar);
  56. // Cross-ratios for S_B
  57. std::complex<double> zB = crossRatio(x1, x5, x6, x4);
  58. std::complex<double> zBbar = crossRatio(x1bar, x5bar, x6bar, x4bar);
  59. // Cross-ratios for S_{A union B}
  60. std::complex<double> z2 = crossRatio(x1, x2, x6, x4);
  61. std::complex<double> z2bar = crossRatio(x1bar, x2bar, x6bar, x4bar);
  62. std::complex<double> z5 = crossRatio(x1, x5, x3, x4);
  63. std::complex<double> z5bar = crossRatio(x1bar, x5bar, x3bar, x4bar);
  64. // Now we calculate entanglement entropies using Fitzpatrick, Kaplan, Walters formula.
  65. double S_A = c/6 * log(Fitzpatrick(zA, zAbar, y < tMinus + tOmega && tMinus + tOmega < y + L ? 2*M_PI : 0));
  66. // std::cerr << "S_A = " << S_A << std::endl;
  67. // double S_A_analytic = c/6 * log(beta/M_PI/epsilon/sinh(M_PI*L/beta)*sinh(M_PI/beta*(y+L-tMinus-tOmega))*sinh(M_PI/beta*(tMinus+tOmega-y)) * sin(M_PI*alpha)/alpha);
  68. // double S_A_analytic = 0;
  69. // std::cerr << "S_A (analytic) = " << S_A_analytic << std::endl;
  70. double S_B = c/6 * log(Fitzpatrick(zB, zBbar, 0));
  71. double S_union = c/6 * log(Fitzpatrick(z2, z2bar, y < tMinus + tOmega ? 2*M_PI : 0) * Fitzpatrick(z5, z5bar, y +L < tMinus + tOmega ? -2*M_PI : 0));
  72. // std::cerr << "S_A+B = " << S_union << std::endl;
  73. double S_union_analytic = c/3 * log(beta/M_PI/epsilon * sin(M_PI*alpha)/alpha) + c/6* log(sinh(M_PI/beta*(tMinus + tOmega - y)) * cosh(M_PI/beta*(tPlus + tOmega - y)) * sinh(M_PI/beta*(tMinus + tOmega - y - L)) * cosh(M_PI/beta*(tPlus + tOmega - y - L)));
  74. // std::cerr << "S_A+B (analytic) = " << S_union_analytic << std::endl;
  75. double S_thermal = 2*c/3 * log(sinh(M_PI*L/beta)/cosh(M_PI/beta*(tMinus-tPlus)));
  76. double I = S_A + S_B - S_union + S_thermal;
  77. // double I_analytic = S_A_analytic - S_union_analytic + S_thermal;
  78. std::cout << tOmega << "\t" << I /*<< "\t" << I_analytic*/ << std::endl;
  79. }
  80. return 0;
  81. }
  82. std::complex<double> crossRatio (std::complex<double> x1, std::complex<double> x2, std::complex<double> x3, std::complex<double> x4)
  83. {
  84. double pb = M_PI/beta;
  85. return sinh(pb*(x1-x2))*sinh(pb*(x3-x4))/sinh(pb*(x1-x3))/sinh(pb*(x2-x4));
  86. }
  87. double Fitzpatrick(std::complex<double> z, std::complex<double> zbar, double phase)
  88. {
  89. // phase is necessary to take into account nontrivial monodromy
  90. std::complex<double> one = 1; // one as a complex number
  91. std::complex<double> i(0,1); // imaginary unit
  92. double alphaExpr = 0.5-alpha/2;
  93. std::complex<double> exponent1 = exp(phase*i*alphaExpr);
  94. std::complex<double> exponent2 = exp(phase*i*alpha);
  95. return real(exponent1 * pow(z, alphaExpr) * pow(zbar, alphaExpr)
  96. * (one - exponent2 * pow(z, alpha)) * (one - pow(zbar, alpha))
  97. / ( alpha*alpha * (one-z) * (one-zbar) ) );
  98. }