# First principles study on the segregation of metallic solutes and non-
metallic impurities in
Cu grain boundary

Vasileios Fotopoulos<EMAIL_ADDRESS>[  Jack Strand  [  Manuel
Petersmann  [  Alexander L. Shluger  [

###### Abstract

Metallic dopants have the potential to increase the mechanical strength of
polycrystalline metals. These elements are expected to aggregate in regions of
lower coordination, such as grain boundaries. At the grain boundaries, they
can have a beneficial (toughening) or detrimental effect (e.g. grain boundary
embrittlement). In this study, we employ Density Functional Theory (DFT) to
compute the segregation energies of various metallic and other non-metallic
elements to determine their effect when introduced in a symmetric Cu grain
boundary. The study results may be used to qualitatively rank the beneficial
effect of certain metallic elements, such as V, Zr, and Ag, as well as the
strong weakening effect of non-metallic impurities like O, S, F and P.
Furthermore, the induced local distortion is found to be proportional to the
weakening effect of the elements.

###### keywords:

American Chemical Society, LaTeX

University College London] Department of Physics and Astronomy, University
College London, Gower Street, London WC1E 6BT, United Kingdom University
College London] Department of Physics and Astronomy, University College
London, Gower Street, London WC1E 6BT, United Kingdom [Nanolayers]Nanolayers
Research Computing Ltd., London, UK KAI] KAI–Kompetenzzentrum Automobil- und
Industrieelektronik GmbH, Europastrasse 8, 9524 Villach, Austria University
College London] Department of Physics and Astronomy, University College
London, Gower Street, London WC1E 6BT, United Kingdom IR,NMR,UV

##  1 Introduction

Degradation phenomena in metals have been reported to initiate at certain
microstructural features, such as grain boundaries (GBs) 1, 2, 3, 4.
Experimental results have illustrated how the level of corrosion and
degradation of metals under high-stress conditions depends on the
concentration of certain elements 5. Impurity-induced embrittlement accounts
for many cases of brittle failure of polycrystalline metals 6, 7, 8. Impurity-
induced embrittlement is commonly attributed to a chemical effect of the
segregation of solutes, which is believed to change the bonding strength at
GBs 6, 9, or to a size (strain) effect that is associated with the size
mismatch between the solute and host atoms 10, 11. The segregation of metallic
dopants and non-metallic impurity elements from GBs in metals plays a key role
in the design of novel materials 12. However, a clear correlation between the
properties of solutes and the strengthening or weakening effect that they can
have when introduced in the GBs of metals is still missing.

In light of recent experimental studies in bicrystalline Fe, solutes that
decorate grain boundaries significantly affect their chemical composition,
charge distribution, structural properties, and therefore their mechanical
stability 13. GBs in polycrystalline metals can be engineered and decorated
with certain elements to increase their resistance against corrosion effects.
Thus, identifying the elements that could potentially be employed as GB
decorative solutes is of great interest. First-principles simulations on 3d
transition metal solutes in Fe 14, 15, Au 12, Mo 16 and Ni 17 illustrated that
their presence at GBs can significantly increase their tensile strength and
resistance against embrittlement phenomena. Cu, a metal used in a wide range
of applications, including electronic devices, showed prominent degradation
effects under thermal cycling 18, 1. Such phenomena often initiate at crystal
grain boundaries and can have devastating effects on the performance of these
devices. Doping of grain boundaries with metallic elements can be an efficient
solution to this problem 19. However, although several properties have been
proposed to play a critical role in the weakening or strengthening effect that
these elements can have, a correlation between the above-mentioned effects and
the elements’ induced local distortion into the lattice hasn’t been
established.

This work complements previous theoretical studies of metallic segregants in
Cu GBs. Density functional theory (DFT) simulations are conducted to determine
the most favorable segregation sites of 3d transition and other metals in a
symmetric Cu GB, along with their local strengthening or weakening effect. We
highlight the beneficial effects of V, Zr, and Ag. Finally, we show how the
computed segregation energies are linked with the solutes’ induced lattice
relaxation. Metallic solutes with low GB segregation energies induce small
relaxation. On the other hand, non-metallic impurities with a high computed
strengthening energy, such as O, S, F, and P, cause severe local displacements
of Cu atoms. Such findings highlight the possibility of decorating Cu GBs with
specific metallic elements to reduce the degradation effects in Cu.

##  2 Methods

DFT calculations are performed using the Vienna Ab Initio Simulation Package
(VASP) 20, 21, 22, with the GGA exchange-correlation functional Perdew-Burke-
Ernzerhof (PBE) 23. For the investigation of grain boundaries and bulk
properties of copper, 76-atom and 108-atom periodic cells are utilized,
respectively. The specific grain boundary chosen for examination is the
(210)[100]$\Sigma 5$ twin boundary (shown in Figure 1(a)), due to its high
symmetry and low-energy characteristics 24. Consistent with previous studies
25, 26, the calculations employed converged 5$\times$4$\times$1 and
4$\times$4$\times$4 k-point grids for the 76 and 108-atom simulation cells,
respectively. An energy cut-off of 450 eV is used, as established in previous
Cu simulations 27.

Segregation energies were calculated using the formula:

$E_{seg}=(E_{GB+X}-E_{GB})-(E_{Bulk+X}-E_{Bulk}),$ (1)

where $E_{GB+X}$ and $E_{GB}$ are the energies of GB cells with and without
segregants (X) and $E_{Bulk+X}$, $E_{Bulk}$ the energies of bulk cells with
and without segregants, respectively. Negative and positive segregation
energies illustrate favorable GB segregation and antisegregation,
respectively. For the determination of the segregation energies related to
interstitial impurities, octahedral and tetrahedral interstitial sites are
considered in the bulk. Figure 1(a) presents the four interstitial and
substitutional positions considered at the GB for each segregant.

Figure 1: (a) Examined substitutional and interstitial sites of
dopants/impurities in (210)[100]$\Sigma 5$ Cu grain boundaries. (b) Schematics
illustrating (top) the segregation energy and (bottom) strengthening energies
of substitutional impurities/dopants in the grain boundaries and triple
points, respectively.

The effect of impurities on the strength of the grain boundary can be derived
from the strengthening energy. This parameter provides insight into whether
the solutes exhibit a preference for locating themselves at grain boundaries
as opposed to surface positions. The calculation of DFT strengthening energies
is carried out according to the following expression:

$E_{str}=(E_{GB+X}-E_{GB})-(E_{Sur+X}-E_{Sur}),$ (2)

where ESur and ESur+X symbolize the energies associated with the pure Cu
surface and the Cu surface simulation cells containing one solute atom,
respectively. For the surface simulations, 108-atom simulations (100) slab
models are used with a thickness of 10.86 Å. The strengthening energy, when
negative, indicates an enhancement of the grain boundary strength due to the
impurity, while a positive value implies a weakening effect. Schematic
representations of the segregation and strengthening energies are shown in
Figure 1(b).

##  3 Results

###  3.1 Segregation Energies

Figure 2: Computed segregation energies for (a) metallic and (b) non-metallic
elements for different sites. Different colored areas correspond to elements
of different groups. Negative energies correspond to favorable segregation at
the GB, whereas positive energies correspond to antisegregation. The int
represents interstitial atoms, whereas the rest of the elements correspond to
substitutional atoms.

To determine the most favorable sites, as shown in Figure 1(a), of the
impurities examined, the segregation energies of all elements are calculated.
Figure 2(a) includes all the computed segregation energies for the various
sites. The plots are divided into colored regions on the basis of the group of
elements on the periodic table. As seen in the figure, most of the metallic
dopants are found to occupy most preferably substitutional sites 1 and 3. Fe
and Co dopants show segregation energies close to 0 eV, indicating that these
elements do not show an energetic preference to be at the GB instead of the
bulk region.

For non-metallic impurities, our results show that all substitutional elements
will occupy site 2. Also, non-metallic elements show a significantly stronger
tendency to segregate at the GB compared to metallic solutes, with F, S, O,
and P resulting in the lowest segregation energies. C is the only element that
shows a strong preference to segregate as an interstitial atom. On the basis
of these results, elements such as C, O, F, P, S, and Se, when introduced in
polycrystalline Cu, are expected to be at considerably higher concentrations
at the GBs of the crystal.

Table 1: Segregation energies comparison between the results from the current work and previous DFT studies for metallic and non-metallic elements. Values with ∗ and ∗∗ correspond to segregation energies computed in (221)[110]$\Sigma$9 and (310)[001]$\Sigma$5 GB, respectively. The rest of the values, including our results, were obtained using (210)[100]$\Sigma$5 GB simulations cells. Element | Reference |  Eseg (eV) |  Eseg (eV) [Our work] |
---|---|---|---|---
C (int) |  Wurmshuber 28 | -0.8 | -0.85 |
C (int) |  Huang (2018) 19 | -0.9 | -0.85 |
O (int) |  Bodlos 29 | -1.2 | -1.27 |
O |  Bodlos 29 | -1.15 | -1.26 |
Si |  Huang (2018) 19 | -0.3 | -0.36 |
Si |  Li 30 |  -0.58∗∗ | -0.36 |
P |  Li 30 | -0.95 | -0.9 |
P |  Lousada 31 |  -0.65∗ | -0.9 |
S |  Bodlos 29 | -1.2 | -1.14 |
S |  Wang 32 | -1 | -1.14 |
S |  Li 30 |  -1.15∗∗ | -1.14 |
S |  Lousada 31 |  -0.75∗ | -1.14 |
Ge |  Razumovskiyy 33 | -0.3 | -0.42 |
Se |  Razumovskiyy 33 | -0.8 | -0.98 |
Mg |  Huang (2020) 34 |  -1∗∗ | -0.57 |
V |  Huang (2020) 34 |  -0.2∗∗ | -0.34 |
Co |  Razumovskiyy 33 | 0 | -0.04 |
Ag |  Razumovskiyy 33 | -0.6 | -0.5 |
Ag |  Huang (2019) 35 | -0.8 | -0.5 |
Zr |  Huang (2018) 19 | -1.65 | -1.23 |

Table 1 includes the comparison between the segregation energies calculated in
the current work and previous studies. For non-metallic impurities, the
computed segregation energies are within 0.2 eV compared to previous studies.
The only exceptions are seen in the case where a different $\Sigma$5 GB
((310)[001]$\Sigma$5 GB 30, 34), or a different $\Sigma$ symmetry
((221)[110]$\Sigma$9 31) is used. For the segregation energies of the metallic
dopants, our results agree well with the results of ref.33, while the study
ref.19 showed lower energies by 0.3 and 0.42 eV for Ag and Zr, respectively.
In general, the table shows that there is a clear trend for certain non-
metallic (C, O, P, S, Se) and metallic dopants (Zr) to strongly segregate at
the grain boundaries. The rest of the elements show small to close to no
preference for being at the grain boundaries instead of the bulk of Cu.

Figure 3: Computed strengthening energies for (a) metallic and (b) non-
metallic elements for different sites. Different colored areas correspond to
elements of different groups. Negative strengthening energies correspond to a
strengthening effect, whereas positive strengthening energies correspond to a
weakening effect.

###  3.2 Strengthening Energies

In Figure 3, the effect of the studied elements on the strength of the GB is
shown when introduced at various sites. Negative strengthening energies
correspond to a strengthening effect, whereas positive strengthening energies
correspond to a weakening effect. On the basis of these results, one can
deduce that most of the metallic solutes show negative strengthening energies
(strengthening effect), whereas all non-metallic elements show positive
strengthening energies (weakening effect). The latter shows that our findings
support previous studies and highlight the differences between metals and non-
metals when it comes to their impact on GBs 33, 19, 35.

Table 2 summarizes the main conclusions of our investigation so far. One can
see that all of the examined non-metallic elements are expected to have a
weakening effect when they are introduced into the GB of Cu. From the metallic
dopants, excluding Mg and Cd, all of the rest showed a strengthening effect.
The weakening effect of Mg and Cd can be attributed to their small atomic
radii. The grain boundary energy is expected to decrease as the atomic radius
of the metallic dopants increases 35. These results agree with the conclusions
drawn in ref.33, where metallic dopants with energies close to zero, in our
case Co and Fe, are expected to strengthen the GB. Also, in agreement with the
aforementioned study, non-metallic impurities relaxing to the site next to the
GB plane (site 2) weaken the GB. Furthermore, as highlighted in previous
theoretical works 34, non-metallic solutes with larger atomic radius segregate
more easily at the grain boundary but cause more significant grain boundary
embrittlement.

