# On the composition dependence of the microscopic structure, thermodynamic,
dynamic and dielectric properties of water-dimethyl formamide model mixtures.
Molecular dynamics simulation results

H. Dominguez, O. Pizio111On sabbatical leave from Instituto de Quimica de la
UNAM, corresponding author<EMAIL_ADDRESS>
(Received August 10, 2017, in final form September 11, 2017)

###### Abstract

Isothermal-isobaric molecular dynamics simulations have been performed to
examine an ample set of properties of the model water-N,N-dimethylformamide
(DMF) mixture as a function of composition. The SPC-E and TIP4P-Ew water
models together with two united atom models for DMF [Chalaris M., Samios J.,
J. Chem. Phys., 2000, 112, 8581; Cordeiro J., Int. J. Quantum Chem., 1997, 65,
709] were used. Our principal analyses concern the behaviour of structural
properties in terms of radial distribution functions, and the number of
hydrogen bonds between molecules of different species as well as thermodynamic
properties. Namely, we explore the density, excess mixing molar volume and
enthalpy, the heat capacity and excess mixing heat capacity. Finally, the
self-diffusion coefficients of species and the dielectric constant of the
system are discussed. In addition, surface tension of water-DMF mixtures has
been calculated and analyzed.

Key words: water models, N,N-dimethylformamide models, thermodynamic
properties, self-diffusion coefficient, dielectric constant, surface tension,
molecular dynamics

PACS: 61.20.-p, 61.20-Gy, 61.20.Ja, 65.20.Jk

###### Abstract

Çäéñíåíî ñèìóëÿö¿ ìåòîäîì ìîëåêóëÿðíî¿ äèíàìêè â çîòåðìíî-çîáàðèíîìó àíñàìáë ç
ìåòîþ âèâåííÿ øèðîêîãî íàáîðó âëàñòèâîñòåé ìîäåë ñóìø âîäà-N,N-äèìåòèëôîðìàìä
(DMF) ÿê ôóíêöé êîíöåíòðàö¿. Âèêîðèñòàíî ìîäåë âîäè SPC-E TIP4P-Ew ðàçîì ç
äâîìà îá’äíàíèìè àòîìíèìè ìîäåëÿìè äëÿ DMF [Chalaris M., Samios J., J. Chem.
Phys., 2000, 112, 8581; Cordeiro J., Int. J. Quantum Chem., 1997, 65, 709].
Íàø îñíîâíèé àíàëç ñòîñóòüñÿ ïîâåäíêè ñòðóêòóðíèõ âëàñòèâîñòåé â òåðìíàõ
ðàäàëüíèõ ôóíêöé ðîçïîäëó èñëà âîäíåâèõ çâ’ÿçêâ ìæ ìîëåêóëàìè ðçíèõ ñîðòâ, à
òàêîæ òåðìîäèíàìíèõ âëàñòèâîñòåé. Çîêðåìà, ìè äîñëäæóìî ãóñòèíó, íàäëèøêîâ
ìîëÿðíèé îá’ì òà åíòàëüïþ çìøóâàííÿ, ïèòîìó òåïëîìíñòü íàäëèøêîâó ïèòîìó
òåïëîìíñòü çìøóâàííÿ. Íàêíåöü, îáãîâîðþþòüñÿ êîåôöíòè ñàìîäèôóç¿ ñîðòâ
äåëåêòðèíà ñòàëà ñèñòåìè. Êðì òîãî, îáèñëþòüñÿ òà àíàëçóòüñÿ ïîâåðõíåâèé íàòÿã
ñóìøåé âîäà-DMF.

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

##  1 Introduction

