# The role of attraction in the phase diagrams and melting scenarios of
generalized 2D Lennard-Jones systems

Elena N. Tsiok  Institute of High Pressure Physics RAS, 108840 Kaluzhskoe
shosse, 14, Troitsk, Moscow, Russia  Yuri D. Fomin  Institute of High Pressure
Physics RAS, 108840 Kaluzhskoe shosse, 14, Troitsk, Moscow, Russia  Moscow
Institute of Physics and Technology (National Research University), 9
Institutskiy Lane, Dolgoprudny, Moscow region, 141701, Russia  Evgenii A.
Gaiduk  Institute of High Pressure Physics RAS, 108840 Kaluzhskoe shosse, 14,
Troitsk, Moscow, Russia  Elena E. Tareyeva  Institute of High Pressure Physics
RAS, 108840 Kaluzhskoe shosse, 14, Troitsk, Moscow, Russia  Valentin N. Ryzhov
<EMAIL_ADDRESS>Institute of High Pressure Physics RAS, 108840
Kaluzhskoe shosse, 14, Troitsk, Moscow, Russia  Pavel A. Libet  Institute of
High Pressure Physics RAS, 108840 Kaluzhskoe shosse, 14, Troitsk, Moscow,
Russia  Bauman Moscow State Technical University, 2nd Baumanskaya street 5,
105005 Moscow, Russia  Nikita A. Dmitryuk  Bauman Moscow State Technical
University, 2nd Baumanskaya street 5, 105005 Moscow, Russia  Nikita P.
Kryuchkov  Bauman Moscow State Technical University, 2nd Baumanskaya street 5,
105005 Moscow, Russia  Stanislav O. Yurchenko<EMAIL_ADDRESS>Bauman
Moscow State Technical University, 2nd Baumanskaya street 5, 105005 Moscow,
Russia

###### Abstract

Monolayer and two-dimensional (2D) systems exhibit rich phase behavior,
compared with 3D systems, in particular, due to the hexatic phase playing a
central role in melting scenarios. The attraction range is known to affect
critical gas-liquid behavior (liquid-liquid in protein and colloidal systems),
but the effect of attraction on melting in 2D systems remains unstudied
systematically. Here, we reveal how the attraction range affects the phase
diagrams and melting scenarios in a 2D system. Using molecular dynamics
simulations we considered the generalized Lennard-Jones system with a fixed
repulsion branch and different power indices of attraction, from long-range
dipolar to short-range sticky-spheres-like. A drop in the attraction range has
been found to reduce the temperature of the gas-liquid critical point,
bringing it closer to the gas-liquid-solid triple point. At high-temperatures,
attraction does not affect the melting scenario that proceeds through the
cascade of solid-hexatic (Berezinski-Kosterlitz-Thouless) and hexatic-liquid
(first-order) phase transitions. In the case of dipolar attraction, we
observed _two triple points_ , inherent in a 2D system: hexatic-liquid-gas and
crystal-hexatic-gas, the temperature of the crystal-hexatic-gas triple point
is _below_ the hexatic-liquid-gas triple point. This observation may have far-
reaching consequences for future studies, since phase diagrams determine
possible routes of self-assembly in molecular, protein, and colloidal systems,
whereas the attraction range can be adjusted with complex solvents and
external electric or magnetic fields. The results obtained may be widely used
in condensed matter, chemical physics, materials science, and soft matter.

##  I Introduction

Understanding phase transitions in 2D systems has prominent importance in a
number of areas, from photonics and electronics, to novel materials and
biotechnologies, since the knowledge on phase behavior opens a way to design
systems with desired properties. Despite many studies, the fundamental
questions here are still related to the effect of particular interaction
between individual particles on their collective behavior. For classical
systems, one of the simplest models capable of reproducing their behavior,
including gas, liquid, and solid phases, is the Lennard-Jones (LJ) system. The
LJ model is widely used for analysis of behavior in molecular, protein,
polymer, emulsion, and colloidal soft matter Fernandez-Nieves and Puertas
(2016). The generalized LJ potential (or the LJn-m-potential, where indices
$n$ and $m$ are responsible for the algebraic branches of repulsion and
attraction, see below) is a suitable model for studies, aimed to reveal the
effects of repulsion and attraction on liquids and solids, and the phase
transitions between them.

Today, it has been established that 2D melting scenarios depend on _repulsion
softness_ , providing the following microscopic scenarios of 2D melting Ryzhov
_et al._ (2017, 2020): (i) the Berezinskii-Kosterlitz-Thouless-Halperin-
Nelson-Young (BKTHNY) theory, according to which melting occurs via two
continuous transitions with an intermediary hexatic phase with quasi-long-
range orientational order and short-range translational order Berezinsky
(1971); Kosterlitz and Thouless (1973); Halperin and Nelson (1978); Nelson and
Halperin (1979); Young (1979); (ii) melting through a first order phase
transition Chui (1983); Ryzhov (1991); (iii) two-step melting, including a
continuous (Berezinskii-Kosterlitz-Thouless, BKT) crystal-hexatic phase
transition and a first-order phase transition between the hexatic phase and
isotropic liquid Bernard and Krauth (2011); Engel  _et al._ (2013); Qi,
Gantapara, and Dijkstra (2014); Kapfer and Krauth (2015); Thorneywork  _et
al._ (2017). The second and third scenarios are inherent in hard-sphere-like
systems, whereas the first one was observed for soft repulsion between
particles Dudalov  _et al._ (2014a, b); Ryzhov and Tareyeva (2002). It has
been established unambiguously that the softness of repulsion affects melting
scenarios, thermodynamics and excitation spectra in monolayer systems Tsiok
_et al._ (2015); Pomirchi, Ryzhov, and Tareyeva (2002); Brazhkin  _et al._
(2012); Kryuchkov, Khrapak, and Yurchenko (2017); Khrapak, Kryuchkov, and
Yurchenko (2018); Kryuchkov  _et al._ (2018, 2019); Kryuchkov, Brazhkin, and
Yurchenko (2019); Yakovlev  _et al._ (2020); Kryuchkov  _et al._ (2020);
Khrapak  _et al._ (2021); Fomin  _et al._ (2018); Tsiok  _et al._ (2020,
2021). However, to the best of our knowledge, the role of _attraction_ in the
melting scenario of monolayer systems remains unstudied systematically.

LJ interactions were among the first systems studied to understand the role of
attraction in melting. However, a lot of reported results on the critical
point and melting scenario for 2D LJ crystals still remain questionable. For
instance, to provide critical temperature depending on the truncation radius,
numerical simulations of a vapor-liquid curve performed in the Gibbs ensemble
were reported in Ref. Smit and Frenkel (1991). For critical temperature and
density, the authors obtained $T_{c}=0.515\pm 0.002$ and $\rho_{c}=0.355\pm
0.003$, respectively, for the full potential; and $T_{c}=0.459\pm 0.001$ and
$\rho_{c}=0.35\pm 0.01$ for the truncated and shifted potential at
$2.5\sigma$. Contradictory melting scenarios of the triangular crystal were
reported in the early works Frenkel and McTague (1979); Toxvaerd (1978);
Abraham (1980); Phillips, Bruch, and Murphy (1981); Bakker, Bruin, and
Hilhorst (1984); Strandburg, Zollweg, and Chester (1984), including two
continuous transitions with an intermediate hexatic phase according to the
BKTHNY theory Frenkel and McTague (1979) and a first-order transition Toxvaerd
(1978); Abraham (1980); Phillips, Bruch, and Murphy (1981); Bakker, Bruin, and
Hilhorst (1984); Strandburg, Zollweg, and Chester (1984).

Thanks to growth of computing capabilities, simulations of large systems
($\gtrsim 10^{5}$ particles) have recently provided new results on 2D melting
of LJ crystals and related systems. Simulations of the systems followed by
analysis of their equation of state and long-range asymptotics of the
translational correlation function (that accurately provides the stability
limit of the crystal) allowed identifying the melting scenarios unambiguously.
For instance, a change in the melting scenario was reported in Ref. Hajibabaei
and Kim (2019), where the authors studied 2D systems of particles interacting
via generalized LJ potential with different repulsive branches ($\propto
1/r^{12}$ and $\propto 1/r^{64}$). The scenario was found to occur through
first-order phase transitions at low temperatures and via two continuous BKT
transitions (according to the BKTHNY theory) at high. An LJ system at high
temperatures is typically assumed to be close to soft repulsive disks
$1/r^{12}$, but such extrapolation to a melting scenario contradicts the
results of Ref. Kapfer and Krauth (2015), wherein the soft disks $1/r^{n}$
with $n>6$ were shown to melt according to the third melting scenario. The
Mayer-Wood loop, inherent in first-order transition, is assumed to disappear
at high temperatures with an increase in the system size. However, the
explanation of the effect by finite-size scaling seems unconvincing: With an
increase in the system size, the loop should flatten and ultimately approach
the plateau Engel  _et al._ (2013); Alonso and Fernández (1999).

A first-order hexatic-liquid transition and a continuous crystal-hexatic BKT
transition at high temperatures and one first-order crystal-liquid transition
at low temperatures were identified in Ref. Li and Ciamarra (2020a) for 2D LJ
particles. The results of numerical simulation of a 2D LJ system and
attractive polygons (squares, pentagons and hexagons) of the same authors in
Ref. Li and Ciamarra (2020b) revealed the role of attraction in a melting
scenario. Thus, at low temperatures where the role of attraction is dominant
all systems melt via first-order transition due to suppression of the hexatic
phase. At high temperatures LJ disks melt in accordance with the third
scenario, the same as soft disks Kapfer and Krauth (2015); Tsiok, Fomin, and
Ryzhov (2018); Gaiduk  _et al._ (2019), whereas hexagons and squares do
according to the BKTHNY theory with participation of the hexatic and tetratic
phases, respectively. The melting scenario of pentagons does not change with
an increase in temperature and is a first-order transition.

LJ crystals compared with the Morse system in Ref. Toledano  _et al._ (2021)
were shown to melt via the third scenario at low temperatures and two
continuous BKT transitions at high. This agrees with Ref. Hajibabaei and Kim
(2019), but contradicts Refs. Kapfer and Krauth (2015); Li and Ciamarra
(2020a). The BKTHNY scenario at high temperatures was questioned because of
seeming disappearance of the Mayer-Wood loop (an analog of the Van-der-Waals
loop in the 3D case). For soft Morse interactions, the third melting scenario
was observed for all temperatures considered in Ref. Toledano  _et al._
(2021), whereas the authors expected to observe the BKTHNY scenario at higher
temperatures. However, at some parameters of potential softness, two
continuous transitions were observed already at low temperatures in the
presence of long-range attraction.