Table 2: Properties of the examined non-metallic impurities and metallic solutes along with the identified most favorable sites and their effect. W and S correspond to weakening and strengthening effects, respectively. Element | Group | Mass (amu) | Radius (pm) | Site | Effect (S/W) |
---|---|---|---|---|---|---
C | 2 | 12.01 | 170 | 3 (int) | W |
O | 2 | 15.99 | 152 | 4 (int) | W |
O | 2 | 15.99 | 152 | 2 | W |
F | 2 | 18.99 | 147 | 2 | W |
Si | 3 | 28.09 | 210 | 2 | W |
P | 3 | 30.97 | 195 | 2 | W |
S | 3 | 32.07 | 180 | 2 | W |
Ge | 4 | 72.64 | 211 | 3 | W |
Se | 4 | 78.96 | 190 | 2 | W |
Mg | 3 | 24.31 | 173 | 1 | W |
V | 4 | 50.94 | 179 | 4 | S |
Fe | 4 | 55.85 | 126 | 4 | S |
Co | 4 | 58.93 | 200 | 4 | S |
Zr | 5 | 91.22 | 206 | 3 | S |
Ag | 5 | 107.87 | 172 | 1/3 | S |
Cd | 5 | 112.41 | 158 | 1 | W |
Pt | 6 | 175 | 195.08 | 1 | S |

###  3.3 Lattice Distortion Effects

In the previous sections, the segregating and strengthening energies were
computed, and the most favorable segregation sites were identified. The
computed energetic parameters follow the pattern already identified by
previous studies 33, 19, 35. We now move on to examine the effect of
impurities/dopants on local lattice distortions. Figure 4(a) includes the
fully relaxed configurations for some of the tested metallic dopants along
with non-metallic elements at their most favorable sites. In addition,
displacements are also included, both in the form of colored atoms and as
displacement vectors. Relaxed configurations for O, Si, P, Zr and Pt are shown
in Figure 4(a). As can be deduced, substitutional metallic solutes cause
minimal relaxation. Non-metallic elements, either substitutional or
intersitial, induce significant local relaxation, which explains the weakening
effect that these elements have. Non-metallic impurities with a small
preference to segregate at the GB, like Si, caused negligible local
distortion. Figure 4(b) includes the results of the bond analysis. Bond
analysis illustrates the distribution of distances between the
dopants/impurities and nearby Cu atoms within a cut-off distance of 3.0 Å. Any
distances beyond the mentioned cut-off are ignored. All metallic dopants
showed distances from the adjacent Cu atoms of approximately 2.6 Å, which is
the same as the Cu-Cu distance in the perfect lattice. The latter highlights
the minimal distortion that the metallic dopants introduced. On the other
hand, non-metallic impurities show peaks at lower than 2.5 Å.

Figure 4: (a) Fully relaxed configurations along with their displacement
magnitudes (bottom) for O, Si, P, Zr, Pt. Yellow arrows correspond to the
displacement vectors. For displacements, the configuration prior to relaxation
is used as a reference. (b) Bond analysis for the aforementioned elements.
Bond analysis illustrates the distribution of distances between the
dopants/impurities and nearby Cu atoms within a cut-off distance of 3 Å. The
y-axis corresponds to the number of Cu atoms at a certain distance from the
dopant/impurity.

##  4 Conclusions

We employed DFT to investigate the energetic properties of metallic dopants
and non-metallic impurities introduced into the grain boundaries of Cu. Among
the metallic dopants, such as Mg, V, Fe, Co, Ag, and Zr, a clear trend towards
segregation within the grain boundaries rather than the bulk was found. The
aforementioned elements also showed negative strengthening energies,
indicating that they are expected to locally increase the resistance of the
crystal against decohesion. Furthermore, the results highlighted a preference
for non-metallic impurities, like O, P, F, S, and Se, to segregate at the
grain boundaries as opposed to the bulk of Cu. This enhanced affinity of
impurities to grain boundaries, along with the predicted local weakening
effect of these elements, could potentially explain experimental findings that
point to the initiation of degradation phenomena in these regions. A link was
determined between the segregation energies and the induced relaxation
effects. Metallic solutes with segregation energies close to zero were shown
to induce minimal distortion, whereas non-metallic elements that strongly
segregate at the GB caused severe displacement of nearby Cu atoms.

Grain boundaries are known to play an important role in the properties of
metals. On the basis of our findings, the incorporation of metallic dopants
such as V, Zr, and Ag at the grain boundaries of Cu has the potential to
improve the mechanical properties of these regions and enhance their
resistance against degradation. Further studies are needed to understand the
behavior of these impurities when introduced into polycrystalline Cu.

## Acknowledgements

A.L.S. acknowledges funding by EPSRC (grant EP/P013503/1). V.F. would like to
acknowledge funding by EPSRC (grant EP/L015862/1) as part of the CDT in
molecular modeling and materials science. Computational resources on ARCHER2
(http://www.archer2.ac.uk) were provided via our membership in the UK’s HPC
Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202,
EP/R029431). V.F. and A.L.S. would like to thank Rishi Bodlos, Lorenz Romaner,
and Ernst Kozeschnik for useful comments and help with calculations.

## References

  * Moser et al. 2021 Moser, S.; Kleinbichler, M.; Kubicek, S.; Zechner, J.; Cordill, M. J. Electropolishing—a practical method for accessing voids in metal films for analyses. _Applied Sciences_ 2021, _11_ , 7009
  * Konishi et al. 2002 Konishi, S.; Moriyama, M.; Murakami, M. Effect of annealing atmosphere on void formation in copper interconnects. _Materials Transactions_ 2002, _43_ , 1624–1628
  * Hallberg et al. 2018 Hallberg, H.; Ås, S. K.; Skallerud, B. Crystal plasticity modeling of microstructure influence on fatigue crack initiation in extruded Al6082-T6 with surface irregularities. _International Journal of Fatigue_ 2018, _111_ , 16–32
  * Stopka and McDowell 2020 Stopka, K. S.; McDowell, D. L. Microstructure-sensitive computational multiaxial fatigue of Al 7075-T6 and duplex Ti-6Al-4V. _International Journal of Fatigue_ 2020, _133_ , 105460
  * Hu et al. 2018 Hu, T.; Yang, S.; Zhou, N.; Zhang, Y.; Luo, J. Role of disordered bipolar complexions on the sulfur embrittlement of nickel general grain boundaries. _Nature communications_ 2018, _9_ , 2764
  * Duscher et al. 2004 Duscher, G.; Chisholm, M. F.; Alber, U.; Rühle, M. Bismuth-induced embrittlement of copper grain boundaries. _Nature materials_ 2004, _3_ , 621–626
  * Ludwig et al. 2005 Ludwig, W.; Pereiro-López, E.; Bellet, D. In situ investigation of liquid Ga penetration in Al bicrystal grain boundaries: grain boundary wetting or liquid metal embrittlement? _Acta Materialia_ 2005, _53_ , 151–162
  * Laporte and Mortensen 2009 Laporte, V.; Mortensen, A. Intermediate temperature embrittlement of copper alloys. _International Materials Reviews_ 2009, _54_ , 94–116
  * Kang et al. 2013 Kang, J.; Glatzmaier, G. C.; Wei, S.-H. Origin of the bismuth-induced decohesion of nickel and copper grain boundaries. _Physical review letters_ 2013, _111_ , 055502
  * Schweinfest et al. 2004 Schweinfest, R.; Paxton, A. T.; Finnis, M. W. Bismuth embrittlement of copper is an atomic size effect. _Nature_ 2004, _432_ , 1008–1011
  * Lozovoi et al. 2006 Lozovoi, A.; Paxton, A.; Finnis, M. Structural and chemical embrittlement of grain boundaries by impurities: A general theory and first-principles calculations for copper. _Physical Review B_ 2006, _74_ , 155416
  * Scheiber 2021 Scheiber, D. Segregation and embrittlement of gold grain boundaries. _Computational Materials Science_ 2021, _187_ , 110110
  * Zhou et al. 2023 Zhou, X.; Ahmadian, A.; Gault, B.; Ophus, C.; Liebscher, C. H.; Dehm, G.; Raabe, D. Atomic motifs govern the decoration of grain boundaries by interstitial solutes. _Nature Communications_ 2023, _14_ , 3535
  * Xu et al. 2019 Xu, Z.; Tanaka, S.; Kohyama, M. Grain-boundary segregation of 3d-transition metal solutes in bcc Fe: ab initio local-energy and d-electron behavior analysis. _Journal of Physics: Condensed Matter_ 2019, _31_ , 115001
  * Mai et al. 2022 Mai, H. L.; Cui, X.-Y.; Scheiber, D.; Romaner, L.; Ringer, S. P. The segregation of transition metals to iron grain boundaries and their effects on cohesion. _Acta materialia_ 2022, _231_ , 117902
  * Scheiber et al. 2018 Scheiber, D.; Romaner, L.; Pippan, R.; Puschnig, P. Impact of solute-solute interactions on grain boundary segregation and cohesion in molybdenum. _Physical Review Materials_ 2018, _2_ , 093609
  * Razumovskiy et al. 2015 Razumovskiy, V. I.; Lozovoi, A.; Razumovskii, I. First-principles-aided design of a new Ni-base superalloy: Influence of transition metal alloying elements on grain boundary and bulk cohesion. _Acta Materialia_ 2015, _82_ , 369–377
  * Moser et al. 2019 Moser, S.; Zernatto, G.; Kleinbichler, M.; Nelhiebel, M.; Zechner, J.; Cordill, M. J.; Pippan, R. A novel setup for in situ monitoring of thermomechanically cycled thin film metallizations. _JOM_ 2019, _71_ , 3399–3406
  * Huang et al. 2018 Huang, Z.; Chen, F.; Shen, Q.; Zhang, L.; Rupert, T. J. Uncovering the influence of common nonmetallic impurities on the stability and strength of a $\Sigma$5 (310) grain boundary in Cu. _Acta Materialia_ 2018, _148_ , 110–122
  * Kresse and Hafner 1993 Kresse, G.; Hafner, J. Ab initio molecular dynamics for open-shell transition metals. _Physical Review B_ 1993, _48_ , 13115
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Physical review B_ 1996, _54_ , 11169
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. _Computational Materials Science_ 1996, _6_ , 15–50
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. _Physical review letters_ 1996, _77_ , 3865
  * Wu et al. 2016 Wu, X.; You, Y.-W.; Kong, X.-S.; Chen, J.-L.; Luo, G.-N.; Lu, G.-H.; Liu, C.; Wang, Z. First-principles determination of grain boundary strengthening in tungsten: Dependence on grain boundary structure and metallic radius of solute. _Acta Materialia_ 2016, _120_ , 315–326
  * Nazarov et al. 2012 Nazarov, R.; Hickel, T.; Neugebauer, J. Vacancy formation energies in fcc metals: influence of exchange-correlation functionals and correction schemes. _Physical Review B_ 2012, _85_ , 144118
  * Nazarov et al. 2014 Nazarov, R.; Hickel, T.; Neugebauer, J. Ab initio study of H-vacancy interactions in fcc metals: Implications for the formation of superabundant vacancies. _Physical Review B_ 2014, _89_ , 144108
  * Fotopoulos et al. 2023 Fotopoulos, V.; Grau-Crespo, R.; Shluger, A. L. Thermodynamic analysis of the interaction between metal vacancies and hydrogen in bulk Cu. _Physical Chemistry Chemical Physics_ 2023, _25_ , 9168–9175
  * Wurmshuber et al. 2022 Wurmshuber, M.; Burtscher, M.; Doppermann, S.; Bodlos, R.; Scheiber, D.; Romaner, L.; Kiener, D. Mechanical performance of doped W–Cu nanocomposites. _Materials Science and Engineering: A_ 2022, _857_ , 144102
  * Bodlos et al. 2023 Bodlos, R.; Scheiber, D.; Spitaler, J.; Romaner, L. Modification of the Cu/W interface cohesion by segregation. _Metals_ 2023, _13_ , 346
  * Li et al. 2017 Li, Y.; Korzhavyi, P. A.; Sandström, R.; Lilja, C. Impurity effects on the grain boundary cohesion in copper. _Physical Review Materials_ 2017, _1_ , 070602
  * Lousada and Korzhavyi 2020 Lousada, C. M.; Korzhavyi, P. A. Segregation of p and s impurities to a $\Sigma$9 grain boundary in Cu. _Metals_ 2020, _10_ , 1362
  * Wang and Benabou 2018 Wang, X.; Benabou, L. First-principles investigation of intergranular fracture in copper by grain boundary segregation of sulfur. _Journal of Engineering Materials and Technology_ 2018, _140_ , 011008
  * Razumovskiy et al. 2018 Razumovskiy, V.; Divinski, S.; Romaner, L. Solute segregation in Cu: DFT vs. Experiment. _Acta Materialia_ 2018, _147_ , 122–132
  * Huang et al. 2020 Huang, Z.; Wang, P.; Chen, F.; Shen, Q.; Zhang, L. Understanding solute effect on grain boundary strength based on atomic size and electronic interaction. _Scientific Reports_ 2020, _10_ , 16856
  * Huang et al. 2019 Huang, Z.; Chen, F.; Shen, Q.; Zhang, L.; Rupert, T. J. Combined effects of nonmetallic impurities and planned metallic dopants on grain boundary energy and strength. _Acta Materialia_ 2019, _166_ , 113–125