N,N-dimethylformamide (DMF) is a multipurpose organic solvent with widespread
applications in chemistry and related areas of much practical importance, see
e.g., [1, 2, 3]. In several processes it is used as a co-solvent. On the other
hand, as one example of amides, it is frequently chosen as a model system for
peptides. Given the importance of DMF and of its mixtures with water and
several organic solvents, investigation of the properties of such complex
molecular systems has been undertaken using a variety of experimental
techniques. Namely, the studies of microscopic structure have been performed
using neutron and electron diffraction, NMR and dielectric spectroscopy; data
resulting from measurements of various properties, like density, excess molar
volume and enthalpy, excess molar heat capacity, self-diffusion coefficients,
dielectric constant have been reported in [4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. Some experimental
observations were supported by quantum chemical calculations [26, 27].

In addition, aiming at rather general and more profound insights into the
behavior of various properties on temperature, pressure and chemical
composition of water-DMF mixtures and at interpretation of experimental data,
computer simulation methods have been applied to these mixtures.
N,N-dimethylformamide is an interesting molecule, hydrophobic methyl groups
can disturb the water structure. Besides, “cross” hydrogen bonds are expected
to form between molecules belonging to two different species. Computer
modelling studies can be classified according to the details of the force
fields describing the internal structure of the DMF molecule. Specifically,
one set of models with either five or six force centers (in each of models the
methyl group has been considered as a single site) was explored in [28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39]. Another set of apparently more
sophisticated models, i.e., at all-atom level, has been designed and then
studied in detail in several publications [40, 41, 42, 43, 44]. A very
comprehensive description of the properties of an ample set of organic liquids
resulting from simulations with various force fields has been discussed in
[45] (see also the file that contains a vast supporting information of [45]
and the http://virtualchemistry.org web site).

In this work we restrict our attention solely to the modelling of water-DMF
mixtures and to the critical evaluation of several properties as a function of
composition resulting from isobaric-isothermal computer simulations. Moreover,
our attention is restricted to the united-atom type models for DMF molecule
developed by Cordeiro and by Chalaris and Samios. The reason of studying these
models, in addition to and in spite of the available knowledge, is that we
believe several even basic issues have not been comprehensively and critically
discussed. Specifically, we calculate density and excess mixing density,
excess mixing volume and enthalpy, heat capacity and excess mixing heat
capacity, self-diffusion coefficients of species of the mixture, dielectric
constant and the excess dielectric constant, average number of hydrogen bonds,
all as functions of mixture composition. The microscopic structure is briefly
discussed in terms of the pair distribution functions. Finally, we performed
calculations of the surface tension switching to the constant volume-constant
temperature canonical ensemble. For the majority of properties calculated, we
perform comparisons with the experimental results. Having in mind our recent
investigation of the properties of water-methanol [46, 47], water-DMSO [48]
and water-1,2-dimethoxyethane (DME) mixtures [49], the results of the present
work are discussed in a wider context in qualitative terms, making comparisons
of the trends observed for different organic liquids mixed with water.

Exploration of the behavior of water-DMF mixtures by using all-atom type
models for DMF will be presented elsewhere in a separate work. Our longer term
objective of the project, however, is in the theoretical exploration of the
behavior of solutes of different complexity in water-DMF and other water-
organic liquid mixtures in the spirit of important and interesting
experimental studies [50, 51, 52, 53].

##  2 Models and simulation details

Our calculations have been performed in the isothermal-isobaric (NPT) ensemble
(at 1 bar, and at a temperature of 298.15 K), unless specified in the surface
tension subsection. We used the GROMACS package [54] version 4.6.5. For water
in the present study, the SPC-E model [55] and the TIP4P-Ew model [56] were
used. As for the DMF molecules, we explore two united-atom type models with
six force sites, H, O, C, N, each methyl group, CH3 is considered as a single
site, C3. The nomenclature of the DMF models is as follows: we use the model
and notation CS2 as proposed by Chalaris and Samios, see table I of [34]. On
the other hand, the model denominated as Cordeiro is the one proposed by this
author and described in detail in table I of [30]. All the parameters for the
intramolecular geometry of the rigid CS2 model are given in table II of [34].
The model of Cordeiro is rigid as well. For the sake of numerical convenience,
in the present work we relaxed rigidity of bonds and angles, taking the force
constants from the OPLS data basis as appropriate. Moreover, to keep the DMF
molecule planar, the improper dihedral angle for H-C-N-C3 and for O-C-N-C3 has
been assumed from the OPLS data basis as well.

The geometric combination rule (rule 3 in GROMACS nomenclature) has been
applied for cross interactions in all our calculations. The nonbonded
interactions were cut off at 1.4 nm, and the long-range electrostatic
interactions were handled by the particle mesh Ewald method implemented in the
GROMACS software package (fourth order, Fourier spacing equal to 0.12). The
van der Waals tail correction terms to the energy and pressure were taken into
account. In order to maintain the geometry of water and DMF molecules, LINCS
algorithm was used.

For each system a periodic cubic simulation box was set up. The GROMACS genbox
tool was employed to randomly place all particles into the simulation box. The
total number of molecules was kept fixed at 3000. The composition of the
mixture is described by the mole fraction of DMF molecules, $X_{\text{dmf}}$,
$X_{\text{dmf}}=N_{\text{dmf}}/(N_{\text{dmf}}+N_{\text{w}})$.

To remove possible overlaps of particles introduced by the procedure of
preparation of the initial configuration, each system underwent energy
minimization using the steepest descent algorithm implemented in the GROMACS
package. Minimization was followed by a 50 ps NPT equilibration run at 298.15
K and 1 bar using a timestep of 0.25 fs. We applied the Berendsen thermostat
and barostat with $\tau_{T}=1$ ps and $\tau_{P}=1$ ps during equilibration.
Constant value of $4.5\times 10^{-5}$ bar-1 for the compressibility of the
mixtures was set up. The V-rescale thermostat and Parrinello-Rahman barostat
with $\tau_{T}=0.5$ ps and $\tau_{P}=2.0$ ps and the time step 2 fs were used
during all production runs. To test the thermostat and barostat, we have
obtained 88.5 J/mol$\cdot$K for the heat capacity of the SPC/E model, this
value being close to 86.6 J/mol$\cdot$K reported by Vega et al. [57]
(experimental result is 75.3 J/mol$\cdot$K). On the other hand, we obtained
155.55 J/mol$\cdot$K for the heat capacity of the DMF for the CS2 model, and
158.78 J/mol$\cdot$K for the Cordeiro model, the experimental values reported
from independent measurements are as follows: 150.16 J/mol$\cdot$K [14], 150.8
J/mol$\cdot$K [11]. Moreover, the DMF molar volume coming out from our
calculations is 76.53 cm3/mol (CS2 model) and 78.23 cm3/mol (Cordeiro model),
the deviation from the experimental result, 77 cm3/mol, is very small in the
case of CS2. Berendsen type control of temperature and pressure is not
satisfactory in this aspect.

Statistics for each mole fraction for some of the properties were collected
over several 10 ns NPT runs, each started from the last configuration of the
preceding run. The total trajectory was not shorter than 60 ns. Actually, the
heat capacity and the dielectric constant are the most time demanding
properties.

##  3 Results and discussion

###  3.1 Density, mixing properties

Figure 1: (Color online) Panel (a): Composition dependence of the density of
water-DMF mixtures from constant pressure-constant temperature simulations
($T=298.15$ K, $P=1$ bar) in comparison with the experimental data from [10,
11]. Panel (b): Excess mixing density on composition.

We begin a discussion of the results by presenting the dependence of density
of water-DMF mixture with united-atom type models for DMF on composition,
panel (a) of figure 1. In spite of a meticulous search, we have not found
previous comparisons of this sort. Two water models and two DMF models are
involved. The experimental data are taken from [10, 11]. Experimental results
show a peculiar behavior of mixture density at low values of $X_{\text{dmf}}$,
namely while DMF is added to water, the mixture density exhibits a weakly
pronounced maximum, at $X_{\text{dmf}}<0.2$, the density preserves values
close to the density of pure water, if $X_{\text{dmf}}$ is low. Actually,
$\rho(X_{\text{dmf}})$ decreases with DMF concentration only for
$X_{\text{dmf}}>0.2$. A combination of SPC-E and Cordeiro model leads to a
pronounced underestimation of density in the entire composition range. On the
other hand, the SPC-E-CS2 as well as TIP4P-Ew-CS2 exhibit a more satisfactory
behavior. Both models underestimate the values of density of the mixture at
small values of $X_{\text{dmf}}$ and slightly overestimate the density for
high values of $X_{\text{dmf}}$.

It is common to characterize the behavior of various properties of binary
mixtures in terms of excess mixing contribution as follows,

$\Delta Y^{\text{mix}}=Y_{\text{m}}-[X_{1}Y_{1}+(1-X_{1})Y_{2}],$ (1)

where $Y_{\text{m}}$ refers to the mixture at a given composition, $Y_{i}$
($i=1,2$) are the values for pure components, $Y$ can be any property of the
mixture, e.g., density or volume or enthalpy or another one.

Here, to begin with, we obtain a finer insight into the decay of density from
pure water to pure DMF. It is provided by the results given in panel (b) of
figure 1. Namely, the deviation of mixture density from an ideal type behavior
is closer to the experimental trends if the SPC-E or TIP4P-Ew water models are
combined with the CS2 model for DMF. On the other hand, predictions of the
SPC-E combined with Cordeiro model are not satisfactory. Positive deviation of
excess mixing density is similar to what was observed in mixtures of alcohols
with water in the work of Wensink [58]. In fact, the behavior of excess mixing
density can indicate the accuracy of theoretical curves, but it does not
provide explicit insights into the mixing trends; methanol or ethanol mix with
water perfectly well, but the $\Delta\rho$ is positive in the entire range of
composition.

Figure 2: (Color online) Excess mixing molar volume and excess mixing enthalpy
of water-DMF mixtures on composition.

The excess mixing volume and enthalpy are given in two panels of figure 2.
From panel (a) of this figure we learn that the excess mixing volume describes
a contraction of the mixture upon composition qualitatively well, although the
magnitude of $\Delta V^{\text{mix}}$ is underestimated. The minimum value from
the simulations is at $X_{\text{dmf}}\approx 0.5$, whereas the experimental
result attains minimum at $X_{\text{dmf}}=0.4$. The closest, to the
experiment, set of theoretical data comes out from the SPC-E-CS2 model.

One can see quite similar trends in the behavior of excess mixing enthalpy,
panel (b) of figure 2. The values for $\Delta H^{\text{mix}}(X_{\text{dmf}})$
are underestimated by all the models considered and the minimum occurs at a
slightly higher $X_{\text{dmf}}$ (between 0.4 and 0.5), in comparison to the
experimental result at $X_{\text{dmf}}\approx 0.35$. Again, the SPC-E-CS2
model is better comparing to, e.g., SPC-E-Cordeiro version of the model. In
general, one can conclude that the united-atom type models of this work
underestimate geometric and energetic trends for mixing of species given by
$\Delta V^{\text{mix}}$ and by $\Delta H^{\text{mix}}$ in the entire interval
of mixture compositions. The agreement with the experimental results is
qualitative for both properties in question.

In addition, it is worth studying the energetic trends of mixing of species in
terms of fluctuations. Namely, in figure 3 we explore the dependence of the
constant pressure heat capacity and the corresponding excess property as
functions of mixture composition. Actually, the situation is rather
satisfactory. From panel (a) of figure 3, we learn that trends of dependence
of $C_{P}$ on $X_{\text{dmf}}$ are well reproduced by all the models studied,
just the theoretical curves are shifted upward in comparison to the
experimental dependence. The shape and magnitude of $\Delta
C_{P}^{\text{mix}}$ is reproduced by the simulated models very well. Moreover,
the maximum is predicted at 0.4 whereas in the experiment it is at $\approx
0.3$. Seemingly, the SPC-E-CS2 combination of models is the best in this
aspect.

Figure 3: (Color online) Panel (a): Composition dependence of the heat
capacity at constant pressure of water-DMF mixtures from constant pressure-
constant temperature simulations ($T=298.15$ K, $P=1$ bar) in comparison with
the experimental data from [14, 11]. Panel (b): Excess mixing constant
pressure heat capacity on composition in comparison with experimental data.

After exploring thermodynamic aspects of mixing, we proceed to the brief
description of the microscopic structure in terms of the pair distribution
functions and hydrogen bonding between molecules.

###  3.2 Pair distribution functions, hydrogen bonding

All the pair distribution functions discussed below result from simulations of
the SPC-E-CS2 model. We have checked, however, that other models, like
TIP4P-Ew-CS2 and SPC-E-Cordeiro, yield a qualitatively similar picture of the
trends observed.

Figure 4: (Color online) Evolution of the pair distribution functions, OD-OD,
with changing solvent composition.

The system is characterized by a large set of the site-site distribution
functions. Our focus is only in some of them, namely in the functions that
describe most essential changes of the structure of the system upon changes of
its chemical composition. Evolution of the distribution of DMF molecules on
changes of $X_{\text{dmf}}$ is shown in terms of the distribution functions
between DMF oxygens, OD-OD (notation OD is used as an abbreviation for
O${}_{\text{dmf}}$), see two panels of figure 4. If the DMF molecules are at
small amount in the medium with predominant number of waters, the DMF
particles seem preferring to stay rather close to one another, as it follows
from the first two maxima of the pair distribution function (pdf) in panel (a)
of figure 4. However, even the first maximum is not high, the second one is
much lower, hence mutual contacts of two and more DMF molecules are not very
probable. If the concentration of DMF species increases, their distribution
becomes even more “homogeneous”, see, e.g., the curves are
$X_{\text{dmf}}=0.2$ and $0.3$ in panel (a). In other words, the first and
second maxima as well as the minimum between them, all slightly decrease in
magnitude and at larger inter-particle separations, the pdf is close to unity.
This rather uniform distribution of DMF species, does not change significantly
in an ample interval of compositions, say from $X_{\text{dmf}}=0.3$ up to
$X_{\text{dmf}}=0.9$. Moreover, we do not observe a shift of the position of
the first maximum of the pdf. It is located at $r\approx 0.43$ nm. Only, the
second maximum of the OD-OD distribution slightly grows and shifts to smaller
inter-particle separations when the composition changes from
$X_{\text{dmf}}=0.4$ up to $X_{\text{dmf}}=0.9$, indicating a weakly enhanced
“crowding” of the DMF particles if their amount increases. At this point we
would like to refer to the visualization of the distribution of a small amount
of DMF molecules in the “sea” of waters, see figure 5, panel (a). This
snapshot corresponds to the final configuration of particles after $\approx
60$ ns. Both, the pdf at a low value of $X_{\text{dmf}}$ and the snapshot at
similar conditions, permit to conclude that there is no clustering of DMF
species.

Figure 5: (Color online) Visualization of the distribution of DMF molecules
(nitrogen atoms — blue spheres) in water medium ($X_{\text{dmf}}=0.03333$) in
the left-hand panel (water oxygens and hydrogens — red and grey spheres,
respectively), and of the distribution of water molecules in the DMF medium
(nitrogens — blue spheres) at $X_{\text{dmf}}=0.9$ in the right-hand panel
(water oxygens and hydrogens — red and green spheres, respectively).

Let us proceed now to the description of the evolution of distribution of
water molecules in terms of the pdf OW-OW, figure 6. The most pronounced trend
is that the first maximum of $g_{\text{OW-OW}}$ strongly increases in
magnitude upon increasing $X_{\text{dmf}}$, i.e., when the concentration of
water species decreases. However, important observation is that such a
behavior is not accompanied by a substantial growth of the second maximum of
this function. To summarize, the contacts between water molecules become more
and more probable whereas the probabilistic aspects of the structure of
organic subsystem are less affected by changes of $X_{\text{dmf}}$. Water
molecules efficiently fill the space between the DMF particles, cf. the OW-OW
first maximum grows at $r\approx 0.28$ nm whereas the OD-OD first maximum
remains almost intact at $r\approx 0.43$ nm. The shape of OW-OW pdf together
with the typical configuration of molecules given in the snapshot, right-hand
panel of figure 5, leads to the conclusion that water molecules are rather
uniformly distributed in the “sea” of DMF particles, closely situated waters
are restricted to pairs at most.

Figure 6: (Color online) Evolution of the pair distribution functions, OW-OW,
with changing solvent composition.

Figure 7:  (Color online) Evolution of the pair distribution functions, OW-OD,
and OD-HW with changing solvent composition.

Next, we would like to explore the trends of cross correlations between
particles belonging to different species in this binary mixture. We analyze
them by using the OW-OD and OD-HW distribution functions. A few representative
examples are shown in figure 7. The shape of the pdf OW-OD changes principally
at small interparticle separations. Namely, the first maximum increases in
magnitude with an increasing $X_{\text{dmf}}$ in the entire range of
composition, starting from low values of $X_{\text{dmf}}$ up to the high
values. Its position on the $r$-axis does not change, however. As we have seen
from figure 4, the shape of OD-OD is weakly affected by $X_{\text{dmf}}$
values, thus changes of the function OW-OD are quite small in the interval
from $r\approx 0.4$ nm up to $r\approx 1$ nm. In contrast to this behavior,
the function OD-HW is sensitive to the value of $X_{\text{dmf}}$ in the region
of the first and of the second maximum. The first maximum of the OD-HW pdf
substantially increases in magnitude with an increasing $X_{\text{dmf}}$.
Concerning the second maximum of this function, it can be seen that its growth
can be related to changes of the OW-OW pdf. Namely, the second maximum of OD-
HW becomes higher than unity for $X_{\text{dmf}}>0.4$. Development of this
maximum can be attributed to an essential growth of OW-OW first maximum at a
high concentration of DMF molecules in the mixture. It is intuitively expected
that hydrogen bonded water molecules occupying the space close to OD atom
predominantly contribute to this second maximum of OD-HW, because if hydrogens
belonging to two different water molecules are close to OD, they should
contribute to the first maximum.

Figure 8:  (Color online) Changes of the average number of water-water and
water-DMF hydrogen bonds (per water molecule) upon composition changes. The
cross water-DMF bonds for SPC-E-CS2 model are given by solid black triangles
up.

Our final concern in this subsection is in the trends of behavior of the
average number of hydrogen bonds (normalization is performed per water
molecule) between water molecules as well as cross water-DMF bonds (figure 8).
The GROMACS software using a default distance-angle criterion was applied to
count the numbers of bonded molecules in each case. The behavior of $\langle
n_{\text{HB}}\rangle$ as a function of $X_{\text{dmf}}$ is qualitatively
similar to what was recently discussed for other water-organic liquid mixtures
[48, 49]. The average number of H-bonds between water molecules decreases
whereas the fraction of the bonds between water molecules and DMF oxygen
steadily increases with an increasing fraction of DMF species. Only at very
low values of $X_{\text{dmf}}$, the value for water-water $\langle
n_{\text{HB}}\rangle$ is high, yielding evidence of the hydrogen bonded
network between water molecules. If DMF species are introduced into water, the
number of cross contacts increases, as it follows from the OW-OD and OD-HW
functions, cf. panels (a) and (d) of figure 7. At high values of
$X_{\text{dmf}}$, the average value of cross bonds becomes higher than the one
between waters. This is explicable straightforwardly, because at high
concentration of DMF species in the mixture there is no room to form fragments
of expanded network of hydrogen bonds and, besides, the majority of water
molecules are surrounded by DMF, either as individual entities or very small
groups such as pairs. It is difficult to establish precisely at what values of
$X_{\text{dmf}}$ the hydrogen bonded network cease to exist as such. That
would require additional and a more profound analysis of a set of topologic
elements characterizing the cooperativity of bonds. All three models explored
in the present work (SPC-E-CS2, TIP4P-Ew-CS2 and SPC-E-Cordeiro), provide a
qualitatively similar picture concerning $\langle
n_{\text{HB}}\rangle(X_{\text{dmf}})$. Fine details depend on water modelling
and on the internal structure of DMF molecule. A model that predicts slightly
more probable water-water bonding in consequence leads to a slightly smaller
fraction of water-DMF bonds. For the moment, we have not explored the life-
time and strength of the bonds formed in the system, however.

Nevertheless, one more remark is pertinent. As we have seen above, maxima on
the experimentally measured excess mixing properties are observed in the
interval $0.3<X_{\text{dmf}}<0.4$ (cf. figures 2 and 3) whereas the simulation
predictions for the extrema are in the interval $0.4<X_{\text{dmf}}<0.5$. This
interval of compositions seems to correspond to the systems in which the
effect of hydrogen bonded network gradually ceases. It is tempting to
attribute this change of the behavior to the range of compositions in which
the curves for $\langle n_{\text{HB}}\rangle(X_{\text{dmf}})$ start to deviate
from linearity. However, this hypothesis should be explored and possibly
confirmed by more sophisticated tools. Now, we would like to see if these
observations of the trends of thermodynamic and structural properties are
related to dynamic properties.

###  3.3 Self-diffusion coefficients

Figure 9:  (Color online) Composition dependence of the self-diffusion
coefficients of species in water-DMF mixtures from constant pressure-constant
temperature simulations ($T=298.15$ K, $P=1$ bar). In panel a the experimental
data are taken from [9] at $T=293$ K.

One of the principal objectives of the present study is to combine a
discussion of an ample set of properties of water-DMF mixtures. Their dynamic
properties undoubtfully deserve a separate work. Here, we focus solely on the
self-diffusion coefficients of water and DMF species. As usual, they are
calculated from the mean-square displacement (MSD) of a particle via Einstein
relation,

$D_{i}=\frac{1}{6}\lim_{t\rightarrow\infty}\frac{\mathrm{d}}{\mathrm{d}t}\langle|{\bf
r}_{i}(\tau+t)-{\bf r}_{i}(\tau)|^{2}\rangle,$ (2)

where $i$ refers to water or DMF molecules and $\tau$ is the time origin.
Default settings of GROMACS were used for the separation of the time origins.
A set on trajectories coming from several consecutive simulations of $10$ ns
were combined to get an entire trajectory not less than 60–70 ns. The fitting
interval then has been chosen from $\approx 10$% to $\approx 50$% of the
entire trajectory, to obtain $D_{\text{dmf}}$ and $D_{\text{w}}$. A set of
results coming from the application of three models is given in two panels of
figure 9. In spite of an extensive search in the literature, we have found
only one set of experimental data describing $D_{\text{dmf}}$ self-diffusion
coefficient at $T=293$ K. Our calculations refer to $T=298.15$ K, however. An
overall shape of $D_{\text{dmf}}(X_{\text{dmf}})$ is similar from experiments
and from simulations of the models, panel (a) of figure 9. At each end of the
curves there is a maximum value separated by a minimum, the minimum in
experiment is at $X_{\text{dmf}}\approx 0.2$ whereas all the simulated models
predict a minimum at $\approx 0.3$. Moreover, the models used in simulations
yield lower values for $D_{\text{dmf}}$ in a wide interval of compositions,
from 0.1 to 1, comparing to experimental data. Only at a rather low
$X_{\text{dmf}}$, less than 0.1, the simulations predict higher values for
$D_{\text{dmf}}$ than from experiment.

The self-diffusion coefficient of water species in binary water-DMF mixtures
is presented in panel (b) of figure 9. Unfortunately, we were unable to find
the experimental data in spite of the efforts. One well documented point is
for pure water [57]. Starting from this point, the values for $D_{\text{w}}$
substantially decrease with an increasing $X_{\text{dmf}}$ till the minimum at
$X_{\text{dmf}}=0.5$. The value of $D_{\text{w}}$ at minimum is approximately
five times lower than for pure water. While $X_{\text{dmf}}$ increases
further, the self-diffusion coefficient, $D_{\text{w}}$, increases in
magnitude to reach the values above unity (twice higher value than at minimum)
for mixtures predominantly composed of DMF species. Thus, for water-rich
mixtures, $D_{\text{w}}$ is much higher than $D_{\text{dmf}}$, in contrast to
DMF-rich mixtures, where $D_{\text{dmf}}$ is higher than $D_{\text{w}}$. Note
that the $y$-scale is different in two panels. The magnitudes of
$D_{\text{dmf}}$ and of $D_{\text{w}}$ at respective minima are close to each
other.

Huge difference between the self-diffusion coefficients of each species at two
opposite sides of the composition axis leads to the conclusion that there
exists only a weak coupling of the dynamics of water and DMF particles. At low
$X_{\text{dmf}}$, the self-diffusion coefficient of water $D_{\text{w}}$ is
approximately twice higher than $D_{\text{dmf}}$. One can suppose that a
molecule of DMF is trapped in a water medium. The existence of the hydrogen
bond network does not suppress the self-diffusion of water molecules but
contributes to diminish self-diffusion of an individual DMF molecule in such a
surrounding. On the other hand, on the opposite side, for DMF-rich mixtures,
the self-diffusion of a water molecule is suppressed much less comparing to
DMF, cf. ratio of $D_{\text{dmf}}$ and of $D_{\text{w}}$ at high
$X_{\text{dmf}}$. The strongest coupling of the dynamics of two species is
observed at intermediate values of composition ($0.3<X_{\text{dmf}}<0.5$),
i.e., in the interval where the deviations from ideality of thermodynamic
characteristics are at maximum. At these composition values, it seems that the
effects of a hydrogen bonding between particles become much less important
since the network cease to exist as such, and just the packing and energetic
aspects of mixing for systems close to equimolarity determine the minima of
the self-diffusion coefficients. It would be of interest to extend the
exploration of the dynamic properties by calculating viscosity and
characteristic relaxation times for mixtures of different composition. A
broader set of properties — stronger conclusions should come out.

###  3.4 Composition dependence of the dielectric constant

Finally, we would like to put attention at one of various manifestations of
the dielectric properties. The long-range, asymptotic behavior of correlations
between molecules possessing, e.g., a dipole moment is determined by the
dielectric constant, $\varepsilon$. It is commonly accepted that long
molecular dynamics runs are necessary to obtain reasonable estimates for this
property, because $\varepsilon$ is calculated from the time-average of the
fluctuations of the total dipole moment of the system,

$\varepsilon=1+\frac{4\piup}{3k_{\text{B}}TV}\big{(}\langle\bf
M^{2}\rangle-\langle\bf M\rangle^{2}\big{)},$ (3)

where $k_{\text{B}}$ is the Boltzmann constant and $V$ is the simulation cell
volume. In this work we have taken care that each of the runs is of sufficient
length, not less than $60$ ns in the majority of cases. The curves coming out
from simulations for three models of water-DMF mixtures are shown in figure 10
(a). A general trend of the behaviour of $\varepsilon(X_{\text{dmf}})$ is that
it decreases with an increasing $X_{\text{dmf}}$, starting from a high value
for pure water to a lower value corresponding to pure DMF. As it follows from
the comparison of the simulation results and experimental data [13], all three
models substantially underestimate the values for $\varepsilon$ in the entire
composition range. The static dielectric constant for pure DMF in the
framework of CS2 model is $\varepsilon\approx$ 27.83, whereas for the model of
Cordeiro it is 28.52, the experiment yields 39.88 [13]. On the water-rich
side, the discrepancy between simulation data and experimental value evidently
results from water modelling. It is highly probable to improve the dependence
of the static dielectric constant on composition by applying the model for
water specifically parametrized to reproduce the dielectric constant [59].
Simultaneously, it would require parametrization of the force field for the
DMF molecule. This task seems to be attainable.

Figure 10:  (Color online) Panel (a): Composition dependence of the static
dielectric constant of water-DMF mixtures from constant pressure-constant
temperature simulations ($T=298.15$ K, $P=1$ bar) in comparison with the
experimental data from [13]. Panel (b): Excess mixing static dielectric
constant on composition in comparison with experimental data.

Another sensitive test is provided by a comparison of the excess dielectric
constant [see equation (1)],
$\Delta\varepsilon^{\text{mix}}=\varepsilon_{\text{m}}-[X_{\text{dmf}}\varepsilon_{\text{dmf}}+(1-X_{\text{dmf}})\varepsilon_{\text{w}}]$,
with the experimental predictions [13]. Experimental points indicate a
negative deviation from ideality in the entire composition range, figure 10
(b). Interestingly, this behaviour is contrary to what follows for water-DMSO
liquid mixtures [48] but similar to what we recently discussed for water-DME
mixtures [49]. Maximal (negative) deviation from the ideal type behaviour
reported from the experimental measurements is at $X_{\text{dme}}\approx
0.45$. The simulations results reproduce the position of a minimum
approximately; namely the minimum is in the interval 0.3–0.4, dependent on the
model employed. These trends are in concordance with our observations
concerning deviations from ideality of thermodynamic properties and self-
diffusion coefficients. However, the magnitude of the excess static dielectric
constant is overestimated from the simulated models. The TIP4P-Ew combined
with CS2 model is the most close to the experimental data. It is worth
mentioning that the excess dielectric constant curve as a function of chemical
composition of the mixture can be related to the excess refractive index
measurements, see, e.g., [15]. This issue has not been studied at the level of
the united atom models for DMF so far.

All the above results have been obtained by using isobaric-isothermal
simulations of water-DMF mixtures. The final configurations of particles,
however, can be conveniently used to explore one of the interfacial properties
by switching to the canonical, $NVT$, ensemble, namely the surface tension on
composition.

###  3.5 Surface tension of water-DMF mixtures

The simulations aiming in surface tension calculations at each point of
composition axis have been performed by taking the final configuration of
particles from the isobaric run. Next, the box edge along $z$-axis was
extended by a factor of 3, generating a box with liquid slab and two liquid-
mixture-vacuum interfaces in the $x$-$y$ plane, in close similarity to the
procedure applied in [45]. The total number of particles, three thousand, is
reasonable to yield an area of the $x$-$y$ face of the liquid slab
sufficiently big. The elongation of the liquid slab along $z$-axis is
sufficient as well. The executable molecular dynamics file was modified by
deleting a fixed pressure condition just preserving V-rescale thermostatting
with the same parameters as in $NPT$ runs. Other corrections have not been
employed.

The values for the surface tension, $\gamma$, follow from the combination of
the time averages for the components of the pressure tensor,

$\gamma=\frac{1}{2}L_{z}\Big{\langle}\Big{[}P_{zz}-\frac{1}{2}(P_{xx}+P_{yy})\Big{]}\Big{\rangle},$
(4)

where $P_{ij}$ ($i,j=x,y,z$) are the components of the pressure tensor, and
$\langle\ldots\rangle$ denotes the time average. We performed a set of runs at
a constant volume, 5–6 each with the time duration of $10$ ns, and obtained
the result for $\gamma$ by taking the block average.

As in the previous subsections, the results concern the SPC-E-CS2, SPC-E-
Cordeiro and TIP4P-EW-CS2 models. It is necessary to emphasize that our
calculations were not focused to precisely reproduce the values for $\gamma$
at a single point corresponding to pure water at a given density, see, e.g.,
[60]. We rather obtain $\gamma$ at a density of each mixture coming out from
the $NPT$ simulation. These values deviate from the experimental data as we
have discussed in the corresponding subsection reporting the results. The
experimental value of $\gamma$ reported in literature for water is 71.73 mN/m
whereas for DMF it is 34.4 mN/m. Our data for pure water models (SPC-E and
TIP4P-Ew) agree with the data given in [60] (see table III and table IV) prior
the application of the tail correction.

Figure 11:  (Color online) Composition dependence of the surface tension of
water-DMF mixtures from constant volume-constant temperature simulations
($T=298.15$ K) in comparison with the experimental data from [61].

A set of data obtained in the present simulation work is given in figure 11,
the experimental data are from the pioneering study [61]. The change of slope
of the surface tension experimental curve around 30 mole percent of DMF
together with the maximum of the viscosity curve at this composition (not
shown here), apparently may be related to the formation of a dihydrate of DMF.
Moreover, peculiarities of the mixing properties seen as a maximal deviation
from an ideal type of behavior and the minimum of $D_{\text{dmf}}$ support
this assumption. We are not aware of any spectroscopy results supporting this
point of view.

However, general trends of the experimental dependence of the surface tension
on $X_{\text{dmf}}$ are reproduced by the models used in simulations. The
values of $\gamma$ are underestimated in the water-rich mixtures. It seems
that a better description of the surface tension of pure water should decrease
this inaccuracy. By contrast, in the DMF-rich mixtures, the values for surface
tension of mixtures are closer to the experimental data. The CS2 model for DMF
satisfactorily describes the surface tension of a pure organic liquid. The
Cordeiro model also yields a reasonable value of the surface tension. On the
other hand, the application of TIP4P-Ew water model does not provide an
improvement of the trends of behavior of $\gamma$, comparing to the SPC-E-CS2
model. Computer simulation results of this set of models do not evidence a
peculiar behavior of surface tension at $X_{\text{dmf}}\approx 0.3$. We
proceed now to the final comments, summary and conclusions.

##  4 Summary and conclusions

The mixtures explored in this work are an example of combined water-organic
liquid solvents. We performed extensive molecular dynamic simulations in the
isobaric-isothermal ensemble to study the density, an ample set of mixing
properties and the microscopic structure of water-DMF mixtures in the entire
range of solvent composition. The self-diffusion coefficients of species and
the dielectric constant were calculated as well. All the simulations were
performed at room temperature and ambient pressure, 1 bar. Two water models
(SPC-E and TIP4P-Ew), combined with the united atom type CS2 and Cordeiro
models for DMF were studied. These insights are complemented by the surface
tension simulations for water-DMF mixtures using constant volume-constant
temperature conditions. The values for volume were chosen from the results of
the isobaric simulations.

From a comparison with the available experimental data for different
properties and with the results of other authors on this and related systems,
we conclude that the predictions obtained are qualitatively correct and give a
physically sound picture of the properties explored. All the properties
investigated are sensitive to the composition of a water-organic liquid
solvent.

The principal conclusions of the present study can be resumed as follows. We
explored the evolution of the microscopic structure in terms of the pair
distribution functions. The simulation results witness that the structure of
the subsystem of DMF species is more inert or less sensitive to the
composition, in comparison to the structure of an aqueous subsystem. The pair
distribution functions evidence a heterogeneous density distribution at local
scale upon adding of DMF molecules. These trends of changes of the microscopic
structure bear similarity to water-DMSO and water-DME co-solvent mixtures,
see, e.g., [48, 49]. Possibly the “associated” species involving water and DMF
molecules can be formed in the system, with or without hydrogen bonds. This
observation is in accordance with the interpretation of experimental data.
Water-water and water-organic co-solvent average number of hydrogen bonds do
not show a peculiar behaviour, in comparison to qualitatively similar mixtures
of water with organic co-solvent [48, 49].

Dynamic properties are discussed in terms of the self-diffusion coefficients
of species, $D_{\text{w}}$ and $D_{\text{dmf}}$. Both of them exhibit a
minimum in the interval of composition that corresponds to most “packed”
structures according to the behaviour of a mixing volume. The values of self-
diffusion coefficients of the species for pure components are reasonably well
described. Concerning the dependence of the static dielectric constant on
composition, we observe that its excess is better described while using the
TIP4P-Ew water model combined with the CS2 model for DMF. It would be of
interest to relate the behaviour of the static dielectric constant on the
composition with refractive index and viscosity data. Finally, we obtained a
reasonable description of the dependence of surface tension on composition, a
better agreement with experimental data is expected from application of other
water models.

At the present stage of the development, there are several missing elements
worthwhile a detailed investigation. Namely, the application of all-atom force
fields for the DMF molecules for the description of water-DMF mixtures will be
reported in a separate publication. A comparison of united-atom type models
with all-atom modelling is of primordial importance. On the other hand,
insights into the behaviour of dynamic and dielectric properties by
exploration of, e.g., the relaxation times, viscosity, hydrogen-bonds life-
time and complex dielectric constant, would be desirable. Further research is
planned in exploration of solutions with complex molecules in water-DMF
mixtures along the lines of the available in literature experimental research.

## Acknowledgements

O.P. is grateful to D. Vazquez and M. Aguilar for technical support of this
work at the Institute of Chemistry of the UNAM.

## References

  * [1] Kroschwitz J.I., Seidel A. (Eds.), Kirk-Othmer Encyclopedia of Chemical Technology, Vol. 1, Wiley-Interscience, Hoboken, 2004.
  * [2] Fiorito A., Larese F., Molinari S., Zanin T., Am. J. Ind. Med., 1997, 32, 255,
doi:10.1002/(SICI)1097-0274(199709)32:3<255::AID-AJIM11>3.0.CO;2-U.

  * [3] Ohara M., Takagaki A., Nishimura S., Ebitani K., Appl. Catal. A, 2010, 383, 149,
doi:10.1016/j.apcata.2010.05.040.

  * [4] Ohtaki H., Itoh S., Yamaguchi T., Ishiguro S., Rode B.M., Bull. Chem. Soc. Jpn., 1983, 56, 3406,
doi:10.1246/bcsj.56.3406.

  * [5] Okada M., Ibuki K., Ueno M., Bull. Chem. Soc. Jpn., 2012, 85, 189, doi:10.1246/bcsj.20110233.
  * [6] Biliškov N., Baranović G., J. Mol. Liq., 2009, 144, 155, doi:10.1016/j.molliq.2008.11.004.
  * [7] Schoester P.C., Zeidler M.D., Radnai T., Bopp P.A., Z. Naturforsch. A: Phys. Sci., 1995, 50, 38,
doi:10.1515/zna-1995-0106.

  * [8] Guarino G., Ortona O., Sartorio R., Vitagliano V., J. Chem. Eng. Data, 1985, 30, 366, doi:10.1021/je00041a039.
  * [9] Volpe C.D., Guarino G., Sartorio R., Vitagliano V., J. Chem. Eng. Data, 1986, 31, 37, doi:10.1021/je00043a012.
  * [10] Bernal-García J.M., Guzmán-López A., Cabrales-Torres A., Estrada-Baltazar A., Iglesias-Silva G.A., J. Chem. Eng. Data, 2008, 53, 1024, doi:10.1021/je700671t.
  * [11] De Visser C., Perron G., Desnoyers J.E., Heuvelsland W.J.M., Somsen G., J. Chem. Eng. Data, 1977, 22, 74, doi:10.1021/je60072a016.
  * [12] Cilense M., Benedetti A.V., Vollet D.R., Thermochim. Acta, 1983, 63, 151,
doi:10.1016/0040-6031(83)80080-X.

  * [13] Kumbharkhane A.C., Puranik S.M., Mehrotra S.C., J. Solution Chem., 1993, 22, 219, doi:10.1007/BF00649245.
  * [14] Checoni R.F., Volpe P.L.O., J. Solution Chem., 2010, 39, 259, doi:10.1007/s10953-010-9500-6.
  * [15] Gofurov Sh., Ismailova O., Makhmanov U., Kokhkharov A., Int. J. Chem. Mol. Nucl. Mater. Metall. Eng., 2017, 11, 330.
  * [16] Ueno M., Mitsui R., Iwahashi H., Tsuchihashi N., Ibuki K., J. Phys. Conf. Ser., 2010, 215, 1,
doi:10.1088/1742-6596/215/1/012074.

  * [17] Scharlin P., Steinby K., Domańska U., J. Chem. Thermodyn., 2002, 34, 927, doi:10.1006/jcht.2002.0946.
  * [18] Han K.-J., Oh J.-H., Park S.-J., Gmehling J., J. Chem. Eng. Data, 2005, 50, 1951, doi:10.1021/je050209y.
  * [19] García B., Alcalde R., Leal J.M., Matos J.S., J. Phys. Chem. B, 1997, 101, 7991, doi:10.1021/jp9626374.
  * [20] Shokouhi M., Jalili A.H., Hosseini-Jenab M., Vahidi M., J. Mol. Liq., 2013, 186, 142,
doi:10.1016/j.molliq.2013.07.005.

  * [21] Chen L., Groß T., Lüdemann H.-D., Z. Phys. Chem., 2000, 214, 239, doi:10.1524/zpch.2000.214.2.239.
  * [22] Jadżyn J., Świergiel J., Phys. Chem. Chem. Phys., 2012, 14, 3170, doi:10.1039/C2CP23960D.
  * [23] Egorov G.I., Makarov D.M., Kolker A.M., Russ. J. Phys. Chem. A, 2007, 81, 528,
doi:10.1134/S003602440704005X.

  * [24] Bai T.-C., Yao J., Han S.-J., J. Chem. Thermodyn., 1998, 30, 1347, doi:10.1006/jcht.1998.0402.
  * [25] Bai T.-C., Yao J., Han S.-J., Fluid Phase Equilib., 1998, 152, 283, doi:10.1016/S0378-3812(98)00402-6.
  * [26] Schmid E.D., Brodbek E., J. Chem. Phys., 1983, 78, 1117, doi:10.1063/1.444895.
  * [27] Schmid E.D., Brodbek E., J. Mol. Struct. THEOCHEM, 1984, 108, 17, doi:10.1016/0166-1280(84)80095-0.
  * [28] Jorgensen W.L., Swenson C.J., J. Am. Chem. Soc., 1985, 107, 569, doi:10.1021/ja00289a008.
  * [29] Yashonath S., Rao C.N.R., Chem. Phys., 1991, 155, 351, doi:10.1016/0301-0104(91)80111-T.
  * [30] Cordeiro J.M.M., Int. J. Quantum Chem., 1997, 65, 709,
doi:10.1002/(SICI)1097-461X(1997)65:5<709::AID-QUA37>3.0.CO;2-U.

  * [31] Cordeiro J.M.M., Freitas L.C.G., Z. Naturforsch. A: Phys. Sci., 1999, 54, 110, doi:10.1515/zna-1999-0204.
  * [32] Cordeiro M.A.M., Santana W.P., Cusinato R., Cordeiro J.M.M., J. Mol. Struct. THEOCHEM, 2006, 759, 159, doi:10.1016/j.theochem.2005.11.016.
  * [33] Radnai T., Bakó I., Jedlovszky P., Pálinkás G., Mol. Simul., 1996, 16, 345, doi:10.1080/08927029608024084.
  * [34] Chalaris M., Samios J., J. Chem. Phys., 2000, 112, 8581, doi:10.1063/1.481460.
  * [35] Chalaris M., Samios J., J. Mol. Liq., 1998, 78, 201, doi:10.1016/S0167-7322(98)00092-0.
  * [36] Chalaris M., Koufou A., Samios J., J. Mol. Liq., 2002, 101, 69, doi:10.1016/S0167-7322(02)00103-4.
  * [37] Zoranić L., Mazighi R., Sokolić F., Perera A., J. Phys. Chem. C, 2007, 111, 15586, doi:10.1021/jp0736894.
  * [38] Zoranić L., Mazighi R., Sokolić F., Perera A., J. Chem. Phys., 2009, 130, 124315, doi:10.1063/1.3093071.
  * [39] Gao J., Pavelites J.J., Habibollazadeh D., J. Phys. Chem., 1996, 100, 2689, doi:10.1021/jp9521969.
  * [40] Biswas S., Mallik B.S., J. Chem. Eng. Data, 2014, 59, 3250, doi:10.1021/je5002544.
  * [41] Vasudevan V., Mushrif S.H., J. Mol. Liq., 2015, 206, 338, doi:10.1016/j.molliq.2015.03.004.
  * [42] Lei Y., Li H., Pan H., Han S., J. Phys. Chem. A, 2003, 107, 1574, doi:10.1021/jp026638+.
  * [43] Jia G.-Z., Huang K.-M., Yang L.-J., Yang X.-Q., Int. J. Mol. Sci., 2009, 10, 1590, doi:10.3390/ijms10041590.
  * [44] Razzokov D., Ismailova O.B., Mamatkulov Sh.I., Trunilina O.V., Kokhkharov A.M., Russ. J. Phys. Chem. A, 2014, 88, 1500, doi:10.1134/S0036024414090271.
  * [45] Fischer N.M., van Maaren P.J., Ditz J.C., Yildirim A., van der Spoel D., J. Chem. Theory Comput., 2015, 11, 2938, doi:10.1021/acs.jctc.5b00190.
  * [46] Galicia-Andrés E., Pusztai L., Temleitner L., Pizio O., J. Mol. Liq., 2015, 209, 586,
doi:10.1016/j.molliq.2015.06.045.

  * [47] Galicia-Andrés E., Dominguez H., Pusztai L., Pizio O., Condens. Matter Phys., 2015, 18, 43602, doi:10.5488/CMP.18.43602.
  * [48] Gujt J., Cázares Vargas E., Pusztai L., Pizio O., J. Mol. Liq., 2017, 228, 71, doi:10.1016/j.molliq.2016.09.024.
  * [49] Gujt J., Dominguez H., Sokolowski S., Pizio O., Condens. Matter Phys., 2017, 20, 33603,
doi:10.5488/CMP.20.33603.

  * [50] Takamuku T., Shimomura T., Tachikawa M., Kanzaki R., Phys. Chem. Chem. Phys., 2011, 13, 11222, doi:10.1039/C0CP00338G.
  * [51] Takamuku T., Yamaguchi A., Matsuo D., Tabata M., Kumamoto M., Nishimoto J., Yoshida K., Yamaguchi T., Nagao M., Otomo T., Adachi T., J. Phys. Chem. B, 2001, 105, 6236, doi:10.1021/jp003011n.
  * [52] Takamuku T., Yamaguchi A., Matsuo D., Tabata M., Yamaguchi T., Otomo T., Adachi T., J. Phys. Chem. B, 2001, 105, 10101, doi:10.1021/jp011692w.
  * [53] Takamuku T., Noguchi Y., Yoshikawa E., Kawaguchi T., Matsugami M., Otomo T., J. Mol. Liq., 2007, 131–132, 131, doi:10.1016/j.molliq.2006.08.048.
  * [54] Van der Spoel D., Lindahl E., Hess B., Groenhof G., Mark A.E., Berendsen H.C., J. Comput. Chem., 2005, 26, 1701, doi:10.1002/jcc.20291.
  * [55] Berendsen H.J.C., Grigera J.R., Straatsma T.P., J. Phys. Chem., 1987, 91, 6269, doi:10.1021/j100308a038.
  * [56] Horn H.W., Swope W.C., Pitera J.W., Madura J.D., Dick T.J., Hura G.L., Head-Gordon T., J. Chem. Phys., 2004, 120, 9665, doi:10.1063/1.1683075.
  * [57] Vega C., Abascal J.L.F., Phys. Chem. Chem. Phys., 2011, 13, 19663, doi:10.1039/C1CP22168J.
  * [58] Wensink E.J.W., Hoffmann A.C., van Maaren P.J., van der Spoel D., J. Chem. Phys., 2003, 119, 7308, doi:10.1063/1.1607918.
  * [59] Alejandre J., Chapela G.A., Saint-Martin H., Mendoza H., Phys. Chem. Chem. Phys., 2011, 13, 19728, doi:10.1039/C1CP20858F.
  * [60] Vega C., de Miguel E., J. Chem. Phys., 2007, 126, 154707, doi:10.1063/1.2715577.
  * [61] Blankenship F., Clampitt B., Proc. Oklahoma Acad. Sci., 1950, 31, 106.

Ukrainian 0 Çàëåæíñòü âä êîíöåíòðàö¿ ìêðîñêîïíî¿ ñòðóêòóðè, òåðìîäèíàìíèõ,
äèíàìíèõ òà åëåêòðèíèõ âëàñòèâîñòåé ìîäåëüíî¿ ñóìø âîäà-äèìåòèë ôîðìàìä.
Ðåçóëüòàòè ñèìóëÿöé ìåòîäîì ìîëåêóëÿðíî¿ äèíàìêè

Ã. Äîìíãåñ, O. Ïçî

Iíñòèòóò ìàòåðiàëîçíàâñòâà, Íàöiîíàëüíèé àâòîíîìíèé óíiâåðñèòåò ì. Ìåõiêî,
Ìåõiêî, Ìåêñèêà


# Active and Passive IRS Jointly Aided Communication: Deployment Design and
Achievable Rate

Min Fu, and Rui Zhang M. Fu is with the Department of Electrical and Computer
Engineering, National University of Singapore, Singapore 117583 (e-mail:
fumin@nus.edu.sg). R. Zhang is with the Chinese University of Hong Kong,
Shenzhen, and Shenzhen Research Institute of Big Data, Shenzhen, China 518172
(e-mail: rzhang@cuhk.edu.cn). He is also with the Department of Electrical and
Computer Engineering, National University of Singapore, Singapore 117583
(e-mail: elezhang@nus.edu.sg).

###### Abstract

In this letter, we study the wireless point-to-point communication from a
transmitter (Tx) to a receiver (Rx), which is jointly aided by an active
intelligent reflecting surface (AIRS) and a passive IRS (PIRS). We consider
two practical transmission schemes by deploying the two IRSs in different
orders, namely, Tx$\rightarrow$PIRS$\rightarrow$AIRS$\rightarrow$Rx (TPAR) and
Tx$\rightarrow$AIRS$\rightarrow$PIRS$\rightarrow$Rx (TAPR). Assuming line-of-
sight channels, we derive the achievable rates for the two schemes by
optimizing the placement of the AIRS with the location of the PIRS fixed. Our
analysis shows that when the number of PIRS elements and/or the AIRS
amplification power is small, the AIRS should be deployed closer to the Rx in
both schemes, and TAPR outperforms TPAR with their respective optimized
AIRS/PIRS placement. Simulation results validate our analysis and show the
considerable performance gain achieved by the jointly optimized AIRS/PIRS
deployment over the baseline double-PIRS system under the same power and IRS
element budgets.

###### Index Terms:

Intelligent reflecting surfaces (IRS), active IRS, double IRSs, IRS
deployment, rate maximization.

##  I Introduction

Intelligent reflecting surface (IRS) [1] has recently received significant
attention from both academia and industry to improve the spectral and energy
efficiency of future wireless networks cost-effectively. An IRS typically
consists of an array of passive reflecting elements, each of which reflects
the incident signal with a desired phase shift and/or amplitude. Unlike
conventional active relays, passive IRSs (PIRSs) do not require costly
transmit/receive radio-frequency (RF) chains and thus incur significantly
lower power consumption [2]. As IRSs can be flexibly deployed in wireless
networks to reconfigure wireless channels dynamically, they have been studied
for achieving various functions such as coverage extension, rate enhancement,
interference mitigation, etc (see, e.g., [1] and the references therein).

Most of the existing works (e.g., [3, 4, 5, 6]) on IRS have considered the
PIRS. Equipped with passive loads (positive resistance), PIRS reflects the
incident signal with desired phase shift and reflection gain no larger than
one. In addition, PIRSs operate in the full-duplex mode without
amplification/processing noise or self-interference [1]. In particular, [3]
showed that deploying a single PIRS with $N$ reflecting elements can lead to a
squared power scaling order (i.e., $\mathcal{O}(N^{2})$) of the reflected
signal, which is even higher than that ($\mathcal{O}(N)$) of the active
arrays. Furthermore, [4] proposed a double-PIRS system where $N$ reflecting
elements are equally allocated over two distributed PIRSs, which yields an
even higher power scaling order ($\mathcal{O}(N^{4})$) as compared to the
single-PIRS system. Recently, wireless systems aided by multiple PIRSs with
multiple signal reflections among them have also been investigated (see [7]
and the references therein). Despite the low cost of PIRSs, their signal
coverage performance is severely limited by the product-distance path-loss of
the cascaded multi-reflection channel[7]. As a result, one PIRS needs to be
allocated massive elements and/or deployed near the transmitter (Tx), the
receiver (Rx), or other PIRSs in practice to make its signal reflection
effective [7].

To address this issue with PIRS, the active IRS (AIRS) has emerged as a
promising solution [8, 9, 10, 11, 12]. Equipped with negative resistance
components connected to an additional power supply, AIRSs enable adjustable
phase shifts like PIRSs as well as reflection amplitude with a gain larger
than one. Compared to PIRSs, although AIRSs offer an amplification gain (that
is typically lower than that of active relays with dedicated RF amplifiers),
they also induce non-negligible amplification noise in their reflected signals
[8]. In particular, given the IRS location and the total power budget, an AIRS
was shown to achieve higher spectral efficiency[8], energy efficiency[11], and
reliability[10] than its PIRS counterpart. Furthermore, with optimized
AIRS/PIRS placements, [9] showed that an AIRS only provides a power scaling of
$\mathcal{O}(N)$ with $N$ active reflecting elements due to the amplification
noise, but it achieves a higher rate than the PIRS when $N$ is small or the
AIRS amplification power is high. Nevertheless, the AIRS system design (e.g.,
AIRS deployment and amplification factor optimization) is more complicated
since it needs to balance signal and noise amplification at the AIRS. The
existing literature on AIRS has mainly compared its performance to PIRS.
However, there has been limited work on exploring their joint use for further
enhancing wireless communication performance.

In this letter, we study a new wireless communication system jointly aided by
a pair of AIRS and PIRS, as illustrated in Fig. 1. To study IRS-enhanced link
performance, we assume a pair of single-antenna Tx and Rx communicating
through a double-reflection link via the two IRSs, while the other links
between them are unavailable due to severe blockages. By deploying the two
IRSs in different orders, we consider two practical transmission schemes,
which are Tx$\rightarrow$PIRS$\rightarrow$AIRS$\rightarrow$Rx (TPAR) and
Tx$\rightarrow$AIRS$\rightarrow$PIRS$\rightarrow$Rx (TAPR) schemes,
respectively (see Fig. 1). Furthermore, we assume that the two IRSs are
properly deployed and thus they form a cascaded line-of-sight (LOS) link from
the Tx to Rx via them successively. Under the above settings, we aim to
characterize the achievable rates for the two schemes with their respective
optimized AIRS/PIRS placement. Our analytic results show that for both
schemes, when the number of PIRS elements is small, the AIRS should be placed
closer to the Rx with decreasing amplification power. Besides, the TAPR scheme
generally outperforms the TPAR scheme when the number of PIRS elements and/or
the AIRS amplification power is small. Simulation results validate our
analysis and demonstrate that the jointly optimized AIRS/PIRS deployment
achieves a much higher rate than the baseline double-PIRS deployment under the
same power and IRS element budgets.

##  II System Model and Problem Formulation

Figure 1: The illustrations of AIRS and PIRS jointly aided wireless
communication.

As illustrated in Fig. 1, a single-antenna Tx transmits data to a single-
antenna Rx that is $D$ meters (m) away. Under a three-dimensional Cartesian
coordinate system, the locations of the Tx and Rx are denoted by
$\bm{u}_{T}=(0,0,0)$ and $\bm{u}_{R}=(D,0,0)$, respectively. A pair of PIRS
and AIRS are deployed between them to assist in their communication. The PIRS
and AIRS are respectively placed at given altitudes, denoted by $H_{P}$ and
$H_{A}$, so that the Tx can communicate with the Rx via a double-reflection
link through the two IRSs over three LOS channels. In particular, to reduce
the product-distance path loss of the double-reflection channel, we assume
that the PIRS is horizontally located above the Tx or Rx while the AIRS can be
flexibly placed between them. We consider two transmission schemes by
deploying the PIRS and AIRS in different orders, i.e., TAPR shown in Fig. 1(a)
and TPAR shown in Fig. 1(b).

Let $N_{p}$ and $N_{a}$ denote the number of passive/active reflecting
elements on the PIRS/AIRS, respectively. We further define $\mathcal{N}_{a}$
and $\mathcal{N}_{p}$ as the sets containing all the elements on the AIRS and
PIRS, respectively. Let
$\bm{\Theta}_{P}=\text{diag}(\eta_{P,1}e^{j\phi_{P,1}},\ldots,\eta_{{P,N_{p}}}e^{j\phi_{P,N_{p}}})$
denote the reflection matrix of the PIRS, where $\eta_{P,n}\in[0,1]$ and
$\phi_{P,n}\in[0,2\pi)$ denote respectively the reflection amplitude and
phase-shift at passive element $n\in\mathcal{N}_{p}$. For maximal reflection
of the PIRS and ease of practical implementation [1], we further set
$\eta_{P,n}=1,\forall n\in\mathcal{N}_{p}$. In addition, let
$\bm{\Phi}_{A}=\text{diag}(\eta_{A,1}e^{j\phi_{A,1}},\ldots,\eta_{{A,N_{a}}}e^{j\phi_{A,N_{a}}})$
denote the reflection matrix of the AIRS, where $\eta_{A,n}\geq 1$ and
$\phi_{A,n}\in[0,2\pi)$ denote respectively the reflection amplitude and
phase-shift at active element $n\in\mathcal{N}_{a}$. Since, under LOS
channels, the transmitted signal arriving at each active element experiences
the same pathloss, we assume that all active elements at the AIRS have the
same amplification gain $\eta$ [9, 8], i.e., $\eta_{A,n}=\eta,\forall
n\in\mathcal{N}_{a}$. Moreover, different from the PIRS that reflects signals
without incurring amplification noise, the AIRS generates non-negligible
amplification noise at all reflecting elements, which is denoted by
$\bm{n}_{\rm a}\in\mathbb{C}^{N_{a}\times 1}$ and assumed to follow the
independent circularly symmetric complex Gaussian distribution, i.e.,
$\bm{n}_{\rm a}\sim\mathcal{CN}(\bm{0}_{N_{a}},\sigma^{2}_{\rm F}{\bf
I}_{N_{a}})$ with power $\sigma^{2}_{\rm F}$.

###  II-A TAPR Transmission Scheme

First, we consider the TAPR scheme with the PIRS fixed above the Rx, as shown
in Fig. 1(a).

####  II-A1 Signal model

Denoting the Tx-AIRS horizontal distance as $x_{TA}$, the locations of the
PIRS/AIRS are represented by $\bm{u}_{P}=(D,0,H_{P})$ and
$\bm{u}_{A}=(x_{TA},0,H_{A})$, respectively. Thus, the distances between the
Rx and PIRS, the Tx and AIRS, and the AIRS and PIRS are respectively given by

$\displaystyle d_{PR}=H_{P},d_{TA}(x_{TA})=\sqrt{x_{TA}^{2}+H_{A}^{2}},$ (1)
$\displaystyle d_{AP}(x_{TA})=\sqrt{(D-x_{TA})^{2}+(H_{P}-H_{A})^{2}},$ (2)

where $x_{TA}$ is considered as the design variable.

Let $\bm{g}_{\rm TA}\in\mathbb{C}^{N_{a}\times 1}$, $\bm{S}_{\rm
AP}\in\mathbb{C}^{N_{p}\times N_{a}}$, and $\bm{h}_{\rm RP}^{\sf
H}\in\mathbb{C}^{1\times N_{p}}$ denote the baseband equivalent LOS channels
for Tx$\rightarrow$AIRS, AIRS$\rightarrow$PIRS, and PIRS$\rightarrow$Rx links,
respectively. Before modeling LOS channels, we first define the one-
dimensional steering vector function for a uniform linear array (ULA) as

$\displaystyle\bm{e}(\theta,M^{\prime})=[1,e^{-j\pi\theta},\ldots,e^{-j\pi(M^{\prime}-1)\theta}]^{\sf
T}\in\mathbb{C}^{M^{\prime}\times 1},$ (3)

where $\theta$ denotes the constant phase-shift difference between the signals
at two adjacent elements and $M^{\prime}$ denotes the size of the ULA. For the
uniform planar array (UPA) model, each array response vector can be expressed
as the Kronecker product of two steering vector functions in the x-axis and
y-axis directions, respectively. For example, the array response vector for
the UPA from the Tx to AIRS is expressed as

$\displaystyle\\!\\!\\!\\!\\!\bm{a}_{TA}(\theta_{TA},\vartheta_{TA},N_{a})=\bm{e}(\frac{2\Delta_{A}}{\lambda}\sin(\theta_{TA})\cos(\vartheta_{TA}),N_{a,x})\otimes$
$\displaystyle\hskip
40.00006pt\bm{e}(\frac{2\Delta_{A}}{\lambda}\sin(\theta_{TA})\sin(\vartheta_{TA}),N_{a,y})\in\mathbb{C}^{N_{a}\times
1},$ (4)

where $\lambda$ denotes the signal wavelength, $\Delta_{A}$ denotes the
element space at the AIRS, $(\theta_{TA},\vartheta_{TA})$ is the angle-of-
arrival pair, and $N_{a,x}$ and $N_{a,y}$ denote the number of
horizontal/vertical elements at the AIRS with $N_{a}=N_{a,x}\times N_{a,y}$.
As such, $\bm{S}_{\rm PA}$, $\bm{g}_{\rm AT}$, and $\bm{h}_{\rm PR}^{\sf H}$
are respectively given by

$\displaystyle\bm{S}_{\rm AP}\\!=\\!\frac{\sqrt{\beta}}{d_{AP}}e^{\frac{-j2\pi
d_{AP}}{\lambda}}\bar{\bm{s}}_{AP}\tilde{\bm{s}}_{PA}^{\sf H},$ (5)
$\displaystyle\bm{g}_{\rm TA}\\!=\\!\frac{\sqrt{\beta}}{d_{TA}}e^{\frac{-j2\pi
d_{TA}}{\lambda}}\bar{\bm{g}}_{TA},\bm{h}_{\rm RP}^{\sf
H}=\frac{\sqrt{\beta}}{H_{P}}e^{\frac{-j2\pi
H_{p}}{\lambda}}\bar{\bm{h}}_{RP}^{\sf H},$ (6)

where $\beta$ denotes the reference channel power at a distance of 1 m,
$\bar{\bm{s}}_{AP}=\bm{a}_{AP}(\theta_{AP},\vartheta_{AP},N_{p})$,
$\tilde{\bm{s}}_{PA}=\bm{a}_{PA}(\theta_{PA},\vartheta_{PA},N_{a})$,
$\bar{\bm{g}}_{TA}=\bm{a}_{TA}(\theta_{TA},\vartheta_{TA},N_{a})$, and
$\bar{\bm{h}}_{RP}=\bm{a}_{RP}(\theta_{RP},\vartheta_{RP},N_{p})$.

Though the Tx$\rightarrow$AIRS$\rightarrow$PIRS$\rightarrow$Rx reflecting
link, the signal received at the Rx is given by

$\displaystyle\\!\\!\\!\\!\\!\\!y_{a}$ $\displaystyle=$
$\displaystyle\bm{h}_{\rm RP}^{\sf{H}}\bm{\Theta}_{P}\bm{S}_{\rm
AP}\eta\bm{\Theta}_{A}\Big{(}\bm{g}_{\rm TA}s+\bm{n}_{\rm a}\Big{)}+e$ (7)
$\displaystyle=$ $\displaystyle\underbrace{\bm{h}_{\rm
RP}^{\sf{H}}\bm{\Theta}_{P}\bm{S}_{\rm AP}\eta\bm{\Theta}_{A}\bm{g}_{\rm
TA}s}_{\text{Double-reflected signal}}+\underbrace{\bm{h}_{\rm
RP}^{\sf{H}}\bm{\Theta}_{P}\bm{S}_{\rm AP}\eta\bm{\Theta}_{A}\bm{n}_{\rm
a}}_{\text{Noise introduced by the AIRS}}+e,$

where $s$ denotes the transmitted signal with power $P_{\rm t}$ of the Tx, and
$e$ denotes the additive white Gaussian noise at the Rx with power
$\sigma^{2}$. As observed from (7), the two IRSs can jointly enhance the
double-reflected signal power whereas the AIRS introduces additional
amplification noise. The received signal-to-noise ratio (SNR) at the Rx is
expressed as

$\displaystyle\\!\\!{\sf
SNR}_{a}(\bm{\Theta}_{A},\bm{\Theta}_{P},\eta,x_{TA})=\frac{\|\bm{h}_{\rm
RP}^{\sf{H}}\bm{\Theta}_{P}\bm{S}_{\rm AP}\eta\bm{\Theta}_{A}\bm{g}_{\rm
TA}\|^{2}P_{\rm t}}{\|\bm{h}_{\rm RP}^{\sf{H}}\bm{\Theta}_{P}\bm{S}_{\rm
AP}\eta\bm{\Theta}_{A}\|^{2}\sigma^{2}_{\rm F}+\sigma^{2}},$ (8)

and the achievable rate (in bits per second per Hertz or bps/Hz) for the TAPR
scheme is given by

$\displaystyle{\sf R}_{a}=\log_{2}\left(1+{\sf
SNR}_{a}(\bm{\Theta}_{A},\bm{\Theta}_{P},\eta,x_{TA})\right).$ (9)

To avoid the impractical case where the AIRS attenuates the signal, the
amplification factor needs to satisfy

$\displaystyle\eta\geq 1.$ (10)

Denoting $P_{\rm F}$ as the given amplification power budget of the AIRS, the
amplification power constraint over the signal reflected by the AIRS is given
by [9, 8]

$\displaystyle\eta^{2}\left(P_{\rm{t}}\|\bm{\Theta}_{A}\bm{g}_{\rm
TA}\|^{2}+\sigma_{\rm F}^{2}\|{\bm{\Theta}_{A}}{\bf
I}_{N_{a}}\|^{2}\right)\leq P_{\rm F}.$ (11)

####  II-A2 Problem formulation

We aim to maximize the achievable rate for the Rx by jointly optimizing the
reflection matrices at the PIRS/AIRS (i.e., $\bm{\Theta}_{P}$ and
$\eta\bm{\Theta}_{A}$) and the AIRS deployment (i.e., $x_{TA}$). Accordingly,
for the TAPR scheme, the optimization problem is formulated as

$\displaystyle\mathop{\text{maximize}}\limits_{\begin{subarray}{c}\bm{\Theta}_{A},\bm{\Theta}_{P},\eta,\\\
x_{TA}\end{subarray}}$ $\displaystyle\log_{2}\left(1+{\sf
SNR_{a}}(\bm{\Theta}_{A},\bm{\Theta}_{P},\eta,x_{TA})\right)$ subject to
$\displaystyle|\bm{\Theta}_{P}(n,n)|=1,\forall n\in{{\mathcal{N}_{p}}},$
$\displaystyle|\bm{\Theta}_{A}(n,n)|=1,\forall n\in{{\mathcal{N}_{a}}},$
$\displaystyle 0\leq x_{TA}\leq D,$
$\displaystyle\text{Constraints~{}}\eqref{Cons:activeIRS1-TAPR},\eqref{Cons:activeIRS2-TAPR}.$

Since the objective function is increasing with $\eta$, the amplification
power constraint in (11) needs to be active at the optimal solution of problem
(12). Furthermore, the amplification power in (11) is independent of the phase
values in $\bm{\Theta}_{A}$ and $\bm{\Theta}_{P}$ while the objective function
of problem (12) is maximized by designing the phases of the AIRS and PIRS to
align in the cascaded Tx-AIRS-PIRS-Rx channel. Therefore, for given $x_{TA}$,
the optimal phase design of the AIRS/PIRS and the optimal amplification factor
are respectively given by

$\displaystyle\bm{\Theta}_{A}(n,n)=e^{-j(\angle[\tilde{\bm{s}}_{PA}^{\sf
H}]_{n}+\angle[\bm{g}_{\rm TA}]_{n})},$ (13)
$\displaystyle\bm{\Theta}_{P}(n,n)=e^{-j(\angle[\bar{\bm{s}}_{AP}]_{n}+\angle[\bm{h}^{\sf
H}_{\rm PR}]_{n})},$ (14) $\displaystyle\eta=\sqrt{P_{\rm
F}d^{2}_{TA}(x_{TA})/(P_{\rm{t}}N_{a}\beta+d^{2}_{TA}(x_{TA})\sigma_{\rm
F}^{2}N_{a}}).$ (15)

With $\eta$, $\bm{\Theta}_{A}$, and $\bm{\Theta}_{P}$ given in the above, the
receiver SNR in (8) is rewritten as

$\displaystyle\\!\\!\\!\\!\\!{\sf SNR}_{a}(x_{TA})$
$\displaystyle\\!\\!\\!\\!\\!=\frac{P_{\rm{t}}P_{\rm
F}\beta^{2}N_{a}N_{p}^{2}}{C_{1}d^{2}_{TA}(x_{TA})+C_{2}d^{2}_{AP}(x_{TA})+C_{3}d^{2}_{TA}(x_{TA})d^{2}_{AP}(x_{TA})},$

where $C_{1}=\sigma_{F}^{2}P_{\rm F}\beta N_{p}^{2}$,
$C_{2}=\sigma^{2}P_{\rm{t}}H_{P}^{2}$, and $C_{3}=H_{P}^{2}\sigma_{\rm
F}^{2}\sigma^{2}/\beta$. After substituting (13), (14), and (15) into (12),
problem (12) is equivalent to the following,

$\displaystyle\mathop{\text{maximize}}\limits_{x_{TA}}$
$\displaystyle\log_{2}(1+{\sf SNR}_{a}(x_{TA}))$ subject to $\displaystyle
x_{a}\leq x_{TA}\leq D,$ (17a)

where constraint (17a) corresponds to (10) in problem (12) and

$\displaystyle x_{a}=\sqrt{\max\\{0,N_{a}\beta P_{\rm{t}}/(P_{\rm
F}-N_{a}\sigma_{\rm F}^{2})-H_{A}^{2}\\}}.$ (18)

Although the optimal solution to problem (17) is difficult to be characterized
in closed form due to the complicated expression of ${\sf SNR}_{a}(x_{TA})$,
its value can be obtained by using a one-dimensional search over
$x_{TA}\in[x_{a},D]$. In the following, we first characterize the effects of
AIRS amplification power and the number of PIRS elements on the optimal AIRS
placement.

######  Lemma 1.

The optimal Tx-AIRS horizontal distance (i.e., $x^{\star}_{TA}$) to problem
(17) is non-increasing with $P_{\rm F}$ and monotonically decreasing with
$N_{p}$.

###### Proof.

First, it can be shown that $x_{a}$ is non-increasing with $P_{\rm F}$. In
addition, it is observed from the first denominator term of ${\sf
SNR}_{a}(x_{TA})$ that $C_{1}$ and $d_{TA}(x_{TA})$ monotonically increases
with both $P_{\rm F}$ and $x_{TA}$. If constraint (17a) is not considered,
$x_{TA}$ should be decreased to maximize the objective function of problem
(17) as $P_{\rm F}$ increases. Thus, the optimal solution $x^{\star}_{TA}$ is
non-increasing with $P_{\rm F}$.

Second, it can be shown that $x_{a}$ is independent of $N_{p}$ while $C_{1}$
monotonically increases with $N_{p}$. If without constraint (17a), $x_{TA}$
should be decreased to maximize the objective function of problem (17) as
$N_{p}$ increases. Thus, the optimal solution $x^{\star}_{TA}$ is decreasing
with $N_{p}$.

This thus completes the proof. ∎

Lemma 1 is expected for the TAPR scheme in order to balance the received
signal power maximization and the AIRS-induced noise minimization. On one
hand, as observed from (15), the AIRS should be deployed farther away from the
Tx as $P_{\rm F}$ decreases, so that the AIRS can provide a higher
amplification factor with $\eta>1$ while the path loss of the cascaded channel
can be reduced to maximize the received signal power. On the other hand, based
on (8), although the received signal power increases with the number of PIRS
reflecting elements (i.e., $N_{p}$), the AIRS-induced noise is also increased.
As such, as $N_{p}$ increases, the AIRS has to be deployed closer to the Tx to
attenuate the AIRS-PIRS channel for suppressing its induced noise power at the
Rx.

###  II-B TPAR Transmission Scheme

Next, we consider the TPAR transmission scheme with the PIRS fixed above the
Tx, as shown in Fig. 1(b).

Similarly to the TAPR case, denoting the horizontal distance between the Tx
and AIRS as $x_{TA}$, the locations of the PIRS/AIRS are represented by
$\bm{u}_{P}=(0,0,H_{P})$ and $\bm{u}_{A}=(x_{TA},0,H_{A})$, respectively. The
signal model of TPAR is also similar to that of TAPR, while the expressions of
the AIRS amplification noise and its power constraint are different since the
order of the AIRS/PIRS placements is reversed. The detailed expressions are
omitted for brevity.

Similarly to the design of $\eta$, $\bm{\Theta}_{A}$, and $\bm{\Theta}_{P}$ in
the TAPR scheme, the expressions of the AIRS amplification factor and receiver
SNR in the case of TPAR are respectively given by

$\displaystyle\eta=\sqrt{P_{\rm
F}d^{2}_{AP}/(P_{\rm{t}}N_{a}\beta^{2}N_{p}^{2}/H_{P}^{2}+d^{2}_{AP}\sigma_{\rm
F}^{2}N_{a})},$ $\displaystyle\\!\\!\\!\\!\\!{\sf SNR}_{b}(x_{TA})$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\frac{N_{a}\beta^{2}N_{p}^{2}P_{\rm t}P_{\rm
F}}{C_{2}/ad^{2}_{AP}(x_{TA})\\!+\\!aC_{1}d^{2}_{AR}(x_{TA})\\!+\\!C_{3}d^{2}_{AP}(x_{TA})d^{2}_{AR}(x_{TA})},$

where $a=\frac{P_{t}\sigma^{2}}{P_{\rm F}\sigma_{\rm F}^{2}}$. As a result,
the optimization problem of the TPAR scheme is formulated as

$\displaystyle\mathop{\text{maximize}}\limits_{\ x_{TA}}$
$\displaystyle\log_{2}(1+{\sf SNR}_{b}(x_{TA}))$ subject to $\displaystyle
x_{b}\leq x_{TA}\leq D,$ (21a)

where constraint (21a) is due to the amplification power constraint at the
AIRS, and

$\displaystyle\\!x_{b}\\!=\\!\sqrt{\\!\max\\{0,N_{a}N_{p}^{2}\beta^{2}P_{\rm{t}}/((P_{\rm
F}\\!-\\!N_{a}\sigma_{\rm
F}^{2})H_{P}^{2})\\!-\\!(H_{P}\\!-\\!H_{A})^{2}\\}}.$

Similarly to problem (17), the optimal solution of problem (21) can also be
obtained by using the one-dimensional search over $x_{TA}\in[x_{b},D]$.
Likewise, we have the following lemma for the TPAR scheme.

######  Lemma 2.

The optimal Tx-AIRS horizontal distance (i.e., $x^{\star}_{TA}$) to problem
(21) is non-increasing with $P_{\rm F}$ and monotonically increasing with
$N_{p}$.

###### Proof.

The proof is similar to that of Lemma 1, and thus omitted for brevity. ∎

Lemma 2 can be intuitively understood for the TPAR scheme. The AIRS-induced
noise is independent of the number of PIRS reflecting elements. As observed
from (II-B), with a smaller $P_{\rm F}$ and/or a larger $N_{p}$, the AIRS
should be deployed closer to the Rx to provide a higher amplification factor
with $\eta>1$ and reduce the path loss of the cascaded channel, thus
increasing the signal power at the Rx.

##  III Low-Complexity Placement Design and Performance Comparison

In this section, we present suboptimal AIRS deployment solutions in closed-
form for the TAPR and TPAR schemes, respectively. Based on these results, we
compare the performance of the two schemes in terms of achievable rate.

###  III-A Suboptimal Solutions

First, to simplify problems (17) and (21), we present a useful lemma as
follows (for which the proof is omitted for brevity).

######  Lemma 3.

Given $D\gg\max\\{H_{A},H_{P}\\}$, we have
$C_{1}d^{2}_{TA}(x_{TA})+C_{2}d^{2}_{AP}(x_{TA})\gg
C_{3}d^{2}_{TA}(x_{TA})d^{2}_{AP}(x_{TA})$ and
$C_{2}/ad^{2}_{AP}(x_{TA})+aC_{1}d^{2}_{AR}(x_{TA})\gg
C_{3}d^{2}_{AP}(x_{TA})d^{2}_{AR}(x_{TA})$, $\forall x_{TA}\in[0,D]$, if

$\displaystyle\min\\{\frac{\sqrt{\beta}N_{p}\sqrt{P_{\rm
t}\beta}}{H_{p}\sigma_{\rm F}}\\!+\\!\frac{\sqrt{P_{\rm
F}\beta}}{\sigma},\frac{\sqrt{\beta}N_{p}\sqrt{P_{\rm
F}\beta}}{H_{p}\sigma}\\!+\\!\frac{\sqrt{P_{\rm t}\beta}}{\sigma_{\rm
F}}\\}\gg D.$

####  III-A1 TAPR Scheme

When condition (3) is satisfied, $\sf{R}_{a}$ in problem (17) can be
approximated as $\overline{\sf{R}}_{a}=\log_{2}(1+\overline{\sf{SNR}}_{a})$,
where
$\overline{\sf{SNR}}_{a}=\frac{P_{\rm{t}}P_{\rm{F}}\beta^{2}N_{a}N_{p}^{2}}{C_{1}x^{2}_{TA}+C_{2}(D-x_{TA})^{2}}$.
Therefore, problem (17) is approximated as

$\displaystyle\mathop{\text{maximize}}\limits_{x_{TA}}$
$\displaystyle\log_{2}\Big{(}1+\frac{P_{\rm{t}}P_{\rm{F}}\beta^{2}N_{a}N_{p}^{2}}{C_{1}x^{2}_{TA}+C_{2}(D-x_{TA})^{2}}\Big{)}$
subject to $\displaystyle\text{constraint~{}}\eqref{Cons:deployment1}.$ (23)

By setting the first derivative of the objective function of (III-A1) with
respect to $x_{TA}$ to zero and then projecting it into the feasible set of
problem (III-A1), the optimal solution to problem (III-A1) can be obtained as

$\displaystyle x^{a}_{TA}=\max\\{\frac{aD}{a+b},x_{a}\\},$ (24)

where $a=\frac{P_{t}\sigma^{2}}{P_{\rm F}\sigma_{\rm F}^{2}}$ and
$b=\frac{\beta N_{p}^{2}}{H_{P}^{2}}$. Note that when $\frac{aD}{a+b}>x_{a}$,
the amplification factor (i.e., $\eta$) in problem (12) is larger than 1.
Given $\eta>1$, the receiver SNR is approximated as

$\displaystyle\overline{\sf{SNR}}_{a}=\frac{P_{\rm F}\beta
N_{a}}{D^{2}\sigma^{2}}(a+b)=\frac{\beta N_{a}}{D^{2}}(\frac{P_{\rm
t}}{\sigma_{\rm F}^{2}}\\!+\\!\frac{P_{\rm
F}N_{p}^{2}\beta}{\sigma^{2}H_{p}}).$ (25)

####  III-A2 TPAR Scheme

Similarly, based on Lemma 3, problem (21) is approximated as

$\displaystyle\mathop{\text{maximize}}\limits_{x_{TA}}$
$\displaystyle\log_{2}\Big{(}1+\frac{P_{\rm t}P_{\rm
F}\beta^{2}N_{a}N_{p}^{2}}{1/aC_{2}x_{TA}^{2}+aC_{1}(D-x_{TA})^{2}}\Big{)}$
subject to $\displaystyle\text{constraint~{}}\eqref{Cons:TPAR-deployment}.$
(26)

Similarly to problem (III-A1), the optimal solution to problem (III-A2) is
obtained as

$\displaystyle x^{b}_{TA}=\max\\{\frac{abD}{ab+1},x_{b}\\}.$ (27)

Given $\eta>1$, the approximated receiver SNR is given by

$\displaystyle\overline{\sf{SNR}}_{b}=\frac{P_{\rm F}\beta
N_{a}}{D^{2}\sigma^{2}}(ab+1)=\frac{\beta N_{a}}{D^{2}}(\frac{P_{\rm
t}N_{p}^{2}\beta}{\sigma_{\rm F}^{2}H_{p}}+\frac{P_{\rm F}}{\sigma^{2}}).$
(28)

###  III-B Comparisons between TPAR and TAPR

With the suboptimal AIRS deployments given in (24) and (27), the achievable
rates for the TPAR and TAPR schemes are compared as follows.

Figure 2: Optimal AIRS placement versus amplification power of the AIRS.

Figure 3: Achievable rate versus amplification power of the AIRS.

######  Proposition 1.

Given $\eta>1$ and $N_{p}\leq\frac{H_{P}}{\sqrt{\beta}}$, the achievable rate
of the TPAR scheme is no less than that of the TAPR scheme, i.e.,
$\overline{\sf{R}}_{a}\geq\overline{\sf{R}}_{b}$, if

$\displaystyle P_{\rm F}\leq\frac{P_{\rm t}\sigma^{2}}{\sigma_{\rm F}^{2}}.$
(29)

Otherwise, $\overline{\sf{R}}_{a}<\overline{\sf{R}}_{b}$. The result is
reversed when $N_{p}\geq\frac{H_{P}}{\sqrt{\beta}}$.

###### Proof.

Comparing $\overline{\sf{R}}_{a}$ and $\overline{\sf{R}}_{b}$ is equivalent to
comparing $\overline{\sf{SNR}}_{a}$ and $\overline{\sf{SNR}}_{b}$. Thus, with
$\eta>1$, based on (25) and (28), we have

$\displaystyle\overline{\sf{SNR}}_{a}-\overline{\sf{SNR}}_{b}$
$\displaystyle=$ $\displaystyle\frac{P_{\rm F}\beta
N_{a}}{D^{2}\sigma^{2}}(a+b-ab-1)$ (30) $\displaystyle=$
$\displaystyle\frac{\beta N_{a}}{D^{2}}\frac{(\frac{P_{\rm t}}{\sigma_{\rm
F}^{2}}-\frac{P_{\rm F}}{\sigma^{2}})(H_{p}^{2}-N_{p}^{2}\beta)}{H_{p}^{2}}.$

Based on (30), we can conclude that
$\overline{\sf{R}}_{a}\geq\overline{\sf{R}}_{b}$, if

$\displaystyle N_{p}\leq\frac{H_{P}}{\sqrt{\beta}},\frac{P_{\rm
F}}{\sigma^{2}}\leq\frac{P_{\rm t}}{\sigma_{\rm
F}^{2}},\text{or~{}}N_{p}\geq\frac{H_{P}}{\sqrt{\beta}},\frac{P_{\rm
F}}{\sigma^{2}}\geq\frac{P_{\rm t}}{\sigma_{\rm F}^{2}}.$ (31)

Otherwise, $\overline{\sf{R}}_{a}<\overline{\sf{R}}_{b}$. This thus completes
the proof. ∎

##  IV Numerical Results

In this section, we present numerical results to compare the rate performance
of the TAPR and TPAR schemes in an AIRS-PIRS jointly aided communication
system. If not specified otherwise, the simulation parameters are set as
follows. The AIRS and PIRS are assumed to be deployed at altitudes of
$H_{A}=6$ m and $H_{P}=5$ m, respectively. The carrier frequency is 3.5 GHz
and thus the reference channel gain is $\beta=(\lambda/4\pi)^{2}=-43$ dB, with
the carrier wavelength $\lambda=0.087$ m. Other parameters are set as $D=30$
m, $P_{\rm t}=20$ dBm, $\sigma^{2}=-80$ dBm, and $\sigma_{\rm
F}^{2}=4\sigma^{2}=-74$ dBm. Under our simulation setup, we consider a limited
budget of reflecting elements with
$N_{a}<N_{p}<\frac{H_{P}}{\sqrt{\beta}}\approx 700$ (see Proposition 1). The
numbers of AIRS and PIRS elements are thus fixed as $N_{a}=450$ and
$N_{p}=600$, respectively.

For performance comparison, we consider a benchmark case with double PIRSs
[4]. Specifically, both PIRSs are equally equipped with $(N_{a}+N_{p})/2$
elements and deployed above the Tx and Rx, respectively. Furthermore, the
transmit power of Tx in the double-PIRS system is $(P_{\rm t}+P_{\rm F})$. As
such, it is fair to compare it with the AIRS-PIRS jointly aided system. For
the latter, we also consider four fixed AIRS/PIRS deployments, where $x_{AP}$
equals $D/2$ or $D$ in both TAPR/TPAR schemes.

Fig. 3 shows the optimal Tx-AIRS horizontal distance versus AIRS amplification
power. We observe that the optimal AIRS placement in both schemes is closer to
the Rx (i.e., larger $x_{TA}$) as $P_{\rm F}$ decreases. This is consistent
with Lemmas 1 and 2. Besides, the suboptimal deployment in both schemes is
almost the same as the optimal one.

Fig. 3 compares the achievable rates versus AIRS amplification power. First,
as $P_{F}$ increases, we observe that the rate performance of all AIRS/PIRS
deployment strategies in the AIRS-PIRS jointly aided system increases due to
the more amplification gain from AIRS. However, the rate improvement for the
Rx-side AIRS deployment in TPAR is marginal with $P_{F}$ since the signal
amplification at the AIRS is gradually counteracted by its induced
amplification noise. Compared to these fixed deployment strategies, the AIRS
deployment design in both TAPR/TPAR schemes dynamically balances the signal
and noise amplification at the AIRS and thus achieves further performance
improvement. Furthermore, with the optimal AIRS placement, TAPR performs
better than TPAR when $P_{\rm F}<\frac{P_{\rm t}\sigma^{2}}{\sigma_{\rm
F}^{2}}$, because it suffers less path loss of the cascaded channel and thus
results in higher received power. These results are consistent with
Proposition 1. In addition, thanks to the optimal AIRS/PIRS placement and the
resulting amplification gain at AIRS, the AIRS-PIRS jointly aided system
achieves a much higher rate than the baseline double-PIRS system under a
limited IRS element budget.

##  V Conclusion

In this letter, we studied a new AIRS-PIRS jointly aided wireless
communication system. Under the LOS channel model, we analyzed the effects of
the AIRS deployment in two transmission schemes (i.e., TAPR and TPAR) on their
achievable rates. Our analysis revealed that the AIRS should be deployed
closer to the Rx in both schemes, and TAPR achieves a higher rate than TPAR
when the number of PIRS elements and/or AIRS amplification power is limited.
Simulation results validated our analysis and demonstrated that our proposed
system considerably outperforms the double-PIRS system in terms of achievable
rate, and the performance gain highly depends on the AIRS/PIRS deployment.

## References

  * [1] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface-aided wireless communications: A tutorial,” _IEEE Trans. Commun._ , vol. 69, no. 5, pp. 3313–3351, May 2021.
  * [2] D. Renzo _et al._ , “Reconfigurable intelligent surfaces vs. relaying: Differences, similarities, and performance comparison,” _IEEE Open J. Commun. Soc._ , vol. 1, pp. 798–807, Jun. 2020.
  * [3] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, p. 5394–5409, Aug. 2019.
  * [4] Y. Han, S. Zhang, L. Duan, and R. Zhang, “Double-IRS aided MIMO communication under LoS channels: Capacity maximization and scaling,” _IEEE Trans. Commun._ , vol. 70, no. 4, pp. 2820–2837, Apr. 2022.
  * [5] B. Zheng, C. You, and R. Zhang, “Double-IRS assisted multi-user MIMO: Cooperative passive beamforming design,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 7, pp. 4513–4526, Jul. 2021.
  * [6] M. Fu, Y. Zhou, Y. Shi, and K. B. Letaief, “Reconfigurable intelligent surface empowered downlink non-orthogonal multiple access,” _IEEE Trans. Commun._ , vol. 69, no. 6, pp. 3802–3817, Jun. 2021.
  * [7] W. Mei, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface-aided wireless networks: From single-reflection to multireflection design and optimization,” _Proc. IEEE_ , 2022, doi:10.1109/JPROC.2022.3170656.
  * [8] R. Long, Y.-C. Liang, Y. Pei, and E. G. Larsson, “Active reconfigurable intelligent surface-aided wireless communications,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 8, pp. 4962–4975, Aug. 2021.
  * [9] C. You and R. Zhang, “Wireless communication aided by intelligent reflecting surface: Active or passive?” _IEEE Wireless Commun. Lett._ , vol. 10, no. 12, pp. 2659–2663, Dec. 2021.
  * [10] M. H. Khoshafa, T. M. N. Ngatched, M. H. Ahmed, and A. R. Ndjiongue, “Active reconfigurable intelligent surfaces-aided wireless communication system,” _IEEE Commun. Lett._ , vol. 25, no. 11, pp. 3699–3703, Nov. 2021.
  * [11] K. Liu, Z. Zhang, L. Dai, S. Xu, and F. Yang, “Active reconfigurable intelligent surface: Fully-connected or sub-connected?” _IEEE Commun. Lett._ , vol. 26, no. 1, pp. 167–171, Jan. 2022.
  * [12] K. Zhi, C. Pan, H. Ren, K. K. Chai, and M. Elkashlan, “Active RIS versus passive RIS: Which is superior with the same power budget?” _IEEE Commun. Lett._ , vol. 26, no. 5, pp. 1150–1154, May 2022.