The role of attraction can be tested experimentally in colloidal systems,
known for a long time as model systems demonstrating a wide range of
“molecular-like” phenomena Fernandez-Nieves and Puertas (2016); Ivlev  _et
al._ (2012); Löwen (1994); Li, Zhou, and Han (2016); Stradner and
Schurtenberger (2020), in particular, crystallization and melting Alsayed
(2005); Taffs  _et al._ (2013); Zahn, Lenke, and Maret (1999); Zahn and Maret
(2000); Sprakel  _et al._ (2017); Paloli  _et al._ (2013); Wang  _et al._
(2012); Yakovlev  _et al._ (2019); Kryuchkov  _et al._ (2021), solid-solid
phase transitions Royall  _et al._ (2006); Yethiraj and van Blaaderen (2003);
Peng  _et al._ (2014), condensation and critical phenomena Royall, Aarts, and
Tanaka (2007); Zhang  _et al._ (2013); Li  _et al._ (2019) caused by
relatively _short-range_ attractive depletion forcesLekkerkerker and Tuinier
(2011). These collective phenomena can be visualized _in real-time_ with
spatial resolution of individual particles. _Long-range_ dipolar attraction
$\propto 1/r^{3}$ in colloidal systems can be induced and controlled _in situ_
with in-plane rotating magnetic Martin and Snezhko (2013); Byrom and Biswal
(2013); Du  _et al._ (2013); Di Du  _et al._ (2017); Hilou  _et al._ (2018) or
electric Snoswell  _et al._ (2006); Elsner  _et al._ (2009); Juarez  _et al._
(2011); Edwards and Bevan (2014); Juárez, Feicht, and Bevan (2012); Yakovlev
_et al._ (2017); Komarov, Yarkov, and Yurchenko (2019); Yakovlev  _et al._
(2021) fields. Using conically-rotating magnetic or electric fields with magic
angles, _Van-der-Waals-like attraction_ can be created $\propto 1/r^{6}$ with
“magic” fields Müller  _et al._ (2014); Osterman  _et al._ (2009); Pham  _et
al._ (2017); Komarov, Kryuchkov, and Yurchenko (2018). Recently, tunable
interactions have been achieved by using spatial hodographs of an external
electric or magnetic field Komarov and Yurchenko (2020), by engineering the
internal structureKomarov, Mantsevich, and Yurchenko (2021) and
geometryKomarov and Yurchenko (2021) of colloidal particles.

In the present paper, we studied the role of the attraction range in the phase
behavior and melting scenarios of a 2D system. Using MD simulations, we
considered a generalized Lennard-Jones system with a fixed repulsion branch
and different attraction, from long-range dipolar attraction to short-range
sticky-spheres-like. A decrease in the attraction range is found to suppress
the temperature of the gas-liquid (liquid-liquid in the case of colloid and
protein systems) critical point. We found that attraction does not affect
melting at high temperatures, occurring through the cascade of solid-hexatic
(BKT) and hexatic-liquid (first-order) phase transitions. Surprisingly, at
moderate temperatures we observed _two triple points_ , hexatic-liquid-gas and
crystal-hexatic-gas. The crystal-hexatic-gas triple point for the case of
isotropic dipolar attraction has been discovered to be _below_ the hexatic-
liquid-gas triple point.

##  II System and methods

###  II.1 MD details

We studied a system of particles interacting via the generalized Lennard-Jones
potential (LJn-m):

$U_{nm}(r)=\frac{\epsilon}{n-m}\left[m\left(\frac{\sigma}{r}\right)^{n}-n\left(\frac{\sigma}{r}\right)^{m}\right],$
(1)

where $n$ and $m$ are the indices of repulsive and attractive branches,
respectively, and $\sigma$ and $\epsilon$ are the characteristic length of
interaction and the depth of the potential well. The potential (1) has minimum
$-\epsilon$ at $r/\sigma=1$. In what follows, we normalized the distances and
energies to $\sigma$ and $\epsilon$, respectively, and considered the
particles of equal mass $m=1$.

We considered cases with fixed repulsion branch $n=12$ and different
attraction, from long-range isotropic dipolar attraction with $m=3$ to short-
range almost sticky-spheres-like attraction with $m=11$. Examples of the
potentials with $m=3$, $6$, $9$, and $11$ are given in Fig. 1. A case similar
to usual LJ interaction corresponds to $n=12$ and $m=6$.

Figure 1: Generalised LJ potentials with different attraction indices: The solid lines correspond to fixed repulsion with $n=12$ and different indices $m=3$, $6$, $9$, and $11$, illustrating the transition from long- to short-range attractive interactions. The values of cutoff radii $r_{c}$ (used for calculations of the equation of state) decrease as shown for each $m$.  Table 1: Simulation parameters for potentials LJ12-m used for calculations in the phase identification method, where $a=\rho^{-1/2}$ LJn-m | $\rho$ | $r_{c}$ | $T_{start}$ | $T_{stop}$ | $T_{step}$
---|---|---|---|---|---
LJ12-3 | 0.3 | 7.5$a$ | 0.7 | 3.05 | 0.06
LJ12-4 | 0.4 | 7.5$a$ | 0.1 | 2.0 | 0.03
LJ12-5 | 0.4 | 7.5$a$ | 0.1 | 1.5 | 0.02
LJ12-6 | 0.4 | 7.5$a$ | 0.1 | 1.5 | 0.02
LJ12-7 | 0.4 | 7.5$a$ | 0.454 | 0.503 | 0.003
LJ12-8 | 0.4 | 7.5$a$ | 0.2 | 0.4 | 0.003
LJ12-9 | 0.4 | 7.5$a$ | 0.14 | 0.3 | 0.002
LJ12-10 | 0.4 | 5.0$a$ | 0.09 | 0.151 | 0.001
LJ12-11 | 0.5 | 5.0$a$ | 0.025 | 0.076 | 0.001

The molecular dynamics (MD) simulations were performed for $N=2\times 10^{4}$
particles using LAMMPS in the canonical ($NVT$) ensemble in a wide range of
densities and temperatures, from $T=0.3$ to $T=50.0$, using $10^{8}$ steps
with time step $dt=10^{-3}$ (in dimensionless units of time normalized to
$\sqrt{m\sigma^{2}/\epsilon}$). The value of potential energy within cutoff
radius $r_{c}$ was $|U(r_{c})|<5\times 10^{-4}$: it was $U(r_{c})=-4\times
10^{-4}$ for $m=3$, $U(r_{c})=-3.4\times 10^{-4}$ for $m=6$,
$U(r_{c})=-2\times 10^{-4}$ for $m=9$, $U(r_{c})=-3.2\times 10^{-4}$ for
$m=11$. To determine accurately the boundaries of two-phase areas, crystal
melting was examined (at some points) for large systems with $N=256^{2}$ and
$512^{2}$ particles. The phase diagram was obtained using the dependence of
pressure on density (the equation of state) along the isotherms, analysis of
the radial distribution functions, orientational and translational order
parameters, and the corresponding correlation functions.

###  II.2 The method of phase identification

The condensate-gas binodals were obtained using the method of phase
identification, proposed in Ref. Ovcharov  _et al._ (2017). The method allows
obtaining a binodal in the coordinates density-temperature using analysis of
Voronoi cells in the system and, in particular, in terms of the following
order parameter Ovcharov  _et al._ (2017); Kryuchkov  _et al._ (2021):

$\begin{split}&\lambda_{i}=\frac{1}{N_{ni}+1}\left(\sigma_{i}+\sum_{j=1}^{N_{ni}}{\sigma_{j}}\right),\\\
&\sigma_{i}=\frac{1}{a_{i}N_{ni}}\sqrt{\sum_{j<k}^{N_{ni}}{(r_{ij}-r_{ik})^{2}}/2},\quad
r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|,\end{split}$ (2)

where $\mathbf{r}_{i}$ is the radius-vector of the $i$-th particle, $N_{ni}$
is the number of its neighboring cells, $a_{i}=\sqrt{S_{i}/\pi}$, $S$ is the
area of a corresponding Voronoi cell. The scalar $\lambda^{2}$-field has small
values for condensed (liquid and solid) phases, and becomes large at the
interface gas-condensate and in gaseous state. Then, analysis of corresponding
statistics of Voronoi cells allows obtaining densities $n_{c}$ and $n_{g}$ of
condensate and gas, respectively (see Ref. Ovcharov  _et al._ (2017) for
details).

Figure 2: The equation of state for the 2D LJ system: (a) The isotherms of the
system with $m=6$ below ($T=0.3$) and above ($T=0.45$) the triple point and
above the critical point at $T=1.0$. (b) Isotherm $T=0.45$ gas-liquid and
isotherm $T=0.3$ gas-crystal at low densities. (c) Isotherm $T=0.45$ gas-
liquid and isotherm $T=0.3$ gas-crystal at high densities.

To apply the method, additional simulations of a small system with
$N=3.6\times 10^{3}$ particles under the same $NVT$ conditions (as we
explained above), from $T_{start}$ to $T_{stop}$, with step $T_{step}$ were
conducted where each step included $6\times 10^{5}$ time-steps with
dimensionless time-step $\Delta t=5\times 10^{-4}$ (see the parameters in
Table 1). At each temperature, we performed phase identification using
$10^{2}$ simulation frames.

Table 2: The temperatures and densities in the critical and triple points obtained using the phase identification method LJn-m | $T_{CP}$ | $\rho_{CP}$ | $T_{TP}$ | $\rho_{TP}$
---|---|---|---|---
LJ12-3 | 1.442 | 0.366 | 0.577 | 1.095
LJ12-4 | 0.748 | 0.453 | 0.455 | 1.042
LJ12-5 | 0.606 | 0.466 | 0.408 | 1.017
LJ12-6 | 0.510 | 0.488 | 0.390 | 1.012
LJ12-7 | 0.480 | 0.486 | 0.399 | 1.015
LJ12-8 | 0.436 | 0.509 | 0.390 | 1.016
LJ12-9 | 0.429 | 0.516 | 0.393 | 1.029
LJ12-10 | 0.410 | 0.544 | 0.386 | 0.996
LJ12-11 | 0.4053 | 0.550 | 0.3844 | 1.0478

This method of phase identification becomes unsuitable near the critical
point, where the condensed and gaseous phases cannot be distinguished.
However, one can obtain the values of critical density $n_{\mathrm{CP}}$ and
temperature $T_{\mathrm{CP}}$ using the following approximation of the binodal
(the results are provided in Table 2):

$n_{c}-n_{g}\simeq A\tau^{\beta},\quad n_{c}+n_{g}\simeq
a\tau+2n_{\mathrm{CP}},$ (3)

where $\tau=T_{\mathrm{CP}}-T-$ is reduced temperature, $\beta$ is a critical
index, $A$ and $a$ are the free parameters of the model. Critical index
$\beta$ is related to the universality class of the system and is determined
by the range of attractionLuijten and Blöte (2002): In 2D systems with dipolar
isotropic attraction, $\beta=1/2$ (mean-field critical behavior), whereas in
the case of shorter-range attraction, at $m>3$, $\beta=1/8$ corresponds to 2D
Ising behaviorLuijten and Blöte (2002).

###  II.3 Order parameters

For analysis of the phase transition scenarios, we used orientational and
translational correlation functions. The orientational correlation function
(OCF) of the global orientational order parameter is calculated as:

$G_{6}(r)=\frac{\left<\Psi_{6}(\mathbf{r})\Psi_{6}^{*}(0)\right>}{g(r)},\quad\Psi_{6}=\frac{1}{N}\left<\left|\sum_{i}\psi_{6}(\mathbf{r}_{i})\right|\right>,$
(4)