# Cluster expansion for the description of condensed state: crystalline cell
approach

G.S. Bokun${}^{\textsf{\footnotesize{\mbox{?}}}}$, M.F.
Holovko${}^{\textsf{\footnotesize{\mbox{?}}}}$

(Received May 2, 2018, in final form July 10, 2018)

###### Abstract

A well-known cluster expansion, which leads to virial expansion for the free
energy of low density systems, is modified in such a way that it becomes
applicable to the description of condensed state of matter. To this end, the
averaging of individual clusters over the states of an ideal gas is replaced
by the averaging over the states of a non-correlated crystal using single-
particle cell potentials. As a result, we arrive at the expansion of the
partition function in correlations on the basis of single-particle functions
corresponding to the multiplicative approximation. The cell potentials
defining these functions are found from the condition of the minimum of the
remainder in the constructed decomposition.

Key words: lattice models, cluster expansions, single-particle cell potential,
free energy

PACS: 5.20.-y, 61.72jd, 64.30.-t, 65.40.-b

###### Abstract

Øèðîêî âäîìå ãðóïîâå ðîçâèíåííÿ, ÿêå ïðèâîäèòü äî âðàëüíîãî ðîçâèíåííÿ äëÿ
âëüíî¿ åíåðã¿ ðîçðäæåíèõ ñèñòåì, ìîäèôêîâàíî òàê, ùîá éîãî ìîæíà áóëî
çàñòîñîâóâàòè äî êîíäåíñîâàíîãî ñòàíó ðåîâèíè. Äëÿ öüîãî óñåðåäíåííÿ îêðåìèõ
êëàñòåðâ ïî ñòàíàõ äåàëüíîãî ãàçó çàìíþòüñÿ óñåðåäíåííÿì ïî ñòàíàõ
íåêîðåëüîâàíîãî êðèñòàëà, âèêîðèñòîâóþè êîìðêîâ îäíîàñòèíêîâ ïîòåíöàëè. Â
ðåçóëüòàò îòðèìàíî ðîçâèíåííÿ ñòàòèñòèíî¿ ñóìè ïî êîðåëÿöÿõ íà áàçèñ
îäíîàñòèíêîâèõ ôóíêöé, ÿê âäïîâäàþòü ìóëüòèïëêàòèâíîìó íàáëèæåííþ. ðàòêîâ
ïîòåíöàëè, ùî âèçíààþòü âêàçàí ôóíêö¿, çíàõîäÿòüñÿ ç óìîâè ìíìçàö¿ çàëèøêó â
ñêîíñòðóéîâàíîìó ðîçâèíåíí.

Ключов слова: ðàòêîâ ìîäåë, ãðóïîâ ðîçâèíåííÿ, îäíîàñòèíêîâèé ðàòêîâèé
ïîòåíöàë, âëüíà åíåðãÿ

##  1 Introduction

The cellular theory, which is based on the structuring of the states of the
system, has made it possible to solve a number of important problems in the
physics of condensed state [1, 2]. The lattice version of this theory has been
used to address fundamental problems of statistical physics as well as to
calculate specific properties of various systems [3, 4]. The development and
introduction of nanomaterials has driven these approaches to be applied to the
description of highly heterogeneous structures, which, in addition to large
gradients of order parameter fields, are characterized by a heterogeneous
state accompanied by phase transitions of various types [5, 6]. Here, as in
the case of homogeneous systems, it seems effective to involve mean-field
representations based on the use of cell potentials [7, 8].

The present paper is devoted to the development of the basis of this approach.
To this end, a system of particles in the mean field of single-particle cell
potentials is used as the reference system. The cell potentials determine the
average forces acting on a particle fixed in a selected cell from the
particles distributed in other molecular cells. Initially, only those states
are taken into consideration where each cell of the system is occupied by one
molecule or is vacant. The heterogeneity in the system is taken into account
by the fact that the cell volumes are not the same. In order to calculate the
partition function of the initial condensed system, the perturbation theory is
used based on the expansions over the states of the reference system using
generalized Mayer functions. These functions contain not only intermolecular
interaction potentials but also the necessary mean potentials. The latter are
determined from the condition for optimizing the deviation of the properties
of the reference and the real systems. Furthermore, the thermodynamic
consistency of the theory is studied at the level of calculation of the first
derivatives with respect to thermodynamic parameters. The developed
microscopic approach is generalized to the case of long-range interaction. A
possibility to take account of a more complete set of occupation number values
is discussed.

The paper is arranged as follows. In section 2, the lattice model and the
description of a corresponding reference system are presented. In order to
take into account the inter-particle correlations in section 3, we formulate
the perturbation theory based on the expansions over the states of the
reference system. In section 4, we formulate the optimal choice of a single-
particle cell potentials needed for the description of the reference system.
In section 5, the verification of thermodynamic self-consistency of the
presented approach is considered at the level of calculation of the first
derivatives with respect to thermodynamic parameters. In section 6, the
generalization of the considered approach for the case of a more complete set
of occupation numbers values and the generalization to the case of the systems
with long-range interaction are discussed. We conclude in section 7.

##  2 The model and the reference system

In this paper we consider a lattice model with the Hamiltonian

$H_{M}=\frac{1}{2}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\Phi\big{(}q_{n_{i}},q_{n_{j}}\big{)}+\sum\limits_{i=1}^{M}\mu_{i}n_{i}\,,$
(2.1)

where $M$ is the total number of lattice sites (i.e., cells in the system),
$\sum_{j\left(i\right)}^{Z}$ denotes consecutive summation over all nodes
surrounding the selected node $i$ taking into account $Z$ counted neighbors,
$\mu_{i}$ is the chemical potential value in the node $i$. The classification
of states considered is similar to that adopted in [9]. The variable
$q_{n_{i}}$ determines the position of the particle ($n_{i}=1_{i}$) or the
vacancy ($n_{i}=0_{i}$) in the cell $i$,
$\Phi\big{(}q_{n_{i}},q_{n_{j}}\big{)}$ is the interaction potential between
two particles with coordinates $q_{n_{i}}$ and $q_{n_{j}}$. We suppose that
this potential is short-ranged enough and can be taken, for instance, in the
Lennard-Jones form [10].

In order to describe the considered model with the Hamiltonian (2.1), we apply
the concept of the reference system widely used in the statistical theory of
various condensed systems. Usually the reference system is a simplified
version of the real condensed system. It should include the main features of
the real model and should be described analytically with sufficient accuracy.
For example, the simplest reference system is the model of the ideal gas. The
application of this model leads to the virial density expansion for
thermodynamic properties of real gases [11]. Another important reference
system is the model of hard spheres, which has been successfully employed in
the modern liquid state theory [10]. In this paper, to study the system with
the Hamiltonian (2.1), we use the reference system that can be described by
the Hamiltonian

$H_{0}=\sum\limits_{i=1}^{M}\mu_{i}n_{i}+\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\phi_{j}\left(q_{n_{i}}\right)$
(2.2)

represented by single-particle cell potentials
$\phi_{j}\left(q_{n_{i}}\right)$, which has the meaning of the potential of an
external field whose source is located formally at the center of the cell $j$.
$\phi_{j}\left(q_{n_{i}}\right)$ depends on the variables that determine the
position of particles or vacancies in the cell $i$ and parametrically is a
function of the quantities characterizing the average distribution of
particles or vacancies in the system and its macroscopic state. In this case,
$n=1$ corresponds to the distribution of particles and $n=0$ corresponds to
that of vacancies.

The connection between the Hamiltonian $H_{M}$ and $H_{0}$ and the calculation
of the single-particle potentials is discussed in the next sections. In this
section we consider only the description of the reference system without
elaboration of the single-particle cell potential
$\phi_{j}\left(q_{n_{i}}\right)$.

We consider the case of an inhomogeneous system with inhomogeneity
characterized by the field of mean occupation numbers

$\rho_{0}=\langle n_{i}\rangle_{0}\,.$ (2.3)

In order to simplify further expressions, we use abbreviations wherever
possible, for example $U(q_{n_{i}})=U_{n_{i}}$ denote the potential acting on
a particle or a vacancy in the position $q_{n_{i}}$. In order to take into
account the variable number of particles and vacancies in the system, we
consider both $q_{1_{i}}$ and $q_{0_{i}}$ as two states of a certain “virtual”
particle. The first one corresponds to the position of a real particle and the
second one to that of a vacancy. In this case, it is possible to replace fixed
particles by arbitrary ones.

This can be seen from the definition of $Z_{M}$

$Z_{M}=\sum\limits_{N=0}^{M}\frac{1}{N!}Q_{N}\text{e}^{-\mu N}.$ (2.4)

Due to indiscernibility of particles, $N!$ is reduced, which makes it possible
to represent $Z_{M}^{(0)}$ of the reference system in the form

$Z_{M}^{(0)}=\sum\limits_{n_{1}=0}^{1}\int\limits_{\omega_{1}}\mathrm{d}q_{n_{1}}...\sum\limits_{n_{i}=0}^{1}\int\limits_{\omega_{i}}\mathrm{d}q_{n_{i}}...\sum\limits_{n_{M}=0}^{1}\int\limits_{\omega_{M}}\mathrm{d}q_{n_{M}}\exp\Bigg{\\{}-\beta\left[\sum\limits_{l=1}^{M}\left(\mu_{l}n_{l}+U_{n_{l}}\right)\right]\Bigg{\\}},$
(2.5)

where $\beta=1/(kT)$, $k$ is the Boltzmann constant, $T$ is the temperature,
$\omega_{i}$ is the volume of the lattice cell $i$,

$U_{n_{l}}=\sum\limits_{j\left(l\right)}^{Z}\phi_{j}\left(q_{n_{l}}\right).$
(2.6)

In accordance with the representation (2.5), the particle distribution
function over the volume of the system turns out to be factorized and can be
written in the form

$D_{M}^{\left(0\right)}=\prod\limits_{i=1}^{M}\rho\left(q_{n_{i}}\right),$
(2.7)

where

$\rho\left(q_{n_{i}}\right)=\exp\left[-\beta\left(\mu_{i}n_{i}+U_{n_{i}}\right)\right].$
(2.8)

The normalization of this function, respectively, is represented by the
expression

$Z_{M}^{\left(0\right)}=\prod\limits_{i=1}^{M}z_{i}^{0}\,,$ (2.9)

where

$z_{i}^{0}=\sum\limits_{n_{i}=0}^{1}\int\limits_{\omega_{i}}\exp\left[-\beta\left(\mu_{i}n_{i}+U_{n_{i}}\right)\right]\mathrm{d}q_{n_{i}}\,,$
(2.10) $z_{i}^{0}=Q_{0_{i}}+\text{e}^{\beta\mu_{i}}Q_{1_{i}}\,,$ (2.11)
$Q_{n_{i}}=\int\limits_{\omega_{i}}\exp\Bigg{[}-\beta\sum\limits_{j\left(i\right)}^{Z}\phi_{j}\left(q_{n_{i}}\right)\Bigg{]}\mathrm{d}q_{n_{i}}\,.$
(2.12)

