# Voltage-dependent first-principles simulation of insertion of chloride ions
into Al/Al2O3 interfaces using the Quantum Continuum Approximation

Quinn Campbell  Sandia National Laboratories, Albuquerque, NM USA
<EMAIL_ADDRESS>
###### Abstract

Experiments have shown that pitting corrosion can develop in aluminum surfaces
at potentials $>-0.5$ V relative to the standard hydrogen electrode (SHE).
Until recently, the onset of pitting corrosion in aluminum has not been
rigorously explored at an atomistic scale because of the difficulty of
incorporating a voltage into density functional theory (DFT) calculations. We
introduce the Quantum Continuum Approximation (QCA) which self-consistently
couples explicit DFT calculations of the metal-insulator and insulator-
solution interfaces to continuum Poisson-Boltzmann electrostatic distributions
describing the bulk of the insulating region. By decreasing the number of
atoms necessary to explicitly simulate with DFT by an order of magnitude, QCA
makes the first-principles prediction of the voltage of realistic
electrochemical interfaces feasible. After developing this technique, we apply
QCA to predict the formation energy of chloride atoms inserting into oxygen
vacancies in Al(111)/$\alpha$-Al2O3 (0001) interfaces as a function of applied
voltage. We predict that chloride insertion is only favorable in systems with
a grain boundary in the Al2O3 for voltages $>-0.2$ V (SHE). Our results
roughly agree with the experimentally demonstrated onset of corrosion,
demonstrating QCA’s utility in modeling realistic electrochemical systems at
reasonable computational cost.

  *

## Introduction

Corrosion of metallic interfaces leads to weakening of their structural
integrity and has been widely estimated to cost billions of dollars
annually.[1] While this has inspired numerous studies on corrosion and
strategies to protect against it, the fundamental atomistic mechanisms behind
corrosion are difficult to identify for any given system. Corrosion in
aluminum provides an excellent example of the issues. In aluminum, it has been
experimentally demonstrated that pitting corrosion begins at potentials
$\Phi>-0.50$ V relative to the standard hydrogen electrode (SHE), dependent on
the salt and pH concentration of the surrounding solution.[2, 3] The exact
mechanism that causes this pitting corrosion is still under debate, however.
Pitting corrosion is often described using the point defect model,[4, 5, 6]
which states that corrosion is initially driven by the formation of charged
defects in the surface of the material, of which chloride ions are
particularly damaging.[7, 8] On aluminum surfaces, a several nm thick aluminum
oxide layer forms [9] and the insertion of chloride into oxygen vacancies in
this layer has been theorized as one of the main mechanisms leading to
corrosion. Experiments have investigated the dependence of pitting corrosion
on chlorine concentration in the surrounding electrolyte and other
factors,[10, 11, 12] but the results have not been fully conclusive as to the
causal mechanism.

While many previous density functional theory (DFT) studies have examined
Al/Al2O3 interfaces at an atomic level,[13, 14, 15, 16, 17] they have
typically not included oxygen vacancies and chloride insertion. For those that
have examined vacancies,[18, 19, 20, 21] they do not include Al metal in the
calculation. Previous DFT studies on chloride insertion into Al/Al2O3
interfaces that have included contributions from the metal, the oxide,
vacancies, and chloride ions [22, 23, 24] have either implicitly limited their
scope to defects with no charge, or tested only a small range of voltages.

This lack of study of charged ions in the full electrochemical environment
stems from the difficulty of needing to include several different length
scales in a simulation to get a realistic picture of corrosion behavior.
Modeling the realistic interaction of point defects with the surface requires
a quantum mechanical treatment, with DFT providing a convenient compromise
between accuracy and speed. DFT is typically limited, however, to simulating
$\sim$1000s of atoms as larger systems are computationally prohibitive. Since
a realistic model of an interface, including the surrounding solution, can be
10s to 100s of nm wide, DFT simulations must necessarily focus on a single
region of interest, typically the surface. The entire surrounding electrode
helps determine the electric field (and therefore voltage) that the region of
interest operates at. Simulations of these isolated regions therefore do not
provide information as to the voltage the electrode is under and, in effect,
operate at “flatband” conditions where the electric field is flat across the
system. This necessarily limits the scope of traditional DFT calculations when
examining these interfaces and defects. Since the defects under examination
are charged, it is crucial to incorporate a realistic model of voltage into
DFT calculations.

Figure 1: A schematic illustration of the Quantum Continuum Approximation
applied to a metal-insulator interface immersed in solution.

In this paper, we develop a Quantum Continuum Approximation (QCA) methodology
that allows for the first-principles simulation of the metal-insulator
interfaces at controlled voltages using limited DFT simulation cell sizes. QCA
self-consistently couples DFT calculations of the metal-insulator interface
and the insulator-solution interface with Poisson-Boltzmann distribution of
charges in the bulk insulator and solution to provide an equilibrium potential
drop for the entire electrode, as schematically illustrated in Fig. 1. By
reducing the number of atoms needed for DFT calculations to include voltage
information by at least an order of magnitude, QCA enables voltage-dependent
first-principles simulation of realistic electrochemical interfaces at
reasonable computational cost. We demonstrate this methodology for the
simulation of chloride insertion into the aluminum oxide layer of an
Al(111)/$\alpha$-Al2O3(0001) system as a function of voltage. We predict that
insertion is not favorable for Al/Al2O3 interfaces that do not include grain
boundaries at any voltage. We do see, however, that in systems with grain
boundaries in the Al2O3, chloride insertion is favorable at $\Phi>-0.2$ V
(SHE), in reasonable agreement with experiment and helping confirm the
importance of grain boundaries in the corrosion process. This work provides a
template for attaching a potentiostat to electrochemical DFT calculations of
interfaces.

## Methods

### Chloride insertion energy calculations

We calculate the formation energy of a Cl- ion inserting into an oxygen
vacancy as

$\Delta E_{\text{Cl-O}}[\Phi]=E_{\text{DFT}}(\text{Cl
ins.})[\Phi]-E_{\text{DFT}}(\text{O vac.})[\Phi]-E_{\text{Cl}^{-}}[\Phi],$ (1)

where $E_{\text{DFT}}(\text{O vac.})[\Phi]$ is the DFT energy of the slab with
the oxygen vacancy as a function of the global potential $\Phi$ of the system,
$E_{\text{DFT}}(\text{Cl ins.})[\Phi]$ is the DFT energy of the slab with a
chlorine ion inserted into the previous location of the oxygen vacancy as a
function of the global potential, and $E_{\text{Cl}^{-}}[\Phi]$ is the energy
of the Cl- ion in solution as a function of voltage.

For the chloride ion, we follow a computational hydrogen electrode-like
approach,[25, 26, 27] where we determine the energy of the chloride ion as

$E_{\text{Cl}^{-}}[\Phi]=\mu(\text{Cl}^{-})+e\Phi^{\text{SHE}}.$ (2)

Here, $\mu(\text{Cl}^{-})$ is the equilibrium chemical potential of the
chloride ion, and $\Phi^{\text{SHE}}$ is the potential that the electrode is
operating at, referenced to the standard hydrogen electrode (SHE). We
determine the chemical potential by referencing the following chemical
reaction

$\text{Cl}_{2}(\text{g})+2e^{-}\rightarrow 2\text{Cl}^{-}(\text{aq}),$ (3)

which is at equilibrium at $\Phi_{\text{Cl}_{2}/\text{Cl}^{-}}=1.36$ V
(SHE).[28] Using this reaction, we can then calculate the chemical potential
of a chloride ion as

$\mu(\text{Cl}^{-})=\frac{1}{2}E_{\text{DFT}}(\text{Cl}_{2}(\text{g}))-e\Phi_{\text{Cl}_{2}/\text{Cl}^{-}}.$
(4)

This energy calculation is inherently dependent on the calculated dissociation
energy of the Cl2 molecule, which we calculate as 3.46 eV with PBE, compared
to a value of 2.47 eV from experiment.[29] However, to maintain consistency
with the PBE calculations of Cl insertion throughout the work, we have not
applied any corrections to the calculated PBE energies as this would interfere
with the same overbinding tendency when using PBE to calculate the energy of
inserting Cl into a slab. With this energetic relationship established, we can
now calculate the formation energy of chloride insertion into the Al/Al2O3
surface as a function of the electrode potential. For this relationship to be
strictly accurate, however, the calculated DFT energies of the slabs need to
be taken at simulation conditions consistent with $\Phi$ for the electrode. We
outline how we can achieve this using the QCA methodology in the next section.

### Quantum Continuum Approximation (QCA) Methodology for Predicting Voltage
at an interface

#### Interface Energetics and Electrostatics

At a fundamental level, our approach to calculations with metal-semiconductor-
solution interfaces works by assigning a global potential to the system. This
potential drop can then be split up between several different voltage drop
regions within the interface. Finding the correct distribution of potential
drops requires a self-consistent process for minimizing the free energy of the
system while taking into account the electrostatic requirements of following
the Poisson equation and maintaining electronic equilibrium throughout the
interface.

To accomplish this task, we employ a new extension of the Quantum Continuum
Approximation (QCA) methodology that was first developed by Campbell and Dabo
for semiconductor-solution interfaces,[30, 31] which allowed for the voltage
of the system to be additionally predicted by DFT slab calculations. It was
later extended to solid-state metal-semiconductor interfaces, providing
predictions of Schottky barriers that corresponded well with nano-ARPES
measurements.[32] Here, we extend this to include metal-insulator or metal-
semiconductor interfaces that are also immersed in solution, which are common
in electrochemical applications. The basic goal is to self-consistently couple
DFT calculations of an interface with analytic Poisson-Boltzmann descriptions
of the bulk semiconductor and solution regions. As shown in Fig. 1, we perform
DFT calculations only on the crucial interface sections, such as the metal-
metal oxide transition and the metal oxide-solution transition, or around
atomic defects in the system. We then apply Poisson-Boltzmann calculations to
the oxide regions in between the DFT calculations and the surrounding
solution.

The free energy of the system is

$\mathcal{F}=\mathcal{F}^{\text{m}|\text{ox.}}+\mathcal{F}^{\text{ox.}}+\mathcal{F}^{\text{ox.}|\text{sol.}}+\mathcal{F}^{\text{sol.}}-\frac{1}{2}\int
dr\epsilon_{0}\epsilon(\mathbf{r})|\nabla(\phi(\mathbf{r}))|^{2},$ (5)

where $\epsilon_{0}$ is the vacuum permittivity, $\epsilon(\mathbf{r})$ is the
dielectric constant as a function of position, $\phi$ is the electrostatic
potential, and $\mathcal{F}^{\text{m}|\text{ox.}}$ and
$\mathcal{F}^{\text{ox.}|\text{sol.}}$ are the free energies of the metal-
oxide and oxide-solution interface, respectively, which can be calculated
using the typical calculations for DFT energy. Similarly,
$\mathcal{F}^{\text{ox.}}$ and $\mathcal{F}^{\text{sol.}}$ are the
contributions of the bulk oxide and solution regions, which can be calculated
as analytic functions of the distribution of charge carriers within each
system. This decomposition of the free energy naturally leads to the following
Poisson equation for the system:

$\nabla(\epsilon_{0}\epsilon(\mathbf{r})\nabla\phi(\mathbf{r}))=p(\mathbf{r})-n(\mathbf{r})-p_{d}(\mathbf{r})+n_{d}(\mathbf{r})+c_{+}(\mathbf{r})-c_{-}(\mathbf{r})+\rho_{+}(\mathbf{r})-\rho_{-}(\mathbf{r}),$
(6)