where the averaging for $G_{6}$ is performed over all particles in the system,
$\psi_{6}=1/n(i)\sum_{j}{e^{6i\theta_{ij}}}$, $\theta_{ij}$ is the angle of
the vector between particles $i$ and $j$ with respect to the reference axis,
and the sum over $j$ is counting the $n(i)$ nearest-neighbors of $j$, obtained
from the Voronoi construction,
$g(r)=\langle\delta(\mathbf{r}_{i})\delta(\mathbf{r}_{j})\rangle$ is an
isotropic pair distribution function (here, $\mathbf{r}_{i}$ is the position
vector of particle $i$, and $r=|\mathbf{r}_{i}-\mathbf{r}_{j}|$). In the
hexatic phase, $G_{6}(r)$ behaves at large distances like $G_{6}(r)\propto
r^{-\eta_{6}}$ with $\eta_{6}\leq 1/4$ Halperin and Nelson (1978); Nelson and
Halperin (1979).

The translational correlation function (TCF) is

$G_{t}(r)=\frac{\langle\exp(i\mathbf{a}(\mathbf{r}_{i}-\mathbf{r}_{j}))\rangle}{g(r)},$
(5)

where $\mathbf{a}$ is the reciprocal-lattice vector of the first shell of the
crystal lattice. In the solid phase, $G_{t}(r)$ behaves at large distances
like $G_{t}(r)\propto r^{-\eta_{T}}$ with $\eta_{T}\leq{1}/{3}$ Halperin and
Nelson (1978); Nelson and Halperin (1979). In the hexatic phase and isotropic
liquid, $G_{t}$ decays exponentially.

Figure 3: The isotherms of LJ12-3 system in $P-V$ coordinates: The gas-liquid
transition boundary is obtained using Maxwell’s construction at different
temperatures and is shown with the black dashdotted line.

##  III Results and discussion

###  III.1 The equations of state and phase diagrams

The equations of state (isotherms) for all considered systems behave in a
similar manner and have peculiarities that, depending on the temperature and
density, attributed to gas-liquid transition or crystal melting, respectively.
This is shown in Fig. 2(a): Here, at isotherm $T=0.45$, one can see a wide
loop peculiar to gas-liquid transition and a narrow loop related to crystal
melting, whereas only one loop related to crystal melting is seen above the
gas-liquid critical point at $T=1.0$. Below the triple point, at $T=0.3$,
there is only one wide loop corresponding to gas-crystal transition.

Figure 4: The effect of the attraction range on the liquid-gas coexistence
area on the phase diagram: (a) The binodals condensate-gas for different
LJ12-m potentials; the circles are the points of the binodal and median
(obtained using the phase identification method), the diamonds are points
obtained from the equation of state, the grey lines are the approximations of
the binodal with Eq. (3), the stars denote the critical points from Table 2.
(b) The dependencies of triple and critical temperatures on attraction index
$m$ for LJ12-m interaction, the ratio $T_{CP}/T_{TP}$ is shown in the inset.

The binodals gas-liquid were calculated using Maxwell’s construction on the
isotherms in the $P-V$ coordinates (an example is shown in Fig. 3 for the
system with $m=3$), as well as with the phase identification method. The
results for the binodal condensate-gas obtained from the phase identification
method and the equation of state are provided in Fig. 4(a). Here, the colored
circles denote the densities of gas, condensate, and their average for each
potential we considered. The grey solid lines show the regions of the data
wherein we used the approximation (3) to obtain data about the critical point.
The grey dashed lines show the extrapolation of the phase diagram to the
critical points gas-liquid, depicted with colored stars (the corresponding
densities and temperatures are given in Table 2).

We observed that the drop in the attraction range reduces the critical
temperature, as well as the ratio between temperatures of the critical and
triple points, as shown in Fig. 4(b) and the corresponding inset (based on the
data in Table 2). With an increase in $m$ (short-range attraction), the two-
phase area becomes narrower towards lower densities and the ratio between the
critical and triple temperatures becomes closer to unity, as shown in Fig.
4(b). For LJ interaction ($m=6$), the critical temperature we obtained is
$T_{c}=(0.51\ldots 0.52)$ (depending on the method of evaluation), which
agrees well with the previous results $T_{c}=0.515\pm 0.002$ for the LJ
potential reported in Ref. Smit and Frenkel (1991).

At high temperatures, the effect of attraction should vanish, and only the
repulsive branch of potential (1) should play a more significant role. Indeed,
analyzing the melting scenarios of triangular crystals at high temperatures
(and high densities), we observed that LJ12-m systems demonstrated the same
melting scenarios as soft disks repulsive potential $1/r^{12}$ (the third
melting scenario Kapfer and Krauth (2015); Li and Ciamarra (2020a); Gaiduk
_et al._ (2019)).

Figure 5: Evolution of correlation functions on an isotherm: (a) The isotherm
$T=3.0$ for LJ system consisting of $N=512^{2}$, obtained by averaging over 20
independent replicas. (b) and (c) The OCF and TCF for the same system at
different densities. Figure 6: The third melting scenario of the triangular
crystals at high temperatures and densities: The solid lines are the
boundaries of the Mayer-Wood loop corresponding to the two-phase hexatic-
liquid area. The dashed lines are the boundaries of BKT transition from the
crystal to hexatic phase. The results are provided for the cases $m=3$, $6$,
$9$, and $11$, as labeled in the figure.

The Mayer-Wood loop was observed in all considered systems at high
temperatures, proving first-order transition. However, to determine an
appropriate melting scenario, we compared the stability limits of the crystal
and the hexatic phases with the boundaries of the Mayer-Wood loop. The
stability limits of the hexatic and crystal phases were derived from
asymptotic behavior of the correlation functions of the orientational and
translational order parameters, respectively. An example of evolution of the
correlation functions on the isotherm of the LJ system at $T=3.0$ is shown in
Fig. 5.

In Fig. 5(a), the equation of state is shown at $T=3.0$ (obtained by averaging
over 20 independent replicas) for the LJ system consisting of $N=512^{2}$
particles. One can see the clearly defined Mayer-Wood loop, inherent in first-
order transitions. Analyzing the correlation functions in Figs. 5(b,c), we
conclude that the loop corresponds to the two-phase hexatic-liquid area,
whereas the melting scenario is the third, with a continuous BKT transition
from crystal to hexatic and a first-order transition from hexatic to liquid.
Here, the hexatic phase becomes unstable at $\rho_{hex-liq}=1.014$ (Fig.
5(b)), inside the loop, whereas the crystal becomes unstable at $\rho_{sol-
hex}=1.037$ (Fig. 5(c)) getting far out of the loop. The obtained results
agree well with the melting scenario for repulsive disks $1/r^{12}$ (reported
in Refs. Kapfer and Krauth (2015); Li and Ciamarra (2020a); Gaiduk  _et al._
(2019)), but do not exhibit two continuous BKT transitions in the LJ system at
high temperatures Engel  _et al._ (2013); Khali, Chakraborty, and Chaudhuri
(2021). The obtained results agree well with the melting scenario for
repulsive disks $1/r^{12}$ (reported in Refs. Kapfer and Krauth (2015); Li and
Ciamarra (2020a); Gaiduk  _et al._ (2019)), but do not exhibit two continuous
BKT transitions in the LJ system at high temperatures Engel  _et al._ (2013);
Khali, Chakraborty, and Chaudhuri (2021).

All triangular crystals, we studied, were found to melt at high temperatures
via the third scenario, irrespective to the attraction index $m$. The melting
lines in this region of the phase diagram are shifted to low densities with an
increase in $m$ and behave similarly to a system of soft disks $1/r^{12}$, as
shown in Fig. 6.

However, with a decrease in melting temperature, the balance of the effects
provided by the repulsive and attractive branches is changed. As a result, we
observed that, at low temperatures, the melting scenario _changes_ from the
third one to a first-order transition (without inclusion of the hexatic phase)
in the systems with $m=6$, $9$, and $11$. The change in the melting scenario
is shown for two isotherms of the LJ12-9 system in Fig. 7.

The temperature of the observed change in the melting scenario decreases with
$m$ (to more soft and long-range attraction). The attraction range affects
essentially the region of gas-liquid transition, as well as temperatures
$T_{c}$ and $T_{m}$, at which the hexatic phase appears.

Figure 7: The change in the melting scenario for the case $m=9$: The melting
scenario changes from a first-order transition ($T=0.6$) to the third scenario
($T=0.7$). The arrows show the stability limits of the hexatic phase, $G_{6}$
are obtained using the OCF, and the stability limits of the crystalline phase,
$G_{t}$ are obtained with TCF, which become caught by the loop with a decrease
in temperature.

###  III.2 The results for isotropic dipolar attraction

Figure 8: The cascades of phase transitions: Isotherm $T=0.6>T_{t}$
illustrates the cascade of a first-order gas-liquid transition, first-order
liquid-hexatic transition, and a continuous BKT crystal-hexatic transition. At
$T=0.57<T_{t}$, only a first-order gas-hexatic transition and a continuous BKT
crystal-hexatic transition are seen. Detailed behavior of the isotherms at
high densities are shown in the inset. At negative pressures and
$T=0.57<T_{t}$, one can see a loop related to the first-order transition gas-
hexatic (in the metastable area), while a continuous BKT crystal-hexatic
transition still occurs at positive pressure. The arrows show the limits of
stability of the hexatic phase $G_{6}$ and of the crystal $G_{t}$, obtained
using OCF and TCF analysis, respectively.

In the case of LJ12-3 potential, the crystal was found to melt at low
temperatures according to the third scenario. However, the hexatic phase was
discovered to be preserved up to triple point temperature $T_{t}$, where it
coexists with both liquid and gaseous phases, as given in Fig. 8.

This unexpected result is explained using the Gibbs rule: Due to the reduced
spatial dimension and since melting represents a two-stage process, we obtain
a couple of points where three phases may exist simultaneously: hexatic-
liquid-gas (at temperature $T_{t1}$) and crystal-hexatic-gas (at lower
temperature $T_{t2}<T_{t1}$). This situation is caused by appearance of an
additional (hexatic) phase, inherent in 2D systems, and stands in contrast to
simple bulk systems, where only one triple point gas-liquid-solid may exist.
In the interval $T_{t}2<T_{t}1$, a continuous BKT transition from the crystal
to hexatic phase and a first-order transition from hexatic to gas occur up to
the crystal-gas sublimation line, at which a first-order transition from
crystal to gas is observed on the isotherms. Simulation of first-order gas-
solid transitions at temperatures below $T_{t2}$ for finite systems inevitably
leads to the appearance of metastable states. This is manifested by the Mayer-
Wood loops, strongly extended due to a large difference of densities on the
isotherms. Moreover, these loops turn out to have a complex structure, as
given in Fig. 9.

Figure 9: The isotherms of the LJ system at temperatures $T<T_{t2}$ in the
high-density region.