Based on equations (2.7)–(2.12), we write the expression for the normalized
distribution function, which is necessary for averaging the magnitudes of the
singlet, binary, and other types. Thus, denoting the normalized functions by a
cap above the notations of respective functions, we write

$\widehat{D}_{M}=\prod\limits_{i=1}^{M}\widehat{\rho}\left(q_{n_{i}}\right),$
(2.13)

where

$\widehat{\rho}\left(q_{n_{i}}\right)=\frac{\exp\left[-\beta\left(\mu_{i}n_{i}+U_{n_{i}}\right)\right]}{z_{i}^{0}}\,.$
(2.14)

We transform the above expressions based on the relation between $\mu_{i}$ and
the mean values of the occupation numbers $\rho_{1_{i}}$ and $\rho_{0_{i}}$.
Integrating (2.14) with respect to $q_{n_{i}}$, we find

$\rho_{n_{i}}=\frac{\text{e}^{\beta\mu_{i}n_{i}}Q_{n_{i}}}{z_{i}^{0}}.$ (2.15)

From (2.15) it follows that

$\frac{\rho_{1_{i}}}{\rho_{0_{i}}}=\frac{\text{e}^{\beta\mu_{i}}Q_{1_{i}}}{Q_{0_{i}}}.$
(2.16)

The relation (2.16) allows us to write

$\text{e}^{\beta\mu_{i}}=\frac{\rho_{1_{i}}Q_{0_{i}}}{\rho_{0_{i}}Q_{1_{i}}}.$
(2.17)

Substituting (2.17) into (2.11), we obtain

$z_{i}^{0}=Q_{0_{i}}+\frac{\rho_{1_{i}}Q_{0_{i}}}{\rho_{0_{i}}}=\frac{Q_{0_{i}}}{\rho_{0_{i}}}.$
(2.18)

The substitution of (2.17) and (2.18) allows us to write

$\widehat{\rho}\left(q_{n_{i}}\right)=\rho_{n_{i}}\frac{\exp\left[-\beta\sum\limits_{j\left(i\right)}^{Z}\phi_{j}\left(q_{n_{i}}\right)\right]}{Q_{n_{i}}}.$
(2.19)

Let us consider the averaging that employs (2.13), for example the
characteristics of a binary type $L\big{(}q_{n_{i}},q_{n_{j}}\big{)}$
according to the definition

${\langle
L_{ij}\rangle_{0}=\sum\limits_{n_{1}=0}^{1}\sum\limits_{n_{2}=0}^{1}...\sum\limits_{n_{m}=0}^{1}\int\limits_{\omega_{1}}\widehat{\rho}\left(q_{n_{1}}\right)...}{\int\limits_{\omega_{j}}\widehat{\rho}\big{(}q_{n_{j}}\big{)}...\int\limits_{\omega_{M}}\widehat{\rho}\big{(}q_{n_{j}}\big{)}...\int\limits_{\omega_{m}}\widehat{\rho}\left(q_{nm}\right)L\big{(}q_{n_{i}}q_{n_{j}}\big{)}}{}.$
(2.20)

Transposing the summation in (2.20) and taking into account the independence
of integration variables and the normalized condition

$\sum\limits_{n_{i}=0}^{1}\int\limits_{\omega_{i}}\widehat{\rho}\left(q_{n_{i}}\right)\mathrm{d}q_{n_{i}}=\sum\limits_{n_{i}=0}^{1}\rho_{n_{i}}=1$
(2.21)

we obtain

$\langle
L_{ij}\rangle_{0}=\sum\limits_{n_{i=0}}^{1}\sum\limits_{n_{j}=0}^{1}\int\limits_{\omega_{i}}\mathrm{d}q_{n_{i}}\int\limits_{\omega_{j}}\mathrm{d}q_{n_{j}}L\big{(}q_{n_{i}},q_{n_{j}}\big{)}\widehat{\rho}\left(q_{n_{i}}\right)\widehat{\rho}\big{(}q_{n_{j}}\big{)}\,,$
(2.22)

or

$\langle
L_{ij}\rangle_{0}=\sum\limits_{n_{i=0}}^{1}\sum\limits_{n_{j}=0}^{1}L_{n_{i,}n_{j}}\rho_{n_{i}}\rho_{n_{j}}\,,$
(2.23)

$L_{n_{i,}n_{j}}=\int\limits_{\omega_{i}}\mathrm{d}q_{n_{i}}\int\limits_{\omega_{j}}\mathrm{d}q_{n_{j}}L\big{(}q_{n_{i}},q_{n_{j}}\big{)}\widehat{F}_{11}\left(q_{n_{i}}\right)\widehat{F}_{11}\big{(}q_{n_{j}}\big{)}.$
(2.24)

The functions $\widehat{F}_{11}\left(q\right)$ in (2.24) correspond to singlet
distribution functions $F_{11}$ in the approximation of the method of
conditional distributions [9]

$\widehat{F}_{11}\left(q_{n_{i}}\right)=\frac{1}{Q_{n_{i}}}\exp\Bigg{[}-\beta\sum\limits_{j\left(i\right)}^{Z}\phi_{j}\left(q_{n_{i}}\right)\Bigg{]}.$
(2.25)

##  3 Perturbation theory

In this section, we consider a perturbation scheme for the treatment of the
remainder

$\Delta
H_{M}=H_{M}-H_{0}=\frac{1}{2}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\Delta\phi\big{(}q_{n_{i}},q_{n_{j}}\big{)}\,,$
(3.1)

where

$\Delta\phi\big{(}q_{n_{i}},q_{n_{j}}\big{)}=\Phi\big{(}q_{n_{i}},q_{n_{j}}\big{)}-\phi_{j}\big{(}q_{n_{i}}\big{)}-\phi_{i}\big{(}q_{n_{j}}\big{)}.$
(3.2)

Now, we represent the partition function of the original system in the form

$Z_{M}=Z_{M}^{\left(0\right)}\langle\text{e}^{-\beta\Delta
H_{M}}\rangle_{0}\,,$ (3.3)

where $\langle\dots\rangle_{0}$ is the averaging represented by the expression

$\langle
L\rangle_{0}=\sum\limits_{n_{1}=0}^{1}...\sum\limits_{n_{M}=0}^{1}\int\limits_{\omega_{1}}\mathrm{d}q_{n_{U}}\widehat{\rho}\left(q_{n_{U}}\right)...\int\limits_{\omega_{M}}\widehat{\rho}\left(q_{n_{M}}\right)\mathrm{d}q_{n_{M}}L\,.$
(3.4)

To calculate (3.3), a cumulant expansion [12, 10] is used, leading to an
expansion in powers of the density if the averaging in (3.3) is performed over
the states corresponding to an ideal gas [11]. In our case, the distribution
characteristic of an ideal crystal is used as the reference system. This
allows us to obtain a suitable description of the properties of a condensed
system. The virial coefficients in this case become density functions. In
other words, the transformation (3.3) is an expansion over the cluster
correlations, although formally it has the form of an expansion over the Mayer
functions. For the latter we use the renormalized Mayer functions of the form

$f\big{(}q_{n_{i}},q_{n_{j}}\big{)}=\exp\Big{[}-\beta\Delta\phi\big{(}q_{n_{i}},q_{n_{j}}\big{)}\Big{]}-1.$
(3.5)

We should note that for vacancies
$\Phi(q_{0_{i}},q_{0_{j}})=\Phi(q_{0_{i}},q_{1_{j}})=\Phi(q_{1_{i}},q_{0_{j}})=0$.
However, the Mayer functions $f(q_{n_{i}},q_{n_{j}})\neq 0$ because in this
case $\bigtriangleup\phi(q_{n_{i}},q_{n_{j}})\neq 0$.

Using the procedure of group expansion in (3.3), we obtain

$\displaystyle Z_{M}$
$\displaystyle=Z_{M}^{\left(0\right)}\exp\Bigg{[}\frac{1}{2}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\sum\limits_{n_{i}=0}^{1}\sum\limits_{n_{j}=0}^{1}\int\limits_{\omega_{i}}\mathrm{d}q_{n_{i}}\int\limits_{\omega_{j}}\mathrm{d}q_{n_{j}}\widehat{\rho}\left(q_{n_{i}}\right)\widehat{\rho}\big{(}q_{n_{j}}\big{)}f\big{(}q_{n_{i}},q_{n_{j}}\big{)}$
$\displaystyle+\frac{1}{6}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\sum\limits_{l\left(i\right)}^{Z}\sum\limits_{n_{i}=0}^{1}\sum\limits_{n_{j}=0}^{1}\sum\limits_{n_{l}=0}^{1}\int\limits_{\omega_{i}}\mathrm{d}q_{n_{i}}\int\limits_{\omega_{j}}\mathrm{d}q_{n_{j}}\int\limits_{\omega_{l}}\mathrm{d}q_{n_{l}}\widehat{\rho}\left(q_{n_{i}}\right)\widehat{\rho}\big{(}q_{n_{j}}\big{)}\widehat{\rho}\left(q_{n_{l}}\right)$
$\displaystyle\times
f\big{(}q_{n_{i}},q_{n_{j}}\big{)}f\big{(}q_{n_{j}},q_{n_{l}}\big{)}f\big{(}q_{n_{i}},q_{n_{l}}\big{)}+...\Bigg{]}.$
(3.6)

Or, due to representations (2.23) and (2.24),

$\displaystyle Z_{M}$
$\displaystyle=Z_{M}^{\left(0\right)}\exp\Bigg{[}\frac{1}{2}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\sum\limits_{n_{i}=0}^{1}\sum\limits_{n_{j}=0}^{1}\rho_{n_{i}}\rho_{n_{j}}f_{n_{i}n_{j}}$
$\displaystyle+\frac{1}{6}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\sum\limits_{l\left(i\right)}^{Z}\sum\limits_{n_{i}=0}^{1}\sum\limits_{n_{j}=0}^{1}\sum\limits_{n_{l}=0}^{1}\rho_{n_{i}}\rho_{n_{j}}\rho_{n_{l}}f_{n_{i}n_{j}n_{l}}+...\Bigg{]}\,,$
(3.7)

where

$\displaystyle f_{n_{i}n_{j}}$
$\displaystyle=\int\limits_{\omega_{i}}\mathrm{d}q_{n_{i}}\int\limits_{\omega_{j}}\mathrm{d}q_{n_{j}}f\big{(}q_{n_{i}},q_{n_{j}}\big{)}\widehat{F}_{11}\big{(}q_{n_{i}}\big{)}\widehat{F}_{11}\big{(}q_{n_{j}}\big{)}\,,$
(3.8) $\displaystyle f_{n_{i}n_{j}n_{l}}$
$\displaystyle=\int\limits_{\omega_{i}}\mathrm{d}q_{n_{i}}\int\limits_{\omega_{j}}\mathrm{d}q_{n_{j}}\int\limits_{\omega_{l}}\mathrm{d}q_{n_{l}}f\big{(}q_{n_{i}},q_{n_{j}}\big{)}f\big{(}q_{n_{j}},q_{n_{l}}\big{)}f\left(q_{n_{i}},q_{n_{l}}\right)$
$\displaystyle\times\widehat{F}_{11}\left(q_{n_{i}}\right)\widehat{F}_{11}\big{(}q_{n_{j}}\big{)}\widehat{F}_{11}\left(q_{n_{l}}\right).$
(3.9)

From (3), we have the cluster expansion for the free energy

$\displaystyle F$ $\displaystyle=-kT\ln Z_{M}=-kT\Bigg{[}\ln
Z_{M}^{0}+\frac{1}{2}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\sum\limits_{n_{i}=0}^{1}\sum\limits_{n_{j}=0}^{1}\rho_{n_{i}}\rho_{n_{j}}f_{n_{i}n_{j}}$
$\displaystyle+\frac{1}{6}\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\sum\limits_{l\left(i\right)}^{Z}\sum\limits_{n_{i}=0}^{1}\sum\limits_{n_{j}=0}^{1}\sum\limits_{n_{l}=0}^{1}\rho_{n_{i}}\rho_{n_{j}}\rho_{n_{l}}f_{n_{i}n_{j}n_{l}}+...\Bigg{]}.$
(3.10)