where and $\rho_{-}(\mathbf{r})$ is the electron density obtained from the
self-consistent Kohn-Sham equation for a given distribution
$\rho_{+}(\mathbf{r})$ of atomic cores, and the source terms can be expressed
as

$n(\mathbf{r})=n_{d}(\mathbf{r})\left[1+\exp\left(\frac{\phi(\mathbf{r})+\epsilon_{c}-\epsilon_{F}}{k_{\text{B}}T}\right)\right]^{-1},$
(7)
$p(\mathbf{r})=p_{d}(\mathbf{r})\left[1+\exp\left(\frac{\epsilon_{F}-\epsilon_{v}-\phi(\mathbf{r})}{k_{\text{B}}T}\right)\right]^{-1},$
(8)
$c_{\pm}(\mathbf{r})=c^{\circ}(\mathbf{r})\exp\left(\mp\frac{\phi(\mathbf{r})}{k_{\text{B}}T}\right).$
(9)

Here, $\epsilon_{F}$ denotes the Fermi energy, $\epsilon_{v}$ is the energy at
the valence band maximum, $\epsilon_{c}$ is the energy at the conduction band
minimum, $n_{d}(\mathbf{r})$ is the concentration of $n$-type carriers in the
semiconducting/insulating region and is 0 outside this region,
$p_{d}(\mathbf{r})$ is similarly the concentration of $p$-type carriers within
the semiconducting regions, and $c^{\circ}(\mathbf{r})$ is the equilibrium
ionic concentration inside the electrolyte region and is 0 outside of the
electrolyte region. In an $n$-type semiconducting region, the charge carriers
can be solved for using the Boltzmann distribution as

$p(\mathbf{r})=p_{d}=0,$ (10)
$n(\mathbf{r})=n_{d}\exp\left(\frac{\phi_{0}-\phi(\mathbf{r})}{k_{\text{B}}T}\right),$
(11)

where $\phi_{0}$ is the asymptotic value of the potential in this region.
Analogous equations can be written for $p$-type semiconducting regions.

#### Process for solving electrostatics

To practically solve Eq. 6 for a given potential, we start by dividing the
system into two portions that we will simulate with DFT separately. The first
is the surface of the interface, i.e. where the passivating oxide layer meets
the solution. We will refer to this region with a superscript ox.$|$sol. in
variables throughout this work. The other DFT region is the interface of the
metal with the passivating oxide, referred to as a superscript m$|$ox. within
this work. For both of these regions we start by doing DFT calculations on
these systems when no charge has been applied. This corresponds to the so-
called ‘flatband’ conditions. We then measure both the planar averaged
potential profile of these systems $\phi_{\text{fb}}^{\text{m}|\text{ox.}}(z)$
and $\phi_{\text{fb}}^{\text{ox.}|\text{sol.}}(z)$, and the Fermi level of
these systems $\varepsilon_{\text{F,fb}}^{\text{m}|\text{ox.}}$ and
$\varepsilon_{\text{F,fb}}^{\text{ox.}|\text{sol.}}$.

Next we assign a total charge $Q$ to the electrode. The electrode is here
defined as the region from the metal to the interface with the solution. The
distribution of this charge $Q$ within the interface mechanically determines
the total potential drop via Gauss’ law. For a perfectly
insulating/semiconducting system, the charge distribution can be well
approximated by the following Poisson–Boltzmann relationship:[33]

$\frac{d^{2}\phi}{dz^{2}}=\frac{n_{\text{d}}}{\epsilon_{o}\epsilon_{\text{sc}}}\left[1-\exp\left(\frac{-\phi}{k_{\text{B}}T}\right)\right].$
(12)

Here, $n_{\text{d}}$ is the concentration of charge carriers (dopants) within
the semiconductor, $\epsilon_{o}$ is the vacuum permittivity,
$\epsilon_{\text{sc}}$ is the dielectric constant of the semiconductor, and
$k_{\text{B}}$ is the Boltzmann constant.

For most real-world systems, however, surface states will form at both the
metal-insulator and insulator-solution interface.[34, 35, 36] These surface
states result from the abrupt nature of the interface disrupting the symmetry
of the material, and from defects induced by interaction. These surface states
trap charge at the interface, thus reducing the charge throughout the
remainder of the insulator and typically reducing the potential drop of the
electrode. Capturing these surface states requires a quantum mechanical
treatment of the interfaces using, e.g., DFT.

Figure 2: A schematic illustration of the potential profile of a metal-oxide-
solution interface within QCA with important variables depicted. Distances and
potentials shown are approximate.

To simulate charge accumulation at these surface interfaces we place Helmholtz
planar counter-charges on either side of a DFT slab, as illustrated in Fig. 2.
We place the Helmholtz planes 4 Å away from the surfaces of the slab to mimic
the Debye length, which is estimated between 3-5 Å for typical ion
concentrations in water. Previous work has found the exact distance of the
Helmholtz plane in this range to not significantly affect results.[37] We obey
charge neutrality for the DFT systems simulated, such that

$q_{\text{h}}^{\text{ox.}|\text{sol.}}+q_{\text{DFT}}^{\text{ox.}|\text{sol.}}+q_{\text{h}}^{\text{sol.}}=0$
(13)

at the oxide-solution interface, where $q_{\text{h}}^{\text{ox.}|\text{sol.}}$
and $q_{\text{h}}^{\text{sol.}}$ is the charge on the Helmholtz planes on the
oxide and solution side of the DFT slab, respectively. Additionally,
$q_{\text{DFT}}^{\text{ox.}|\text{sol.}}$ is the charge on the DFT slab
representing the interface between the oxide and solution. Since the electrode
charge $Q$ is defined as including the metal-oxide interface, the continuum
oxide, and the oxide-solution interface, we can assign the charge on the
Helmholtz plane representing the solution as the opposite of this electrode
charge, i.e. $q_{\text{h}}^{\text{sol.}}=-Q$, maintaining charge neutrality
within our larger system.

This arrangement of charge creates a potential gradient across the system. We
measure this as
$\bar{\phi}^{\text{ox.}|\text{sol.}}(z)=\phi^{\text{ox.}|\text{sol.}}_{q}-\phi^{\text{ox.}|\text{sol.}}_{\text{fb}}$,
where $\phi^{\text{ox.}|\text{sol.}}_{q}$ is the planar averaged potential of
the charged system and $\bar{\phi}^{\text{ox.}|\text{sol.}}(z)$ is the
difference in the planar averaged potential due to the charge. We can finally
couple this DFT potential to a bulk charge distribution by extracting the
derivative of the subtracted potential $\frac{d\bar{\phi}}{dz}$ at a specific
cutoff point $z_{\text{edge}}^{\text{ox.}|\text{sol.}}$. In DFT, our oxide
will necessarily only extend for a few layers, but in reality it should extend
for at least several nanometers. This nonphysical oxide surface creates
spurious surface effects on the potential. By setting the cutoff value for the
potential $z_{\text{edge}}^{\text{ox.}|\text{sol.}}$ a few layers within the
material, we avoid these surface impacts changing the potential away from the
correct bulk-like potential.

We can now use this information to inform the potential distribution of the
bulk oxide region using the Poisson-Boltzmann conditions outlined in Eq. 12.
The potential within the continuum oxide region can then be approximated as:

$\bar{\phi}^{\text{ox.}}(z)=\bar{\phi}(z_{\text{edge}}^{\text{ox.}|\text{sol.}})+\frac{\epsilon_{o}\epsilon_{\text{sc}}}{2n_{\text{d}}}\left[\left(\frac{d\bar{\phi}^{\text{ox.}}}{dz}(z)\right)^{2}-\left(\frac{d\bar{\phi}}{dz}(z_{\text{edge}}^{\text{ox.}|\text{sol.}})\right)^{2}\right]-k_{\text{B}}T.$
(14)

This relationship applies throughout the bulk oxide region, which lasts from
$z_{\text{edge}}^{\text{ox.}|\text{sol.}}$ to
$z_{\text{edge}}^{\text{m}|\text{ox.}}$. We will specify the potential drop
across this region as
$\Phi_{\text{ox.}}=\bar{\phi}(z_{\text{edge}}^{\text{m}|\text{ox.}})-\bar{\phi}(z_{\text{edge}}^{\text{ox.}|\text{sol.}})$.

Similarly, for the metal-oxide interface we can say that

$q_{\text{DFT}}^{\text{m}|\text{ox.}}+q_{\text{h}}^{\text{m}|\text{ox.}}=0,$
(15)

where $q_{\text{h}}^{\text{m}|\text{ox.}}$ is the charge on the Helmholtz
planes on the oxide sides of the DFT slab, and
$q_{\text{DFT}}^{\text{m}|\text{ox.}}$ is the charge on the DFT metal-oxide
interface. We note that for metals, all excess charge can be expected to
accumulate at the interface, leading to a flat potential within the bulk
metal. This means that it is unnecessary to include a symmetric Helmholtz
plane representing charge within the rest of the metal electrode. Here, the
value of $q_{\text{h}}^{\text{m}|\text{ox.}}$ should be set equal to the
amount of charge that has accumulated in the system up to that point. Using
Gauss’ law, this would mean that

$q_{\text{h}}^{\text{m}|\text{ox.}}=\epsilon_{o}\epsilon_{\text{sc}}A\frac{d\bar{\phi}^{\text{ox.}}}{dz}\left(z_{\text{edge}}^{\text{m}|\text{ox.}}\right),$
(16)

where $A$ is the lateral surface area of the metal-oxide DFT interface. Since
this charge will be localized to the metal-oxide oxide interface in the
subsequent DFT calculation, the area of the metal-oxide DFT interface is
therefore the relevant one for Gauss’ law.

While this setup ensures continuity between the two different DFT equations,
it does not yet provide guidance on what the correct distribution of charges
in Eq. 13 should be. The final condition to find the correct distribution of
charges is to specify that a system at equilibrium must have the Fermi level
be the same throughout the system:

$\varepsilon_{\text{F}}^{\text{ox.}|\text{sol.}}=\varepsilon_{\text{F}}^{\text{m}|\text{ox.}}+\Phi_{\text{ox.}},$
(17)

where $\varepsilon_{\text{F}}^{\text{ox.}|\text{sol.}}$ and
$\varepsilon_{\text{F}}^{\text{m}|\text{ox.}}$ are the Fermi level of the DFT
oxide-solution interface and DFT metal-oxide interface, respectively. By
adding in the potential shift from the oxide $\Phi_{\text{ox.}}$, we make sure
these are evaluated from the same reference point. We can relate these DFT
calculated Fermi levels to the vacuum energy referenced work function by
adding on a potential shift from the solution $\Phi_{\text{sol.}}$:

$W=\Phi_{\text{sol.}}+\varepsilon_{\text{F}}^{\text{ox.}|\text{sol.}}.$ (18)

This solution potential shift can be determined directly by another Poisson-
Boltzmann shift,[38, 39] or by a constant solvation shift. In real world
scenarios, this potential drop would further be determined by the salt
concentration and pH of the solution. These effects can be partially included
when calculating the formation energy of ions as in Eqn. 2, and can also be
somewhat included in the implicit Poisson-Boltzmann solvers mentioned above.
For the sake of simplicity, within this work we will run our simulations in
vacuum and do not add a solution contribution since we already including water
molecules and Helmholtz planes as a first-order approximation for the voltage-
charge response. Future work should focus on integrating more detailed
knowledge of the surrounding solution into the voltage drop within the
solution.