At high densities in the vicinity of the crystalline phase, additional
modulations of the loops may be interpreted as two-stage melting of the
metastable crystal appearing in the equation of state. In this case, phase
identification was performed using analysis of the radial distribution
functions, behavior of the orientational and translational order parameters,
and their correlation functions. We observed a transition of the metastable
crystal to the hexatic phase through a continuous BKT transition with a
subsequent first-order transition from hexatic to gas. With a decrease in
temperature, the region of the hexatic phase narrows on the phase diagram, and
disappears completely at $T=0.45$, which leads to one first-order crystal-gas
transition, as shown in Fig. 10.

Figure 10: The change of the transition scenario from the third scenario
(isotherm $T=0.5$) to a first-order transition (isotherm $T=0.45$) in the area
of negative pressures at $T<T_{t2}$, which is related to suppression of the
hexatic phase at $T=0.45$. Figure 11: The phase diagram of the LJ12-3 system:
(a) The phase diagram in the ($P-T$) coordinates. (b) The phase diagram in the
($\rho-T$) coordinates. $T_{CP}$ is a gas-liquid critical point, $T_{t1}$ is a
hexatic-liquid-gas triple point, and $T_{t2}$ is a crystal-hexatic-gas triple
point. Gas, liquid, hexatic phase and crystal are denoted as G, L, Hex, and S,
respectively. The regions of their coexistence are marked as G+L, G+Hex, and
G+S. The solid and dashed lines correspond to a first-order transition and to
continuous BKT crystal-hexatic phase transition, respectively. The two-stage
crystal-gas transition in the metastable area of the hexatic phase and crystal
at $T<T_{t2}$ are denoted with open symbols. (c) The enlarged region of the
phase diagram in the area of triple points $T_{t1}$ and $T_{t2}$.

For the case of isotropic dipolar attraction $m=3$, the phase diagrams in
different coordinates are presented in Figs. 11(a,b). Here, the lines of phase
coexistence, critical, and triple points are shown in Fig. 11(a). In the
density-temperature plane, a two-stage crystal-gas transition occurs in the
metastable area of the hexatic phase and crystal at $T<T_{t2}$, see Fig.
11(b). In Fig. 11(c) an enlargement of the phase diagram is represented the
behavior of different domain boundaries on the phase diagram near triple
points $T_{t1}$ and $T_{t2}$. We note in conclusion that, with an increase in
the number of particles in the system, the loops corresponding to gas-crystal
transition should flatten. This leads to the disappearance of the metastable
area of the crystal and to the first-order gas-crystal transition, without the
modulation of the equation of state mentioned above.

##  IV Conclusions

In this work, we studied the evolution of the phase diagrams and melting
scenarios of two dimensional systems of particles interacting via generalized
LJ potential with a different attraction range, whereas the repulsion branch
was fixed. Gas-liquid transition was studied using analysis of the equation of
state and the phase identification method. The results obtained by the two
methods stand in good agreement indicating their consistency. The drop in the
attraction range reduces the gas-liquid coexistence region and temperatures of
the critical and triple points.

Melting at high temperatures (and high densities) is found to behave like a
system of soft disks $1/r^{12}$, via the third scenario. However, at low
melting temperatures, a change in the melting scenario was identified from the
third scenario to a first-order transition (without the hexatic phase) in the
systems with $m=6$, $9$, and $11$. The temperature of the change in the
scenarios is shifted toward lower temperatures with an increase in the
attraction range (corresponding to the decrease in $m$). Analysis of the case
$m=9$ (LJ12-9), very close to that studied in Ref. Li  _et al._ (2019),
demonstrated that for short-range attraction in complex solvents we should
witness the third melting scenario.

The largest region of gas-liquid coexistence was observed in the phase diagram
in the $\rho-T$ coordinates in the case of $m=3$. Due to soft long-range
isotropic dipolar attraction, we unexpectedly saw two triple points at
$T_{t1}$ and $T_{t2}$, corresponding to hexatic-liquid-gas and to crystal-
hexatic-gas equilibrium, respectively. On the isotherms in the interval
between two triple points $T_{t2}$ and $T_{t1}$, we observed a continuous BKT
transition from crystal to hexatic phase and a first-order transition from
hexatic to gas. At the crystal-gas sublimation line, a first-order transition
from crystal to gas has been identified.

The main conclusions related to the scenario of phase transitions in a system
with dipolar attraction can be tested in future experiments with colloidal
systems, wherein the tunable long-range dipolar attraction can be created with
rotating magnetic or electric fields. A nontrivial question here is related to
the role of three-body forces, inherent in atomic materials and tunable
colloids Komarov, Yarkov, and Yurchenko (2019); Yakovlev  _et al._ (2021). LJ-
like interactions can be created with magic hodographs of rotating electric or
magnetic fields Komarov and Yurchenko (2020). However, three-body interactions
in this case behave as $\propto 1/r^{6}$, with the same asymptotics as
pairwise potential. Due to this, the scenario of melting can change, and we
leave the corresponding study to future work.

###### Acknowledgements.

The study was supported by the Russian Science Foundation, Grant No.
19-12-00092 (MD simulations to obtain phase diagrams) and 20-12-00356
(analysis with phase identification method). This work was carried out using
computing resources of the federal collective usage center “Complex for
simulation and data processing for mega-science facilities” at NRC “Kurchatov
Institute”, http://ckp.nrcki.ru, and supercomputers at Joint Supercomputer
Center of the Russian Academy of Sciences (JSCC RAS). E.N.T., Y.D.F., and
E.A.G. performed MD simulations to obtain phase diagrams, N.A.D., P.A.L., and
N.P.K. obtained condensate-gas binodals with phase identification method,
E.N.T., V.N.R., E.E.T., N.P.K., and S.O.Y. analysed and discussed the results,
E.N.T., V.N.R., and S.O.Y. wrote the manuscript, V.N.R. and S.O.Y. conceived
and directed the study.