##  4 An optimal choice of a single-particle potential

In order to apply the considered theory, we should specify the single particle
potential. In this section we propose an optimal choice for this potential.
This approach is in some sense similar to the problem of the connection
between models with soft and hard core repulsions [13] which are successfully
used in modern liquid state theory [14]. The equation defining the single-
particle potentials of the reference system is determined from the self-
consistent condition, which has several different formulations. One of them is
connected with the extremum of the remainder in the expansion (3) due to the
fact that the sum of all terms contained in (3) does not depend on the choice
of single-particle potentials. Since a change of the potentials leads to
redistribution of contributions of individual terms, the best choice would be
the one with the maximum contribution of the terms that are taken into
account. This is analogous to the requirement of the minimum susceptibility of
the system to a virtual external field and leads to the condition

$\frac{\delta\ln Z_{M}}{\delta\phi_{k}(q_{n_{m}})}=0.$ (4.1)

Namely, the condition of the extremum of the part of the functional written in
(3) corresponds simultaneously to the condition that the sum in brackets in
(3) tends to zero. Let us consider the proof of the foregoing. We vary (3)
over all the potentials in (2.6) assuming that they can all be independent of
each other, which makes it possible to substantially simplify the procedure of
transformations. Performing the variation with respect to an individual
$\phi_{k}\left(q_{n_{m}}\right)$ and using (2.8) and (2.18), for the reference
system part we obtain

$\displaystyle\delta\ln z_{m}^{0}$
$\displaystyle=\frac{1}{z_{m}^{0}}\Bigg{\\{}\int\limits_{\omega_{m}}\delta\phi_{k}\big{(}q_{0_{m}}\big{)}\exp\bigg{[}-\beta\sum_{j(m)}^{Z}\phi_{j}\big{(}q_{0_{m}}\big{)}\bigg{]}\mathrm{d}q_{0_{m}}$
$\displaystyle+\int\limits_{\omega_{m}}\delta\phi\big{(}q_{1_{m}}\big{)}\exp\bigg{[}-\beta\sum_{j(m)}^{Z}\phi_{j}\big{(}q_{1_{m}}\big{)}\bigg{]}{\rm
e}^{\beta\mu_{m}}\mathrm{d}q_{1_{m}}\Bigg{\\}}.$ (4.2)

Since the variation in $\delta\ln Z_{M}$ is satisfied for fixed $\mu_{i}$, we
transform the second term in $Z_{M}$ into a form that contains explicit
$\mu_{i}$.

$\displaystyle\mathcal{I}_{m}=\frac{1}{2}\sum\limits_{j\left(m\right)}^{Z}\sum\limits_{n_{m}=0}^{1}\sum\limits_{n_{j}=0}^{1}\rho_{n_{m}}\rho_{n_{j}}f_{n_{m}n_{j}}=\frac{1}{2}\sum\limits_{j\left(m\right)}^{Z}\sum\limits_{n_{m}=0}^{1}\sum\limits_{n_{j}=0}^{1}\int\limits_{\omega_{m}}\int\limits_{\omega_{j}}\exp\bigg{(}-\beta\Big{\\{}\Phi\big{(}q_{n_{m}},q_{n_{j}}\big{)}$
$\displaystyle+\sum\limits_{s\neq
m,j}\left[\phi_{s}\big{(}q_{n_{m}}\big{)}+\phi_{s}\big{(}q_{n_{j}}\big{)}\right]\Big{\\}}\bigg{)}{\text{e}^{\beta\mu_{m}n_{m}}\text{e}^{\beta\mu_{j}n_{j}}\mathrm{d}q_{n_{m}}\mathrm{d}q_{n_{j}}\frac{1}{z_{m}^{0}z_{j}^{0}}}-\frac{1}{2}\,.$
(4.3)

Let us write the expression for the variation of (4) with respect to
$\phi_{k}\left(q_{n_{m}}\right)$

$\displaystyle\delta\mathcal{I}_{m,k}$ $\displaystyle=\sum\limits_{j\neq
k,m}\sum\limits_{n_{m}=0}^{1}\sum\limits_{n_{j}=0}^{1}\int\limits_{\omega_{m}}\int\limits_{\omega_{j}}\delta\phi_{k}\left(q_{n_{m}}\right)\exp\left[-\beta\Phi\big{(}q_{n_{m}},q_{n_{j}}\big{)}\right]$
$\displaystyle\times\exp\Big{\\{}-\beta\sum\limits_{s\neq
m,j}\left[\phi_{s}\big{(}q_{{}_{m}}\big{)}+\phi_{s}\big{(}q_{n_{j}}\big{)}\right]\Big{\\}}\text{e}^{\beta\mu_{m}n_{m}}\text{e}^{\beta\mu_{j}n_{j}}\mathrm{d}q_{n_{m}}\mathrm{d}q_{n_{j}}\displaystyle\frac{1}{z_{m}^{0}z_{j}^{0}}$
$\displaystyle-\sum\limits_{j\left(m\right)}^{Z}\sum\limits_{n_{m}=0}^{1}\sum\limits_{n_{j}=0}^{1}\int\limits_{\omega_{m}}\int\limits_{\omega_{j}}\exp\Big{\\{}-\beta\Phi\big{(}q_{n_{m}},q_{n_{j}}\big{)}+\sum\limits_{s\neq
m,j}\left[\phi_{s}\left(q_{n_{m}}\right)+\phi_{s}\big{(}q_{n_{j}}\big{)}\right]\Big{\\}}$
$\displaystyle\times\text{e}^{\beta\mu_{m}n_{m}}\text{e}^{\beta\mu_{j}n_{j}}\mathrm{d}q_{n_{m}}\mathrm{d}q_{n_{j}}\displaystyle\frac{1}{z_{m}^{0}z_{{}_{j}}^{0}}\delta\ln
z_{m}^{0}\,.$ (4.4)

The first sum on the right-hand side of equation (4) does not contain a term
with $j=m$, which is convenient to add and subtract, which allows the sum (4)
and (4) to be represented in a form that allows the separation of variables by
cell numbers. Namely,

$\displaystyle\delta\ln z_{m}^{0}+\delta\mathcal{I}_{m,k}$
$\displaystyle=\frac{1}{z^{0}_{m}}\sum\limits_{n_{m}=0}^{1}\int\limits_{\omega_{m}}\delta\phi_{k}\left(q_{n_{m}}\right)\exp\bigg{[}-\beta\sum\limits_{j\left(m\right)}^{Z}\phi_{j}\left(q_{n_{m}}\right)\bigg{]}\text{e}^{\beta\mu_{m}n_{m}}\mathrm{d}q_{n_{m}}$
$\displaystyle-\frac{1}{z^{0}_{m}z^{0}_{k}}\sum\limits_{n_{m}=0}^{1}\sum\limits_{n_{k}=0}^{1}\int\limits_{\omega_{k}}\delta\phi_{k}\left(q_{n_{m}}\right)\mathrm{d}q_{n_{m}}\int\limits_{\omega_{k}}\mathrm{d}q_{n_{k}}\exp\left[-\beta\Phi\left(q_{n_{m}},q_{n_{k}}\right)\right]$
$\displaystyle\times\exp\Big{\\{}-\beta\sum\limits_{s\neq
m,k}\left[\phi_{s}\left(q_{n_{m}}\right)+\phi_{s}\left(q_{n_{k}}\right)\right]\Big{\\}}\exp\left[\beta\left(\mu_{m}n_{m}+\mu_{k}n_{k}\right)\right]$
$\displaystyle+\sum\limits_{j\neq m}A_{m,j}=0.$ (4.5)

In the relation (4), the term $A_{m,j}$ is the symmetrized part of
$\delta\mathcal{I}_{m,k}$ determined from equation (4). So, for $A_{m,j}$ we
can write

$\displaystyle A_{m,j}$
$\displaystyle=\displaystyle\frac{1}{z^{0}_{m}z^{0}_{j}}\sum\limits_{n_{m}=0}^{1}\sum\limits_{n_{j}=0}^{1}\text{e}^{\beta\left(\mu_{m}n_{m}+\mu_{j}n_{j}\right)}\int\limits_{\omega_{m}}\int\limits_{\omega_{j}}\left[\delta\phi_{k}\left(q_{n_{m}}\right)-\delta\ln
z^{0}_{m}\right]$
$\displaystyle\times\exp\left[-\beta\Phi\big{(}q_{n_{m}},q_{n_{j}}\big{)}\right]\exp\bigg{\\{}-\beta\sum\limits_{s\neq
m,j}\left[\phi_{s}\left(q_{n_{m}}\right)+\phi_{s}\big{(}q_{n_{j}}\big{)}\right]\bigg{\\}}\mathrm{d}q_{n_{m}}\mathrm{d}q_{n_{j}}\,.$
(4.6)

Subsequent separation of variables makes it possible to obtain an equation for
the required potentials $\phi_{j}\left(q_{n_{i}}\right)$. Since both
$\delta\phi_{k}\left(q_{1_{m}}\right)$ and
$\delta\phi_{k}\left(q_{0_{m}}\right)$ are independent, after simplifications
of equations (4), we obtain a system of defining equations

$\displaystyle\exp\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]=\frac{1}{z_{k}^{0}}\sum\limits_{n_{k}=0}^{1}\text{e}^{\beta\mu_{k}n_{k}}\int\limits_{\omega_{k}}\mathrm{d}q_{n_{k}}\exp\Bigg{\\{}-\beta\bigg{[}\Phi\left(q_{n_{m}},q_{n_{k}}\right)+\sum\limits_{s\neq
m,k}^{Z}\phi_{s}\left(q_{n_{k}}\right)\bigg{]}\Bigg{\\}}.$ (4.7)

Equations (4.7) can be rewritten in another form, namely when the density
field is used as a variable that defines the system. Substituting (2.15) into
(4.7) we arrive at a description

$\exp\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]=\sum\limits_{n_{k}=0}^{1}\displaystyle\frac{\rho_{n_{k}}}{Q_{n_{k}}}\int\limits_{\omega_{k}}\mathrm{d}q_{n_{k}}\exp\Bigg{\\{}-\beta\bigg{[}\Phi\left(q_{n_{m}},q_{n_{k}}\right)+\sum\limits_{s\neq
m,k}^{Z}\phi_{s}\big{(}q_{n_{k}}\big{)}\bigg{]}\Bigg{\\}}.$ (4.8)

Equations (4.7) and (4.8) define the single-particle cell potentials under the
condition that the two-vertex diagrams in the cluster expansion of the
original partition function (3) are equal to zero. As was shown in [15], this
condition leads to the results that are equivalent to the quasichemical (or
Bether-Peierls) approximation. The inclusion of three-vertex diagrams in the
equations (4.7) or (4.8) is also possible but it would lead to more
complicated equations and, therefore, we will not consider it here.

Now, we show that when equation (4.7) is satisfied, each term ${A_{j}}$ in (4)
turns out to be zero, which indicates that the separation constant of the
variables in (4) chosen to be zero is correct.

Thus, substituting (4.7) into (4), we obtain

$A_{m,j}=\frac{1}{z^{0}_{m}}\sum\limits_{n_{m=0}}^{1}\text{e}^{\beta\mu_{m}n_{m}}\int\limits_{\omega_{m}}\left[\delta\phi_{k}\left(q_{n_{m}}\right)-\delta\ln
z^{0}_{m}\right]\exp\Big{[}-\beta\sum\limits_{s\neq
m}\phi_{s}\left(q_{n_{m}}\right)\Big{]}\mathrm{d}q_{n_{m}}\,.$ (4.9)

Using the definitions (2.11) and (2.12) in (4.9), we obtain

$A_{m,j}=\frac{1}{z^{0}_{m}}\sum\limits_{n_{m=0}}^{1}\text{e}^{\beta\mu_{m}n_{m}}\int\limits_{\omega_{m}}\delta\phi_{k}\left(q_{n_{m}}\right)\exp\Big{[}-\beta\sum\limits_{s\neq
m}\phi_{s}\left(q_{n_{m}}\right)\Big{]}\mathrm{d}q_{n_{m}}-\delta\ln
z_{m}^{0}\,.$ (4.10)