To relate the calculated work function from DFT to the equilibrium voltage of
a DFT calculation, we will use the Trasatti convention,[40] where we define
the potential of a given DFT slab as (versus SHE) as

$\Phi=W/|e|-4.44\text{ V},$ (19)

where $W$ is the work function of the system, defined as the opposite of the
DFT calculated Fermi level, i.e. $W=-\varepsilon_{\text{F}}$. (This relation
is only strictly true for pure water, but we will accept the approximation for
this system). Thus, by using this methodology, we can specify a charge $Q$ on
the electrode, and then test different distributions of
$q_{\text{h}}^{\text{ox.}|\text{sol.}}$ and
$q_{\text{DFT}}^{\text{ox.}|\text{sol.}}$ to find the equilibrium charge state
and voltage of the system. This allows standard DFT slab calculations of
metal-insulator-solution interfaces to fully include voltage effects– without
the need to simulate 10s of nm of material.

By reducing the number of atoms needed to simulate within DFT by at least an
order of magnitude, this work significantly speeds up calculations of
realistic voltages for the system, to the point of making $>$ 5 nm systems
computationally feasible. Nevertheless, several smaller DFT calculations are
needed to optimize the free energy within the system, finding the charge
distribution which leads to a Fermi level being in equilibrium throughout the
interface. In this work, we manually test different combinations of charge to
find the correct equilibrium distribution, which leads to roughly five to ten
DFT slab calculations needed for each data point. We have previously
implemented an automated Newton-Raphson optimization algorithm for finding
equilibrium free energy within semiconductor-solution systems.[31] A similar
algorithm could also be extended to these metal-insulator-solution systems for
automated voltage processing.

### DFT Details

We use structures based on the interfaces used in Leung’s work on Al
surfaces.[24] Leung investigated chloride insertion into
Al(111)/$\alpha$-Al2O3(0001) surfaces, finding that chloride insertion was
only favorable in situations with grain boundaries. By using these atomic
structures, we can make more direct comparisons with the Leung work. The
$\alpha$ phase of Al2O3 is selected as an approximation of the amorphous
alumina that would develop on a surface in contact with water molecules. The
Al/$\alpha$-Al2O3 structure is based on the “FCC” interface structure,[17] but
with double lattice cell dimensions. We work with two interfaces throughout
the work, one without a grain boundary, shown in Fig. 3, and one with a grain
boundary, shown in Fig. 4. For each interface, we create a structure both at
the metal-oxide interface and the oxide-solution interface. All the oxide-
solution structures in question include a monolayer of water on the surface
representing a buildup from humidity. In Fig. 3b and d, this represents 6
water molecules, and in Fig. 4b and d this represents 12 water molecules. The
arrangement of water molecules in this work was taken from Leung’s work,[24]
which was from a snapshot from ab-initio molecular dynamics (AIMD) chosen to
be typical of a system at equilibrium. The exact arrangement of water
molecules has been shown to have significant effects on the total work
function of DFT calculations.[41, 42, 43] At non-zero temperatures, the water
molecules will necessarily be dynamic, and this represents a reasonable
snapshot. By keeping the arrangement of water molecules consistent between the
structures with and without a grain boundary, moreover, we minimize the
impacts this has on the relative error between the two structures. The
structures used in this work have been included in the supplementary files.

All electronic structure calculations are done using the quantum espresso
package [44]. We use norm-conserving pseudopotentials from the PseudoDojo
repository [45] and the Perdew-Burke-Ernzerhof exchange-correlation
functional.[46] The Al atom uses a Neon core with a valence configuration of
3s23p1 (Z=3), the O atom uses a Helium core with a 2s22p4 valence
configuration (Z=6), and the Cl atom uses a Neon core with a valence
configuration of 3s23p5 (Z=7). We use kinetic energy cutoffs of 50 Ry and 200
Ry for the plane wave basis sets used to describe the Kohn-Sham orbitals and
charge density, respectively. We use a 2$\times$2$\times$1 Monkhorst-Pack grid
[47] to sample the Brillioun zone in our calculations. Ionic relaxations are
done until forces are converged below 50 meV/Å and the electronic calculations
are done with a 10-10 Ry energy threshold. All atoms were relaxed and none
frozen.

We use the Environ package [48] to calculate the parabolic corrections to the
surface dipole. For simplicity, all calculations are done in vacuum, with a
vacuum region of $\geq$ 20 Å.

In the results calculated from this work, we use the following
parameterizations for the QCA equations outlined above:
$\epsilon_{\text{sc}}=9.1$ to match Al2O3, $p_{d}=0$, $n_{d}=10^{15}$ cm-3,
$T=300$ K, and $c^{\circ}=0$ since we are simplifying the problem as in vacuum
and thus not including Poisson-Boltzmann contributions from the solution. We
instead use the Helmholtz distribution of ionic charge in the solution region
combined with an explicit water monolayer, which his often used as a first
approximation of full ionic distributions.[33]

## Results

Figure 3: The structures used for the Al(0001)/$\alpha$-Al2O3-water interface
calculations with (a) the Al(0001)/$\alpha$-Al2O3 interface with Cl inserted,
(b) the $\alpha$-Al2O3/water interface used for the Cl calculations, (c) the
Al(0001)/$\alpha$-Al2O3 interface with an oxygen vacancy, highlighted by the
red circle, and (d) the $\alpha$-Al2O3/water interface used for the O-vacancy
calculations.

We first consider chloride ion insertion in an interface without any grain
boundaries. To start, we calculate the flatband potential of the system before
any oxygen vacancies or chlorine atoms are inserted. For the oxide-solution
system we find a Fermi level of
$\varepsilon_{\text{F,fb}}^{\text{ox.}|\text{sol.}}=-6.2$ eV and for the
metal-oxide system we find
$\varepsilon_{\text{F,fb}}^{\text{m.}|\text{ox.}}=-6.1$ eV. To find the
flatband potential of the total system it is necessary to find the arrangement
of charge that causes the Fermi level to be in equilibrium throughout the
system, and the total charge on the electrode must be neutral, i.e. $Q=0$ .
This constrains the problem sufficiently that the flatband potential of the
entire system can be calculated as $\Phi_{\text{fb}}\approx 1.7$ V (SHE).

We then investigate placing an oxygen vacancy into the structure. Following
Leung’s work, we place the oxygen vacancy 4.3 Å away from the metal surface,
where the zero $z$-location is here defined as the first oxygen atoms of the
oxide structure moving away from the metal. We allow the system to relax
around the vacancy, which is uncharged. We then place a chlorine atom into
this oxygen vacancy and again allow the system to relax. We assume that oxygen
vacancies are not mobile within the timeframe of chloride insertion. For the
sake of this demonstration, we do not examine the kinetics of chloride
insertion, but the QCA techniques presented here do provide an opportunity for
future work calculating kinetic barriers while the surface is under applied
voltage. When we double the lateral size of the system, we observe a -0.16 V
in the calculated flatband potential for the insertion of the Cl ion,
indicating that the size of the system is reasonably converged so that
concentration effects will not overwhelm the calculation. We examine a number
of changes to the DFT system and their impact on the Cl insertion flatband
potential in the Appendix. Using Bader charge analysis,[49] we confirm that
the Cl atom carries $\approx$ 8 electrons, and thus is, in fact, acting as a
Cl- ion within the system.

Figure 4: The structures used for the Al(0001)/$\alpha$-Al2O3-water interface
calculations including a grain boundary for (a) the Al(0001)/$\alpha$-Al2O3
interface with Cl inserted, (b) the $\alpha$-Al2O3/water interface used for
the Cl calculations, (c) the Al(0001)/$\alpha$-Al2O3 interface with an oxygen
vacancy, highlighted by the red circle, and (d) the $\alpha$-Al2O3/water
interface used for the O-vacancy calculations.

We then apply QCA as outlined in above over a variety of voltages, and
calculate the resulting formation energy for chloride insertion into the
oxygen vacancy. Here, the DFT energies for the chloride insertion and oxygen
vacancy, $E_{\text{DFT}}(\text{Cl ins.})$ and $E_{\text{DFT}}(\text{O vac.})$
in Eq. 1 need to be calculated including the impact of voltage as found in the
equilibrium QCA configuration. As seen in Fig. 5, the formation energy for
chloride insertion is highly positive across a wide range of voltages,
implying that chloride is highly unlikely to form at these potentials. The
chloride insertion energy only dips below 0 eV at $\gtrapprox$ 1.2 V (SHE).
This matches previous work by Leung who reports that chloride ions are not
stable inserted into crystalline $\alpha$-Al2O3.[24]

We predict that the relationship between the applied voltage and the chloride
insertion energy is largely linear. This indicates that the dominant term in
the energy is in fact the chloride ion energy introduced in Eqn. 2. This can
be attributed to a majority of the charge in the system accumulating at the
oxide-solution interface, leaving only a small potential drop in the bulk
oxide region, where the parabolic relationships for voltage would become more
dominant. This also highlights the importance of pairing QCA with a
computational standard hydrogen electrode or similar approach which allows for
inclusion of the voltage in the calculation of the ionic energy.

We next consider an interface that includes grain boundaries. We again first
measure the Fermi level of the oxide-solution system as
$\varepsilon_{\text{F,fb}}^{\text{ox.}|\text{sol.}}=-5.54$ eV and the metal-
oxide system as $\varepsilon_{\text{F,fb}}^{\text{m.}|\text{ox.}}=-6.1$ eV.
Using the same procedure as before we can thus find a flatband potential for
the grain boundary system of $\Phi_{\text{fb}}\approx 1.4$ V (SHE). We then
insert oxygen vacancies or chlorine atoms into the structure and allow the
system to relax. It should be noted that due to periodic boundary conditions,
the density of grain boundaries is much higher in this system than would be
realistically expected. We again confirm with Bader charge analysis that the
Cl atom is acting as a Cl- ion within the system.

Figure 5: A stability diagram indicating the calculated insertion energy of a
chloride ion as a function of voltage. When the chloride ion is inserted into
a grain boundary it becomes stable (i.e. $<0$ eV, indicated by the horizontal
black dashed line) at voltages $>-0.2$ V (SHE). This compares favorably with
the experimentally measured onset of chloride pitting at -0.5 V (SHE)
indicated by the vertical red dashed line. In contrast, in systems without
grain boundaries, chloride ion insertion is not favorable until much higher
voltages.

As seen in Fig. 5, we then test a variety of voltages across the system with
grain boundaries. We predict that chloride ion insertion becomes stable at
voltages $>-0.2$ V (SHE). This correlates relatively well with the previous
experiments showing pitting corrosion at voltages $>-0.5$ V (SHE) [2, 3] (see
below for a discussion on comparing theoretical and experimental voltages). It
also matches previous results from Leung, who similarly finds that chloride
ion insertion is only favorable at grain boundaries within the aluminum
oxide.[24] This work emphasizes the role that grain boundaries likely play in
pitting corrosion according to the point defect model, serving as locations
for chloride insertion into the aluminum oxide and serving as corrosion
centers. This result also follows intuition: Cl- is a larger anion than oxygen
and since the grain boundary is less dense than the surrounding Al2O3
inserting the Cl in this grain boundary causes less distortion to the
underlying lattice and is thus more favorable.

Figure 6: The equilibrium charge density at the interface with a Cl- ion
inserted both with and without grain boundaries. Here, negative values are
taken to indicate excessive electrons.

