smash::ChemicalPotentialSolver Class Reference

#include <chemicalpotential.h>

A class which encapsulates a GSL algorithm for finding the effective chemical potential and supporting functions.

Definition at line 24 of file chemicalpotential.h.

Collaboration diagram for smash::ChemicalPotentialSolver:

Classes
struct	ParametersForChemPotSolver
	Struct, root equations, and procedure for finding the effective chemical potential for a given particle species. More...

Public Member Functions
double	effective_chemical_potential (double degeneracy, double mass, double number_density, double temperature, double statistics, double solution_precision)
	Convenience wrapper for finding the effective chemical potential for a given particle species and performing sanity checks on the result. More...

Static Public Member Functions
static double	density_one_species (double degeneracy, double mass, double temperature, double effective_chemical_potential, double statistics, Integrator *integrator)
	Vector number density of one particle species, obtained through integrating the Juttner distribution function over the momentum space. More...
static double	root_equation_effective_chemical_potential (double degeneracy, double mass, double number_density, double temperature, double effective_chemical_potential, double statistics, Integrator *integrator)
	Root equation for finding the value of the effective chemical potential for one particle species. More...
static int	root_equation_effective_chemical_potential_for_GSL (const gsl_vector *roots_array, void *parameters, gsl_vector *function)
	Root equation for finding the value of the effective chemical potential for one particle species, suited for the GSL root finding procedure. More...
static void	print_state_effective_chemical_potential (unsigned int iter, gsl_multiroot_fsolver *solver)
	A GSL utility which allows for printing out the status of the solver during the root finding procedure. More...
static int	find_effective_chemical_potential (double degeneracy, double mass, double number_density, double temperature, double statistics, double mu_initial_guess, double solution_precision, Integrator *integrator, double *effective_chemical_potential)
	A GSL root solver for finding the effective chemical potential. More...

Private Attributes
Integrator	integrator_ = Integrator(1E6)
	A wrapper for gsl numerical integration. More...

Member Function Documentation

density_one_species()

double smash::ChemicalPotentialSolver::density_one_species ( double  degeneracy, double  mass, double  temperature, double  effective_chemical_potential, double  statistics, Integrator *  integrator  )

static

Vector number density of one particle species, obtained through integrating the Juttner distribution function over the momentum space.

Vector number density of one particle species.

Parameters

[in]	degeneracy	degeneracy g of the particle species
[in]	mass	(pole) mass m of the particle species [GeV]
[in]	temperature	temperature T of the system [GeV]
[in]	effective_chemical_potential	effective chemical potential of the system [GeV]
[in]	statistics	quantum statistics of the particles species (+1 for Fermi, -1 for Bose, 0 for Boltzmann)
[in]	integrator	wrapper for gsl numerical integration

Returns

the vector number density for a given species of particles [GeV^3]

Definition at line 29 of file chemicalpotential.cc.

                            {
  const auto integral = (*integrator)(0.0, 1.0, [&](double x) {
    return ((1.0 - x) * (1.0 - x) / (x * x * x * x)) *
           juttner_distribution_func((1.0 - x) / x, mass, temperature,
                                     effective_chemical_potential, statistics);
  });
  const double integrated_number_density =
      (degeneracy / (2.0 * M_PI * M_PI)) * integral;
 
  return integrated_number_density;
}

Here is the call graph for this function:

Here is the caller graph for this function:

root_equation_effective_chemical_potential()

double smash::ChemicalPotentialSolver::root_equation_effective_chemical_potential ( double  degeneracy, double  mass, double  number_density, double  temperature, double  effective_chemical_potential, double  statistics, Integrator *  integrator  )

static

Root equation for finding the value of the effective chemical potential for one particle species.

Parameters

[in]	degeneracy	degeneracy g of the particle species
[in]	mass	(pole) mass m of the particle species
[in]	number_density	number density n of the particle species [GeV^3]
[in]	temperature	temperature T of the system [GeV]
[in]	effective_chemical_potential	effective chemical potential of the system [GeV]
[in]	statistics	quantum statistics of the particles species (+1 for Fermi, -1 for Bose, 0 for Boltzmann)
[in]	integrator	wrapper for gsl numerical integration

Returns

the root equation

Definition at line 49 of file chemicalpotential.cc.

                            {
  const double integrated_number_density =
      density_one_species(degeneracy, mass, temperature,
                          effective_chemical_potential, statistics, integrator);
 
  return number_density - integrated_number_density;
}

Here is the call graph for this function:

Here is the caller graph for this function:

root_equation_effective_chemical_potential_for_GSL()

int smash::ChemicalPotentialSolver::root_equation_effective_chemical_potential_for_GSL ( const gsl_vector *  roots_array, void *  parameters, gsl_vector *  function  )

static

Root equation for finding the value of the effective chemical potential for one particle species, suited for the GSL root finding procedure.

Parameters

[in]	roots_array	an array holding the current best estimate of the roots of the solved equation
[in]	parameters	refers to the parameters as provided in the GSL root solving procedure
[in]	function	refers to the root equation(s) as provided in the GSL root solving procedure (where it's called "function")

Returns

the root equation suited for GSL root finding procedure

Definition at line 60 of file chemicalpotential.cc.

                                                                           {
  struct ParametersForChemPotSolver *par =
      static_cast<ParametersForChemPotSolver *>(parameters);
 
  const double degeneracy = (par->degeneracy);
  const double mass = (par->mass);
  const double number_density = (par->number_density);
  const double temperature = (par->temperature);
  const double statistics = (par->statistics);
  Integrator *integrator = (par->integrator);
 
  const double effective_chemical_potential = gsl_vector_get(roots_array, 0);
 
  gsl_vector_set(function, 0,
                 root_equation_effective_chemical_potential(
                     degeneracy, mass, number_density, temperature,
                     effective_chemical_potential, statistics, integrator));
 
  return GSL_SUCCESS;
}

Here is the call graph for this function:

Here is the caller graph for this function:

print_state_effective_chemical_potential()

void smash::ChemicalPotentialSolver::print_state_effective_chemical_potential ( unsigned int  iter, gsl_multiroot_fsolver *  solver  )

static

A GSL utility which allows for printing out the status of the solver during the root finding procedure.

Parameters

[in]	iter	variable keeping track of how many steps in the root solving procedure have been taken
[in]	solver	GSL solver object, which has access to the current best estimate of the roots and the corresponding function values

Returns

message about the current state of the solver

Definition at line 82 of file chemicalpotential.cc.

                                                      {
  printf(
      "\n***\nfind_effective_chemical_potential(): iter = %3u \t"
      "x = % .3f \t"
      "f(x) = % .3e \n",
      iter, gsl_vector_get(solver->x, 0), gsl_vector_get(solver->f, 0));
}

find_effective_chemical_potential()

int smash::ChemicalPotentialSolver::find_effective_chemical_potential ( double  degeneracy, double  mass, double  number_density, double  temperature, double  statistics, double  mu_initial_guess, double  solution_precision, Integrator *  integrator, double *  effective_chemical_potential  )

static

A GSL root solver for finding the effective chemical potential.

In practice one should use the convenience wrapper effective_chemical_potential(), which also performs sanity checks on the obtained solution.

Solver of the equation for chemical potential \( \mu \)

\[ n = \frac{g}{2 \pi^2} \int \frac{p^2 dp}{ e^{(\sqrt{p^2 + m^2} - \mu)/T} \pm 1} \]

given the density \( n \), and temperature \( T \).

Parameters

[in]	degeneracy	degeneracy g of the particle species
[in]	mass	m (pole) mass m of the particle species [GeV]
[in]	number_density	number density n of the particle species n [GeV^3]
[in]	temperature	temperature T of the system in GeV
[in]	statistics	quantum statistics of the particles species (+1 for Fermi, -1 for Bose, 0 for Boltzmann)
[in]	mu_initial_guess	the initial guess for the value of the effective chemical potential [GeV]
[in]	solution_precision	precision with which the solution is found; note that this precision also goes into the precision of integrals calculated in the process
[in]	integrator	wrapper for gsl numerical integration
[out]	effective_chemical_potential	the solution stored in an array object (we use an array in anticipation of generalizing to multidimensional root finding, needed for example when scalar interactions are present and effective mass is calculated at the same time as the effective chemical potential)

Definition at line 91 of file chemicalpotential.cc.

                                                                  {
  const gsl_multiroot_fsolver_type *Solver_name;
  gsl_multiroot_fsolver *Root_finder;
 
  int status;
  size_t iter = 0;
 
  const size_t problem_dimension = 1;
 
  struct ParametersForChemPotSolver parameters = {
      degeneracy, mass, number_density, temperature, statistics, integrator};
 
  gsl_multiroot_function EffectiveChemicalPotential = {
      &(ChemicalPotentialSolver::
            root_equation_effective_chemical_potential_for_GSL),
      problem_dimension, &parameters};
 
  double roots_array_initial[problem_dimension] = {mu_initial_guess};
 
  gsl_vector *roots_array = gsl_vector_alloc(problem_dimension);
  gsl_vector_set(roots_array, 0, roots_array_initial[0]);
 
  Solver_name = gsl_multiroot_fsolver_hybrids;
  Root_finder = gsl_multiroot_fsolver_alloc(Solver_name, problem_dimension);
  gsl_multiroot_fsolver_set(Root_finder, &EffectiveChemicalPotential,
                            roots_array);
 
  // print_state_effective_chemical_potential (iter, Root_finder);
 
  do {
    iter++;
 
    /*
     * gsl_multiroot_fsolver_iterate returns either 0 for a correct behavior,
     * or an error code (a positive integer) when the solver is stuck
     */
    status = gsl_multiroot_fsolver_iterate(Root_finder);
 
    // print_state_effective_chemical_potential (iter, Root_finder);
 
    /*
     * Check whether the solver is stuck
     */
    if (status) {
      std::cout << "\nGSL error message:\n"
                << gsl_strerror(status) << "\n\n"
                << std::endl;
      const double conversion_factor =
          smash::hbarc * smash::hbarc * smash::hbarc;
      logg[LogArea::Main::id].warn(
          "\n\nThe GSL solver\nfind_effective_chemical_potential\nis stuck!"
          "\n\nInput parameters:"
          "\n            degeneracy = ",
          degeneracy, "\n            mass [GeV] = ", mass,
          "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
          "\n     temperature [GeV] = ", temperature,
          "\n            statistics = ", statistics,
          "\n    solution_precision = ", solution_precision,
          "\n\n"
          "Initialization cannot continue without calculating the "
          "chemical potential. \nTry adjusting the initial "
          "guess or the solution precision.\n\n\n");
      throw std::runtime_error(
          "Chemical potential calculation returned no result.\n\n");
      continue;
    }
 
    status = gsl_multiroot_test_residual(Root_finder->f, solution_precision);
 
    if (status == GSL_SUCCESS) {
      effective_chemical_potential[0] = gsl_vector_get(Root_finder->x, 0);
    }
  } while (status == GSL_CONTINUE && iter < 100000);
 
  gsl_multiroot_fsolver_free(Root_finder);
  gsl_vector_free(roots_array);
 
  return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

effective_chemical_potential()

double smash::ChemicalPotentialSolver::effective_chemical_potential	(	double	degeneracy,
		double	mass,
		double	number_density,
		double	temperature,
		double	statistics,
		double	solution_precision
	)

Convenience wrapper for finding the effective chemical potential for a given particle species and performing sanity checks on the result.

Parameters

[in]	degeneracy	degeneracy g of the particle species
[in]	mass	(pole) mass m of the particle species [GeV]
[in]	number_density	number density n of the particle species [GeV^3]
[in]	temperature	temperature T of the system [GeV]
[in]	statistics	quantum statistics of the particles species (+1 for Fermi, -1 for Bose, 0 for Boltzmann)
[in]	solution_precision	precision with which the solution is found; note that this precision also goes into the precision of integrals calculated in the process

Returns

effective chemical potential mu, the solution [GeV]

Integration precision should be better than solving precision

value of small_p in GeV; value is chosen to be small (arbitrary choice)

Definition at line 179 of file chemicalpotential.cc.

                                                  {
  assert(temperature > 0);  // zero temperature case is not considered here
  assert(degeneracy > 0);
  assert(mass > 0);
  /*
   * The initial guess for the GSL chemical potential finding
   * procedure; the guess has to be different for Bose and Fermi
   * particles.
   */
  double initial_guess = 0.0;
  if (statistics == -1.0) {
    /*
     * Such initial guess allows one to sample really close to the Bose
     * condensate region; at the same time, it works fine in other regions.
     */
    initial_guess = 0.999999 * mass;
  } else {
    initial_guess = 1.05 * mass;
  }
 
  const double integration_precision = 0.1 * solution_precision;
  integrator_.set_precision(integration_precision, integration_precision);
 
  /*
   * Array for the GSL chemical potential finder.
   * Use of array is motivated by anticipating a generalization or overloading
   * of this solver to include effective mass calculation, which will lead to
   * multidimensional root finding, in which case the result is best passed as
   * an array.
   */
  double chemical_potential[1];
  chemical_potential[0] = 0.0;
 
  /*
   * Call the GSL effective chemical potential finding procedure.
   */
  find_effective_chemical_potential(
      degeneracy, mass, number_density, temperature, statistics, initial_guess,
      solution_precision, &integrator_, chemical_potential);
 
  /*
   * Sanity checks are performed for the obtained value of the chemical
   * potential .
   */
 
  /*
   * Check that the calculated chemical potential reproduces the input number
   * density within acceptable error.
   */
  // precision for number density in units of fm^-3
  const double precision_fm = 1e-3;
  const double conversion_factor = smash::hbarc * smash::hbarc * smash::hbarc;
  double nDensityCheck =
      density_one_species(degeneracy, mass, temperature, chemical_potential[0],
                          statistics, &integrator_);
 
  if (abs((number_density - nDensityCheck) / conversion_factor) >
      precision_fm) {
    logg[LogArea::Main::id].warn(
        "\n\nThe calculated chemical potential = ", chemical_potential[0],
        " [GeV] does not reproduce the input number density to within ",
        precision_fm,
        "!"
        "\nnumber_density [GeV^3] = ",
        number_density, "\n nDensityCheck [GeV^3] = ", nDensityCheck,
        "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
        "\n nDensityCheck [fm^-3] = ", nDensityCheck / (conversion_factor),
        "\n\n"
        "Input parameters:"
        "\n            degeneracy = ",
        degeneracy, "\n            mass [GeV] = ", mass,
        "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
        "\n     temperature [GeV] = ", temperature,
        "\n            statistics = ", statistics,
        "\n    solution_precision = ", solution_precision,
        "\n\n"
        "Initialization cannot continue without a correct "
        "calculation of the chemical potential."
        "\nTry adjusting the initial guess or the solution precision.\n\n\n");
    throw std::runtime_error(
        "Chemical potential calculation does not reproduce the "
        "input number density.\n\n");
  }
 
  /*
   * Check if the obtained value of the effective chemical potential indicates
   * that the number density is such that for bosons one has encountered the
   * Bose-Einstein condensate. We need mu < mass in order for the system NOT to
   * be in the Bose-Einstein condensate region.
   */
  if ((statistics == -1) && (chemical_potential[0] > mass)) {
    logg[LogArea::Main::id].warn(
        "\n\nThe calculated chemical potential indicates the "
        "input number density is such that Bose-Einstein "
        "condensate is encountered."
        "Input parameters:"
        "\n            degeneracy = ",
        degeneracy, "\n            mass [GeV] = ", mass,
        "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
        "\n     temperature [GeV] = ", temperature,
        "\n            statistics = ", statistics,
        "\n    solution_precision = ", solution_precision,
        "\n\n"
        "Initialization cannot proceed with a pathological "
        "distribution function leading to the Bose-Einstein "
        "condensate. \nTry increasing the temperature or "
        "decreasing the number density.\n\n\n");
    throw std::runtime_error(
        "Chemical potential calculation indicates that the Bose-Einstein "
        "condensate is produced.");
  }
 
  /*
   * Check if the distribution is pathological, i.e., has negative values
   * for low momenta (Bose statistics, in the Bose-Einstein condensate regime,
   *  may produce pathological distributions for which f(p) < 0 ).
   */
  const double small_p = 0.005;
  const double distribution_at_small_p = juttner_distribution_func(
      small_p, mass, temperature, chemical_potential[0], statistics);
  if (distribution_at_small_p < 0.0) {
    logg[LogArea::Main::id].warn(
        "\n\nNegative values of the distribution function at "
        "small momenta indicate that Bose-Einstein "
        "condensate is encountered."
        "Input parameters:"
        "\n            degeneracy = ",
        degeneracy, "\n            mass [GeV] = ", mass,
        "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
        "\n     temperature [GeV] = ", temperature,
        "\n            statistics = ", statistics,
        "\n    solution_precision = ", solution_precision, "\n\nf(p=", small_p,
        ") = ", distribution_at_small_p,
        "\n\n"
        "Initialization cannot proceed with a pathological "
        "distribution function leading to the Bose-Einstein "
        "condensate. \nTry increasing the temperature or "
        "decreasing the number density.\n\n\n");
    throw std::runtime_error(
        "Negative values of the distribution function at low momenta "
        "indicate that the Bose-Einstein condensate is produced.");
  }
 
  return chemical_potential[0];
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

integrator_

Integrator smash::ChemicalPotentialSolver::integrator_ = Integrator(1E6)

private

A wrapper for gsl numerical integration.

Definition at line 184 of file chemicalpotential.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/chemicalpotential.h
• src/chemicalpotential.cc