In turn, it is clear from (4) that the first sum in (4.10) is identical to the
variation of the second term, which proves that $A_{j}=0$ when choosing
single-particle potentials satisfying equation (4.7). As already noted for
specific calculations, (4.7) is preferable in the form (4.8) because this
option allows for an implicit relationship between the chemical potential and
the density to be replaced by an explicit one. Namely, for the chosen
$\rho_{n_{i}}$ the solution (4.8) is found with (2.12) taken into account.
Then, due to (2.17), the chemical potential and the free energy of the system
are determined. To calculate the latter, it is convenient to use the relation

$-\beta F=\sum\limits_{i=1}^{M}\big{(}\ln
z^{0}_{i}-\rho_{1_{i}}\mu_{i}\big{)}.$ (4.11)

Taking into account equation (2.17), expression (4.11) is represented in the
form

$-\beta F=\sum\limits_{i=1}^{M}\left[\ln
z^{0}_{i}-\rho_{1_{i}}\ln\left(\frac{\rho_{1_{i}}Q_{0_{i}}}{\rho_{0_{i}}Q_{1_{i}}}\right)\right].$
(4.12)

As a result, it follows from equation (4.7) that

$\sum\limits_{n_{j}=0}^{1}\rho_{n_{j}}f_{n_{i}n_{j}}=0.$ (4.13)

In order to prove (4.13), we multiply equation (4.8) by

$\frac{\rho_{n_{m}}}{Q_{n_{m}}}\exp\Big{[}-\beta\sum\limits_{s\neq
m}\phi_{s}\left(q_{n_{m}}\right)\Big{]}.$ (4.14)

After integration, we obtain

$1=\sum\limits_{n_{m}=0}^{1}\sum\limits_{n_{k}=0}^{1}\rho_{n_{m}}\rho_{n_{k}}\int\limits_{\omega_{m}}\int\limits_{\omega_{k}}\exp\left[-\beta\Delta\phi\left(q_{n_{m}},q_{n_{k}}\right)\right]\widehat{F}_{11}\left(q_{n_{m}}\right)\widehat{F}_{11}\left(q_{n_{k}}\right)\mathrm{d}q_{n_{m}}\mathrm{d}q_{n_{k}}\,.$
(4.15)

The identity obtained with the definition of (3.8) proves the validity of
equation (4.13). Hence, it follows that when the potentials are determined
from (4.8), the results of all three approaches are the same. Due to (4.13),
$\ln z_{i}$ coincides with $\ln z_{i}^{0}$ which is defined by the relation
(2.18).

Substituting (2.18) into the formula (4.12), we obtain an expression for the
free energy of the system in the form

$-\beta F=\sum\limits_{i=1}^{M}\left(\ln
Q_{0_{i}}-\ln\rho_{0_{i}}+\rho_{1_{i}}\ln\rho_{0_{i}}+\rho_{1_{i}}\ln
Q_{1_{i}}-\rho_{1_{i}}\ln\rho_{1_{i}}-\rho_{1_{i}}\ln Q_{0_{i}}\right)$

or

$-\beta F=\sum\limits_{i=1}^{M}\left(\rho_{0_{i}}\ln Q_{0_{i}}+\rho_{1_{i}}\ln
Q_{1_{i}}-\rho_{0_{i}}\ln\rho_{0_{i}}-\rho_{1_{i}}\ln\rho_{1_{i}}\right).$
(4.16)

Such a representation of the right-hand side of equation (4.16) corresponds to
a configurational integral in the form

$Q_{N}=Q_{N}^{0}=\sum\limits_{i=1}^{M}\frac{Q_{1_{i}}^{\rho_{1_{i}}}Q_{0_{i}}^{\rho_{0_{i}}}}{\rho_{1_{i}}^{\rho_{1_{i}}}\rho_{0_{i}}^{\rho_{0_{i}}}}\,.$
(4.17)

##  5 Verification of thermodynamic self-consistency of the theory

The expression for $Q_{N}^{0}$ in the form (4.17) can be obtained using the
Hamiltonian

$H_{N}^{0}=\sum\limits_{i=1}^{M}\sum\limits_{j\left(i\right)}^{Z}\phi_{j}\big{(}q_{n_{i}}\big{)}.$
(5.1)

To this end, it is necessary to consider the states corresponding to (5.1)
based on the methods for forming local equilibrium distributions [16]. We show
that the conditions (4.8), when the contributions of the third and the
subsequent virial coefficients are not taken into account, give the identical
$\rho_{i}$ determined by formulae

$\mu_{i}=\frac{\partial F}{\partial\rho_{i}}$ (5.2)

and

$\rho_{1_{i}}=\frac{\partial\ln
Z_{M}^{0}}{\partial\left(\beta\mu_{i}\right)}.$ (5.3)

In order to verify the thermodynamic consistency of the theory, we first
consider the validity of equation (5.3). Using (2.18), we obtain

$\displaystyle\rho_{1_{i}}$
$\displaystyle=\frac{\partial\ln\left(Q_{0_{i}}+\text{e}^{\beta\mu_{i}}Q_{1_{i}}\right)}{\partial\left(\beta\mu_{i}\right)}+\sum\limits_{k\neq
i}^{M}\frac{\partial\ln z_{k}^{0}}{\partial\left(\beta\mu_{i}\right)}$
$\displaystyle=\frac{1}{\left(Q_{0_{i}}+\text{e}^{\beta\mu_{i}}Q_{1_{i}}\right)}\left[\text{e}^{\beta\mu_{i}}Q_{1_{i}}+\frac{\partial
Q_{0_{i}}}{\partial\left(\beta\mu_{i}\right)}+\frac{\partial
Q_{1_{i}}}{\partial\left(\beta\mu_{i}\right)}\right]+\sum\limits_{k\neq
1}^{M}\frac{\partial\ln z_{k}^{0}}{\partial\left(\beta\mu_{i}\right)}$
$\displaystyle=\rho_{1_{i}}+\frac{1}{\left(Q_{0_{i}}+\text{e}^{\beta\mu_{i}}Q_{1_{i}}\right)}\left[\frac{\partial}{\partial\left(\beta\mu_{i}\right)}Q_{0_{i}}+\text{e}^{\beta\mu_{i}}\frac{\partial
Q_{1_{i}}}{\partial\left(\beta\mu_{i}\right)}\right]+\sum\limits_{k\neq
i}^{M}\frac{\partial\ln z_{k}^{0}}{\partial\left(\beta\mu_{i}\right)}\,.$
(5.4)

It follows from (5) that (4.7) must satisfy the additional condition

$\frac{\partial}{\partial\left(\beta\mu_{i}\right)}Q_{0_{i}}+\text{e}^{\beta\mu_{i}}\frac{\partial
Q_{1_{i}}}{\partial\left(\beta\mu_{i}\right)}=0.$ (5.5)

Then, $\rho_{1_{i}}$ calculated from the formulae (5.3) and (2.16) will
coincide. To prove (5.5) we differentiate the condition (4.7) with respect to
$\beta{\mu_{m}}$

$\displaystyle\frac{\partial}{\partial\left(\beta\mu_{m}\right)}\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]\exp\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]=-\frac{1}{z_{k}^{0}}\frac{\partial\ln
z_{k}^{0}}{\partial\left(\beta\mu_{m}\right)}\sum\limits_{n_{k}=0}^{1}\exp\left(\beta\mu_{k}n_{k}\right)$
(5.6)
$\displaystyle\times\int\limits_{\omega_{k}}\mathrm{d}q_{n_{k}}\exp\bigg{\\{}-\beta\bigg{[}\Phi\left(q_{n_{m}},q_{n_{k}}\right)+\sum\limits_{s\neq
m,k}^{Z}\phi_{s}\left(q_{n_{k}}\right)\bigg{]}\bigg{\\}}+\frac{1}{z_{k}^{0}}\sum\limits_{n_{k}=0}^{1}\exp\left(\beta\mu_{k}n_{k}\right)$
$\displaystyle\times\int\limits_{\omega_{k}}\left[\sum\limits_{s=m,k}^{Z}\frac{\partial\phi_{s}\left(q_{n_{k}}\right)}{\partial\left(\beta\mu_{m}\right)}\right]\exp\bigg{\\{}-\beta\bigg{[}\Phi\left(q_{n_{m}},q_{n_{k}}\right)+\sum\limits_{s\neq
m,k}^{Z}\phi_{s}\left(q_{n_{k}}\right)\bigg{]}\bigg{\\}}\mathrm{d}q_{n_{k}}.$

Multiplying (4.10) by
$\displaystyle\frac{1}{z_{m}^{0}}\exp\bigg{[}-\beta\sum\limits_{s\neq
k,m}\phi_{s}\left(q_{n_{m}}\right)\bigg{]}$, after integration over
$q_{n_{m}}$ we obtain

$\displaystyle\frac{1}{z_{m}^{0}}\int\limits_{\omega_{m}}\frac{\partial}{\partial\left(\beta\mu_{m}\right)}\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]\exp\left[-\beta\sum\limits_{s\neq
k}^{Z}\phi_{s}\left(q_{n_{m}}\right)\right]\mathrm{d}q_{n_{m}}$ (5.7)
$\displaystyle=-\frac{1}{z_{m}^{0}}\frac{\partial\ln
z_{k}^{0}}{\partial\left(\beta\mu_{m}\right)}\int\limits_{\omega_{m}}\exp\left[-\beta\sum\limits_{s\neq
k}^{Z}\phi_{s}\left(q_{n_{m}}\right)\right]\mathrm{d}q_{n_{m}}$
$\displaystyle+\frac{1}{z_{k}^{0}z_{m}^{0}}\sum\limits_{n_{k}=0}^{1}\exp\left(\beta\mu_{k}n_{k}\right)\int\limits_{\omega_{k}}\int\limits_{\omega_{m}}\mathrm{d}q_{n_{m}}\mathrm{d}q_{n_{k}}\exp\left[-\beta\Phi\left(q_{n_{k}},q_{n_{m}}\right)\right]$
$\displaystyle\times\left[-\beta\sum\limits_{s\neq
m,k}^{Z}\displaystyle\frac{\partial\phi_{s}\left(q_{n_{k}}\right)}{\partial\left(\beta\mu_{m}\right)}\right]\exp\left[-\beta\sum\limits_{s\neq
k}^{Z}\phi_{s}\left(q_{n_{m}}\right)\right].$

Multiplying (5.7) by $\sum\limits_{n_{m}}^{1}\text{e}^{-\beta\mu_{m}n_{m}}$
and taking into account (2.11) and (4.7), we have

$\displaystyle\frac{2}{z_{m}^{0}}\sum\limits_{n_{m}=0}^{1}\text{e}^{-\beta\mu_{m}n_{m}}\int\limits_{\omega_{m}}\frac{\partial}{\partial\left(\beta\mu_{n}\right)}\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]\exp\left[-\beta\sum\limits_{s\neq
k}^{Z}\phi_{s}\left(q_{n_{m}}\right)\right]\mathrm{d}q_{m}$ (5.8)
$\displaystyle=-\frac{\partial\ln
z_{k}^{0}}{\partial\left(\beta\mu_{m}\right)}+\frac{1}{z_{k}^{0}}\sum\limits_{n_{k}=0}^{1}\exp\left(\beta\mu_{k}n_{k}\right)\int\limits_{\omega_{k}}\left[-\beta\sum\limits_{S\neq
k}^{Z}\frac{\partial\phi_{s}\left(q_{n_{k}}\right)}{\partial\left(\beta\mu_{m}\right)}\right]\exp\left[-\beta\sum\limits_{s\neq
k}^{Z}\phi_{s}\left(q_{n_{m}}\right)\right]\mathrm{d}q_{n_{k}}.$

As a result,

$\displaystyle\frac{2}{z_{m}^{0}}\sum\limits_{n_{m}=0}^{1}\text{e}^{-\beta\mu_{m}n_{m}}\int\limits_{\omega_{m}}\frac{\partial}{\partial\left(\beta\mu_{m}\right)}\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]\exp\left[-\beta\sum\limits_{s\neq
k}^{Z}\phi_{s}\left(q_{n_{m}}\right)\right]\mathrm{d}q_{n_{m}}$ (5.9)
$\displaystyle=\frac{\partial\ln
z_{k}^{0}}{\partial\left(\beta\mu_{m}\right)}+\frac{1}{z_{k}^{0}}\sum\limits_{n_{k}=0}^{1}\exp\left(\beta\mu_{k}n_{k}\right)\frac{\partial
Q_{n_{k}}}{\partial}\left(\beta\mu_{m}\right)=0.$