The voltage dependence of the charge density, shown in Fig. 6, also matches
the intuitive picture of chloride stability in grain boundaries we have
established. For interfaces with and without a grain boundary, the charge
density is largely linear as a function of voltage. This reflects the fact
that, as previously noted, the majority of the charge in our simulations of
the system is caught at the oxide/solution interface, causing the interface to
behave in a more metallic manner, particularly away from the flatband
potential. Furthermore, the interface with grain boundaries consistently has
lower charge densities at voltages where the Cl- ion is predicted to be
stable, somewhat matching the intuition that more charge is allowed in the
system. (In this notation, “lower” charge densities indicate a greater density
of electrons induced at the interface.) This is only violated at
$\Phi\lessapprox$ -0.2 V (SHE), corresponding to the location where a Cl- ion
is no longer stable in the system. Below this voltage, the charge density for
the two systems becomes nearly equal.

## Discussion

While QCA provides a useful tool for attaching relative voltages to DFT
calculations of interfaces, it does not solve a variety of outstanding
problems in defining the electronic voltage of a neutral DFT system relative
to experiment. In particular, the arrangement of water molecules within the
first monolayer of coverage has been shown to have a large impact on the
estimated flatband potential, with previous studies on this system showing the
impact of different H2O configurations being 0.2-0.5 eV.[50, 51] This is a
difficult problem to overcome as there is no clearly “correct” configuration
of the water molecules on the surface. One possible solution would be to
average over several possible H2O configurations sampled from an AIMD
calculation and attempt to form some averaging.

As indicated previously, another concern is the relative concentration of
charged defects being examined. The periodic nature of these DFT calculations
ensures that the 2D density of any defects will be relatively high.
Furthermore, realistic systems likely have larger macrostructure, with grain
boundaries and defects that are not actually periodic in nature. This can be
partially overcome by creating larger lateral system sizes and performing
convergence tests. Increasing the lateral size of the interfaces typically
necessitates a a tradeoff of the length of the slab to keep calculations
computationally feasible. This shorter slab can still be accommodated to
provide useful results with QCA provided that a few layers have been included
around the defect and interface.

These twin uncertainties make exact comparison of our predicted chloride
insertion onset voltage to macroscopic potentials extracted from experiment
problematic. Our QCA predictions of exact potentials should then be taken as
semi-quantitative about the system, not exact. For instance, the finding that
grain boundaries within the aluminum oxide make chloride insertion much more
favorable are consistent regardless of the exact voltage at which they become
favorable. Similarly, if our results predicted that chloride insertion was
only favorable only at voltages $>$ 3 V (SHE), it would be likely that we were
not correctly representing the interface within DFT (or that QCA itself was an
incorrect framework for the problem). When predicting relative voltage
differences in the same system (i.e. the difference between two voltages
predicted using the same surface termination), however, QCA should be more
rigorous as the previous drawbacks to exact work function determination are
held constant.

QCA nevertheless presents an advantage over other computational approaches to
controlling the voltage of an electrode by coupling DFT to continuum and
Poisson-Boltzmann descriptions of the solution. This is namely the ability to
include an extended insulating (oxide) region within the simulation. Previous
work has needed to focus on metals with a linear relationship between the
charge and voltage of the system, consistent with a constant capacitance
model.[52, 37, 53, 54] While these approaches can be quite useful, their
utility can be limited in comparison to experiment, especially if large
passivating layers are expected to form over the metal electrode and the
capacitance can no longer be expected to be constant as a function of
potential. Given that the contact layer between metal and oxide has been shown
to sometimes have barriers as high as 2-3 V,[55, 56] the impact of the
extended oxide layer is vital to fully understanding the voltage behavior of
an electrode. QCA thus presents a computationally scalable method for
addressing voltages in systems with extended insulating regions.

## Conclusion and Outlook

In this work, we have extended the QCA methodology to encompass metal-oxide
interfaces that are also in contact with a solution. This extension allows us
to provide first principles calculations of atomic mechanisms behind corrosion
and other electrochemical processes as a function of applied voltage, which
has previously been difficult to manage within a DFT framework. We demonstrate
this methodology on an Al/$\alpha$-Al2O3 interface, calculating that chloride
insertion is only favorable above -0.2 V (SHE) when a grain boundary is
included at the interface. This conclusion is broadly consistent with both
previous first principles work, and experimental measurements. This provides
helpful context for parameterizing the point defect model with first-
principles results.

While QCA provides a useful method for examining interfacial problems that
includes voltage, it does not resolve some of the traditional difficulties of
DFT. Band gaps will be underestimated in semiconductors and insulators when
using semi-local functionals. As long as the potential drop across the
electrode is less than the band gap of the DFT material, however, this should
not affect the relative voltage results. Similarly, the exact atomic
structures of both the metal-oxide and oxide-solution interfaces exert a large
influence on the assigned flatband potential of the interfacial structure.
Numerous DFT methods have been developed to address this issue, with high
temperature ab-initio molecular dynamics providing one option for generating
structures. Larger surveys of possible structures should still be made for QCA
analysis of any given system. It nevertheless presents advantages over
previous work in being able to examine systems with large passivation layers.

In this work, we demonstrated QCA in the context of corrosion of aluminum, but
it can clearly also be used in numerous applications of interest. In battery
electrode research, for instance, the exact atomistic interaction of a solid-
electrolyte interface (SEI) with Li ions shuttling through the structure is
often not well understood as a function of voltage.[57] QCA provides one
potential route for exploring these systems while including a realistic
voltage to the system. Another application may be photocatalysis, where a
semiconductor is placed in contact with solution to catalyze a reaction, often
attached to a metal electrode to provide some level of voltage to the
system.[58] Beyond electrochemistry, a similar methodology can be applied to
electronic interfaces of interest such as p-n junctions and Schottky
barriers.[32] Thus, QCA provides a new tool for computational electrochemistry
across a broad variety of interfaces and applications, including voltage
effects at reasonable computational cost.

We would like to thank Kevin Leung for providing the initial Al/Al2O3 system
and a critical reading of the draft manuscript, as well as Andrew Baczewski
for useful feedback. This work was supported by the Laboratory Directed
Research and Development program at Sandia National Laboratories under project
226354 and 229741. Sandia National Laboratories is a multi-mission laboratory
managed and operated by National Technology and Engineering Solutions of
Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for
DOE’s National Nuclear Security Administration under contract DE-NA0003525.
This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government.

##  Appendix A Dependence of the flatband voltage of the Cl- insertion on DFT
system

The results presented in the main text are heavily dependent on the flatband
potential calculated for Cl- insertion. It is thus natural to investigate how
this flatband potential is dependent on the exact configuration used for our
DFT calculations. We investigate a number of alternatives in the no grain
boundary case, such as doubling the lateral area of the system, moving the
position of the Cl insertion, and increasing the number of oxide layers
modeled, with the results displayed in Table 1. We find that while the exact
number changes, the flatband voltage is largely robust to these changes, with
the exception of moving the Cl ion closer or farther from the metal-oxide
interface. This matches the previous work from Leung.[24] We additionally
tested increasing the number of oxide layers and found a change of -0.01 V in
the flatband potential. These results highlight the importance of converging
system size with respect to flatband potential. Once the system is
sufficiently converged, the length of the DFT region included within the QCA
framework outlined in the main text should not impact the results.

Table 1: Flatband voltage of different system configurations for Cl insertion System | Flatband Voltage (V, SHE)
---|---
Cl inserted cell (Fig. 3a) | 1.66
Lateral doubled Cell | 1.50
Cl- position changed, same distance from m$|$ox. interface | 1.60
Cl- position changed, larger distance from m$|$ox. interface | 0.80
Increase in number of oxide layers modeled | 1.65

##  Appendix B Supporting Files

We have added supporting files which include the scripts used to generate QCA
calculations and input structures for both the calculations with and without
grain boundaries.

## References