## References

  * Fernandez-Nieves and Puertas (2016) A. Fernandez-Nieves and A. M. Puertas,  _Fluids, colloids, and soft materials: an introduction to soft matter physics_ (Wiley, 2016).
  * Ryzhov  _et al._ (2017) V. N. Ryzhov, E. E. Tareyeva, Y. D. Fomin, and E. N. Tsiok, “Berezinskii – kosterlitz – thouless transition and two-dimensional melting,” Physics-Uspekhi 60, 857–885 (2017).
  * Ryzhov  _et al._ (2020) V. N. Ryzhov, E. E. Tareyeva, Y. D. Fomin, and E. N. Tsiok, “Complex phase diagrams of systems with isotropic potentials: results of computer simulations,” Physics-Uspekhi 63, 417–439 (2020).
  * Berezinsky (1971) V. L. Berezinsky, “Destruction of long range order in one-dimensional and two-dimensional systems having a continuous symmetry group. i. classical systems,” Sov. Phys. JETP 32, 493–500 (1971).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, “Ordering, metastability and phase transitions in two-dimensional systems,” Journal of Physics C: Solid State Physics 6, 1181–1203 (1973).
  * Halperin and Nelson (1978) B. I. Halperin and D. R. Nelson, “Theory of two-dimensional melting,” Physical Review Letters 41, 121–124 (1978).
  * Nelson and Halperin (1979) D. R. Nelson and B. I. Halperin, “Dislocation-mediated melting in two dimensions,” Physical Review B 19, 2457–2484 (1979).
  * Young (1979) A. P. Young, “Melting and the vector coulomb gas in two dimensions,” Physical Review B 19, 1855–1866 (1979).
  * Chui (1983) S. T. Chui, “Grain-boundary theory of melting in two dimensions,” Physical Review B 28, 178–194 (1983).
  * Ryzhov (1991) V. Ryzhov, “Dislocation-disclination melting of two-dimensional lattices,” JETP 73, 899 (1991).
  * Bernard and Krauth (2011) E. P. Bernard and W. Krauth, “Two-step melting in two dimensions: First-order liquid-hexatic transition,” Physical Review Letters 107, 155704 (2011).
  * Engel  _et al._ (2013) M. Engel, J. A. Anderson, S. C. Glotzer, M. Isobe, E. P. Bernard, and W. Krauth, “Hard-disk equation of state: First-order liquid-hexatic transition in two dimensions with three simulation methods,” Physical Review E 87, 042134 (2013).
  * Qi, Gantapara, and Dijkstra (2014) W. Qi, A. P. Gantapara, and M. Dijkstra, “Two-stage melting induced by dislocations and grain boundaries in monolayers of hard spheres,” Soft Matter 10, 5449 (2014).
  * Kapfer and Krauth (2015) S. C. Kapfer and W. Krauth, “Two-dimensional melting: From liquid-hexatic coexistence to continuous transitions,” Physical Review Letters 114, 035702 (2015).
  * Thorneywork  _et al._ (2017) A. L. Thorneywork, J. L. Abbott, D. G. Aarts, and R. P. Dullens, “Two-dimensional melting of colloidal hard spheres,” Physical Review Letters 118, 158001 (2017).
  * Dudalov  _et al._ (2014a) D. E. Dudalov, E. N. Tsiok, Y. D. Fomin, and V. N. Ryzhov, “Effect of a potential softness on the solid-liquid transition in a two-dimensional core-softened potential system,” The Journal of Chemical Physics 141 (2014a), 10.1063/1.4896825.
  * Dudalov  _et al._ (2014b) D. E. Dudalov, Y. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, “How dimensionality changes the anomalous behavior and melting scenario of a core-softened potential system?” Soft Matter 10, 4966–4976 (2014b).
  * Ryzhov and Tareyeva (2002) V. Ryzhov and E. Tareyeva, “Melting in two dimensions: first-order versus continuous transition,” Physica A-Statistical Mechanics and Its Applications 314, 396–404 (2002), Meeting on Horizons in Complex Systems, MESSINA, ITALY, DEC, 2001.
  * Tsiok  _et al._ (2015) E. N. Tsiok, D. E. Dudalov, Y. D. Fomin, and V. N. Ryzhov, “Random pinning changes the melting scenario of a two-dimensional core-softened potential system,” Physical Review E 92 (2015), 10.1103/PhysRevE.92.032110.
  * Pomirchi, Ryzhov, and Tareyeva (2002) L. Pomirchi, V. Ryzhov, and E. Tareyeva, “Melting of two-dimensional systems: Dependence of the type of transition on the radius of the potential,” Theoretical and Mathematical Physics 130, 101–110 (2002).
  * Brazhkin  _et al._ (2012) V. V. Brazhkin, A. G. Lyapin, V. N. Ryzhov, K. Trachenko, Y. D. Fomin, and E. N. Tsiok, “Where is the supercritical fluid on the phase diagram?” Physics-Uspekhi 55, 1061–1079 (2012).
  * Kryuchkov, Khrapak, and Yurchenko (2017) N. P. Kryuchkov, S. A. Khrapak, and S. O. Yurchenko, “Thermodynamics of two-dimensional yukawa systems across coupling regimes,” The Journal of Chemical Physics 146, 134702 (2017).
  * Khrapak, Kryuchkov, and Yurchenko (2018) S. A. Khrapak, N. P. Kryuchkov, and S. O. Yurchenko, “Thermodynamics and dynamics of two-dimensional systems with dipolelike repulsive interactions,” Physical Review E 97, 022616 (2018).
  * Kryuchkov  _et al._ (2018) N. P. Kryuchkov, E. V. Yakovlev, E. A. Gorbunov, L. Couedel, A. M. Lipaev, and S. O. Yurchenko, “Thermoacoustic instability in two-dimensional fluid complex plasmas,” Physical Review Letters 121, 075003 (2018).
  * Kryuchkov  _et al._ (2019) N. P. Kryuchkov, L. A. Mistryukova, V. V. Brazhkin, and S. O. Yurchenko, “Excitation spectra in fluids: How to analyze them properly,” Scientific Reports 9, 10483 (2019).
  * Kryuchkov, Brazhkin, and Yurchenko (2019) N. P. Kryuchkov, V. V. Brazhkin, and S. O. Yurchenko, “Anticrossing of longitudinal and transverse modes in simple fluids,” The Journal of Physical Chemistry Letters 10, 4470–4475 (2019).
  * Yakovlev  _et al._ (2020) E. V. Yakovlev, N. P. Kryuchkov, P. V. Ovcharov, A. V. Sapelkin, V. V. Brazhkin, and S. O. Yurchenko, “Direct experimental evidence of longitudinal and transverse mode hybridization and anticrossing in simple model fluids,” The Journal of Physical Chemistry Letters 11, 1370–1376 (2020).
  * Kryuchkov  _et al._ (2020) N. P. Kryuchkov, L. A. Mistryukova, A. V. Sapelkin, V. V. Brazhkin, and S. O. Yurchenko, “Universal effect of excitation dispersion on the heat capacity and gapped states in fluids,” Physical Review Letters 125, 125501 (2020).
  * Khrapak  _et al._ (2021) S. Khrapak, N. P. Kryuchkov, L. A. Mistryukova, and S. O. Yurchenko, “From soft- to hard-sphere fluids: Crossover evidenced by high-frequency elastic moduli,” Physical Review E 103, 052117 (2021).
  * Fomin  _et al._ (2018) Y. D. Fomin, E. A. Gaiduk, E. N. Tsiok, and V. N. Ryzhov, “The phase diagram and melting scenarios of two-dimensional Hertzian spheres,” Molecular Physics 116, 3258–3270 (2018).
  * Tsiok  _et al._ (2020) E. N. Tsiok, E. A. Gaiduk, Y. D. Fomin, and V. N. Ryzhov, “Melting scenarios of two-dimensional Hertzian spheres with a single triangular lattice,” Soft Matter 16, 3962–3972 (2020).
  * Tsiok  _et al._ (2021) E. N. Tsiok, Y. D. Fomin, E. A. Gaiduk, and V. N. Ryzhov, “Structural transition in two-dimensional Hertzian spheres in the presence of random pinning,” Physical Review E 103 (2021), 10.1103/PhysRevE.103.062612.
  * Smit and Frenkel (1991) B. Smit and D. Frenkel, “Vapor–liquid equilibria of the two-dimensional lennard-jones fluid(s),” The Journal of Chemical Physics 94, 5663–5668 (1991).
  * Frenkel and McTague (1979) D. Frenkel and J. P. McTague, “Evidence for an orientationally ordered two-dimensional fluid phase from molecular-dynamics calculations,” Physical Review Letters 42, 1632–1635 (1979).
  * Toxvaerd (1978) S. Toxvaerd, “Melting in a two-dimensional lennard-jones system,” The Journal of Chemical Physics 69, 4750–4752 (1978).
  * Abraham (1980) F. F. Abraham, “Melting in two dimensions is first order: An isothermal-isobaric monte carlo study,” Physical Review Letters 44, 463–466 (1980).
  * Phillips, Bruch, and Murphy (1981) J. M. Phillips, L. W. Bruch, and R. D. Murphy, “The two-dimensional lennard-jones system: Sublimation, vaporization, and melting,” The Journal of Chemical Physics 75, 5097–5109 (1981).
  * Bakker, Bruin, and Hilhorst (1984) A. F. Bakker, C. Bruin, and H. J. Hilhorst, “Orientational order at the two-dimensional melting transition,” Physical Review Letters 52, 449–452 (1984).
  * Strandburg, Zollweg, and Chester (1984) K. J. Strandburg, J. A. Zollweg, and G. V. Chester, “Bond-angular order in two-dimensional lennard-jones and hard-disk systems,” Physical Review B 30, 2755–2759 (1984).
  * Hajibabaei and Kim (2019) A. Hajibabaei and K. S. Kim, “First-order and continuous melting transitions in two-dimensional lennard-jones systems and repulsive disks,” Physical Review E 99, 022145 (2019).
  * Alonso and Fernández (1999) J. J. Alonso and J. F. Fernández, “van der waals loops and the melting transition in two dimensions,” Physical Review E 59, 2659–2663 (1999).
  * Li and Ciamarra (2020a) Y.-W. Li and M. P. Ciamarra, “Phase behavior of lennard-jones particles in two dimensions,” Physical Review E 102, 062101 (2020a).
  * Li and Ciamarra (2020b) Y.-W. Li and M. P. Ciamarra, “Attraction tames two-dimensional melting: From continuous to discontinuous transitions,” Physical Review Letters 124, 218002 (2020b).
  * Tsiok, Fomin, and Ryzhov (2018) E. N. Tsiok, Y. D. Fomin, and V. N. Ryzhov, “Random pinning elucidates the nature of melting transition in two-dimensional core-softened potential system,” Physica A-Statistical Mechanics and Its Applications 490, 819–827 (2018).
  * Gaiduk  _et al._ (2019) E. A. Gaiduk, Y. Fomin, E. N. Tsiok, and V. N. Ryzhov, “The influence of random pinning on the melting scenario of two-dimensional soft-disk systems,” Molecular Physics 117, 2910–2919 (2019).
  * Toledano  _et al._ (2021) Ó. Toledano, M. Pancorbo, J. E. Alvarellos, and Ó. Gálvez, “Melting in two-dimensional systems: Characterizing continuous and first-order transitions,” Physical Review B 103, 094107 (2021).
  * Ivlev  _et al._ (2012) A. Ivlev, H. Löwen, G. Morfill, and C. P. Royall,  _Complex plasmas and Colloidal dispersions: particle-resolved studies of classical liquids and solids (Series in soft condensed matter)_ (Singapore: Word Scientific, 2012).
  * Löwen (1994) H. Löwen, “Melting, freezing and colloidal suspensions,” Phys. Rep. 237, 249 – 324 (1994).
  * Li, Zhou, and Han (2016) B. Li, D. Zhou, and Y. Han, “Assembly and phase transitions of colloidal crystals,” Nature Reviews Materials 1, 15011 (2016).
  * Stradner and Schurtenberger (2020) A. Stradner and P. Schurtenberger, “Potential and limits of a colloid approach to protein solutions,” Soft Matter 16, 307–323 (2020).
  * Alsayed (2005) A. M. Alsayed, “Premelting at defects within bulk colloidal crystals,” Science 309, 1207–1210 (2005).
  * Taffs  _et al._ (2013) J. Taffs, S. R. Williams, H. Tanaka, and C. P. Royall, “Structure and kinetics in the freezing of nearly hard spheres,” Soft Matter 9, 297–305 (2013).
  * Zahn, Lenke, and Maret (1999) K. Zahn, R. Lenke, and G. Maret, “Two-stage melting of paramagnetic colloidal crystals in two dimensions,” Physical Review Letters 82, 2721–2724 (1999).
  * Zahn and Maret (2000) K. Zahn and G. Maret, “Dynamic criteria for melting in two dimensions,” Physical Review Letters 85, 3656–3659 (2000).
  * Sprakel  _et al._ (2017) J. Sprakel, A. Zaccone, F. Spaepen, P. Schall, and D. A. Weitz, “Direct observation of entropic stabilization of bcc crystals near melting,” Physical Review Letters 118, 088003 (2017).
  * Paloli  _et al._ (2013) D. Paloli, P. S. Mohanty, J. J. Crassous, E. Zaccarelli, and P. Schurtenberger, “Fluid-solid transitions in soft-repulsive colloids,” Soft Matter 9, 3000 (2013).
  * Wang  _et al._ (2012) Z. Wang, F. Wang, Y. Peng, Z. Zheng, and Y. Han, “Imaging the homogeneous nucleation during the melting of superheated colloidal crystals,” Science 338, 87–90 (2012).
  * Yakovlev  _et al._ (2019) E. V. Yakovlev, M. Chaudhuri, N. P. Kryuchkov, P. V. Ovcharov, A. V. Sapelkin, and S. O. Yurchenko, “Experimental validation of interpolation method for pair correlations in model crystals,” The Journal of Chemical Physics 151, 114502 (2019).
  * Kryuchkov  _et al._ (2021) N. P. Kryuchkov, N. A. Dmitryuk, W. Li, P. V. Ovcharov, Y. Han, A. V. Sapelkin, and S. O. Yurchenko, “Mean-field model of melting in superheated crystals based on a single experimentally measurable order parameter,” Scientific Reports 11, 17963 (2021).
  * Royall  _et al._ (2006) C. P. Royall, M. E. Leunissen, A.-P. Hynninen, M. Dijkstra, and A. van Blaaderen, “Re-entrant melting and freezing in a model system of charged colloids,” The Journal of Chemical Physics 124, 244706 (2006).
  * Yethiraj and van Blaaderen (2003) A. Yethiraj and A. van Blaaderen, “A colloidal model system with an interaction tunable from hard sphere to soft and dipolar,” Nature 421, 513–517 (2003).
  * Peng  _et al._ (2014) Y. Peng, F. Wang, Z. Wang, A. M. Alsayed, Z. Zhang, A. G. Yodh, and Y. Han, “Two-step nucleation mechanism in solid-solid phase transitions,” Nature Materials 14, 101–108 (2014).
  * Royall, Aarts, and Tanaka (2007) C. P. Royall, D. G. A. L. Aarts, and H. Tanaka, “Bridging length scales in colloidal liquids and interfaces from near-critical divergence to single particles,” Nature Physics 3, 636–640 (2007).
  * Zhang  _et al._ (2013) I. Zhang, C. P. Royall, M. A. Faers, and P. Bartlett, “Phase separation dynamics in colloid–polymer mixtures: the effect of interaction range,” Soft Matter 9, 2076 (2013).
  * Li  _et al._ (2019) B. Li, X. Xiao, S. Wang, W. Wen, and Z. Wang, “Real-space mapping of the two-dimensional phase diagrams in attractive colloidal systems,” Physical Review X 9 (2019), 10.1103/physrevx.9.031032.
  * Lekkerkerker and Tuinier (2011) H. N. Lekkerkerker and R. Tuinier,  _Colloids and the Depletion Interaction_ (Springer Netherlands, 2011).
  * Martin and Snezhko (2013) J. E. Martin and A. Snezhko, “Driving self-assembly and emergent dynamics in colloidal suspensions by time-dependent magnetic fields,” Reports on Progress in Physics 76, 126601 (2013).
  * Byrom and Biswal (2013) J. Byrom and S. L. Biswal, “Magnetic field directed assembly of two-dimensional fractal colloidal aggregates,” Soft Matter 9, 9167 (2013).
  * Du  _et al._ (2013) D. Du, D. Li, M. Thakur, and S. L. Biswal, “Generating an in situ tunable interaction potential for probing 2-d colloidal phase behavior,” Soft Matter 9, 6867 (2013).
  * Di Du  _et al._ (2017) Di Du, M. Doxastakis, E. Hilou, and S. L. Biswal, “Two-dimensional melting of colloids with long-range attractive interactions,” Soft Matter 13, 1548–1553 (2017).
  * Hilou  _et al._ (2018) E. Hilou, D. Du, S. Kuei, and S. L. Biswal, “Interfacial energetics of two-dimensional colloidal clusters generated with a tunable anharmonic interaction potential,” Physical Review Materials 2, 025602 (2018).
  * Snoswell  _et al._ (2006) D. R. E. Snoswell, C. L. Bower, P. Ivanov, M. J. Cryan, J. G. Rarity, and B. Vincent, “Dynamic control of lattice spacing within colloidal crystals,” New Journal of Physics 8, 267–267 (2006).
  * Elsner  _et al._ (2009) N. Elsner, C. P. Royall, B. Vincent, and D. R. E. Snoswell, “Simple models for two-dimensional tunable colloidal crystals in rotating ac electric fields,” The Journal of Chemical Physics 130, 154901 (2009).
  * Juarez  _et al._ (2011) J. J. Juarez, B. G. Liu, J.-Q. Cui, and M. A. Bevan, “kT-scale colloidal interactions in high-frequency inhomogeneous AC electric fields. II. concentrated ensembles,” Langmuir 27, 9219–9226 (2011).
  * Edwards and Bevan (2014) T. D. Edwards and M. A. Bevan, “Controlling colloidal particles with electric fields,” Langmuir 30, 10793–10803 (2014).
  * Juárez, Feicht, and Bevan (2012) J. J. Juárez, S. E. Feicht, and M. A. Bevan, “Electric field mediated assembly of three dimensional equilibrium colloidal crystals,” Soft Matter 8, 94–103 (2012).
  * Yakovlev  _et al._ (2017) E. V. Yakovlev, K. A. Komarov, K. I. Zaytsev, N. P. Kryuchkov, K. I. Koshelev, A. K. Zotov, D. A. Shelestov, V. L. Tolstoguzov, V. N. Kurlov, A. V. Ivlev, and S. O. Yurchenko, “Tunable two-dimensional assembly of colloidal particles in rotating electric fields,” Scientific Reports 7, 13727 (2017).
  * Komarov, Yarkov, and Yurchenko (2019) K. A. Komarov, A. V. Yarkov, and S. O. Yurchenko, “Diagrammatic method for tunable interactions in colloidal suspensions in rotating electric or magnetic fields,” The Journal of Chemical Physics 151, 244103 (2019).
  * Yakovlev  _et al._ (2021) E. V. Yakovlev, N. P. Kryuchkov, S. A. Korsakova, N. A. Dmitryuk, P. V. Ovcharov, M. M. Andronic, I. A. Rodionov, A. V. Sapelkin, and S. O. Yurchenko, “2d colloids in rotating electric fields: A laboratory of strong tunable three-body interactions,” Journal of Colloid and Interface Science (2021), 10.1016/j.jcis.2021.09.116.
  * Müller  _et al._ (2014) K. Müller, N. Osterman, D. Babič, C. N. Likos, J. Dobnikar, and A. Nikoubashman, “Pattern formation and coarse-graining in two-dimensional colloids driven by multiaxial magnetic fields,” Langmuir 30, 5088–5096 (2014).
  * Osterman  _et al._ (2009) N. Osterman, I. Poberaj, J. Dobnikar, D. Frenkel, P. Ziherl, and D. Babić, “Field-induced self-assembly of suspended colloidal membranes,” Physical Review Letters 103, 228301 (2009).
  * Pham  _et al._ (2017) A. T. Pham, Y. Zhuang, P. Detwiler, J. E. S. Socolar, P. Charbonneau, and B. B. Yellen, “Phase diagram and aggregation dynamics of a monolayer of paramagnetic colloids,” Physical Review E 95, 052607 (2017).
  * Komarov, Kryuchkov, and Yurchenko (2018) K. A. Komarov, N. P. Kryuchkov, and S. O. Yurchenko, “Tunable interactions between particles in conically rotating electric fields,” Soft Matter 14, 9657–9674 (2018).
  * Komarov and Yurchenko (2020) K. A. Komarov and S. O. Yurchenko, “Colloids in rotating electric and magnetic fields: designing tunable interactions with spatial field hodographs,” Soft Matter 16, 8155–8168 (2020).
  * Komarov, Mantsevich, and Yurchenko (2021) K. A. Komarov, V. N. Mantsevich, and S. O. Yurchenko, “Core–shell particles in rotating electric and magnetic fields: Designing tunable interactions via particle engineering,” The Journal of Chemical Physics 155, 084903 (2021).
  * Komarov and Yurchenko (2021) K. A. Komarov and S. O. Yurchenko, “Diagrammatics of tunable interactions in anisotropic colloids in rotating electric or magnetic fields: New kind of dipole-like interactions,” The Journal of Chemical Physics 155, 114107 (2021).
  * Ovcharov  _et al._ (2017) P. V. Ovcharov, N. P. Kryuchkov, K. I. Zaytsev, and S. O. Yurchenko, “Particle-resolved phase identification in two-dimensional condensable systems,” The Journal of Physical Chemistry C 121, 26860–26868 (2017).
  * Luijten and Blöte (2002) E. Luijten and H. W. J. Blöte, “Boundary between long-range and short-range critical behavior in systems with algebraic interactions,” Physical Review Letters 89, 025703 (2002).
  * Khali, Chakraborty, and Chaudhuri (2021) S. S. Khali, D. Chakraborty, and D. Chaudhuri, “Two-step melting of the weeks–chandler–anderson system in two dimensions,” Soft Matter 17, 3473–3485 (2021).