Performing the summation in (5.9) with respect to $k\neq m$, we find

$\sum\limits_{n_{m}=0}^{1}\text{e}^{-\beta\mu_{m}n_{m}}\frac{\partial
Q_{n_{m}}}{\partial\left(\beta\mu_{m}\right)}=0.$ (5.10)

Likewise, we show that although in accordance with (4.7) each
$\phi_{s}\left(q_{n_{k}}\right)$ depends on the chemical potentials in the
entire range of the $k$-th node, it turns out that

$\frac{\partial z^{0}_{m}}{\partial\left(\beta\mu_{l}\right)}=0$ (5.11)

for $l\neq m$.

In combination with the condition (5.10), this will ultimately prove the
thermodynamic consistency of the theory when calculating the thermodynamic
potentials $F$ and $\Omega$ and their first derivatives, determined by
formulae (2.17) and (5.3). The proof of the condition (5.11) is carried out
similarly to (5.5).

In order to reduce the transformations and make them more transparent, we
write a relation analogous to (5.6) in the form

$\displaystyle\bigg{\\{}\displaystyle\frac{\partial}{\partial\left(\beta\mu_{s}\right)}\left[-2\beta\phi_{k}\left(q_{n_{m}}\right)\right]\bigg{\\}}\exp\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]=\displaystyle-\frac{\partial\ln
z_{k}^{0}}{\partial\left(\beta\mu_{s}\right)}\exp\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]$
(5.12)
$\displaystyle\displaystyle+\frac{1}{z_{k}^{0}}\sum\limits_{n_{k}=0}^{1}\exp\left(\beta\mu_{k}n_{k}\right)\int\limits_{\omega_{k}}\left[\sum\limits_{s\neq
m}^{Z}\displaystyle\frac{\partial\beta\phi_{s}\left(q_{n_{k}}\right)}{\partial\left(\beta\mu_{s}\right)}\right]\exp\Bigg{\\{}-\beta\bigg{[}\Phi\left(q_{n_{m}},q_{n_{k}}\right)$
$\displaystyle+\sum\limits_{s\neq
m,k}^{Z}\phi_{s}\left(q_{n_{k}}\right)\bigg{]}\Bigg{\\}}\mathrm{d}q_{n_{k}}.$

Let us apply the following operation to relation (5.12)

$\widehat{\pi}=\sum\limits_{k\neq
m}^{Z}\frac{1}{z_{m}^{0}}\sum\limits_{n_{m}}^{1}\exp\left(-\beta\mu_{m}n_{m}\right)\int\limits_{\omega_{m}}\exp\left[-\beta\sum\limits_{s\neq
k}^{Z}\phi_{s}\left(q_{n_{m}}\right)\right]\mathrm{d}q_{n_{m}}...$ (5.13)

writing equation (5.12) in a symbolic form

$A=B+C.$ (5.14)

When converting the left-hand side of expression (5.14), the summation sign
with respect to $k\neq m$ is put under the integration sign with respect to
$\omega_{m}$, and taking into account (2.12), we write the result in the form

$\widehat{\pi}A=2\frac{1}{z_{m}^{0}}\sum\limits_{n_{m}=0}^{1}\exp\left(-\beta\mu_{m}n_{m}\right)\frac{\partial
Q_{n_{m}}}{\partial\left(\beta\mu_{s}\right)}.$ (5.15)

Since in (5.15) $s\neq m$ and $\mu_{m}$ and $\mu_{s}$ are independent
variables, we switch the summation over $n_{m}$ and differentiation with
respect to $\mu_{s}$, which due to (2.11) allows us to write

$\widehat{\pi}A=2\frac{\partial\ln
z_{m}^{0}}{\partial\left(\beta\mu_{s}\right)}.$ (5.16)

The result of the transformation of the expression $B$ is equal to 1. It can
be obtained if the operation (5.3) is applied to
$\exp\left[-\beta\phi_{k}\left(q_{n_{m}}\right)\right]$ taking the constant as
the sign of operation (5.13). As a result, direct application of
$\widehat{\pi}$ gives

$\widehat{\pi}B=-\frac{\partial\ln
z_{k}^{0}}{\partial\left(\beta\mu_{s}\right)}.$ (5.17)

To transform the expression $C$, we switch the summation and integration
operations in (5.13), place the operator $\widehat{\pi}$ under the integration
sign with respect to $\omega_{k}$ and take into account that according to
(4.7)

$\widehat{\pi}\exp\big{\\{}-\beta\left[\Phi\left(q_{n_{k}},q_{n_{m}}\right)\right]\big{\\}}=\exp\left[-\beta\phi_{m}\left(q_{n_{k}}\right)\right].$
(5.18)

Then, $\widehat{\pi}C$ becomes of the form

$\widehat{\pi}C=\frac{1}{z_{k}^{0}}\sum\limits_{n_{k}=0}^{1}\exp\left(\beta\mu_{k}n_{k}\right)\frac{\partial\ln
Q_{k}}{\partial\left(\beta\mu_{s}\right)}.$ (5.19)

Taking (2.11) into account, condition (5.19) gives

$\widehat{\pi}C=\frac{\partial\ln
z_{k}^{0}}{\partial\left(\beta\mu_{s}\right)}.$ (5.20)

The results (5.17) and (5.20) show that according to (5.14), the expression
(5.16) is identically equal to zero. This, in turn, proves the validity of
(5.12).

Similarly to the consistency of (5.3) and (2.15) proved by derivations
(5)–(5.20), one can show that (5.2) is consistent with (2.16). To this end, it
is necessary to perform transformations similar to those performed for (4.7)
over (4.8), differentiating each of the equations (5)–(5.20) with respect to
thermodynamic variables $\rho_{n_{i}}$ but not as it was before with respect
to $\mu_{s}$. Further analysis of the consistency problem involves comparing
the results of calculating the second derivatives of thermodynamic potentials,
and will be considered separately.

##  6 Discussion

In this paper, for a lattice model with the Hamiltonian (2.1) and a pair
interaction potential $\Phi(q_{n_{i}},q_{n_{j}})$ we formulate a reference
system with the Hamiltonian (2.2) and a single-particle cell potential
$\phi_{j}(q_{n_{i}})$. The potential $\phi_{j}(q_{n_{i}})$ can be interpreted
as the mean potential exerted by a particle in the lattice cell $\omega_{j}$
on a particle in the lattice cell $\omega_{i}$. It is shown that the system
with the single-cell potential $\phi_{j}(q_{n_{i}})$ reduces to the
description of a Fermi-like lattice model in an external field. Using this
system as the reference system and renormalized Mayer functions in the form
(3.5) with $\Delta\phi(q_{n_{i}},q_{n_{j}})$ in the form (3.2), the
generalized cluster expansion for the free energy of the considered system is
obtained. The cell potentials are calculated from the condition of the minimum
difference of thermodynamic properties of the systems with the Hamiltonians
(2.1) and (2.2). Such a procedure is considered under the condition that the
two-vertex diagrams in the cluster expansion of the partition function for the
system with the Hamiltonian (2.1) are equal to zero. As a result, for the
single-particle cell potential, a system of equations (4.7) was obtained. This
system can also be presented in the form of equation (4.8). As it was noted
previously in the reference [15], such a description is equivalent to the
quasichemical approximation. We should note that the descriptions of the
considered system in the framework of the Hamiltonians (2.1) and (2.2) are
exactly equivalent only in the case when all the terms in the cluster
expansion (3) are taken into account at the calculation of the single-particle
potentials. There are two principal differences between the considered
equation (4.7) and the corresponding equations in traditional approaches such
as the mean field approximation formulated in the framework of the field
theoretical approach [17, 18] or the density functional approach [19] well
developed for non-lattice fluid systems. The first difference is connected
with the presence of the interparticle potential $\Phi(q_{n_{i}},q_{n_{j}})$
in equation (4.7) in the exponential form. The second one is connected with
the inequality $s\neq m,k$ in the exponent of equation (4.7). It means that
instead of the singlet distribution function $\hat{F}_{11}(q_{n_{k}})$, which
usually appears in the mean field approximation, in the approach considered,
the function
$\hat{F}_{11}(q_{n_{k}})\exp\left[\beta\phi_{m}(q_{n_{k}})\right]$ appears.
This is the result of the peculiarity of the Mayer functions in the form (3.5)
with $\Delta\phi(q_{n_{i}},q_{n_{j}})$ in the form (3.2). If we neglect the
condition $s\neq m,k$ due to the renormalization condition (2.21), the
equation (4.7) can be rewritten in the form

$\text{e}^{-\beta\phi_{k}(q_{n_{m}})}=1+\sum\limits_{n_{k}=0}^{1}\rho_{n_{k}}\int\mathrm{d}q_{n_{k}}\big{\\{}\exp\left[-\beta\Phi\left(q_{n_{m}},q_{n_{k}}\right)\right]-1\big{\\}}\hat{F}_{11}(q_{n_{k}})\,,$
(6.1)

which after linearization of the exponents
$\exp\left[{-\beta\phi_{k}(q_{n_{m}})}\right]$ and
$\exp\left[-\beta\Phi\left(q_{n_{m}},q_{n_{k}}\right)\right]$ leads to a
traditional form for the mean field approximation

$\phi_{k}(q_{n_{m}})=\sum\limits_{n_{k}=0}^{1}\rho_{n_{k}}\int\mathrm{d}q_{n_{k}}\Phi\left(q_{n_{m}},q_{n_{k}}\right)\hat{F}_{11}\left(q_{n_{k}}\right).$
(6.2)

The equation (4.7) describes the single-particle cell potentials
$\phi_{k}(q_{1_{m}})$ for real particles, but for vacancies it is more of a
problem. This problem is similar to the description of solvophobic interaction
in the theory of solutions [20, 21] and for a correct description of single-
particle cell potential $\phi_{k}(q_{0_{m}})$ for vacancies, at least three-
vertex diagrams in the equation (4.7) should be included. In the two-vertex
diagram approximation for $\phi_{k}(q_{0_{m}})$ due to inequality $s\neq m,k$
there appears only some constant corresponding to the change of the chemical
potential due to the creation of a vacancy.

The theory presented here is easily generalized when it is necessary to take
into account a larger number of possible states. This is achieved by expanding
the possible values of the occupation numbers, when $n_{i}=0,1,2\dots$. In
this case, for every microconfiguration given by the set of values
$\\{n_{1},\dots,n_{m}\\}$, the expression for the Hamiltonian is completely
preserved in the form (2.1). So, for example, for $n=0,1,2$ in all formulae
(2.2)–(5.20) one should remember that