## References

  * [1] G. H. Koch, M. P. Brongers, N. G. Thompson, Y. P. Virmani, J. H. Payer, et al.  Corrosion cost and preventive strategies in the United States.  Technical report, United States. Federal Highway Administration, (2002).
  * [2] G. Dibari and H. Read.  Electrochemical behavior of high purity aluminum in chloride containing solutions.  Corrosion, 27(11), 483–494, (1971).
  * [3] J. Bessone, D. Salinas, C. Mayer, M. Ebert, and W. Lorenz.  An EIS study of aluminium barrier-type oxide films formed in different media.  Electrochimica Acta, 37(12), 2283–2290, (1992).
  * [4] L. Lin, C. Chao, and D. Macdonald.  A point defect model for anodic passive films: II. Chemical breakdown and pit initiation.  Journal of the Electrochemical Society, 128(6), 1194, (1981).
  * [5] C. Chao, L. Lin, and D. Macdonald.  A point defect model for anodic passive films: I. Film growth kinetics.  Journal of the Electrochemical Society, 128(6), 1187, (1981).
  * [6] G. Engelhardt and D. D. Macdonald.  Unification of the deterministic and statistical approaches for predicting localized corrosion damage. I. Theoretical foundation.  Corrosion Science, 46(11), 2755–2780, (2004).
  * [7] G. Frankel.  Pitting corrosion of metals: a review of the critical factors.  Journal of the Electrochemical Society, 145(6), 2186, (1998).
  * [8] P. Natishan.  2017 wr whitney award: Perspectives on chloride interactions with passive oxides and oxide film breakdown.  Corrosion, 74(3), 263–275, (2018).
  * [9] J. Evertsson, F. Bertram, F. Zhang, L. Rullik, L. Merte, M. Shipilin, M. Soldemo, S. Ahmadi, N. Vinogradov, F. Carlà, et al.  The thickness of native oxides on aluminum alloys and single crystals.  Applied Surface Science, 349, 826–832, (2015).
  * [10] F. Brown and W. Mackintosh.  The use of Rutherford backscattering to study the behavior of ion-implanted atoms during anodic oxidation of Aluminum: Ar, Kr, Xe, K, Rb, Cs, Cl, Br, and l.  Journal of The Electrochemical Society, 120(8), 1096, (1973).
  * [11] K. Shimizu, G. Brown, K. Kobayashi, P. Skeldon, G. Thompson, and G. Wood.  A novel approach for the study of the migration of Cl- ions in anodic alumina.  Corrosion Science, 41(9), 1835–1847, (1999).
  * [12] L. Serna, K. Zavadil, C. Johnson, F. Wall, and J. Barbour.  A critical implanted cl concentration for pit initiation on aluminum thin films.  Journal of The Electrochemical Society, 153(8), B289, (2006).
  * [13] E.-G. Kim and J.-L. Brédas.  The nature of the aluminum–aluminum oxide interface: A nanoscale picture of the interfacial structure and energy-level alignment.  Organic Electronics, 14(2), 569–574, (2013).
  * [14] I. G. Batyrev and L. Kleinman.  In-plane relaxation of Cu (111) and Al (111)/$\alpha$\- Al 2 O 3 (0001) interfaces.  Physical Review B, 64(3), 033410, (2001).
  * [15] M. Koberidze, M. Puska, and R. Nieminen.  Structural details of Al/Al 2 O 3 junctions and their role in the formation of electron tunnel barriers.  Physical Review B, 97(19), 195406, (2018).
  * [16] D. Costa, T. Ribeiro, F. Mercuri, G. Pacchioni, and P. Marcus.  Atomistic modeling of corrosion resistance: a first principles study of O2 reduction on the Al (111) surface covered with a thin hydroxylated alumina film.  Advanced Materials Interfaces, 1(3), 1300072, (2014).
  * [17] D. J. Siegel, L. G. Hector Jr, and J. B. Adams.  Adhesion, atomic structure, and bonding at the Al (111)/$\alpha$\- Al 2 O 3 (0001) interface: A first principles study.  Physical Review B, 65(8), 085415, (2002).
  * [18] J. Weber, A. Janotti, and C. Van de Walle.  Point defects in Al2O3 and their impact on gate stacks.  Microelectronic Engineering, 86(7-9), 1756–1759, (2009).
  * [19] J. Carrasco, J. R. Gomes, and F. Illas.  Theoretical study of bulk and surface oxygen and aluminum vacancies in $\alpha$\- Al 2 O 3.  Physical Review B, 69(6), 064116, (2004).
  * [20] N. Hine, K. Frensch, W. Foulkes, and M. Finnis.  Supercell size scaling of density functional theory formation energies of charged defects.  Physical Review B, 79(2), 024112, (2009).
  * [21] F. Janetzko, R. A. Evarestov, T. Bredow, and K. Jug.  First-principles periodic and semiempirical cyclic cluster calculations for single oxygen vacancies in crystalline Al2O3.  physica status solidi (b), 241(5), 1032–1040, (2004).
  * [22] M. Liu, Y. Jin, C. Leygraf, and J. Pan.  A DFT-Study of Cl Ingress into $\alpha$-Al2O3 (0001) and Al (111) and its possible influence on localized corrosion of Al.  Journal of the Electrochemical Society, 166(11), C3124, (2019).
  * [23] M. Liu, Y. Jin, B. Chen, C. Leygraf, L. Wang, and J. Pan.  Density Functional Theory Study of Influence of Oxide Thickness and Surface Alloying on Cl Migration within $\alpha$-Al2O3.  Journal of The Electrochemical Society, 168(8), 081508, (2021).
  * [24] K. Leung.  First principles, explicit interface studies of oxygen vacancy and chloride in alumina films for corrosion applications.  Journal of The Electrochemical Society, 168(3), 031511, (2021).
  * [25] J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard, and H. Jonsson.  Origin of the overpotential for oxygen reduction at a fuel-cell cathode.  The Journal of Physical Chemistry B, 108(46), 17886–17892, (2004).
  * [26] J. Rossmeisl, Z.-W. Qu, H. Zhu, G.-J. Kroes, and J. K. Nørskov.  Electrolysis of water on oxide surfaces.  Journal of Electroanalytical Chemistry, 607(1-2), 83–89, (2007).
  * [27] I. C. Man, H.-Y. Su, F. Calle-Vallejo, H. A. Hansen, J. I. Martínez, N. G. Inoglu, J. Kitchin, T. F. Jaramillo, J. K. Nørskov, and J. Rossmeisl.  Universality in oxygen evolution electrocatalysis on oxide surfaces.  ChemCatChem, 3(7), 1159–1165, (2011).
  * [28] P. Atkins and J. De Paula.  Physical Chemistry, volume 8.  W. H. Freeman and Company, (2006).
  * [29] R. J. Le Roy and R. B. Bernstein.  Dissociation energies and long-range potentials of diatomic molecules from vibrational spacings: The halogens.  Journal of Molecular Spectroscopy, 37(1), 109–130, (1971).
  * [30] Q. Campbell and I. Dabo.  Quantum-continuum calculation of the surface states and electrical response of silicon in solution.  Physical Review B, 95(20), 205308, (2017).
  * [31] Q. Campbell, D. Fisher, and I. Dabo.  Voltage-dependent reconstruction of layered Bi 2 WO 6 and Bi 2 MoO 6 photocatalysts and its influence on charge separation for water splitting.  Physical Review Materials, 3(1), 015404, (2019).
  * [32] S. Subramanian, Q. T. Campbell, S. K. Moser, J. Kiemle, P. Zimmermann, P. Seifert, F. Sigger, D. Sharma, H. Al-Sadeg, M. Labella III, et al.  Photophysics and electronic structure of lateral graphene/MoS2 and metal/MoS2 junctions.  ACS Nano, 14(12), 16663–16671, (2020).
  * [33] W. Schmickler and E. Santos.  Interfacial Electrochemistry.  Springer Science & Business Media, (2010).
  * [34] J. Bardeen.  Surface states and rectification at a metal semi-conductor contact.  Physical Review, 71(10), 717, (1947).
  * [35] A. Cowley and S. Sze.  Surface states and barrier height of metal-semiconductor systems.  Journal of Applied Physics, 36(10), 3212–3220, (1965).
  * [36] A. J. Bard, F.-R. F. Fan, A. S. Gioda, G. Nagasubramanian, and H. S. White.  On the role of surface states in semiconductor electrode photoelectrochemical cells.  Faraday Discussions of the Chemical Society, 70, 19–31, (1980).
  * [37] N. Keilbart, Y. Okada, A. Feehan, S. Higai, and I. Dabo.  Quantum-continuum simulation of the electrochemical response of pseudocapacitor electrodes under realistic conditions.  Physical Review B, 95(11), 115423, (2017).
  * [38] G. Fisicaro, L. Genovese, O. Andreussi, N. Marzari, and S. Goedecker.  A generalized Poisson and Poisson-Boltzmann solver for electrostatic environments.  The Journal of Chemical Physics, 144(1), 014103, (2016).
  * [39] K. Mathew, V. C. Kolluru, S. Mula, S. N. Steinmann, and R. G. Hennig.  Implicit self-consistent electrolyte model in plane-wave density-functional theory.  The Journal of Chemical Physics, 151(23), 234101, (2019).
  * [40] S. Trasatti.  The absolute electrode potential: an explanatory note (Recommendations 1986).  Pure and Applied Chemistry, 58(7), 955–966, (1986).
  * [41] J. Cheng and M. Sprik.  Alignment of electronic energy levels at electrochemical interfaces.  Physical Chemistry Chemical Physics, 14(32), 11245–11267, (2012).
  * [42] N. Kharche, J. T. Muckerman, and M. S. Hybertsen.  First-principles approach to calculating energy level alignment at aqueous semiconductor interfaces.  Physical Review Letters, 113(17), 176802, (2014).
  * [43] N. G. Hörmann, Z. Guo, F. Ambrosio, O. Andreussi, A. Pasquarello, and N. Marzari.  Absolute band alignment at semiconductor-water interfaces using explicit and implicit descriptions for liquid water.  npj Computational Materials, 5(1), 1–6, (2019).
  * [44] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al.  QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials.  Journal of Physics: Condensed Matter, 21(39), 395502, (2009).
  * [45] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G.-M. Rignanese.  The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table.  Computer Physics Communications, 226, 39–54, (2018).
  * [46] J. P. Perdew, K. Burke, and M. Ernzerhof.  Generalized gradient approximation made simple.  Physical Review Letters, 77(18), 3865, (1996).
  * [47] H. J. Monkhorst and J. D. Pack.  Special points for Brillouin-zone integrations.  Physical Review B, 13(12), 5188, (1976).
  * [48] O. Andreussi, I. Dabo, and N. Marzari.  Revised self-consistent continuum solvation in electronic-structure calculations.  The Journal of Chemical Physics, 136(6), 064102, (2012).
  * [49] G. Henkelman, A. Arnaldsson, and H. Jónsson.  A fast and robust algorithm for Bader decomposition of charge density.  Computational Materials Science, 36(3), 354–360, (2006).
  * [50] S. Sakong, K. Forster-Tonigold, and A. Groß.  The structure of water at a Pt (111) electrode and the potential of zero charge studied from first principles.  The Journal of Chemical Physics, 144(19), 194701, (2016).
  * [51] S. Sakong and A. Groß.  The electric double layer at metal-water interfaces revisited based on a charge polarization scheme.  The Journal of Chemical Physics, 149(8), 084705, (2018).
  * [52] S. Schnur and A. Groß.  Challenges in the first-principles description of reactions in electrocatalysis.  Catalysis Today, 165(1), 129–137, (2011).
  * [53] J. Haruyama, T. Ikeshoji, and M. Otani.  Electrode potential from density functional theory calculations combined with implicit solvation theory.  Physical Review Materials, 2(9), 095801, (2018).
  * [54] K. Kano, S. Hagiwara, T. Igarashi, and M. Otani.  Study on the free corrosion potential at an interface between an Al electrode and an acidic aqueous NaCl solution through density functional theory combined with the reference interaction site model.  Electrochimica Acta, 377, 138121, (2021).
  * [55] X.-x. Yu and L. D. Marks.  Combining the physics of metal/oxide heterostructure, interface dipole, band bending, crystallography, and surface state to understand heterogeneity contrast in oxidation and corrosion.  Corrosion, 75(2), 152–166, (2019).
  * [56] P. Cornette, D. Costa, and P. Marcus.  Relation between Surface Composition and Electronic Properties of Native Oxide Films on an Aluminium-Copper Alloy Studied by DFT.  Journal of The Electrochemical Society, 167(16), 161501, (2020).
  * [57] K. Leung.  DFT modelling of explicit solid–solid interfaces in batteries: Methods and challenges.  Physical Chemistry Chemical Physics, 22(19), 10412–10425, (2020).
  * [58] K. Wenderich and G. Mul.  Methods, mechanism, and applications of photodeposition in photocatalysis: a review.  Chemical Reviews, 116(23), 14587–14619, (2016).


# Airborne virus transmission under different weather conditions

Santosh K. Das<EMAIL_ADDRESS>School of Physical Sciences, Indian
Institute of Technology Goa, Ponda-403401, Goa, India  Jan-e Alam
<EMAIL_ADDRESS>Variable Energy Cyclotron Centre, 1/AF Bidhan Nagar, Kolkata-
700064, India  Homi Bhabha National Institute, Training School Complex, Mumbai
- 400085, India  Salvatore Plumari<EMAIL_ADDRESS>Department
of Physics and Astronomy, University of Catania, Via S. Sofia 64, I-95125
Catania, Italy  Laboratori Nazionali del Sud, INFN-LNS, Via S. Sofia 62,
I-95123 Catania, Italy  Vincenzo Greco<EMAIL_ADDRESS>Department of
Physics and Astronomy, University of Catania, Via S. Sofia 64, I-95125
Catania, Italy  Laboratori Nazionali del Sud, INFN-LNS, Via S. Sofia 62,
I-95123 Catania, Italy

###### Abstract

The COVID19 infection is known to be disseminated through the droplets ejected
by infected individual during coughing, sneezing, speaking and breathing. The
spread of the infection and hence its menace will depend on how the virus-
loaded droplets evolve in space and time with changing environmental
conditions. In view of this, we investigate here the evolution of the droplets
under the purview of the Brownian motion of the evaporating droplets in the
air under the action of gravity with varying weather conditions of relative
humidity, temperature and wind flow. We track the movement of the droplets
till either they gravitationally settle on the ground or evaporate to aerosols
of size 2$\mu$m or below and remain suspended in the air. The correction to
drag force required for flow of droplets beyond the laminar regime has been
taken into consideration. We find that on the one hand the larger evaporating
droplets falls on the ground under the action of gravity and on the other hand
the smaller evaporating droplets generate aerosols. The intermediate size
droplets are very sensitive to the weather conditions which decides whether
the intermediate size droplets will settle in the ground gravitationally or
evaporate to aerosol. Therefore, it is clear that the weather condition play
crucial in deciding the intensity of spreading of the infection. Since all the
possible interactions of the droplets with the surroundings are taken into
consideration the results obtained here can be applied for both big as well as
small droplets. The effects of relative humidity and temperature on the
evaporation of the droplet are found to be significant. It is also noted that
under strong wind flow the droplet can travel large distances. These weather
dependence effects indicate that the maintenance of physical separation to
evade the virus is not corroborated and hence the use of face mask becomes
indispensable.

###### pacs:

12.38.Mh, 12.39.-x, 11.30.Rd, 11.30.Er

##  I Introduction

It is well known that droplets released by infected persons through coughing,
sneezing, speaking or breathing contain microorganism (bacteria, virus, fungi,
etc) causing a large number of diseases cflugge ; Wells . Several decades ago
it was thought-out that the diseases which are spread by droplets are airborne
Gelfand  and recently ten scientific reasons have been provided in support of
SARS-COV-2 being an airborne infection 10reasons  (see also review  for
review). The transmission routes kutter  of the virus crucially depends on how
the droplets evolve is space with the progression of time. For example, the
big droplets settle gravitationally on different surfaces like, hand rail,
door handle, tables, etc and eventually infect through direct contact. The
small droplets either directly discharged or formed due to evaporation of the
big droplets will remain suspended in the air for long time due to the
dominance of the diffusive force over the gravitational force. It still
remains as a big challenge to understand the density of virus in the air mpan
, their capability to infect bynumbers , their survivability on different
surfaces doremalen  and comparison. The survivability of the virus on solid
surfaces increases under humid conditions AA . Therefore, these factors limits
our ability to evaluate the risk Nazhu . The release of virus via respiration
process richard ; herfst  and speaking is well known kkwang ; asadi ;
stadnytskyl . The saliva of the infected persons contains the Coronavirus
kkwang ; lazzi . The increase in emission of pathogens with the loudness in
speaking which may depend on some unknown physiological factors varying among
individuals has also been reported asadi . Recently it has been reported that
the droplets and aerosols (droplets with size $<5\mu$m) can travel a distance
much larger distance than the prescribed 6 feet and remain suspended in the
air for hours and in fact SARS-COV2 RNA is recovered in air sample jama20 .

The propagation and aerosolisation of the droplets have been investigated
within the purview of statistical mechanics and fluid dynamics dbouk ; mittal
. Euler-Lagrange equation has been applied to study the effects of the size,
ejection velocity and angle of emission of the droplets on the motion of the
droplets during sneezing. The application of the stochastic statistical
mechanics Reif ; pathria  is particularly crucial to study the motion of
aerosols (droplets with diameter $<5\mu$m seta ) for which the airborne
transmission turns out to be very vital. The Brownian motion of the aerosols
in the air-bath can be studied by solving Langevin differential equation which
takes care of the stochastic motion SKDas  normally ignored in Eulerian-
Lagrangian approach dbouk .

We investigate the spread of these ejected droplets in the neighbourhood of
the infected individual under varying weather conditions which will help in
planning the preventive strategies at different climatic conditions. The
droplets ejected will interact with the molecules of the still/flowing air at
temperature ($T$), relative humidity (RH) in presence of gravity. The
interaction force between the droplets and the air molecules is changing
continuously as the molecules are changing their coordinates continuously with
time. This makes the problem very complex, if not exactly unsolvable. The
spread of the droplets, however, will depend on these interactions. Under such
circumstances the air molecules can be regarded as forming a thermal bath
where the droplets are executing Brownian motion with its changing mass due to
evaporation. The interaction of the evaporating droplets with the bath can
then be grouped into drag and diffusive forces quantified through the drag and
diffusion coefficients under the action of gravity. Therefore, there are three
types of force acting on a droplet which are (i) drag and (ii) diffusive
forces between the droplet and the air molecules and the (iii) Newtonian
gravitational force on the droplets due to their non-zero mass. However, it is
crucial to note that droplets are undergoing loss of mass due to evaporation
and hence all these forces are changing with time.

In the present work we investigate the propagation of the virus-containing
droplets subject to continuous evaporation by solving the Langevin stochastic
differential equation of statistical mechanics Reif ; pathria  coupled with
the equation that governs the evaporation of the droplets. It is crucially
important to mention here that the correction to the drag force required for
flow beyond the laminar regime has been taken care of by using the
experimental parameterization the Reynold’s number Holterman . The Langevin
equation is applicable in the present context as the mass of the droplets are
much higher than the mass of oxygen and nitrogen molecules present in the air.
The inclusion of all the forces mentioned above enable us to study the
trajectories of droplets with a wide range of sizes. It will be interesting to
study how the evaporating droplets evolve in space and time under the
influence of gravitation which will act to pull the droplet on the ground in
contrast to the diffusive and drag forces which will prevent it to do so. The
big (hence massive) droplets is expected to settle gravitationally quickly and
the smaller one is expected to remain suspended in the air for longer time.
However, under evaporation the droplet will suffer continuous loss of mass and
consequently a droplet which will otherwise fall on the ground due gravitation
may not do so but may remain suspended as smaller droplets/aerosol or isolated
virus for longer time before complete decomposition.

The paper is organized as follows. In the next section the solution of the
Langevin equation coupled with the equation controlling the evaporation has
been discussed. Section III contains the results and section IV is devoted to
summary and discussions.

##  II Monte-Carlo solution of the Langevin equation in the presence of
evaporation

The Langevin equation governing the motion of the droplet of mass ($M$) in the
still air in the presence of gravitational field Reif  is given by:

$\frac{dr_{i}}{dt}=v_{i}$ (1) $M\frac{dv_{i}}{dt}=-\lambda v_{i}+\xi(t)+F^{G}$
(2)

The mass of the droplet is related its diameter ($D$) by $M=\pi
D^{3}\rho_{L}/6$ which is bound to change due to evaporation. The rate of
decrease of the diameter $D$ of a spherical liquid drop due to evaporation is
given by Kukkonen ; Holterman :

$\frac{dD}{dt\texttimes}=-\frac{4M_{L}D_{v}}{D\rho_{L}RT_{f}\texttimes}\Delta
p(1+0.276Re^{1/2}Sc^{1/3})$ (3)

These coupled set of equations, i.e. Eqs. 1, 2 and 3 are solved simultaneously
to understand the spread of the virus through droplets.

In Eqs. 1 and 2 the $dr_{i}$ and $dv_{i}$ are the shifts of the coordinate and
velocity in each discrete time step $dt$, $i(=x,y,z)$ stands for the Cartesian
components of the position and velocity vectors. The $\lambda$ in Eq. 2 is the
drag coefficient which will be fixed later. The first term in the right hand
side of Eq. 2 represents the dissipative force and the second term stands for
the diffusive (stochastic) force where $\xi(t)$ is regulated by the diffusion
coefficient $\kappa$. $\xi(t)$ is also called noise due to its stochastic
nature. We study the evolution with a white noise ansatz for $\xi(t)$, i.e
$\langle{\bf\xi}(t)\rangle=0$ and
$\langle{\bf\xi}(t){\bf\xi}(t^{\prime})\rangle=\kappa\delta(t-t^{\prime})$.
White noise describes a fluctuating field without memory, whose correlations
have an instantaneous decay called $\delta$ correlation. The third term in Eq.
2, $F^{G}$ represents the gravitational force $(=Mg$, $g=9.8$ m/s2) acting on
a droplet of mass $M$ which will changing with time due to the evaporation.

Eq. 3 describes the evaporation of the droplets with $M_{L}$ ($=0.018$ kg/mol)
and $\rho_{L}$ are the molecular weight and density of the evaporating liquid
respectively, $D_{v}$ is the diffusion coefficient of the vapor molecules in
the saturated film around the droplets and $T_{f}$ is the average temperature
of the film formed around the droplets due to evaporation. $R=8.3144$ J/(mol
K) is the gas constant. $Re$ is the Reynold’s number Landau :

$Re=\frac{\rho_{a}D_{v}v}{\eta_{a}\texttimes}$ (4)

where $v$ is the relative velocity of the droplet to the surrounding air,
$\rho_{a}$ is the density and $\eta_{a}$ is the viscosity of the air at
temperature $T_{f}$. $Sc$ is the Schmidt’s number, given by:

$Sc=\frac{\eta_{a}}{\rho_{a}D_{v}\texttimes}$ (5)

$\Delta p$ is the difference between the vapor pressure near the droplet and
in the atmosphere which acts as the driving force for the transport of vapor
away from the droplet surface. $\Delta p$ can be related to the saturated
vapor pressures at ambient temperature $(T)$ and wet-bulb temperature
($T_{w}$) as:

$\Delta p=p_{sat}-p=\gamma(T-T_{w})$ (6)

where $p_{sat}$ is the vapor pressure near the surface of the droplet, $p$ is
the vapor pressure in the ambient air and $\gamma$ is approximately constant (
67 pa/K). $T_{w}$ can be expressed in terms of $T$ and $RH$ as:

$T_{w}=T-[(a_{0}+a_{1}T)+(b_{0}+b_{1}T)RH+(c_{0}+c_{1}T)RH^{2}]$ (7)

with $a_{0}$=5.1055, $a_{1}$=0.4295, $b_{0}$=-0.04703, $b_{1}$=-0.005951,
$c_{0}$=-0.00004005 and $c_{1}$=0.0000166. For details we refer to Holterman .

The trajectory of the droplets will depend on the initial ejected velocity and
on the flow velocity of the ambience resulting in the path difference in still
and flowing air conditions [such as in a air conditioned (AC) room or even is
open air]. The results for both the scenarios of still and flowing air will be
derived here. The air flow velocity has been taken into account through the
Galilean transformation of the Langevin equation. Since the virus carrying
droplets follow different trajectories in still and flowing air the preventive
strategies for indoor and out door conditions should take care of this fact.
Here we consider a velocity profile for the air flow as:
$u(x)=u_{0}(1-\frac{x}{x_{\text{max}}})$ with its upward and downward
components as zero to serve this purpose, where $x$ is the running coordinate,
$u_{0}$ is the peak value of $u(x)$ at $x=0$ and $x_{\text{max}}$ is the
maximum value of $x$, which may be restricted to the size of an AC room. It is
obvious that any non-zero upward (downward) component will influence the time
of suspension of droplets in the air. We will also consider constant (no
dependence on spatial coordinates) air flow velocity to calculate the distance
that an ejected drop will travel an open air.

Eqs. 1, 2 and 3 have been simultaneously solved by using Monte-Carlo
techniques SKDas ; mc1 ; mc2  with the following inputs. The initial spatial
coordinate of the droplet is , $x=y=0$ and $z=H_{0}$, where $H_{0}$ is the
height (1.7 meter) at which the droplet is released (nose/mouth), i.e.
$(x,y,z)=(0,0,1.7{\text{meter}})$, is the point of ejection. The initial
velocity (at $t=0$) is uniformly distributed in the $x-y$ plane with
$v_{z}=0$. The calculation is done for droplet radii ($R$) varying from 10
$\mu$m to 200$\mu$m RR  and the ejection velocity, $V_{0}=21$ m/s VE . The
value of the drag coefficients, $\lambda$ is estimated by using the relation,

$\lambda=\frac{1}{2}C_{D}\rho_{a}\,S\,v$ (8)

where $S$ is the projected cross sectional area of the spherical droplet of
diameter $D$, $v$ is the velocity and the following expression is used for
$C_{D}$ hongping  to extend the validity of the model beyond the regime of
laminar flow.

$C_{D}=\frac{24}{R_{e}}+\frac{6}{1+\sqrt{R_{e}}}+0.4$ (9)

The diffusion coefficient is obtained by using the Einstein relation pathria ,
$D=K_{B}T\lambda$, where $K_{B}=1.38\times 10^{-23}J/^{\circ}K$, is the
Boltzmann constant. It may be noted that for $C_{D}=24/Re$ the Stokes law is
recovered.