# Evidence of quasi equilibrium in pressure-gradient turbulent boundary layers

Victor Baxerres1  Ricardo Vinuesa2  Hassan Nagib1<EMAIL_ADDRESS>1ILLINOIS
TECH (IIT), Chicago, IL, 60616, USA 2FLOW, Engineering Mechanics, KTH Royal
Institute of Technology, Stockholm, Sweden

###### Abstract

Two sets of measurements utilizing hot-wire anemometry and oil film
interferometry for flat-plate turbulent boundary layers, exposed to various
controlled adverse and favorable pressure gradients, are used to evaluate
history effects of the imposed and varying freestream gradients. The results
are from the NDF wind tunnel at ILLINOIS TECH (IIT) and the MTL wind tunnel at
KTH, over the range $800<Re_{\tau}<22,000$ (where $Re_{\tau}$ is the friction
Reynolds number). The streamwise pressure-gradient parameter
$\beta\equiv(-\ell/\tau_{w})\cdot(\partial P_{e}/\partial x)$ varied between
$-2<\beta<7$, where $\ell$ is an outer length scale for boundary layers
equivalent to the half height of channel flow and the radius of pipe flow, and
is estimated for each boundary-layer profile. Extracting from each profile the
three parameters of the overlap region, following the recent work of Monkewitz
& Nagib (2023) that led to an overlap region of combined logarithmic and
linear parts, we find minimum history effects in the overlap region. Thus, the
overlap region in this range of pressure-gradient boundary layers appears to
be in “quasi equilibrium”.

###### keywords:

wall-bounded turbulence, boundary layers, pressure gradient, history effects

##  1 Introduction

A topic within the study of wall-bounded turbulence that has received
extensive attention for many decades, due to its significant implications on
flow representation and modeling, is the overlap region between the inner or
wall region and the outer region of the flow. According to the classical
literature, the mean velocity profile in this overlap region follows the well-
known logarithmic law (1), which contains $\kappa$, the von Kármán constant.

$\overline{U}^{+}_{x\rm OL}(y^{+}\gg 1~{}\&~{}Y\ll 1)=\frac{1}{\kappa}\ln
y^{+}+B.$ (1)

In this equation, $\overline{U}_{x}$ is the streamwise mean velocity, the
“${\rm+}$” superscript indicates that we are using the viscous scaling, and
the “${\rm OL}$” subscript indicates validity in overlap region. The
independent variable $y^{+}$ is the inner-scaled wall-normal coordinate,
$y^{+}=yu_{\tau}/\nu$, where $u_{\tau}$ is the friction velocity and $\nu$ is
the fluid kinematic viscosity. The variable $Y$ is the wall-normal coordinate
made dimensionless by the “outer scale of the flow” $\ell$. Since the work of
Nagib & Chauhan (2008) (CN08 from now on) the universality of the log law and
the von Kármán “coefficient” have been occasionally challenged or reaffirmed,
e.g., see references by George & Castillo (2006); Monkewitz et al. (2008);
Vinuesa et al. (2014); Luchini (2017). More recently, Monkewitz & Nagib (2023)
(MN23 from now on) shed additional light on this topic, challenging the
accumulated knowledge during the last century. The main point of MN23 is that
a more thorough matched asymptotic expansion (MAE) of wall-bounded flows
yields an overlap region containing a combination of logarithmic and linear
terms of the distance to the wall in (2). They arrived at this result by
considering in the inner asymptotic expansion a term proportional to the wall-
normal coordinate, $\mathcal{O}(Re_{\tau}^{-1})$. Here we will reference this
overlap region as “log+lin”, and use $Re_{\tau}=u_{\tau}\ell/\nu$, which is
the friction Reynolds number with $\ell$ a characteristic outer length. For
consistency with MN23, this $\ell$ can be thought of as the radius $R$ in
pipes, the semi-height in channels, $h$, or a length scale related to the 99%
boundary-layer thickness $\delta_{99}$ in boundary layers.

$\overline{U}^{+}_{x\rm OL}(y^{+}\gg 1~{}\&~{}Y\ll 1)=\kappa^{-1}\ln
y^{+}+S_{0}y^{+}/Re_{\tau}+B_{0}+B_{1}/Re_{\tau}.$ (2)

Earlier publications by Yajnik (1970); Afzal & Yajnik (1973); Lee & Moser
(2015); Luchini (2017) have considered a linear term in the overlap region of
channel and pipe flows. However, unlike MN23, none of them revealed the non-
universality of the Kármán constant, or the dependence of the overlap
coefficients $\kappa$, $S_{0}$ and $B_{0}$ on the pressure gradient along the
flow. In most literature on boundary layers, including for turbulent
conditions, the commonly used pressure-gradient parameter is
$\beta^{*}\equiv(-\delta^{*}/\tau_{w})\cdot(\partial P_{e}/\partial x)$, where
$\delta^{*}$ is the displacement thickness, and $P_{e}$ is the freestream
static pressure. In order to compare results from developing boundary layers
to fully developed flow in channels and pipes, at corresponding pressure
gradients, we will use the parameter
$\beta\equiv(-\ell/\tau_{w})\cdot(\partial P_{e}/\partial x)$. The outer
length scale of all boundary-layer profiles used here, $\ell$, is established
in the next section. Effects of pressure gradient on boundary layers have
often been considered to have “history effects” on the flow. In particular,
several studies have documented flow-history effects on the Reynolds stresses
and the wake region of the mean velocity profile (Harun et al., 2013; Bobke et
al., 2017), while other studies have also reported an effect of flow history
on the overlap region and its coefficients $\kappa$ and $B$ (Vishwanathan et
al., 2023).