$\displaystyle\Phi\big{(}q_{n_{i}},q_{n_{j}}\big{)}=\left\\{\begin{array}[]{l}{0,\,n_{i},n_{j}=0}\,,\\\
{h\left(q_{i},q_{j}\right),\,n_{i},n_{j}=1}\,,\\\
{h\big{(}q_{i},q_{j}\big{)}+h\big{(}q_{i},q^{\prime}_{j}\big{)},\,n_{i},n_{j}=2}\,,\\\
{h\left(q_{i},q_{j}\right)+h\big{(}q_{i},q^{\prime}_{j}\big{)}+h\big{(}q^{\prime}_{i},q_{j}\big{)}+h\big{(}q^{\prime}_{i},q^{\prime}_{j}\big{)},\,n_{i},n_{j}=4}\,.\end{array}\right.$
(6.7)

Here, $h(q_{i},q_{j})$ is the intermolecular interaction potential of two
particles in positions $q_{i}$ and $q_{j}$, when
$(q_{i},q^{\prime}_{i})\in\omega_{i}$, $(q_{j},q^{\prime}_{j})\in\omega_{j}$.

In addition to taking into account (6.7) in all expressions (2.2)–(5.20), the
summation over $n=0,1$ must be extended to the case $n=0,1,2$, additionally
taking into account that

$\int\limits_{\omega_{i}}\mathrm{d}q_{0_{i}}=\frac{1}{\omega_{i}}\int\limits_{\omega_{i}}\mathrm{d}q_{i},\quad\int\limits_{\omega_{i}}\mathrm{d}q_{1_{i}}=\int\limits_{\omega_{i}}\mathrm{d}q_{i},\quad\int\limits_{\omega_{i}}\mathrm{d}q_{2_{i}}=\int\limits_{\omega_{i}}\mathrm{d}q_{i}\int\limits_{\omega_{i}}\mathrm{d}q^{\prime}_{i}\,.$
(6.8)

Thus, it is shown that the properties of a condensed system can be described
by combining the model of an ideal crystal with a group expansion over the
modified Mayer functions that impose correlations on the properties of an
ideal crystal.

It seems justified to extend the developed approach to take into account long-
range effects in essentially inhomogeneous media. Let us demonstrate the
possibility of such propagation using the example when the inhomogeneity of
the medium is described by different sizes of microcells, provided that in
each of them there is one particle. Representing the energy of the system by
short-range $\Phi(i,j)$ and long-range $V(i,j)$ potentials of pair
interactions, respectively, for particles in positions $q_{i}$, $q_{j}$, we
present the configurational integral $Q_{N}$ of the system in the form

$Q_{N}=Q_{N}^{0}\bigg{\langle}\exp\bigg{[}-\beta\sum\limits_{i<j}^{N}V(i,j)\bigg{]}\prod\limits_{i<j}^{N}\big{[}1+f(i,j)\big{]}\bigg{\rangle}_{0}\,,$
(6.9)

where

$f(i,j)=\exp\big{\\{}-\beta\left[\Phi\left(i,j\right)\right]-\phi_{j}\left(i\right)-\phi_{i}\left(j\right)\big{\\}}-1\,,$
(6.10)
$Q_{N}^{0}=\prod\limits_{i=1}^{N}Q_{i}\,,{\rm\;\;\;\;\;\;\;\;}Q_{i}=\int\limits_{v_{i}}\exp\Big{[}-\beta\sum\limits_{k\neq
i}\phi_{k}\left(i\right)\Big{]}\mathrm{d}q_{i}\,.$ (6.11)

$Q_{N}^{0}$ is the configurational integral of an ideal crystal, expressed
through single-particle cell potentials of mean forces
$\phi_{j}\left(i\right)$, $f(i,j)$ is a renormalized Mayer function, the angle
brackets $\left\langle...\right\rangle_{0}$ denote averaging over the
equilibrium states of the reference system. The subsequent cumulant expansion
of the expression (2.2) with respect to the functions (2.3) allows one to
write in the approximation of the second virial coefficient

$\ln Q_{N}=\ln Q_{L}^{0}+\ln Q_{N}^{0}+\sum\limits_{i,j}\left\langle
f(i,j)g(i,j)\right\rangle_{0}+...\,,$ (6.12)

$Q_{L}^{0}=\Big{\langle}\exp\Big{[}-\beta\sum\limits_{i<j}V(i,j)\Big{]}\Big{\rangle}_{0}.$
(6.13)

The averaging in relations (6.11) and (6.13) is realized by multiplying the
unary distribution functions $F_{11}(i)$ and $F_{11}(j)$

$F_{11}\left(i\right)=\frac{1}{Q_{i}}\exp\Big{[}-\beta\sum\limits_{k\neq
1}\phi_{k}(i)\Big{]}.$ (6.14)

A special feature of relation (6.13) is that the renormalized Mayer function
$f(i,j)$ is modulated here by a binary function, for which a consistent
calculation scheme was developed in [10]. Applying the procedure of self-
consistent calculation of the potentials $\phi_{j}\left(i\right)$ in
accordance with the foregoing, we arrive at a closed system of equations of
the form

$\exp\left[-\beta\phi_{j}\left(i\right)\right]=\frac{1}{Q_{j}}\int\limits_{v_{j}}g(i,j)\exp\left[-\beta\Phi\left(i,j\right)\right]\times\exp\Big{[}-\beta\sum\limits_{k\neq
i,j}\phi_{k}\left(j\right)\Big{]}\mathrm{d}j.$ (6.15)

The system of equations (6.15) differs from the one used earlier because its
kernel, in addition to the point short-range potential, contains a binary
function for a system of particles with a Coulomb interaction, the expression
for which has the form [10]

$\displaystyle F_{2}(i,j)$ $\displaystyle=F_{11}(i)F_{11}(j)g(i,j)\,,$
$\displaystyle g(i,j)$ $\displaystyle=\exp[-\beta u(i,j)]\,,$ (6.16)

where $u(i,j)$ is a screening potential. Such a form solves the problem of the
divergence of the integrals when calculating the free energy (6.11). As a
result, long-range effects apply also to the renormalization of single-
particle cell potentials. In expression (6.16), only the initial terms of the
series are written corresponding to the Debye description of ion systems. A
more complete representation for relation (6.16) follows when calculating
additional terms of the distribution function $g(i,j)$ using collective
variables [10]. A generalization of Debye screening to the system of mobile
ions in lattice models was recently discussed in [4]. It leads to the change
of the traditional inverse Debye length $\kappa$ to the new one
$\kappa=\left[\beta e^{2}c\left(1-c\right)/\left(\varepsilon_{0}\varepsilon
h^{3}\right)\right]^{1/2}$, where $e$ is the charge of a mobile ion,
$\varepsilon_{0}$ is the dielectric permittivity of the vacuum, $\varepsilon$
is the relative dielectric permittivity of the medium, $c$ is the
concentration of mobile ions, $h$ is the length of the cell.

The peculiarity of the equations in the case under consideration is due to the
fact that the function $g(i,j)$ is determined by averaging not over the states
of an ideal gas, but over the states of an ideal crystal. Accordingly, for
example for ${i}=1$, ${j}=2$, the function $g(i,j)$ is defined by an
expression of the form

$\displaystyle
g(1,2)=\int\limits_{v_{3}}...\int\limits_{v_{N}}\exp\bigg{[}-\beta\sum\limits_{l<m=1}^{N}V\left(l,m\right)\bigg{]}F_{11}\left(3\right)F_{11}\left(4\right)...F_{11}\left(N\right)\mathrm{d}3\mathrm{d}4...\mathrm{d}N.$
(6.17)

There are reasons to believe that the appearance of additional Gaussian in
(6.17) introduced by the functions $F_{11}$ will lead to an improvement in the
convergence of the series (6.16) obtained in the framework of the collective
variables method [10] since a similar procedure was successful in constructing
the description using an effective potential [22].

##  7 Conclusions

A method for modifying the cluster virial expansion, which makes it possible
to describe condensed state, is proposed. The method is based on replacing the
averaging over the states of an ideal gas by the averaging over the states of
an ideal crystal. The approach outlined differs from the well-known
perturbation theory because the basis with respect to which the expansion is
performed is not a part of the original Hamiltonian but is introduced
independently. The basis distribution conveys the main features of the solid
state of matter, where the motion of molecules is of an oscillatory nature
with respect to the lattice sites. As a result, the Hamiltonian of the basic
reference system is represented by the sum of single-particle cell mean force
potentials. It is shown how the Mayer functions can be modified so that they
would act as a small parameter for the subsequent expansion of the
thermodynamic potential by cumulant expansions. From the condition of
independence of the initial partition function on the introduced potentials of
the mean forces, a system of integral equations determining the above
potentials is obtained. The expression for the free energy functional is
obtained with the two first terms of its expansion in correlations taken into
account. It is shown that the obtained equations satisfy several optimization
conditions for the parameters characterizing the properties of the reference
system and the initial system. The thermodynamic consistency of various
methods for calculating the thermodynamic characteristics of a condensed
medium both in the canonical and in the grand canonical ensembles is proved. A
possibility of using the developed approach to take into account not only
short-range but also long-range interactions is shown.

## Acknowledgements

This project has received funding from European Unions Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie (grant agreement No
734276), the Belarusian Republican Foundation for Fundamental Research (grant
No $\Phi$16K-061) and the State Fund for Fundamental Research of Ukraine
(grant No $\Phi$73/26-2017). We thank Ivan Kravtsiv and Dung di Caprio for the
careful reading of the manuscript and useful comments.

## References

  * [1] Prigogine I., Bellemans A., Mathot V., The Molecular Theory of Solutions, North-Holland Publishing Company, Amsterdam, 1957.
  * [2] Frenkel J.I., Kineticheskaya Teoriya Zhidkostei, AN SSSR, Moskva–Leningrad, 1945 (in Russian), [Kinetic Theory of Liquids, Dover Publications, Inc., New York, 1955].
  * [3] Yukhnovskii I.R., Phase Transitions of the Second Order: Collective Variables Method, Word Scientific, Singapore, 1987.
  * [4] Bokun G., di Caprio D., Holovko M., Vikhrenko V., J. Mol. Liq., 2018, 270, 183,
doi:10.1016/j.molliq.2018.03.123.

  * [5] Bisquert J., Vikhrenko V.S., J. Phys. Chem. B, 2004, 108, 2313, doi:10.1021/jp035395y.
  * [6] Ciach A., Góźdź W.T., Condens. Matter Phys., 2010, 13, 23603, doi:10.5488/CMP.13.23603.
  * [7] Narkevich I.I., Physica A, 1982, 112, 167, doi:10.1016/0378-4371(82)90213-8.
  * [8] Di Caprio D., Badiali J.P., Holovko M., J. Phys. A: Math. Theor., 2009, 42, 214038,
doi:10.1088/1751-8113/42/21/214038.

  * [9] Rott L.A., Statistical Theory of Molecular Systems, Nauka, Moscow, 1979 (in Russian).
  * [10] Yukhnovskii I.R., Holovko M.F., Statistical Theory of Classical Equilibrium Systems, Naukova Dumka, Kiev, 1980 (in Russian).
  * [11] Van Kampen N.G., Physica, 1961, 27, 783, doi:10.1016/0031-8914(61)90097-0.
  * [12] Kubo R., J. Phys. Soc. Jpn., 1962, 17, 1100, doi:10.1143/JPSJ.17.1100.
  * [13] Andersen H.C., Weeks J.D., Chandler D., Phys. Rev. A, 1971, 4, 1597, doi:10.1103/PhysRevA.4.1597.
  * [14] Barker J.A., Henderson D., Rev. Mod. Phys., 1976, 48, 587, doi:10.1103/RevModPhys.48.587.
  * [15] Argyrakis P., Groda Y.G., Bokun G.S., Vikhrenko V.S., Phys. Rev. E, 2001, 64, 066108,
doi:10.1103/PhysRevE.64.066108.

  * [16] Zubarev D.N, Nonequilibrium Statistical Thermodynamics, Nauka, Moscow, 1971 (in Russian).
  * [17] Di Caprio D., Stafiej J., Holovko M., Kravtsiv I., Mol. Phys., 2011, 109, 695,
doi:10.1080/00268976.2010.547524.

  * [18] Kravtsiv I., Patsahan T., Holovko M., di Caprio D., J. Chem. Phys., 2015, 142, 194708, doi:10.1063/1.4921242.
  * [19] Evans R., In: Fundamentals of Inhomogeneous Fluids, Henderson D. (Ed.), Marcel Dekker, New York, 1992, 85–175.
  * [20] Ronis D., Martina E., Deutch J.M., Chem. Phys. Lett., 1977, 46, 53, doi:10.1016/0009-2614(77)85161-0.
  * [21] Bandura A.V., Holovko M.F., Lvov S.N., J. Mol. Liq., 2018, 270, 52, doi:10.1016/j.molliq.2018.01.015.
  * [22] Ma S.-K., Modern Theory of Critical Phenomena, Westview Press, New York, 1976.

Ukrainian 0

Ãðóïîâå ðîçâèíåííÿ äëÿ îïèñó êîíäåíñîâàíèõ ñèñòåì: ïäõä êðèñòàëíèõ êîìðîê Ã.Ñ.
Áîêóí${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ì.Ô.
Ãîëîâêî${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ Áëîðóñüêèé äåðæàâíèé òåõíîëîãíèé óíâåðñèòåò,
âóë. Ñâåðäëîâà, 13à, 220006 Ìíñüê, Áëîðóñü ${}^{\textsf{\footnotesize 2}}$
íñòèòóò ôçèêè êîíäåíñîâàíèõ ñèñòåì ÍÀÍ Óêðà¿íè, âóë. Ñâíööüêîãî, 1, 79011
Ëüââ, Óêðà¿íà