By solving the set of Eqs. 1, 2 and 3 we get the trajectories of the droplets
which enables us to calculate the horizontal distance ($L(t)$) traveled by the
droplets from the point of ejection as a function of time
$L(t)=\sqrt{x(t)^{2}+y(t)^{2}}$. Its maximum value of $L(=L_{\text{max}})$
dictates the stationary distance that to be maintained between infected and
healthy persons to prevent the virus. The solution of these equations can also
be used to estimate the maximum time ($t_{\text{max}}$) of suspension of the
droplets. In the following results for both $L_{\text{max}}$ and
$t_{\text{max}}$ have been presented.

##  III Results

The virus carrying droplets are ejected with different sizes and initial
velocity in the ambience with widely varying climatic conditions. Therefore,
results with different size of the droplets under different meteorological
conditions are exhibited here to understand the prevention measures to be
adopted to avoid the infection. We assume that the ejection velocity is 21 m/s
(unless stated otherwise) which is close to its highest possible value as
found experimentally VE ; zhu . For given wind flow velocity, this will give
the maximum distance that a droplet can travel which in turn will provide the
maximum risk involved and accordingly helps in deciding the safest possible
measures to be taken. In the present work the weather effects through
temperature, relative humidity, wind flow have been taken care of. The effects
of evaporation and correction to drag force for flow beyond the laminar regime
have been included in all the results displayed below unless explicitly
mentioned. The value of $T=20^{\circ}$C whenever the value of $T$ is not
mentioned.

Figure 1: The variation of the height ($H(t)$) with the horizontal distance
($L(t)=\sqrt{x^{2}+y^{2}}$) is shown for droplets of various radii for the
initial ejection velocity, $V_{0}=21$ m/s and the peak value of the wind flow
velocity, $u_{0}$=0.1 m/s. The results contain the effects of evaporation at
ambience temperature $20^{\circ}$ C.

In Fig. 1 the variation of the height ($H$) with longitudinal distance ($L$)
for droplets of different radii released at a height of 1.7 meter is depicted.
It is seen that a large droplets of radius 200 $\mu$m will propagate upto a
distance of 1.4 meter horizontally in still air condition and 1.45 meter in
wind flow condition with the peak velocity of the wind, $u_{0}=0.1$ m/s. The
interplay between the ejection velocity and the wind velocity can be
understood from the following discussions. It may be mentioned here that the
terminal velocity ($v_{t}$) for a droplet of diameter $D$ in air is given by

$v_{t}=\sqrt{\frac{4}{3}\,g\frac{\rho-\rho_{a}}{\rho_{a}}\frac{D}{C_{D}}}$
(10)

which reduces to the well known expression for $v_{t}$ in laminar flow regime
(with $C_{D}=24/R_{e}$) as:

$v_{t}=\sqrt{\frac{g}{18}\,R_{e}\frac{\rho-\rho_{a}}{\rho_{a}}}$ (11)

The effects of wind flow will be insignificant if the terminal velocity of a
droplet is larger than the wind flow velocity. For a large droplet of size 200
$\mu$m the $v_{t}$($\sim\sqrt{D}$) is larger than $u_{0}$ which makes the
effect of wind flow marginal as observed in the results displayed in Fig 1.

A droplet of size 100 $\mu$m will travel a distance 0.55 meter and 0.7 meter
respectively in still and wind flow conditions implying that the wind flow has
larger effects on smaller droplet as their terminal velocity is smaller. For
smaller droplets of size 50 $\mu$m the effect of wind flow is larger as shown
in Fig. 1. We have observed that the droplets of size 25 $\mu$m evaporates to
aerosols with mean radius of 2 $\mu$m. The motion of these aerosols are
controlled by the diffusive force and therefore, these droplets have the
ability float in the air for longer time (see also SChatterjee ). The results
displayed in Fig. 1 illustrate that the social distance to be maintained is
about 1.5 meter. For high wind speed the effects of terminal velocity will be
small. For example, a droplet of radius $100\mu$m can travel about 6.6 meter
for wind speed of 2 m/s HLi .