Our objective here is to examine whether such history effects are present in
the overlap region, or a behavior more like quasi equilibrium exists when the
log+lin overlap is used to analyse well-documented boundary-layer data under
conditions of adverse pressure gradient (APG) or favorable pressure gradient
(FPG), including strong favorable pressure gradient (SFPG), and for reference,
zero pressure gradient (ZPG).

##  2 Boundary-layer databases

Two sets of boundary-layer data were selected and used because the pressure
gradients imposed were well controlled and documented, and independent wall-
shear stress was measured directly by oil film interferometry. One was carried
out in the NDF wind tunnel at ILLINOIS TECH (I.I.T., Chicago) (Nagib et al.
(2004)) and the other was obtained in the MTL wind tunnel of KTH (Stockholm)
(Sanmiguel Vila et al. (2020)). Both sets of data were measured by hot-wire
anemometry and the free stream velocity was measured independently using Pitot
probe measurements at reference points above the boundary layers developing on
a suspended flat plate. In the case of the MTL data we did not use the
”strongly increasing” pressure gradient data, and only the 30 m/s cases to
avoid any low $Re_{\tau}$ effects and insure fully turbulent conditions even
near the start of the boundary layer development. More cases were studied by
Sanmiguel Vila et al. (2020).

The velocity profiles of all the data used by us are shown in figure 1. Since
the work by Chauhan et al. (2009), a composite mean velocity profile has been
used to determine an outer scale of ZPG boundary layers instead of various
other measures of the boundary-layer thickness, including integral
thicknesses, e.g., Samie et al. (2018) and MN23. This outer scale $\delta$ was
systematically compared with the large data base used by Chauhan et al.
(2009), and we concluded that for ZPG, $\delta\approx 1.25\delta_{99}$. In
both the NDF and MTL data sets we used, we also found that approximate measure
of an outer scale to the edge of the boundary layers is valid under APG, FPG
and SFPG conditions in addition to ZPG. Therefore we utilized the outer scale
$\ell\equiv\Delta_{1.25}=1.25\delta_{99}$, and
$\beta\equiv(-\Delta_{1.25}/\tau_{w})\cdot(\partial P_{e}/\partial x)$. We
also define $Y\equiv y^{+}/Re_{\tau}$ to obtain the following from (2):

${U_{x}}^{+}_{\rm OL}(Y\ll 1)\sim\kappa^{-1}\ln Y+\kappa^{-1}\ln
Re_{\tau}+S_{0}Y+B_{0}+B_{1}/Re_{\tau}.$ (3)

Simplifying, we obtain:

${U_{x}}^{+}_{\rm OL}(Y\ll 1)=\kappa^{-1}\ln Y+\kappa^{-1}\ln
Re_{\tau}+S_{0}Y+B_{0}+{\rm H.O.T.}$ (4)

