eos_quark_cfl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file eos_quark_cfl.h
24  \brief File defining \ref o2scl::eos_quark_cfl
25 */
26 #ifndef CFL_NJL_EOS_H
27 #define CFL_NJL_EOS_H
28 
29 #include <iostream>
30 #include <complex>
31 
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/vector_proxy.hpp>
34 #include <boost/numeric/ublas/matrix.hpp>
35 #include <boost/numeric/ublas/matrix_proxy.hpp>
36 
37 #include <gsl/gsl_math.h>
38 #include <gsl/gsl_eigen.h>
39 #include <gsl/gsl_complex_math.h>
40 #include <gsl/gsl_complex.h>
41 #include <gsl/gsl_sort_vector.h>
42 #include <gsl/gsl_poly.h>
43 
44 #include <o2scl/constants.h>
45 #include <o2scl/err_hnd.h>
46 #include <o2scl/mm_funct.h>
47 #include <o2scl/inte_qng_gsl.h>
48 #include <o2scl/poly.h>
49 #include <o2scl/test_mgr.h>
50 #include <o2scl/columnify.h>
51 
52 #include <o2scl/part.h>
53 #include <o2scl/eos_quark_njl.h>
54 
55 #ifndef DOXYGEN_NO_O2NS
56 namespace o2scl {
57 #endif
58 
59  /** \brief Nambu Jona-Lasinio model with a schematic CFL
60  di-quark interaction at finite temperature
61 
62  The variable B0 must be set before use.
63 
64  The original Lagrangian is
65 
66  \f[
67  {\cal L} =
68  {\cal L}_{\mathrm{Dirac}} +
69  {\cal L}_{\mathrm{4-fermion}} +
70  {\cal L}_{\mathrm{6-fermion}} +
71  {\cal L}_{CSC1} +
72  {\cal L}_{CSC2}
73  \f]
74 
75  \f[
76  {\cal L}_{\mathrm{Dirac}} = \bar{q}_{i \alpha}
77  \left( i {\partial} \delta_{i j} \delta_{\alpha \beta} -
78  m_{i j} \delta_{\alpha \beta} - \mu_{i j,~\alpha \beta} \gamma^0
79  \right) q_{j \beta}
80  \f]
81 
82  \f[
83  {\cal L}_{\mathrm{4-fermion}} =
84  G_S \sum_{a=0}^8 \left[
85  \left( \bar{q} \lambda^a_f q \right)^2 +
86  \left( \bar{q} i \gamma_5 \lambda^a_f q \right)^2
87  \right]
88  \f]
89 
90  \f[
91  {\cal L}_{\mathrm{6-fermion}} =
92  - G_{D} \left[
93  {\mathrm det}_{i j} \, \bar{q}_{i \alpha} \left( 1 + i \gamma_5 \right)
94  q_{j \beta} +
95  {\mathrm det}_{i j} \, \bar{q}_{i \alpha} \left( 1 - i \gamma_5 \right)
96  q_{j \beta} \right] \delta_{\alpha \beta}
97  \f]
98 
99  \f[
100  {\cal L}_{CSC1} =
101  G_{DIQ} \sum_k \sum_{\gamma} \left[
102  \left(\bar{q}_{i \alpha} \epsilon_{i j k}
103  \epsilon_{\alpha \beta \gamma} q^C_{j \beta}\right)
104  \left(\bar{q}_{i^{\prime} \alpha^{\prime}}^C
105  \epsilon_{i^{\prime} j^{\prime} k} \epsilon_{\alpha^{\prime}
106  \beta^{\prime} \gamma} q_{j^{\prime} \beta^{\prime}}\right)\right]
107  \f]
108 
109  \f[
110  {\cal L}_{CSC2} =
111  G_{DIQ} \sum_k \sum_{\gamma} \left[
112  \left(\bar{q}_{i \alpha} i \gamma_5 \epsilon_{i j k}
113  \epsilon_{\alpha \beta \gamma} q^C_{j \beta}\right)
114  \left(\bar{q}_{i^{\prime} \alpha^{\prime}}^C i \gamma_5
115  \epsilon_{i^{\prime} j^{\prime} k} \epsilon_{\alpha^{\prime}
116  \beta^{\prime} \gamma} q_{j^{\prime} \beta^{\prime}}\right) \right] \,,
117  \f]
118 
119  where \f$ \mu \f$ is the \quark number chemical potential.
120  couplings \f$ G_S \f$, \f$ G_D \f$, and \f$ G_{DIQ} \f$
121  ultra-violet three-momentum cutoff, \f$ \Lambda \f$
122 
123  The thermodynamic potential is
124  \f[
125  \Omega(\mu_i,\left<\bar{q} q\right>_i,\left< q q\right>_i,T)
126  = \Omega_{\mathrm{vac}}+\Omega_{\mathrm{stat}} + \Omega_0
127  \f]
128  where \f$ i \f$ runs over all nine (three colors times three
129  flavors) quarks. We assume that the condensates are independent
130  of color and
131  that the \quark chemical potentials
132  are of the form
133  \f$ \mu_Q=\mu_{\mathrm{Flavor(Q)}}+\mu_{\mathrm{Color(Q)}} \f$
134  with
135  \f[
136  \mu_{\mathrm{red}} = \mu_3 + \mu_8/\sqrt{3}
137  \quad
138  \mu_{\mathrm{green}} = -\mu_3 + \mu_8/\sqrt{3}
139  \quad
140  \mu_{\mathrm{blue}} = -2 \mu_8 /\sqrt{3}
141  \f]
142  With these assumptions, the thermodynamic potential as given
143  by the function thd_potential(), is a function of 12 variables
144  \f[
145  \Omega(\mu_u, \mu_d, \mu_s, \mu_3, \mu_8, \left<\bar{u} u\right>,
146  \left<\bar{d} d\right>, \left<\bar{s} s\right>,
147  \left< u d\right>, \left< u s\right>, \left< d s\right>,
148  T)
149  \f]
150 
151  The individual terms are
152  \f[
153  \Omega_{\mathrm{stat}} = - \frac{1}{2} \int \frac{d^3
154  p}{\left(2 \pi\right)^3} \, \sum_{i=1}^{72} \left[ \frac{\lambda_i}{2} +
155  T \ln{\left(1 + e^{-\lambda_i/T} \right)} \right]
156  \f]
157 
158  \f[
159  \Omega_{\mathrm{vac}} =
160  - 2 G_S \sum_{i=u,d,s} \langle {\bar q_i} q_i \rangle^2
161  +4 G_D \left<{\bar u} u\right> \left<{\bar d} d \right>
162  \left<{\bar s} s\right>
163  + \sum_k \sum_{\gamma} \frac{\left|\Delta^{k \gamma}\right|^2}{4 G_{DIQ}}
164  \f]
165 
166  where \f$ \lambda_i \f$ are the eigenvalues of the (72 by 72) matrix
167  (calculated by the function eigenvalues())
168  \f[
169  D = \left[
170  \begin{array}{cc}
171  - \gamma^0 \vec{\gamma} \cdot \vec{p} - M_{i} \gamma^0 + \mu_{i \alpha}
172  & \Delta i \gamma^0 \gamma_5 C \\
173  i \Delta \gamma^0 C \gamma_5
174  & - \gamma^0 \vec{\gamma}^T \cdot \vec{p} + M_{i} \gamma^0
175  - \mu_{i \alpha}\end{array}
176  \right]
177  \f]
178  and \f$ C \f$ is the charge conjugation matrix (in the Dirac
179  representation).
180 
181  The values of the various condensates are usually determined by
182  the condition
183  \f[
184  \frac{\partial \Omega}{\left<\bar{q} q\right>_i} = 0
185  \quad
186  \frac{\partial \Omega}{\left<q q\right>_i} = 0
187  \f]
188 
189  Note that setting fixed_mass to \c true and setting all of the
190  gaps to zero when \c gap_limit is less than zero will reproduce
191  an analog of the bag model with a momentum cutoff.
192 
193  The variable eos_quark_njl::fromqq is automatically set to true in
194  the constructor, as computations with \c fromqq=false are not
195  implemented.
196 
197  \future This class internally mixes ubvector, ubmatrix, gsl_vector
198  and gsl_matrix objects in a confusing and non-optimal way. Fix this.
199  \future Allow user to change derivative object? This isn't possible
200  right now because the stepsize parameter of the derivative
201  object is used.
202 
203  \hline
204  <b>References:</b>
205 
206  Created for \ref Steiner02.
207  */
208  class eos_quark_cfl : public eos_quark_njl {
209  public:
210 
211  eos_quark_cfl();
212 
213  virtual ~eos_quark_cfl();
214 
215  /** \brief Set the parameters and the bag constant 'B0'
216 
217  This function allows the user to specify the momentum cutoff,
218  \c lambda, the four-fermion coupling \c fourferm, the
219  six-fermion coupling from the 't Hooft interaction \c sixferm,
220  and the color-superconducting coupling, \c fourgap. If 0.0 is
221  given for any of the values, then the default is used (\f$
222  \Lambda=602.3/(\hbar c), G=1.835/\Lambda^2, K=12.36/\Lambda^5
223  \f$).
224 
225  If the four-fermion coupling that produces a gap is not
226  specified, it is automatically set to 3/4 G, which is the
227  value obtained from the Fierz transformation.
228 
229  The value of the shift in the bag constant eos_quark_njl::B0 is
230  automatically calculated to ensure that the vacuum has zero
231  energy density and zero pressure. The functions set_quarks()
232  and set_thermo() must be used before hand to specify the \ref
233  quark and \ref thermo objects.
234 
235  */
236  virtual int set_parameters(double lambda=0.0, double fourferm=0.0,
237  double sixferm=0.0, double fourgap=0.0);
238 
239  /** \brief Calculate the EOS
240 
241  Calculate the EOS from the quark condensates in \c u.qq, \c
242  d.qq and \c s.qq. Return the mass gap equations in \c qq1, \c
243  qq2, \c qq3, and the normal gap equations in \c gap1, \c gap2,
244  and \c gap3.
245 
246  Using \c fromqq=false as in eos_quark_njl and eos_quark_njl does not
247  work here and will return an error. Also, the quarks must be
248  set through eos_quark::quark_set() before use.
249 
250  If all of the gaps are less than gap_limit, then the
251  eos_quark_njl::calc_temp_p() is used, and \c gap1, \c gap2, and
252  \c gap3 are set to equal \c u.del, \c d.del, and \c s.del,
253  respectively.
254 
255  \todo It surprises me that n3 is not -res[11]. Is there
256  a sign error in the color densities?
257  */
258  virtual int calc_eq_temp_p(quark &u, quark &d, quark &s, double &qq1,
259  double &qq2, double &qq3, double &gap1,
260  double &gap2, double &gap3, double mu3,
261  double mu8, double &n3, double &n8,
262  thermo &qb, double temper);
263 
264  /// Check the derivatives specified by eigenvalues()
265  virtual int test_derivatives(double lmom, double mu3, double mu8,
266  test_mgr &t);
267 
268  /** \brief Calculate the energy eigenvalues as a function of the
269  momentum
270 
271  Given the momentum \c mom, and the chemical potentials
272  associated with the third and eighth gluons (\c mu3 and \c
273  mu8), the energy eigenvalues are computed in egv[0]
274  ... egv[35].
275  */
276  virtual int eigenvalues(double lmom, double mu3, double mu8,
277  double egv[36], double dedmuu[36],
278  double dedmud[36], double dedmus[36],
279  double dedmu[36], double dedmd[36],
280  double dedms[36], double dedu[36],
281  double dedd[36], double deds[36],
282  double dedmu3[36], double dedmu8[36]);
283 
284  /// Set the routine for solving quartics
286  quartic=&q;
287  return 0;
288  }
289 
290  /// The equal mass threshold
291  double eq_limit;
292 
293  /// Set to true to test the integration (default false)
295 
296  /// Test the integration routines
297  int test_integration(test_mgr &t);
298 
299  //
300  //gsl_quadratic_real_coeff quad;
301 
302  /** \brief The default quartic routine
303 
304  Slightly better accuracy (with slower execution times) can be
305  achieved using \ref poly_real_coeff_gsl which polishes the
306  roots of the quartics. For example
307 
308  \code
309  eos_quark_cfl cfl;
310  poly_real_coeff_gsl gp;
311  cfl.set_quartic(gp);
312  \endcode
313  */
315  //@}
316 
317  /** \brief Test the routine to compute the eigenvalues of
318  non-superfluid fermions
319  */
321 
322  /** \brief Test the routine to compute the eigenvalues of
323  superfluid fermions
324  */
326 
327  /** \brief Smallest allowable gap (default 0.0)
328 
329  If any of the gaps are below this value, then it is assumed
330  that they are zero and the equation of state is simplified
331  accordingly. If all of the gaps are less than gap_limit, then
332  the results from eos_quark_njl are used in
333  calc_eq_temp_p(), calc_temp_p() and thd_potential().
334  */
335  double gap_limit;
336 
337  /** \brief If this is true, then finite temperature
338  corrections are ignored (default false)
339 
340  This implements some simplifications in the
341  momentum integration that are not possible at finite temperature.
342  */
343  bool zerot;
344 
345  /** \brief Use a fixed quark mass and ignore the quark
346  condensates
347  */
349 
350  /// If true, then ensure color neutrality
352 
353  /** \brief Diquark coupling constant (default 3 G/4)
354 
355  The default value is the one derived from a Fierz
356  transformation. (\ref Buballa04)
357  */
358  double GD;
359 
360  /** \brief The absolute precision for the integration
361  (default \f$ 10^{-4} \f$ )
362 
363  This is analogous to gsl_inte::epsabs
364  */
365  double inte_epsabs;
366 
367  /** \brief The relative precision for the integration
368  (default \f$ 10^{-4} \f$ )
369 
370  This is analogous to gsl_inte::epsrel
371  */
372  double inte_epsrel;
373 
374  /** \brief The number of points used in the last integration
375  (default 0)
376 
377  This returns 21, 43, or 87 depending on the number of function
378  evaluations needed to obtain the desired precision. If it
379  the routine failes to obtain the desired precision, then
380  this variable is set to 88.
381  */
382  size_t inte_npoints;
383 
384  /// Return string denoting type ("eos_quark_cfl")
385  virtual const char *type() { return "eos_quark_cfl"; };
386 
387 #ifndef DOXYGEN_INTERNAL
388 
389  protected:
390 
391  /** \brief The integrands
392 
393  - res[0] is the thermodynamic potential, \f$ \Omega \f$
394  - res[1] is \f$ d -\Omega / d T \f$
395  - res[2] is \f$ d \Omega / d \mu_u \f$
396  - res[3] is \f$ d \Omega / d \mu_d \f$
397  - res[4] is \f$ d \Omega / d \mu_s \f$
398  - res[5] is \f$ d \Omega / d m_u \f$
399  - res[6] is \f$ d \Omega / d m_d \f$
400  - res[7] is \f$ d \Omega / d m_s \f$
401  - res[8] is \f$ d \Omega / d \Delta_{ds} \f$
402  - res[9] is \f$ d \Omega / d \Delta_{us} \f$
403  - res[10] is \f$ d \Omega / d \Delta_{ud} \f$
404  - res[11] is \f$ d \Omega / d \mu_3 \f$
405  - res[12] is \f$ d \Omega / d \mu_8 \f$
406  */
407  virtual int integrands(double p, double res[]);
408 
409  /// Compute ungapped eigenvalues and the appropriate derivatives
410  int normal_eigenvalues(double m, double lmom, double mu,
411  double lam[2], double dldmu[2],
412  double dldm[2]);
413 
414  /** \brief Treat the simply gapped quarks in all cases gracefully
415 
416  This function uses the quarks \c q1 and \c q2 to construct the
417  eigenvalues of the inverse propagator, properly handling the
418  either zero or finite quark mass and either zero or finite
419  quark gaps. In the case of finite quark mass and finite quark
420  gaps, the quartic solver is used.
421 
422  The chemical potentials are separated so we can add the
423  color chemical potentials to the quark chemical potentials
424  if necessary.
425 
426  This function is used by eigenvalues(). It does not work for
427  the "ur-dg-sb" set of quarks which are paired in a non-trivial
428  way.
429 
430  \todo In the code, the equal mass case seems to be commented
431  out. Why?
432  */
433  int gapped_eigenvalues(double m1, double m2, double lmom,
434  double mu1, double mu2, double tdelta,
435  double lam[4], double dldmu1[4],
436  double dldmu2[4], double dldm1[4],
437  double dldm2[4], double dldg[4]);
438 
439  /// Temperature
440  double temper;
441 
442  /// 3rd gluon chemical potential
443  double smu3;
444 
445  /// 8th gluon chemical potential
446  double smu8;
447 
448  /// \name Numerical methods
449  //@{
450  /// The routine to solve quartics
452 
458 
459  /// \name For computing eigenvalues
460  //@{
461  /// Inverse propagator matrix
462  gsl_matrix_complex *iprop;
463 
464  /// The eigenvectors
465  gsl_matrix_complex *eivec;
466 
467  /// The derivative of the inverse propagator wrt the ds gap
468  ubmatrix_complex dipdgapu;
469  /// The derivative of the inverse propagator wrt the us gap
470  ubmatrix_complex dipdgapd;
471  /// The derivative of the inverse propagator wrt the ud gap
472  ubmatrix_complex dipdgaps;
473 
474  /// The eigenvalues
475  gsl_vector *eval;
476 
477  /// Workspace for eigenvalue computation
478  gsl_eigen_hermv_workspace *w;
479  //@}
480 
481  /// \name For the integration
482  //@{
483  /// The error scaling function for integ_err
484  double rescale_error(double err, double result_abs,
485  double result_asc);
486 
487  /** \brief A new version of inte_qng_gsl to integrate several
488  functions at the same time
489  */
490  int integ_err(double a, double b, const size_t nr,
491  ubvector &res, double &err2);
492  //@}
493 
494  private:
495 
496  eos_quark_cfl(const eos_quark_cfl &);
497  eos_quark_cfl& operator=(const eos_quark_cfl&);
498 
499 #endif
500 
501  };
502 
503 #ifndef DOXYGEN_NO_O2NS
504 }
505 #endif
506 
507 #endif
gsl_matrix_complex * iprop
Inverse propagator matrix.
gsl_vector * eval
The eigenvalues.
quartic_real_coeff_cern def_quartic
The default quartic routine.
Nambu Jona-Lasinio model with a schematic CFL di-quark interaction at finite temperature.
bool fixed_mass
Use a fixed quark mass and ignore the quark condensates.
int normal_eigenvalues(double m, double lmom, double mu, double lam[2], double dldmu[2], double dldm[2])
Compute ungapped eigenvalues and the appropriate derivatives.
int test_gapped_eigenvalues(test_mgr &t)
Test the routine to compute the eigenvalues of superfluid fermions.
virtual int integrands(double p, double res[])
The integrands.
int gapped_eigenvalues(double m1, double m2, double lmom, double mu1, double mu2, double tdelta, double lam[4], double dldmu1[4], double dldmu2[4], double dldm1[4], double dldm2[4], double dldg[4])
Treat the simply gapped quarks in all cases gracefully.
double gap_limit
Smallest allowable gap (default 0.0)
quartic_real_coeff * quartic
The routine to solve quartics.
virtual const char * type()
Return string denoting type ("eos_quark_cfl")
virtual int eigenvalues(double lmom, double mu3, double mu8, double egv[36], double dedmuu[36], double dedmud[36], double dedmus[36], double dedmu[36], double dedmd[36], double dedms[36], double dedu[36], double dedd[36], double deds[36], double dedmu3[36], double dedmu8[36])
Calculate the energy eigenvalues as a function of the momentum.
double GD
Diquark coupling constant (default 3 G/4)
virtual int set_parameters(double lambda=0.0, double fourferm=0.0, double sixferm=0.0, double fourgap=0.0)
Set the parameters and the bag constant &#39;B0&#39;.
ubmatrix_complex dipdgaps
The derivative of the inverse propagator wrt the ud gap.
int test_integration(test_mgr &t)
Test the integration routines.
gsl_matrix_complex * eivec
The eigenvectors.
virtual int test_derivatives(double lmom, double mu3, double mu8, test_mgr &t)
Check the derivatives specified by eigenvalues()
double rescale_error(double err, double result_abs, double result_asc)
The error scaling function for integ_err.
Nambu Jona-Lasinio EOS at zero temperature.
size_t inte_npoints
The number of points used in the last integration (default 0)
ubmatrix_complex dipdgapu
The derivative of the inverse propagator wrt the ds gap.
bool integ_test
Set to true to test the integration (default false)
int test_normal_eigenvalues(test_mgr &t)
Test the routine to compute the eigenvalues of non-superfluid fermions.
int integ_err(double a, double b, const size_t nr, ubvector &res, double &err2)
A new version of inte_qng_gsl to integrate several functions at the same time.
double temper
Temperature.
int set_quartic(quartic_real_coeff &q)
Set the routine for solving quartics.
ubmatrix_complex dipdgapd
The derivative of the inverse propagator wrt the us gap.
double smu8
8th gluon chemical potential
gsl_eigen_hermv_workspace * w
Workspace for eigenvalue computation.
double smu3
3rd gluon chemical potential
double inte_epsrel
The relative precision for the integration (default )
double eq_limit
The equal mass threshold.
bool color_neut
If true, then ensure color neutrality.
virtual int calc_eq_temp_p(quark &u, quark &d, quark &s, double &qq1, double &qq2, double &qq3, double &gap1, double &gap2, double &gap3, double mu3, double mu8, double &n3, double &n8, thermo &qb, double temper)
Calculate the EOS.
bool zerot
If this is true, then finite temperature corrections are ignored (default false)
double inte_epsabs
The absolute precision for the integration (default )

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).