Figure 2: (Color online) Left panel:The change in the height, $H(t)$, as a
function of time for droplets of different sizes have been depicted. Here the
initial ejection velocity, $V_{0}=21$ m/s and the peak value of the wind flow
velocity, $u_{0}$=0.1 m/s at $T=20^{\circ}$ C. Right panel: Same as left panel
showing the sensitivity of results on the ambient temperature ($T=20^{\circ}$
C and $35^{\circ}$ C. The results are derived with the inclusion of
evaporation.

Fig. 2 illustrates the time that droplets of different sizes will take to
settle gravitationally under different conditions of relative humidity and
wind flow. We find that the droplets at smaller $RH$ will evaporate faster and
make the effect of gravity weaker due to loss of mass. In such cases the time
of suspension in the air will be prolonged before settling gravitationally on
the ground. For larger value of $RH$ the evaporation process will be slowed
down resulting in smaller loss of mass and hence forcing it to fall on the
ground under the action of gravitation. Comparison of results for droplet of
initial radius $50\mu$m displayed in Fig. 2 for $RH=40\%,60\%$ and $80\%$
substantiate this fact with the indication that the effects of $RH$ is
significant. The smaller droplets with initial radius 25 $\mu$m, however,
remain suspended in the air after evaporating to droplet nuclei of mean radius
2 $\mu$m which do not settle under the action of gravitation but continue to
diffuse in air making the use of mask mandatory to prevent these aerosol. At
higher temperature ($35^{\circ}$ C) the evaporation becomes faster.

It is observed that at high $RH$ and low $T$ the effect of evaporation is
small and hence the virus may survive for longer time in such cases. It is
important to note that smaller droplets at higher $RH$ and low $T$ will
survive longer in the air. In other words the evaporation is lesser effective
at high $RH$ and low $T$. The relation between the climatic condition and
index of airborne infection rate and concentration rate of particle in saliva
is found in Ref.dbouk3  (see also dbouk2 ).

Figure 3: (Color online) Left panel:The change in the height, $H$, as a
function of $L$ for a droplet of radius 100 $\mu$m for different peak values
($u_{0}$) of the wind velocity profile at $T=20^{\circ}$ and $RH=60\%$. The
droplet is ejected in a room of size 5 metre with initial ejection velocity is
21 m/s. Right panel: Same as left panel for different wind velocity of
constant values (independent of space coordinate) as indicated in the figure.

The effects of wind flow in indoor (left panel) and outdoor (right panel)
conditions have been demonstrated in Fig. 3. The left panel shows the results
for flow profile described above for various values of $u_{0}$, say, in a room
of size 5 metre. The right panel indicates the results for different wind flow
velocities of constant values (independent of space coordinate), say, an open
air. It is clearly observed that the wind velocity for given temperature,
relative humidity and ejection velocity affects the results substantially when
the wind velocity is larger than the terminal velocity discussed above. The
droplet released by an infected individual can travel a distance much larger
than the prescribed 2 meters separation if the wind velocity is large.

Figure 4: The variation of the maximum horizontal distance ($L_{max}$)
traveled by the droplets as a function of the droplet radius for different
relative humidity have been shown. Here initial ejection velocity, $V_{0}=21$
m/s and the peak value of the wind flow velocity, $u_{0}$=0.1 m/s.

In Fig 4, we have displayed the variation of the maximum horizontal distance
traveled by droplets as a function of its radius for different relative
humidity. We have checked that in still air bigger droplets travel larger
distance for the same initial ejection velocity due to the effect of inertia.
The droplets of small and intermediate sizes are strongly influenced by the
$RH$. We notice mild effects of flow and RH on larger droplets. However, the
droplets with intermediate radii are influenced the most due to variation in
RH and flow mainly due to their longer life time in air than the smaller
droplets which are subjected to quick evaporation at lower $RH$. It is also
clear that the prevention of mass loss in absence to evaporation allows small
droplets to travel longer distance.

Figure 5: The variation of the maximum time the droplets take to settle on the
ground under the action of gravity and the time the smaller droplets take to
evaporate to generate aerosols of average radii $2\mu$m are displayed here.
The sensitivity of the results on the relative humidity has also been shown
for $u_{0}=0.1$ m/s. Figure 6: Change in the running radius, $R(t)$ normalized
to the initial radius have been depicted as a function of time for different
relative humidity. Here $R_{0}$ is the initial radius of the droplets. The
value of $V_{0}$ and $u_{0}$ are taken as $21$ m/s and $0.1$ m/s respectively.
Figure 7: Same as Fig. 6 for running mass to initial mass of the droplets.

Now we would like to estimate the maximum time that evaporating droplets of
different radii will remain suspended in the air before settling on the ground
under different weather conditions. Relevant results to address this issue are
displayed in Fig. 5. The effects of evaporation is found to be significant for
small droplets. The time of suspension is more without evaporation. At higher
RH the droplets have the higher survival probability and hence to travel
larger distance in the air as indicated by the results shown here. However,
the time of suspension of intermediate and smaller sized droplets are very
sensitive to both the evaporation and RH. Through the process of evaporation
smaller droplets may generate isolated virus which may survive for more than
an hour swxong . However, droplets having larger radius fall on the ground
quickly due to stronger gravitational force compare to the diffusive force.
Droplets having intermediate radius stay in the air for longer time due to the
balancing of gravitational force by the diffusive forces.

In Fig 6, we have plotted the change in the running radius $(R)$ of the
droplets normalized to their initial values ($R_{0}$) for different humidity.
The reduction of mass due to evaporation at higher RH allows the droplet to
survive longer. It is noticed that smaller droplets are very sensitive to
evaporation and they generate droplet nuclei before reaching the ground. In
Fig 7, the change in the running masses ($M$) of the droplet normalized to
their initial masses ($M_{0}$) is illustrated as a function of time for
different humidity. The results are consistent with those displayed in Fig. 6.
The results clearly indicate the effects of RH on the evolution of the
droplets. For example a droplet of radius $50\mu$m will approximately loss
$60\%$ and $20\%$ of its mass due to evaporation at $RH=60\%$ and $80\%$
respectively.

##  IV Summary and Discussions

We have investigated the evolution of the droplets ejected during coughing,
sneezing or speaking by solving the Langevin equation with the inclusion of
evaporation of the droplets. The drag, diffusive and the gravitational forces
are included in this study. The correction to the drag force in the non-
laminar flow regime due to large Reynolds number has been implemented through
appropriate parametrization of experimental data. The droplets of various
sizes have been considered with large ejection velocity (21 m/s) to get the
idea on the upper limit of the distance cross over by the droplets. The
effects of different weather conditions have been taken into consideration
through temperature, relative humidity and wind flow. It is found that the
maximum distance that a large droplet of size 200 $\mu$m travels is about 1.5
meter ($u_{0}=0.1$ m/s) which is not affected significantly by low wind
velocity. Tinier droplets after traversing a spatial interval evaporate which
can diffuse through air to go far away from the point of ejection Bouroaiba
due to weaker gravitational influence. It is found that Liu  in two Wuhan
hospitals the micrometer and sub-micrometer droplets of Sars-COV-2 were found
at a distance of about 3 meters away from the bed of an infected person.

An infected person emits droplets of varying sizes. In the present work we
found that the smaller droplets can also be created by the mechanism of
evaporation (see also FanLia ). The smaller droplets can also be created by by
fragmentation xuwang . It is found that droplets of smaller size evaporates
and create aerosols, smaller the sizes quicker they evaporate. Isolated virus
generated by the process of evaporation can survive in the air for more than
an hour swxong . The evolution of these droplets will be controlled by the
diffusive forces and such droplets may remain suspended in the air for hours.
However, it has also been reported jama  that a multi-phase turbulent gas
cloud is produced along with the droplets at the time of coughing and
sneezing. The droplets within the envelope of this gas cloud may evade
evaporation and prolong their survivability as isolated droplets. Droplets of
larger size (say, 100 $\mu$m) can travel a distance as large as 3.6 meter in
indoor condition and even more in outdoor condition if the wind flow velocity
is large. It is appropriate to mention at this juncture that the intermediate
size droplets are very sensitive to the weather conditions. The climatic
conditions determines whether the intermediate size droplets will settle in
the ground under the action of gravity or evaporate to aerosol. This indicates
that the weather condition play crucial role in deciding the severity of
spreading of the infection. Therefore, the use of mask and face shield should
be used to prevent the virus Howard ; pnas ; verma . The maintenance of
ventilation in indoor situation is also very crucial to avoid the infection
NSen ; SBathula . The maintenance of only a social distance of 2 meter as norm
to avoid the virus is not corroborated by the present investigation because
the aerosol created by the evaporation may remained suspended in the air for
more than an hour. The impact of diffusion, represented by the term $\xi(t)$
appearing in the right hand side of Eq. 2 on droplets of smaller of size, say
$5\mu$m or less is more important. However, the impact of diffusion on larger
droplets is insignificant.

## Data Availability Statement

The data that supports the findings of this study are available within the
article.

## Acknowledgement

SKD would like to acknowledge IIT Goa for internal funding (No.
2020/IP/SKD/005) and Professor Barada Kanta Mishra for useful discussions and
encouragement. Fruitful discussions with Raja Mitra is thankfully
acknowledged.

## References

  * (1) C. Flugge, Uber Luftinfection. Z. Hyg. Infektionskr 25, 179 (1897).
  * (2) Wells, W. F. and Wells, M. W., 1938, Measurement of sanitary ventilation, Am. J. Public Health Nations Health, 28, 343 (1938).
  * (3) H. M. Gelfand and J. Posh, ”The recent outbreak of smallpox in Meschede, West Germany”, Am. J. Epidem. 93, 234 (1971).
  * (4) T. Greenhalgh, J. L. Jimenez, K. A. Prather, Z. Tufekci, D. Fisman and Robert Schooley, ”Ten scientific reasons in support of airborne transmission of SARS-CoV-2”, The Lancet, 397, 1603 (2021); https://doi.org/10.1016/S0140-6736(21)00869-2.
  * (5) C. C. Wang, K. A. Prather, J. Sznitman, J. L. Jimenez, S. S. Lakdawala, Z. Tufekci, L. C. Marr, ”Airborne transmission of respiratory viruses”, Science 373, 981 (2021).
  * (6) J. S. Kutter, M. I. Spronken, P. L. Fraaij, R. AM. Fouchier and S. Herfst, Current Opinion in Virology, 28, 142 (2018).
  * (7) M. Pan, J. A. Lednicky and C.-Y. Wu, J. Applied Microbiology, 127, 1596 (2019).
  * (8) Y. M. Bar-On, A. Flamholz, R. Phillips, R. Milo, https://elifesciences.org/articles/57309.
  * (9) N. van Doremalen, T. Bushmaker, D. H. Morris, M. G. Holbrook, A. Gamble, B. N.Williamson, A. Tamin, J. L. Harcourt, N. J. Thornburg, S. I. Gerber, J. O. Lloyd-Smith, E. de Wit,and V. J. Munster, The New England Journal of Medicine, 382, 16 (2020).
  * (10) R. Bhardwaj and Amit Agrawal, ”Likelihood of survival of coronavirus in a respiratory droplet deposited on a solid surface”, Phys. Fluids 32, 061704 (2020); doi: 10.1063/5.0012009.
  * (11) N. Zhu, D. Zhang, W. Wang, X. Li, B. Yang, J. Song, X. Zhao, B. Huang, W. Shi, R. Lu, P. Niu, F. Zhan, X. Ma, D. Wang, W. Xu, G. Wu, G. F. Gao, and W. Tan, The New England Journal of Medicine, 382, 8 (2020).
  * (12) M. Richard, J. M. A. van den Brand, T. M. Bestebroer, P. Lexmond, D. de Meulder, R. A. M. Fouchier, A. C. Lowen and S. Herfst, Nature Communication, 11, 766 (2020).
  * (13) S. Herfst1, M. Boḧringer, B. Karo, P. Lawrence, N. S. Lewis, M. J. Mina, C. J. Russell, J. Steel, R. L. de Swart and C. Menge, Current Opinion in Virology, 22, 22 (2017).
  * (14) S. Asadi, A. S. Wexley, C. D. Cappa, S. Barreda, N. M. Bouvier and W. D. Ristenpart, Nature Scientific Report 9, 2348 (2019).
  * (15) V. Stadnytskyi, C. E. Bax, A. Bax and P. Anfinrud, PNAS, 117, 11875 (2020).
  * (16) K. K.-W. To et al., ”Consistent Detection of 2019 Novel Coronavirus in Saliva”, Clin. Infect. Diseases, 71, 841 (2020).
  * (17) L. Azzi et al., ”Saliva is a reliable tool to detect SARS-CoV-2”, Journal of Infection, 81, e45 (2020).
  * (18) M. Klompas, M. A. Baker and C. Rhee, ”Airborne transmission of SARS-COV2 - theoretical considerations and available evidences”, JAMA, 324, 441 (2020).
  * (19) S. K. Das, J. Alam, S. Plumari and V. Greco, Phys. Fluids 32, 097102 (2020); https://doi.org/10.1063/5.0022859.
  * (20) T. Dbouk and D. Drikakis, Phys. Fluid 32, 053310 (2020)
  * (21) R. Mittal, R. Ni and J.-H. Seo, arXiv:2004.09354[physics-fluid.dyn].
  * (22) M.-R. Pendar and J. C. Páscoa, ”Numerical modeling of the distribution of virus carrying saliva droplets during sneeze and cough”, Phys. Fluids 32, 083305 (2020); doi: 10.1063/5.0018432.
  * (23) Reif, F., Fundamentals of Statistical and Thermal Physics, Mcgraw-Hill International Editions, Singapore, 1985.
  * (24) Pathria, R. K., Statistical Mechanics, Butterworth-Heinemann, Oxford, 1996.
  * (25) W. H. Seta, J. Hospital Infection, 89, 225 (2015).
  * (26) Kukkonen, J., Vesala, T., and Kulmala, M., “The interdependence of evaporation and settling for airborne freely falling droplets,” J. Aerosol Sci. 20, 749–763 (1989)
  * (27) Holterman, H., Kinetics and Evaporation of Water Drops in Air (Citeseer, 2003), Vol. 2012.
  * (28) L. D.Landau and E. M. Lifshitz, Fluid Mechanics, Butterworth-Heinemann, Oxford, UK (2005).
  * (29) S. K. Das, F. Scardina, S. Plumari and V. Greco, Phys. Rev. C 90, 044901 (2014).
  * (30) S. K. Das, S. Plumari, S. Chatterjee, J. Alam, F. Scardina and V. Greco, Phys. Lett. B 768, 260 (2017).
  * (31) Vansciver, M., Miller, S. and Hertzberg, J., Particle image velocimetry of human cough, Aerosol Sci. Technol. 45, 415 (2011).
  * (32) S. Zhu, S. Kato, and J.-H. Yang, “Study on transport characteristics of saliva droplets produced by coughing in a calm indoor environment,” Build. Environ. 41, 1691 (2006).
  * (33) Han, Z. Y., Weng, W. G. and Huang, Q. Y., Characterizations of particle size distribution of the droplets exhaled by sneeze, Jour. Roy. Soc. Interface 10, 20130560 (2013).
  * (34) H. Wang, Z. Li, X. Zhang, L. Zhu, Y. Liu and S. Wang, ”The motion of respiratory droplets produced by coughing” Phys. Fluids bf32, 125102 (2020); doi: 10.1063/5.0033849.
  * (35) S. Chatterjee, J. S. Murallidharan, A. Agrawal and R. Bhardwaj, ”How coronavirus survives for hours in aerosols”, Phys. Fluids 33, 081708 (2021); doi: 10.1063/5.0059908.
  * (36) H. Li, F. Y. Leong, G. Xu, Z. Ge, C. W. Kang and K. H. Lim, ”Dispersion of evaporating cough droplets in tropical outdoor environment”, Phys. Fluids 32, 113301 (2020); doi: 10.1063/5.0026360.
  * (37) T. Dbouk and D. Drikakis, ”Fluid dynamics and epidemiology: Seasonality and transmission dynamics”, Phys. Fluids 33, 021901 (2021); doi: 10.1063/5.0037640.
  * (38) T. Dbouk and D. Drikakis, ”Weather impact on airborne coronavirus survival” Phys. Fluids 32, 093312 (2020); https://doi.org/10.1063/5.0024272.
  * (39) S. W. X. Ong, Y. K. Tan, P. Y. Chia, T. H. Lee, O. T. Ng, M. S. Y. Wong, and K. Marimuthu, Air, surface environmental, and personal protective equipment contamination by severe acute respiratory syndrome coronavirus 2 (SARS- CoV-2) from a symptomatic patient, JAMA March 4 (2020).
  * (40) L. Bourouiba, E. Dehandschoewercker and J. W. M. Bush, J. Fluid Mechanics, 745, 537 (2014).
  * (41) Liu, Yuan and Ning, Zhi and Chen, Yu and Guo, Ming and Liu, Yingle and Gali, Nirmal Kumar and Sun, Li and Duan, Yusen and Cai, Jing and Westerdahl, Dane and Liu, Xinjin and Xu, Ke and Ho, Kin-fai and Kan, Haidong and Fu, Qingyan and Lan, Ke, Aerodynamic analysis of SARS-CoV-2 in two Wuhan hospitals, Nature, 27 April 2020; doi.org/10.1038/s41586-020-2271-3.
  * (42) F. Liu et. al., IOP Conf. Series: Materials Science and Engineering 609, 042059IOP (2019). Doi:10.1088/1757-899X/609/4/042059.
  * (43) Z. Xu, T. Wang and Z. Che, Physics Fluid 32, 052109 (2020).
  * (44) J.Howard, A. Huang, Z. Li, Z. Tufekci, V. Zdimal, H. van der Westhuizen, A. von Delft, A. Price, L. Fridman, L. Tang, V. Tang, G. L. Watson, C. E. Bax, R. Shaikh, F. Questier, D. Hernandez, L. F. Chu, C. M. Ramirez, A. W. Rimoin, Face Masks Against COVID-19: An Evidence Review. Preprints 2020, 2020040203 (doi: 10.20944/preprints202004.0203.v1).
  * (45) http://blog.pnas.org/2020/04/fluid-dynamics-work-hints-at-whether-spoken-word-can-spread-covid-19/.
  * (46) S. Verma, M. Dhanak, J. Frankenfield, Physics of Fluid, 32, 061708 (2020).
  * (47) L. Bourouiba, JAMA, 323, 1837 (2020).
  * (48) N. Sen, ”Transmission and evaporation of cough droplets in an elevator: Numerical simulations of some possible scenarios”, Phys. Fluids 33, 033311 (2021); doi: 10.1063/5.0039559.
  * (49) S. Bathula, S. Anand, T. Thajudeen, Y. S. Mayya, P. Chaudhury and S. Chaturvedi, ”Survival of Expiratory Aerosols in a Room: Study Using a Bi-compartment and Bi-component Indoor Air Model”, Aerosol and Air Quality Research, 21, 200547 (2020); https://aaqr.org/articles/aaqr-20-09-covid-0547.