We present a new algorithm to evaluate the grand potential at high and finite temperatures using many-body perturbation theory. This algorithm enables us to calculate the contribution of any Hugenholtz or Feynman vacuum diagrams and formulate the results as a sum of divided differences. Additionally, the proposed method is applicable to any interaction in any dimension, allowing us to calculate thermodynamic quantities efficiently at any given temperature, particularly at high temperatures.Furthermore, we apply this algorithm to the Heisenberg spin-1/2 XXZ chain. We obtain all coefficients of the high-temperature expansion of the free energy and susceptibility per site of this model up to the sixth order.

How to run the code:

1. Copy all files in the same directory.
2. Open the Mathematica file "The free energy of the Heisenberg XXZ spin chain at high temperature.nb".
3. Execute the code by pressing "Shift+Enter".
4. The methods and algorithms applied here are explained in detail in my paper below: https://www.researchgate.net/publication/380518708_Finite_and_high-temperature_series_expansion_via_many-body_perturbation_theory_application_to_Heisenberg_spin-12_XXZ_chain