Figure 1: Mean velocity profiles from boundary layers of data sets used (Nagib
et al. (2004); Sanmiguel Vila et al. (2020), normalized by free stream
velocity $U_{\infty}$ and $\delta_{99}$; (Left) FPG in magenta and SFPG in
black from NDF; (Right) APG in red from MTL and APG in blue from NDF in blue;
dashed red lines at $y/\delta_{99}=1.25$, corresponding to $\Delta_{1.25}$.

With increasing Reynolds number, the higher-order terms (H.O.T.) starting with
$B_{1}/Re_{\tau}$ may be ignored. The pressure-gradient parameter defined
using the outer scale $\Delta_{1.25}$ calculated for all profiles of boundary
layers from the NDF and MTL data sets are displayed in figure 2. Note that in
figure 3 of Sanmiguel Vila et al. (2020) the pressure-gradient parameter used
was $\beta^{*}$ utilizing $\delta^{*}$ as a length scale and yielding
conditions of “approximately constant” or “mildly increasing” APG.

Figure 2: Reynolds-number dependence of pressure-gradient parameter, $\beta$,
in boundary layers of data sets used; APG data in red dots with second-order
best fit for two configurations in MTL by Sanmiguel Vila et al. (2020),
“mildly increasing” APG (dotted red line) and “approximately constant” APG
(red line); data from NDF by Nagib et al. (2004) for three configurations at
six streamwise locations and three freestream velocities, FPG in magenta, SFPG
in black and APG in blue; black lines depict ZPG (solid) fully developed
channel (dotted) and fully developed pipe (dashed) conditions.

Finally, the commonly used indicator function based on the mean velocity
profile for wall-bounded turbulence, $\Xi$, can be obtained from:

$\Xi=y^{+}\frac{{\rm d}{U}_{x}^{+}}{{\rm d}y^{+}}=Y\frac{{\rm
d}{U}_{x}^{+}}{{\rm d}Y}.$ (5)

From equations (4) and (5), one can obtain the equation for $\kappa$ and
$S_{0}$:

$\Xi_{\rm OL}=\kappa^{-1}+S_{0}y^{+}/Re_{\tau}=\kappa^{-1}+S_{0}Y.$ (6)

##  3 Determining overlap coefficients

Applying equation (5) to mean-velocity-profile data of any wall-bounded flow
in inner or outer distances from the wall, $y^{+}$ or $Y$ respectively, can
result in determining $\kappa$ and $S_{0}$ by a linear fit corresponding to
equation (6). One needs to only select the range in $y^{+}$ or $Y$ to fit the
linear relations. In figure 3b of Lee & Moser (2015), for $Re_{\tau}\approx
5200$, they applied a linear fit over around $300<y^{+}<700$ ($0.06<Y<0.14$)
and obtained $\kappa=0.384$ and $S_{0}=0$. Instead, if they had applied a
linear fit over $0.21<Y<0.45$ ($1100<y^{+}<2300$) they would have obtained
$\kappa=0.41$ and $S_{0}=1.11$, which compare very well with values for
channel flows obtained by MN23, i.e., $\kappa=0.417$ and $S_{0}=1.1$.
Recently, Hoyas et al. (2023) obtained values of $\kappa=0.51$ and $S_{0}=2.6$
also for channel flows computed with two different resolutions using this
approach over $0.3<Y<0.6$ for $Re_{\tau}=550$. The value of $B_{0}=6.9$ was
then obtained using the already determined values of $\kappa$ and $S_{0}$ by
minimizing the differences between the velocity profile data to equation (4).

This approach requires highly accurate and closely spaced data for the
velocity profiles to apply the derivatives of the indicator function $\Xi$.
This is often the case from direct-numerical-simulation (DNS) data, and to
some degree in the boundary-layer data from the NDF, where over $200$ data
points with equal physical spacing were measured. Typical wall-bounded-flow
measurements, including in boundary layers, collect the measurements using
points with logarithmic physical spacing, i.e., with increasing physical
separation away from the wall. This has been a tradition based on an expected
“log law” in the overlap region. The MTL data set was obtained in this way
with only $45$ to $75$ measurement points in the wall-normal direction. The
profiles of figure 3 display the measured points and demonstrate the contrast
in spacing between MTL (red) and NDF (blue, green, magenta and black). While
determining the indicator function $\Xi$ from the MTL data was attempted using
several ways, it was found to be not reliable and produced noisy results. On
the other hand, two examples of $\Xi$ obtained from NDF profiles are shown in
figure 4. To process both data sets using the same approach and algorithms a
different method was required. However, the $\Xi$ profiles from all of the NDF
data were very important to select the range of $Y$ to use for the overlap
region. Based on the experience from MN223 and Hoyas et al. (2023), it is
found that fitting in the outer-scaled distances $Y$ from the wall is more
representative of the overlap region than by using inner-scaled distances
$y^{+}$, over wider ranges of $Re_{\tau}$. It also helps in avoiding the wake
region especially for APG conditions, and other pressure-gradient boundary
layers, including ZPG.

Figure 3: (Left) Inner-scaled mean velocity profiles plotted on linear scales
against $y/\Delta_{1.25}$ for boundary layers of data sets used, FPG
(magenta), SFPG (black), ZPG (green), and APG (blue) from NDF, and APG (red)
from MTL; dashed green lines at $y/\Delta_{1.25}=0.05~{}\&~{}0.2$ correspond
to dashed green lines in Figure 4; (Right) smoke-wire visualization of ZPG
boundary layer in NDF at $Re_{\theta}$ around $1,800$ with arrows
corresponding to approximate locations of $y/\Delta_{1.25}=0.05~{}\&~{}0.2$.
Figure 4: (Top) Inner-scaled mean velocity profile data plotted against
logarithmic-scaled $y/\Delta_{1.25}$ for two boundary layers from NDF, for FPG
(magenta) and APG (blue); green lines represent the log+lin fit of the overlap
region between two vertical dashed green lines at
$y/\Delta_{1.25}=0.05~{}\&~{}0.2$, and vertical dashed red lines correspond to
$y/\Delta_{1.25}=0.05~{}\&~{}0.11$, which is consistent with MN23; (Bottom)
Indicator function, $\Xi$ plotted against same scale of top part.

The approach used here is different from the one for channel and pipe flows
described above as we are retrieving all three parameters of equation (4) at
the same time, by utilizing all profiles in $U^{+}$ versus Y (=
$y/\Delta_{1.25}$). After choosing the limits of the overlap region in Y, we
compute the best combination of parameters using nonlinear least-squares curve
fitting. We fit our function to the data by adjusting the parameters of the
function to minimize the sum of squared differences between the observed and
predicted values. This optimization involves minimizing the sum of squared
residuals.

Given the set of data points ($Y_{i}$, $U^{+}_{i}$) and the model function
$f(\kappa,~{}S_{0},~{}B_{0},~{}Y_{i})$, represented by equation (4), we
approximate the combination of parameters that minimizes the difference
between the equation and the data for each mean velocity profile: $\forall
N\in\mathbb{N}^{*}$, $(\kappa_{\rm opt},~{}S_{0_{\rm opt}},~{}B_{0_{\rm
opt}})=\underset{\kappa,~{}S_{0},~{}B_{0}}{\rm
argmin}\sum_{n=1}^{N}(U^{+}_{i}-f(\kappa,~{}S_{0},~{}B_{0},~{}Y_{i}))^{2}$.
This optimization algorithm, with an accuracy of the order of $10^{-8}$, leads
to the final values of the three parameters of the log+lin overlap for each
velocity profile, using $Y_{\rm min}=0.05$ and $Y_{\rm max}=0.2$.

##  4 Discussion and conclusion

With the approach described in the last part of the previous section, a
parametric evaluation of the range of $Y$ to extract the parameters of the
log+lin overlap from the mean velocity profile data was conducted using the
NDF data first and then extended to the MTL profiles to arrive at a single
range for all data. The experience from MN23, especially for ZPG, was a good
starting point. For the lower limit $Y_{\rm OL,min}$ we arrived at a value of
0.05 to accommodate the lower range of $Re_{\tau}$’s in the MTL data. For the
upper limit of the range $Y_{\rm OL,max}$, we tested several values in the
range $0.11<Y<0.4$ (compare to figure 7 of MN23). Based on both sets of data,
and in order to avoid including parts of the wake, which varies with pressure
gradient, we again converged on $Y_{\rm OL,max}=0.2$. Two ranges are indicated
in figure 4 with dashed lines representing $Y_{\rm OL,min}=0.05$ and $Y_{\rm
OL,max}=0.11$ in red, and $Y_{\rm OL,min}=0.05$ and $Y_{\rm OL,max}=0.2$ in
green. The resulting log+lin overlap fit based on equation (6) is shown in
solid curves, also in green, with the mean velocity profiles corresponding to
an FPG case and an APG case. Note that for the FPG case with the smaller wake,
as indicated by the $\Xi$ curves, the profile data are above the log+lin
overlap fit beyond $Y>0.3$. Further systematic study of the best choice may
lead to an ultimate best value of $Y_{\rm OL,max}$ and its possible dependence
on pressure gradient, but we are satisfied with the $0.2$ value especially in
light of the results presented in figures 5 and 6.

Boundary layers with small pressure gradients $-0.5<\beta<0.5$, including
ZPGs, are the most sensitive to the choice of $Y_{\rm OL,max}$ for the
accurate extraction of the log+lin overlap parameters and in particular
$S_{0}$, because of the influence of the wake part. These profiles come from a
mixture of NDF cases and the MTL cases, with the fewer data points in the
profile, and the resulting parameters are included in figure 5 with grey
crosses. The discrimination between the log+lin overlap region and the start
of the outer wake region requires closely spaced measurement points.

On the left part of figure 3 we included the selected values $Y_{\rm
OL,min}=0.05$ and $Y_{\rm OL,max}=0.2$ on the velocity profiles also in green
dashed lines. The arrows starting from those values of $Y=y/\Delta_{1.25}$ to
the upstream edge of the flow-visualization image in the right part of the
figure are intended to point the approximate corresponding distances from the
wall at the bottom of the image.

Figure 5: Best-fit overlap parameters, $\kappa$, $S_{0}$ and $B_{0}$ based on equation (4) plotted against pressure-gradient parameter, $\beta$, for all data sets of Figure 3; NDF in open blue circles and MTL in closed red circles, with ZPG from NDF and Samie et al. (2018) in large black open circle; parameter values from Monkewitz & Nagib (2023) in closed cyan circles; parameter values from DNS of Hoyas et al. (2023) for channel (large brown open circles) and pipe (large red open square); recent experiments in CICLoPE Pipe in magenta triangle; green lines are second-order best fits of all figure data with dashed red line representing 20% deviations from fit (see Table 1); faint grey crosses represent profiles with high fit uncertainty.  Favorable Pressure Gradient | Adverse Pressure Gradient
---|---
$\kappa_{-}=-0.044~{}\beta^{2}-0.142~{}\beta+0.384$ | $\kappa_{+}=-0.006~{}\beta^{2}+0.056~{}\beta+0.384$
$S_{0_{-}}=-2.516~{}\beta^{2}-5.652~{}\beta$ | $S_{0_{+}}=-0.569~{}\beta^{2}+5.068~{}\beta$
$B_{0_{-}}=-0.974~{}\beta^{2}-2.261~{}\beta+4.17$ | $B_{0_{+}}=-0.275~{}\beta^{2}+2.114~{}\beta+4.17$
Table 1: Equations for second-order best fits of overlap parameters, $\kappa$,
$S_{0}$ and $B_{0}$,
for NDF and MTL Data

All the data for best-fit values of $\kappa$, $S_{0}$ and $B_{0}$ for each
boundary-layer velocity profile from both datasets are included in figure 5.
Furthermore, data from MN23, Samie et al. (2018) and Hoyas et al. (2023) are
added to the figure. It is important to point out that MN23 defined their
pressure-gradient parameter $\beta$ without a minus sign so pipe and channel
flows were represented by positive values of $\beta$. All the data in each
part of figure 5 were used to fit second-order polynomials, using an algorithm
to minimize the differences with an accuracy of the order of $10^{-3}$,
resulting in the summary of Table 1.

A correlation to test the universality of the Kármán coefficient was developed
by Nagib & Chauhan (2008) and is shown by their equation (13), which was based
on a log-only overlap region:

$\kappa B=1.6~{}[e^{0.1663B}-1].$ (7)

For a constant or universal Kármán constant, $\kappa$, a linear relation
between $\kappa B$ and $B$ needs to exist. In figure 6 we test the
universality of $\kappa$ based on the log+lin overlap region and plot all the
data from both sets. Lines of various values of constant $\kappa$ are included
for reference and none of them represent the entire data sets or segments of
them with any reasonable accuracy. Fitting the data extracted using the
log+lin overlap by a similar relation to equation (4) we find, using an
algorithm to minimize the differences with an accuracy of the order of
$10^{-3}$, slightly different coefficients in the relation:

$\kappa B_{0}=4.05~{}[e^{0.087B_{0}}-1].$ (8)

Figure 6: Product of $\kappa B_{0}$ plotted against $B_{0}$ to test
universality of Kármán coefficient against lines of constant $\kappa$ based on
best-fit overlap parameters from equation (4); NDF data in blue dots and MTL
data in red dots; magenta line is best fit according to equation (8) with
$10\%$ deviation indicated by dashed lines; green curve is similar fit of a
wide range of data using composite profiles by Nagib & Chauhan (2008) for a
log-only overlap region, shown in equation (7).

Based on the results of figure 5 we conclude that when using the log+lin
overlap given by (4), no history effects are observed in the overlap region in
either of the two data sets used here or between them, over a very wide range
of Reynolds numbers and adverse- and favorable-pressure gradients. We
interpret this behaviour as quasi equilibrium of the overlap between outer and
inner parts when the new understanding of the overlap region of Monkewitz &
Nagib (2023) is adopted. In the literature, some history effects have been
observed in pressure-gradient boundary layers but the well-documented ones
(Harun et al., 2013; Bobke et al., 2017) are reflected primarily in Reynolds
stresses or the wake part of the mean flow. We only used the approximately
constant and mildly increasing pressure gradient cases from MTL. Sanmiguel
Vila et al. (2020) data included a “strongly increasing” case with $\beta^{*}$
up to 2.38. The validity of our conclusions on the history effects will need
to be examined for a wider range of $\beta$ as defined here with an outer
length scale.

## Acknowledgements

VB and HN acknowledge the ongoing support of the John T. Rettaliata Chair of
Mechanical and Aerospace Engineering. RV acknowledges the financial support
from ERC grant no. ‘2021-CoG-101043998, DEEPCONTROL’. Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of European Union or European Research Council. Neither the
European Union nor granting authority can be held responsible for them.

## Declaration of interests

The authors report no conflict of interest.

## References

  * Afzal & Yajnik (1973) Afzal, N. & Yajnik, K. 1973 Analysis of turbulent pipe and channel flows in moderately large reynolds numbers. Journal of Fluid Mechanics 61, 23–31.
  * Bobke et al. (2017) Bobke, A., Vinuesa, R., Örlü, R. & Schlatter, P. 2017 History effects and near equilibrium in adverse-pressure-gradient turbulent boundary layers. J. Fluid Mech. 820, 667–692.
  * Chauhan et al. (2009) Chauhan, K., Monkewitz, P. & Nagib, H. 2009 Criteria for assessing experiments in zero pressure gradient boundary layers. Fluid Dynamics Research 41, 021404.
  * George & Castillo (2006) George, W. K. & Castillo, L. 2006 Recent advancements toward the understanding of turbulent boundary layers. AIAA J. 44, 2435–2449.
  * Harun et al. (2013) Harun, Z., Monty, J. P., Mathis, R. & Marusic, I. 2013 Pressure gradient effects on the large-scale structure of turbulent boundary layers. J. Fluid Mech. 715, 477–498.
  * Hoyas et al. (2023) Hoyas, S., Vinuesa, R., Baxerras, V. & Nagib, H. 2023 Resolution and convergence requirements for extended overlap region in wall-bounded turbulence. arXiv:2311.05204v1 [physics.flu-dyn] pp. 1–6.
  * Lee & Moser (2015) Lee, M. & Moser, R. 2015 Direct numerical simulation of turbulent channel flow up to ${R}e_{\tau}\approx 5200$. J. Fluid Mech. 774, 395–415.
  * Luchini (2017) Luchini, P. 2017 Universality of the turbulent velocity profile. Phys. Rev. Lett. 118, 224501.
  * Monkewitz & Nagib (2023) Monkewitz, P.A. & Nagib, H.M. 2023 The hunt for the Kármán ‘constant’ revisited. J. Fluid Mech. 967, A15.
  * Monkewitz et al. (2008) Monkewitz, P. A., Chauhan, K. A. & Nagib, H. M. 2008 Comparison of mean flow similarity laws in zero pressure gradient turbulent boundary layers. Phys. Fluids 20, 105102.
  * Nagib et al. (2004) Nagib, H., Christophorou, C. & Monkewitz, P. 2004 High Reynolds number turbulent boundary layers subjected to various pressure-gradient conditions. Proceedings of IUTAM Symposium on One Hundred Years of Boundary Layer Research; Göttingen, Germany pp. 383–394.
  * Nagib & Chauhan (2008) Nagib, H. M. & Chauhan, K. A. 2008 Variations of von Kármán coefficient in canonical flows. Phys. Fluids 20, 101518\.
  * Samie et al. (2018) Samie, M., Marusic, I., Hutchins, N., Fu, M.K., Fan, Y., Hultmark, M. & Smits, A.J. 2018 Fully resolved measurements of turbulent boundary layer flows up to $Re_{\tau}=20000$. J. Fluid Mech. 851, 391–415.
  * Sanmiguel Vila et al. (2020) Sanmiguel Vila, C., Vinuesa, R., Discetti, S., Ianiro, A., Schlatter, P. & Örlü, R. 2020 Resolution and convergence requirements for extended overlap region in wall-bounded turbulence. Experimental Thermal and Fluid Sciences] 112, 109975\.
  * Vinuesa et al. (2014) Vinuesa, R., Schlatter, P. & Nagib, H. M. 2014 Role of data uncertainties in identifying the logarithmic region of turbulent boundary layers. Exp. Fluids 55, 1751.
  * Vishwanathan et al. (2023) Vishwanathan, V., Fritsch, D. J., Lowe, K. T. & Devenport, W. J. 2023 History effects and wall-similarity of non-equilibrium turbulent boundary layers in varying pressure gradient over rough and smooth surfaces. Int. J. Heat Fluid Flow 102, 109145\.
  * Yajnik (1970) Yajnik, K. 1970 Asymptotic theory of turbulent shear flows. Journal of Fluid Mechanics 42, 411–427.