# Spin/Phonon Dynamics in Single Molecular Magnets: II. spin/phonon
entanglement

Nosheen Younas1,2,3, Yu Zhang2, Andrei Piryatinski2, Eric R. Bittner1,3 1
Department of Physics, University of Houston, Houston, Texas 77204, USA
2Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico
87545, United States
3Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New
Mexico 87545, United States

###### Abstract

We introduce a new quantum embedding method to explore spin-phonon
interactions in molecular magnets. This technique consolidates various
spin/phonon couplings into a limited number of collective degrees of freedom,
allowing for a fully quantum mechanical treatment. By precisely factorizing
the entire system into ”system” and ”bath” sub-ensembles, our approach
simplifies a previously intractable problem, making it solvable on modest-
scale computers. We demonstrate the effectiveness of this method by studying
the spin relaxation and dephasing times of the single-molecule qubit
VOPc(OH)8, which features a lone unpaired electron on the central vanadium
atom. By using this mode projection method, we are able to perform numerical
exact quantum dynamical calculation on this system which allows us to follow
the flow of quantum information from the single spin qubit into the projected
phonon degrees of freedom. Our results demonstrate both the utility of the
method and suggest how one can engineer the environment as to further optimize
the quantum properties of a qubit system.

##  I Introduction

Single molecular magnets (SMMs) and metal-organic frameworks (MOFs) have
garnered considerable interest because of their potential in quantum
information processing and their ability to maintain extended coherence times.
However, these systems experience spin dephasing because of interactions and
couplings with the molecular framework’s vibrational motions.

Amongst the various MOF architectures, vanadyl phthalocyanine (VOPc) is a
prominent candidate for room-temperature quantum computing using the unpaired
$d_{xy}$ as a prototypical qubit owing to its long $T_{1}$ and $T_{2}$
relaxation times, which are largely determined by weak interaction between the
spin and vibrational degrees of freedom. [1, 2, 3, 4] While first-principles
calculations, particularly density functional theory (DFT), have become
essential in studying the electronic and magnetic properties of molecular
magnets [5, 6] challenges remain in terms of accurately modeling the dynamical
processes due to the high computational cost and complexity of spin-phonon
interactions, which involve a large number of degrees of freedom. Even though
the spin/phonon coupling is weak and can be treated as a perturbation, it is
interesting to ascertain the mechanism by which energy and quantum information
flows irreversibly from the central spin into the environment.

In our accompanying paper, Ref.7, we introduced a systematic
projection/embedding scheme to analyze spin-phonon dynamics in molecular
magnets. By consolidating all spin/phonon couplings into a few collective
degrees of freedom that can be treated fully quantum mechanically, this
approach transforms a numerically intractable problem into one manageable one
on modest-scale computers using an exact factorization of the complete system
into “system” and ”bath” sub-ensembles in which very few collective modes of
the bath couple directly to the quantum system. The projected modes themselves
carry physical importance, since all of the spin/phonon coupling is
concentrated in these collective degrees of freedom. Using parameters derived
from ab initio methods, we applied this approach to calculate the electronic
spin relaxation and dephasing times for the single-molecule qubit VOPc(OH)8,
which contains a single unpaired electron on the central vanadium atom. This
development provides a useful tool for simulating spin relaxation in complex
environments with significantly reduced computational complexity.[8, 9, 10,
11]

Although the focus of Ref. 7 was in the identification of the optical modes
and the computation of the $T_{1}$ and $T_{2}$ relaxation and dephasing times
for VoPC, here we extend the “system” to include the bath modes identified by
the projection method within a unified spin + phonon dynamical scheme in order
to map the flow of quantum information from the spin qubit into the
environment. We find for this system that the spin rapidly becomes entangled
with at least one of the collective modes identified by our approach and that
this mode provides a conduit for energy and quantum information loss from the
qubit. We speculate that if spin/phonon coupling can be manipulated by
chemical modification of the surrounding molecular scaffold, it may be
possible to enhance the $T_{1}$ and $T_{2}$ times of other SMM systems.

##  II Theory

###  II.1 System/bath decomposition

We begin with a brief summary of the approach developed in Ref. 7. Starting
from the Zeeman Hamiltonian describing the spin-orbit coupling for an electron
in a molecular orbital

$\displaystyle
H=\beta\hat{\mathbf{S}}\cdot\mathbf{\stackrel{{\scriptstyle\leftrightarrow}}{{g}}}\cdot\mathbf{B}+\frac{1}{2}\sum_{n}^{3N-6}p_{n}^{2}+\omega_{n}^{2}x_{n}^{2}$
(1)

The first term is the usual Zeeman interaction between an external magnetic
field $\bf B$ and electronic spin $\bf\hat{S}$ and $\beta=\mu_{B}/\hbar$. The
two are coupled via the g-tensor,
$\stackrel{{\scriptstyle\leftrightarrow}}{{g}}$, defined at some optimized
molecular geometry of the molecular system. Likewise, the $3N-6$ normal modes
$\\{x_{n}\\}$ and associated frequencies $\\{\omega_{n}\\}$ are determined at
the optimized molecular geometry. The coupling between the spin and nuclear
motions enters by expanding the g-tensor in terms of nuclear displacements
along the normal modes. Although nuclear displacements do not appear directly
in this expression, $\bm{g}$ depends upon the local molecular geometry via
second-order perturbation theory

$\displaystyle
g_{\alpha\beta}=g_{o}+\sum_{e}\frac{\langle\psi_{g}|\hat{L}_{\alpha}|\psi_{e}\rangle\langle\psi_{e}|\hat{L}_{\beta}|\rangle}{E_{e}-E_{g}},$
(2)

where $g_{o}$ is the g-value for a free electron and $\hat{L}_{\alpha}$ are
the orbital angular momentum operators coupling the ground electronic state
$|\psi_{g}\rangle$ to the unoccupied electronic states $|\psi_{e}\rangle$
taken at the equilibrium nuclear geometry. Of the terms entering this
expression, the energy gap is most sensitive to fluctuations in equilibrium
geometry. Expanding about the equilibrium geometry along the normal
coordinates gives

$\displaystyle
H=\underbrace{\left(\beta\hat{\mathbf{S}}\cdot\mathbf{\stackrel{{\scriptstyle\leftrightarrow}}{{g}}}\cdot\mathbf{B}\right)_{o}}_{\rm
Zeeman\,term}+\underbrace{\sum_{\alpha,n}g^{\prime}_{\alpha
n}\hat{\sigma}_{\alpha}x_{n}}_{\rm
Spin/phonon\,interaction}+\underbrace{\frac{1}{2}\sum_{n}^{3N-6}p_{n}^{2}+\omega_{n}^{2}x_{n}^{2}}_{\rm
normal\,modes}$ (3)

where $g^{\prime}_{\alpha n}$ denotes the derivative coupling between the
$\alpha=x,y,z$ spin operator and the $n$-th normal model.

In our previous work[7], we further partition the normal modes into “primary”
and residual modes by finding collective degrees of freedom that optimize the
spin / phonon coupling by performing a single value decomposition in
$g^{\prime}$ viz.

$\bm{g}^{\prime}=[g^{\prime}_{\alpha n}]_{n_{s}\times
n_{q}}=\bm{U}_{n_{s}\times n_{s}}\cdot\bm{\Sigma}_{n_{s}\times
n_{q}}\cdot\bm{V}^{\dagger}_{n_{q}\times n_{q}}$ (4)

where $\bm{\Sigma}$ is a matrix of singular values with at most
$\min(n_{s},n_{q})$ nonzero singular values. The ${\bm{U}}$ and $\bm{V}$ allow
us define a series of canonical transformations which brings the full
spin+phonon Hamiltonian into the following form

$\displaystyle H\approx$
$\displaystyle\underbrace{\left(\beta\hat{\mathbf{S}}\cdot\mathbf{\stackrel{{\scriptstyle\leftrightarrow}}{{g}}}\cdot\mathbf{B}\right)_{o}}_{\rm
Zeeman\,term}+\underbrace{\sum_{\alpha,k\in{\cal P}}g^{\prime}_{\alpha
k}\hat{\sigma}_{\alpha}X_{k}}_{\rm Spin/phonon\,interaction}$
$\displaystyle+\overbrace{\underbrace{\sum_{k\in{\cal
P}}\frac{1}{2}\left(P_{k}^{2}+\omega_{k}^{2}X_{k}^{2}\right)}_{\rm
primary\,modes}+\underbrace{\sum_{q\in{\cal
Q}}\frac{1}{2}\left({P_{q}^{2}}+\omega_{q}^{2}Y_{q}^{2}\right)}_{\rm
residual\,modes}+\underbrace{\sum_{k\in{\cal P},q\in{\cal
Q}}\gamma_{kq}\hat{X}_{k}\hat{Y}_{q}}_{linear\,couplings}}^{\rm
normal\,modes}$ (5)

The last three terms in this expression define the Hamiltonian for all the
vibrational degrees of freedom, written in terms of the projected “primary”
and “residual” modes in which the “primary” modes are coupled directly to the
spin and the “residual” modes are only coupled to the “primary” modes, both of
which we write in terms of quantum position and momentum operators.

Formally, Eq.3 and Eq. 5 are identical. However, a full quantum treatment of
the system remains intractable. Consequently, we will restrict the quantum
evolution to the spin + primary $\\{X_{k}\\}$ modes and tracing over the
residual degrees of freedom and defining the $1_{k}\to 0_{k}$ relaxation rate
for the $k$-th mode as[12]

$\displaystyle\frac{1}{T_{vib,k}}$ $\displaystyle=\frac{1}{\hbar^{2}}|\langle
0_{k}|X_{k}|1_{k}\rangle|^{2}\sum_{q}\int_{-\infty}^{+\infty}d\tau
e^{i\omega_{k}\tau}\gamma_{kq}^{2}\langle Y_{q}(\tau)Y_{q}(0)\rangle$ (6)
$\displaystyle=\frac{1}{2\hbar\omega_{k}}\sum_{q}\int_{-\infty}^{+\infty}e^{i\omega_{k}\tau}\gamma_{kq}^{2}\langle
Y_{q}(\tau)Y_{q}(0)\rangle d\tau$ (7)

where

$\displaystyle\langle Y_{q}(t)Y_{q}(0)\rangle$
$\displaystyle=\frac{1}{Z_{Q}}{\rm Tr}_{Q}\left[e^{-\beta
H_{Q}}e^{+iH_{Q}t/\hbar}Y_{q}e^{-iH_{Q}t/\hbar}Y_{q}\right]$
$\displaystyle=\frac{\hbar}{2\omega_{q}}((1+n(\omega_{q}))e^{-i\omega_{q}t}+n(\omega_{q})e^{+i\omega_{q}t}),$
(8)

where $n(\omega)$ is the Bose-Einstein population for an oscillator with
frequency $\omega$. Integrating over time

$\displaystyle\frac{1}{T_{vib,k}}=\frac{\pi}{\hbar\omega_{k}}\sum_{q}\frac{\hbar\gamma_{kq}^{2}}{2\omega_{q}}((1+n(\omega_{q}))\delta(\omega_{k}-\omega_{q})+n(\omega_{q})\delta(\omega_{q}-\omega_{k}))$
(9)

The relaxation rates ($k_{i}$) for all three primary modes are presented in
Fig. 1 for a temperature range of 10 to 300 K. The results indicate that all
three modes have nearly identical decay rates at low temperatures. With
increasing temperature, all decay rates increase. At intermediate temperature
(around 100 K), modes 2 and 3 start to exhibit a higher gain in the decay
rate, and their decay rate surpasses that of mode one. This trend continues
until room temperature, with mode two showing the highest decay rate (hence,
the shortest $T_{1}$). This is due to populating bath modes that interact with
modes 2 and 3 relatively more strongly than mode 1.

Figure 1: Vibrational relaxation rates, $k_{i}$, for all three primary modes,
as given by Eq. 7. All three decay rates are nearly identical at low
temperatures. At intermediate temperature, modes 2 and 3 show increased gain
in decay rate that continues till 300K. Consequently, at 300 K, mode two shows
the highest decay rate, followed by modes 3 and 1.

###  II.2 Dynamics

We now consider the full quantum dynamics of a single spin coupled directly to
the projected modes, which in turn are coupled to a completely dissipative
environment. For this we, employ the Lindblad approach[13]

$\displaystyle\dot{\rho}(t)$
$\displaystyle=-\frac{i}{\hbar}[H_{s},\rho(t)]+\sum_{k,\alpha=\pm}\frac{1}{2}\left[2C_{k\alpha}\rho(t)C_{k\alpha}^{\dagger}-\\{\rho(t),C_{k\alpha}^{\dagger}C_{k\alpha}\\}\right]$
(10)

in which $H_{s}$ contains the Zeeman term for the spin and its coupling to the
three primary modes

$\displaystyle
H_{s}=\beta\sum_{\alpha=x,y,z}B_{\alpha}g_{\alpha}\hat{\sigma}_{\alpha}+\sum_{\alpha
k}g^{\prime}_{\alpha
k}\sigma_{\alpha}(a_{k}^{\dagger}+a_{k})+\sum_{k}\hbar\omega_{k}(a^{\dagger}_{k}a_{k}+1/2).$
(11)

Coupling to the residual modes is captured by the Lindblad relaxation
operators as given by

$\displaystyle
C_{k,-}=\sqrt{\frac{n(\omega_{k})+1}{T_{vib,k}}}a\,\,\,\&\,\,\,C_{k,+}=\sqrt{\frac{n(\omega_{k})}{T_{vib,k}}}a^{\dagger}$
(12)

corresponding to vibrational de-excitation and excitation of the $k$-th
primary modes due to the thermal fluctuations of the residual modes.

###  II.3 Isolated system dynamics

We first consider the dynamics in the absence of the residual bath. The
transition frequency between the two possible spin states is set by the
external magnetic field, typically between 1-10T Under these conditions, the
transition frequency between spin levels is much smaller compared to the
frequencies of the system phonons. Consequently, the spin evolves on a much
slower timescale than the phonons. This is advantageous for molecular qubits
since this detuning leads much longer coherence and relaxation times.
Additionally, the couplings ($g_{ik}$) are directly proportional to the
applied magnetic field; therefore, a lower field (e.g. 1 T) will result in a
significantly smaller spin-phonon interaction and quantum exchange. For
purposes of our model, we can increase the applied magnetic field to 200
Tesla, thus bringing the spin transition frequency to 185.51$cm^{-1}$. This is
close to the first system phonon frequency and allows us to examine the
quantum dynamics of this system and the role of spin/phonon entanglement.

###  II.4 Spin with Isolated Projected Modes

We construct a system with spin and three projected modes without the addition
of any interaction with the bath. Initially, the spin is excited, while all
projected phonon modes start in the ground state. In number state
representation,

$|m;\\{\ n_{1},\ n_{2},\
n_{3}\\}\rangle=|+\frac{1}{2};\\{0_{1},0_{2},0_{3}\\}\rangle.$ (13)

This setup allows for a controlled exploration of the dynamics and
interactions of the isolated system. Since the spin-phonon coupling is weak,
the populations of these states should remain close to their initial values.
Fig. 2 shows the results for this simulation with populations given in terms
of the deviation from their initial value, i.e. $\delta\rho=\rho(t)-\rho(0)$.
The population of the excited state for spin is shown in panel (a), while that
of the ground state of mode 1 is shown in panel (b). Both show nearly
identical dynamics, a consequence of the similar energies (or frequencies) of
these 2-level systems. The dynamics of modes 2 and 3, shown in panels (c)
and(d),) is markedly different. Furthermore, mode 1 shows beats (similar to
spin), mode 2 exhibits simple oscillations, and mode 3 undergoes high-
frequency oscillations on top of a low-frequency wave. Relatively, the
amplitude of oscillations for mode 1 is larger than for mode 2, with mode 3
displaying the smallest amplitude.

Figure 2: Population dynamics of spin and three projected modes represented
using $\delta\rho=\rho(t)-\rho(0)$. The initial state is given by
$|spin,phonon_{1},phonon_{2},phonon_{3}\rangle=|1,0,0,0\rangle$. Panel (a)
shows the excited state population of the spin. Panels (b), (c), and (d) show
ground state populations for modes 1, 2, and 3, respectively with
$\delta\rho=\rho_{00}-1$. Since the interaction between spin and phonons is
weak, these populations are very close to their initial value throughout the
evolution.

###  II.5 Thermalized Projected Modes

####  II.5.1 Including Spin

Finally, we simulate the evolution of the spin with all three projected modes,
including the effect of the thermal bath calculated previously. Initially, the
spin is excited, while the phonons are in their ground state. The combined
system evolved at 65 K, a typical temperature for liquid-nitrogen-cooled
experiments. The lifetimes of system modes at this temperature are $43.3$,
$127.1$, and $2879.1$ ps for mode 1, 2, and 3, respectively. The results are
presented in Fig. reffig: simulation4b, which shows the $rho_{1}1$ element of
spin and $rho_{0}0$ elements of modes 1 to 3 (marked by $p1,p2$, and $p3$).
Since modes 1 and 2 exhibit relatively short lifetimes compared to mode 3,
their populations reach thermodynamic equilibrium much faster than mode 3’s.
On the other hand, the spin and mode 3 are still decaying by the end of the
simulation.

Figure 3: Population dynamics for the spin and three system-phonons at 65 K.
The total system starts in $|1,0,0,0\rangle$ state. Consequently, excited
state population is shown for spin ($\rho_{11}^{spin}$) and ground state
populations are shown for phonons ($\rho_{00}^{pi}$). All populations are
represented by elements of density matrices for corresponding subsystems.
Notice that the populations of modes 1 and 2 have reached their thermodynamic
equilibrium but those of mode 3 and spin are still relaxing by the end of
simulation.  Figure 4: Population dynamics of the ground state for (a) mode 1,
(b) mode 2, and (c) mode 3, at 65 K starting from $|1,0,0,0\rangle$. The
exponential thermal decay for each mode has been removed individually for
clarity such that $\delta\rho=\rho(t)-\rho^{th}(t)$. We notice that these
dynamics are identical to those observed in Test Case 1, despite an overall
thermal decay.

At first glance, all three modes seem to be decaying thermodynamically,
following $\rho^{th}(t)$. However, intricate dynamics can be resolved by
removing this exponential decay profile induced by thermalization. We remove
the exponential decay profile for mode 1 by introducing,

$\delta\rho=\rho(t)-\rho^{th}(t),$

The results are presented in Fig. 4(a). This enhanced perspective highlights
that the first system phonon undergoes coherent oscillations primarily driven
by its interaction with the spin. These oscillations also indicate beats with
tens of picoseconds time period.

A similar evolution for mode 2 is presented in Fig. 4(b) indicating population
oscillations with a time period of about 200 ps. Likewise, Fig. 4(c) shows the
ground state evolution for mode 3, displaying fast oscillations at a period of
$\approx 20$ps. The dynamics of all three modes are qualitatively similar to
the case of an isolated system (Fig. 2), although thermalization is due to
interaction with the bath.

###  II.6 Spin-phonon Correlation

To quantify the degree of entanglement between spin and three projected modes,
we analyze mutual information between the spin and the phonon degrees of
freedom as the total system evolves. Mutual information measures correlation
between parts of a composite quantum mechanical system. Consider a bipartite
quantum system ($A\otimes B$), mutual information $I(A:B)$ reflects the extent
of shared information between the two subsystems $A$ and $B$. It is
mathematically defined as

$I(A:B)=S(A)+S(B)-S(A\otimes B)$ (14)

where, $S(A)$, $S(B)$, and $S(A\otimes B)$ represent von Neumann entropies of
the sub-systems $A$, $B$, and the composite system $A\otimes B$, respectively.
$I(A:B)$ serves as a measure of the correlations between components $A$ and
$B$, capturing both classical and quantum aspects. This measure is useful in
quantum information theory, as it indicates the level of entanglement and the
degree to which the state of one subsystem informs about the state of the
other. High mutual information signifies strong correlations, which can be
harnessed in quantum computing and quantum communication protocols to improve
performance and security. In the present context, we use $I(A:B)$ to measure
the degree of entanglement between the single spin qubit and each of the three
possible projected modes.

To investigate the correlation between the spin and a specific system mode, we
first trace over the remaining two modes and then apply the mutual information
formulation (Eq. 14) using von Neumann entropy. The results are shown in Fig.
5. Initially, the spin and phonon degrees of freedom are not entangled.
However, we can see that the first mode shows the most correlation with spin,
followed by the second and third modes. This means that the quantum
information “flows” from the spin state through mode 1 more readily than
through the other two modes. This relatively higher degree of correlation
between the spin and the first mode is reasonable since both are nearly
resonant, even though the first mode does not show the highest coupling with
spin via any Pauli operators.

Figure 5: Mutual information between the central spin and each of the three
phonon degrees of freedom. Mode 1 displays fast oscillations with slow beats,
which are absent in modes 2 and 3.

##  III Summary

The reduction of phonon degrees of freedom via SVD-mode projection has enabled
the full quantum treatment of the spin-phonon system. We investigate the spin-
phonon interaction by building models with and without a thermal bath. These
population oscillations preserve their character even in a thermal bath, a
signature of their quantum nature.

This work introduced a novel quantum embedding approach to study spin-phonon
dynamics in molecular magnets. Our method efficiently explores the spin-
lattice relaxation process with high accuracy. We applied this approach to
investigate the spin dynamics in single molecular magnets. Leveraging this
method, we have deepened our understanding of the physics underlying spin
relaxation in molecular magnets, which is essential for enhancing their
performance across various applications. Furthermore, this study lays the
groundwork for rationalising new materials with customized magnetic properties
suitable for quantum information processing and other applications.

As a future direction, we plan to extend the quantum embedding method to
encompass spin-spin interactions in molecular magnets. This expansion will
provide insights into the interplay between spin-spin and spin-phonon
interactions and elucidate the roles these interactions play in shaping the
overall relaxation dynamics of molecular qubits. Exploring spin-spin
interactions will deepen our understanding of these intricate systems and aid
in creating superior materials for quantum computing and other applications
involving magnetism. Such an extension promises to significantly widen the
quantum embedding approach’s scope and utility, enhancing its value as a
research tool in the molecular magnet domain.

## Acknowlegments

The research presented in this article was supported by the LANL LDRD program
(number 20220047DR). LANL is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the U.S. Department of Energy
(contract no. 89233218CNA000001). We thank the LANL Institutional Computing
(IC) program for access to HPC resources. The work at the University of
Houston was funded in part by the National Science Foundation (CHE-2404788)
and the Robert A. Welch Foundation (E-1337). This work was performed, in part,
at the Center for Integrated Nanotechnologies, an Office of Science User
Facility operated for the U.S. Department of Energy (DOE) Office of Science by
Los Alamos National Laboratory (Contract 89233218CNA000001) and Sandia
National Laboratories (Contract DE-NA-0003525).

## Data Availability Statement

The data supporting this study’s findings are available from the corresponding
author upon reasonable request.

## Author Contribution Statement

The authors confirm their contribution to the paper as follows: study
conception and design: ERB, YZ, and AP; data collection and simulations: NY;
analysis and interpretation of results: NY, YZ, AP, ERB; draft manuscript
preparation: NY, ERB. All authors reviewed the results and approved the final
version of the manuscript.

## References

  * [1] Atzori M, Tesi L, Morra E, Chiesa M, Sorace L, Sessoli R. 2016 Room-Temperature Quantum Coherence and Rabi Oscillations in Vanadyl Phthalocyanine: Toward Multifunctional Molecular Spin Qubits. J Am Chem Soc 138, 2154–2157. (10.1021/jacs.5b13408)
  * [2] Bonizzoni C, Ghirri A, Atzori M, Sorace L, Sessoli R, Affronte M. 2017 Coherent coupling between Vanadyl Phthalocyanine spin ensemble and microwave photons: towards integration of molecular spin qubits into quantum circuits. Scientific Reports 7, 13096. (10.1038/s41598-017-13271-w)
  * [3] Malavolti L, Briganti M, Hänze M, Serrano G, Cimatti I, McMurtrie G, Otero E, Ohresser P, Totti F, Mannini M, Sessoli R, Loth S. 2018 Tunable Spin-Superconductor Coupling of Spin 1/2 Vanadyl Phthalocyanine Molecules. Nano Lett. 18, 7955–7961. (10.1021/acs.nanolett.8b03921)
  * [4] Nadeem S, Khan MN, Muhammad N, Ahmad S. 2019 Mathematical analysis of bio-convective micropolar nanofluid. Journal of Computational Design and Engineering 6, 233–242. (10.1016/j.jcde.2019.04.001)
  * [5] Neese F. 2006 Importance of direct spin spin coupling and spin-flip excitations for the zero-field splittings of transition metal complexes: a case study. J Am Chem Soc 128, 10213–10222. (10.1021/ja061798a)
  * [6] Timco GA, Carretta S, Troiani F, Tuna F, Pritchard RJ, Muryn CA, McInnes EJL, Ghirri A, Candini A, Santini P, Amoretti G, Affronte M, Winpenny REP. 2009 Engineering the coupling between molecular spin qubits by coordination chemistry. Nature Nanotechnology 4, 173–178. (10.1038/nnano.2008.404)
  * [7] Younas N, Zhang Y, Piryatinski A, Bittner ER. submitted Spin/Phonon Dynamics in Single Molecular Magnets: I. quantum embedding. Proc. Roy. Soc. A.
  * [8] Yang X, Bittner ER. 2015 Computing intramolecular charge and energy transfer rates using optimal modes. The Journal of Chemical Physics 142, 244114. (10.1063/1.4923191)
  * [9] Pereverzev A, Bittner ER, Burghardt I. 2009 Energy and charge-transfer dynamics using projected modes. The Journal of Chemical Physics 131, 034104. (10.1063/1.3174447)
  * [10] Thouin F, Srimath Kandada AR, Valverde-Chávez DA, Cortecchia D, Bargigia I, Petrozza A, Yang X, Bittner ER, Silva C. 2019 Electron–Phonon Couplings Inherent in Polarons Drive Exciton Dynamics in Two-Dimensional Metal-Halide Perovskites. Chemistry of Materials 31, 7085–7091. (10.1021/acs.chemmater.9b02267)
  * [11] Yang X, Keane T, Delor M, Meijer AJHM, Weinstein J, Bittner ER. 2017 Identifying electron transfer coordinates in donor-bridge-acceptor systems using mode projection analysis. Nature Communications 8, 14554\. (10.1038/ncomms14554)
  * [12] Kenkre VM, Tokmakoff A, Fayer MD. 1994 Theory of vibrational relaxation of polyatomic molecules in liquids. The Journal of Chemical Physics 101, 10618–10629. (10.1063/1.467876)
  * [13] Lindblad G. 1976 On the generators of quantum dynamical semigroups. Communications in Mathematical Physics 48, 119–130. (10.1007/BF01608499)


# A Fast, Semi-analytical Model for the Venusian Binary Cloud System

L.-K. Dai1, X. Zhang2, and J. Cui1,3
1Planetary Environmental and Astrobiological Research Laboratory (PEARL),
School of Atmospheric Sciences, Sun Yat-sen University, Zhuhai, Guangdong,
People’s Republic of China
2Department of Earth and Planetary Sciences, University of California Santa
Cruz, Santa Cruz, CA, USA
3Center for Excellence in Comparative Planetology, Chinese Academy of
Sciences, Hefei, Anhui, People’s Republic of China E-mail:
<EMAIL_ADDRESS>
###### Abstract

The Venusian clouds originate from the binary condensation of H2SO4 and H2O.
The two components strongly interact with each other via chemistry and cloud
formation. Previous works adopted sophisticated microphysical approaches to
understand the clouds. Here we show that the observed vapor and cloud
distributions on Venus can be well explained by a semi-analytical model. Our
model assumes local thermodynamical equilibrium for water vapor but not for
sulfuric acid vapor, and includes the feedback of cloud condensation and
acidity to vapor distributions. The model predicts strong supersaturation of
the H2SO4 vapor above 60 km, consistent with our recent cloud condensation
model. The semi-analytical model is 100 times faster than the condensation
model and 1000 times faster than the microphysical models. This allows us to
quickly explore a large parameter space of the sulfuric acid gas-cloud system.
We found that the cloud mass loading in the upper clouds has an opposite
response of that in the lower clouds to the vapor mixing ratios in the lower
atmosphere. The transport of water vapor influences the cloud acidity in all
cloud layers while the transport of sulfuric acid vapor only dominates in the
lower clouds. This cloud model is fast enough to be coupled with the climate
models and chemistry models to understand the cloudy atmospheres of Venus and
Venus-like extra-solar planets.

###### keywords:

Planets and satellites: atmospheres – Planets and satellites: terrestrial
planets

††pubyear: 2022††pagerange: A Fast, Semi-analytical Model for the Venusian
Binary Cloud System–References

##  1 Introduction

The globally distributed clouds on Venus are mainly composed of concentrated
sulfuric acid (H2SO4, Esposito et al., 1983; Krasnopolsky, 1985). Their
formation is governed by complex microphysical processes (Gao et al., 2014;
Imamura & Hashimoto, 2001; James et al., 1997; McGouldrick, 2017; McGouldrick
& Toon, 2007; Parkinson et al., 2015; Rossow, 1977; Rossow & Gierasch, 1977),
such as condensation, evaporation, nucleation, coagulation, as well as
coalescence. Moreover, the concentration of sulfuric acid and water inside the
cloud droplets determines the cloud acidity, which in turn strongly impacts
the saturation vapor pressure of both species. Thus, all four components (the
vapor and condensed phases of both sulfur acid and water) need to be solved
simultaneously, further complicating the calculation. A recent model (Dai et
al., 2022, hereafter D22) found that the condensation and evaporation
processes are sufficient to explain the observed vertical distributions of
sulfuric acid and water vapor mixing ratios, cloud mass loading of the Mode-2
particles, and cloud acidity on Venus. The success of the condensation model
implies that the vapor and cloud distribution simulation in previous
microphysical models can be greatly simplified. However, all the
aforementioned models, including D22, adopted the time-stepping method to
solve the cloud formation process to achieve the steady state. With reasonable
assumptions, we present here a semi-analytical model to directly obtain the
steady-state solution of the cloud system on Venus without time stepping,
allowing the orders-of-magnitude increase of the computational efficiency.

Our work is built on a seminal study in Krasnopolsky & Pollack (1994) that was
later updated by Krasnopolsky (2015). We call their model KP94 hereafter. The
KP94 model is a semi-analytical model without time stepping. It mainly focused
on H2SO4 and H2O vapors in the middle and lower cloud region (below 60 km).
The aerosol mass can be derived after the vapor distributions are obtained.
There are two important assumptions in KP94: hydrogen flux conservation and
local thermodynamic equilibrium (LTE). The sulfuric acid vapor is formed via
cloud-top chemistry by sulfuric oxides and water vapor that are diffused from
below. The first assumption thus states that the net upward hydrogen flux in
the water (including both vapor and droplets) should be balanced by the net
downward hydrogen flux in the sulfuric acid at any given altitude. The second
assumption states that both sulfuric acid and water vapor follow the
saturation vapor curves that are determined by the temperature and local cloud
acidity. With these two assumptions, the KP94 model can directly solve for the
cloud acidity and the vapor distribution given the temperature profile. The
model results are consistent with the observed vapor and acidity below 60 km.

The KP94 model can only simulate the region below the sulfuric acid production
(<60 km) where the LTE assumptions are generally valid. Krasnopolsky (2011)
tried to justify the LTE assumption up to 100 km by the comparison between
timescales of condensation and atmospheric mixing, but they did not take into
account the H2SO4 chemical production. In fact, the LTE assumption breaks down
for H2SO4 above 60 km. Using the production rates from various chemical models
(Krasnopolsky, 2012; Zhang et al., 2012; Shao et al., 2020), Dai et al. (2022)
simulated the cloud condensation and found that H2SO4 is highly supersaturated
above 60 km because the condensation is too slow to compete with the H2SO4
vapor production. Since the H2SO4 vapor is not in LTE, one must consider the
condensation process that critically depends on the cloud mass, the droplet
size, the vapor supersaturation, etc. But the KP94 model did not
simultaneously simulate the clouds and their feedback to the condensation of
H2SO4 vapor. On the other hand, Dai et al. (2022) confirmed that the water
vapor generally follows its saturation vapor profile in the upper cloud layer.
Thus, the LTE assumption remains valid for water.

In this study, we developed a semi-analytical model that extends to 80 km to
include the H2SO4 vapor production region. We kept the LTE assumption for the
water vapor as in KP94 but did not assume the LTE condition for the H2SO4
vapor. Instead, we explicitly included the cloud formation and its feedback to
the H2SO4 vapor abundances. Our model is more complicated than KP94 but we can
still obtain a steady-state solution without time stepping. Since both the
complicated microphysics and the time integration are avoided, the semi-
analytical model is significantly faster than the D22 model and previous full
microphysics models (Gao et al., 2014; McGouldrick, 2017; McGouldrick & Toon,
2007; Parkinson et al., 2015). In Section 2 we describe the equations and
numerical procedure of the model. In Section 3 we show the results of our
nominal simulation and explore a large parameter space of the H2SO4-H2O
system. We conclude this study in Section 4.

##  2 Model

###  2.1 Hydrogen Flux Conservation

Following KP94, the hydrogen is transported by vapor diffusion as well as
sedimentation and diffusion in both H2SO4 and H2O condensates. The vertical
fluxes are written as

$\displaystyle\Phi_{1}=-K_{zz}n_{atm}\frac{dq^{g}_{1}}{dz}+$ $\displaystyle
q^{c}_{1}n_{atm}v-$ (1)
$\displaystyle\frac{q^{c}_{1}}{q^{c}_{1}+q^{c}_{2}}K_{zz}n_{atm}\frac{d(q^{c}_{1}+q^{c}_{2})}{dz},$

$\displaystyle\Phi_{2}=-K_{zz}n_{atm}\frac{dq^{g}_{2}}{dz}+$ $\displaystyle
q^{c}_{2}n_{atm}v-$ (2)
$\displaystyle\frac{q^{c}_{2}}{q^{c}_{1}+q^{c}_{2}}K_{zz}n_{atm}\frac{d(q^{c}_{1}+q^{c}_{2})}{dz},$

where the subscripts $1$ and $2$ represent H2SO4 and H2O, respectively, the
superscripts $g$ and $c$ represent gas and condensate, respectively, $q$
represents the molar mixing ratio (with respect to the atmospheric density
$n_{atm}$), $z$ represents the altitude, $\Phi_{i}$ is the vertical net flux
of species $i$, $K_{zz}$ is the eddy diffusivity, $v$ is the sedimentation
velocity. Since Dai et al. (2022) suggested that the condensation of mode-2
particles was able to explain the main observed cloud properties including its
acidity and abundance, we neglect the particle size distribution in our simple
model here. At any altitude, the single-mode cloud particles share the same
acidity and diameter, and consequently the same sedimentation velocity. The
standard vertical profiles of $K_{zz}$ and $n_{atm}$ for Venus are shown in
Fig. 1.

Since the chemical reaction SO3+2H2O=H2SO4+H2O is the main source of H2SO4 and
the main sink of H2O (Bierson & Zhang, 2020; Krasnopolsky, 2012; Krasnopolsky
& Parshev, 1981; Mills, 1998; Shao et al., 2020; Winick & Stewart, 1980; Yung
& Demore, 1982; Parkinson et al., 2015), the vertical transports of the
hydrogen is conserved. We have:

$\Phi_{1}=-\int^{\infty}_{z}P_{1}(z)dz=-\Phi_{2},$ (3)

where $P_{1}$ is the chemical net production rate of H2SO4 adopted from the
chemical model in Krasnopolsky (2012). Both the standard vertical profiles of
$P_{1}$ and $-\Phi_{1}$ are shown in Fig. 1.

Figure 1: The standard vertical profiles of the background atmospheric
properties in our model: (a) the temperature profile (red) from Seiff et al.
(1985) and the eddy diffusivity profile (blue solid line) from Dai et al.
(2022). The eddy diffusivity profiles tested in Section 3.2 are represented as
the blue dashed lines. We increase (or decrease) the standard diffusivity
profile below 54 km by varying the diffusivity at 40 km to 2000, 4000, 6000,
and 8000 cm2 s-1, respectively, while no variation is made above 54 km; (b)
the atmospheric number density profile (red) from Krasnopolsky (2012) and the
cloud condensation nuclei (CCN) number density profile (blue) from Gao et al.
(2014). The observations of CCN number density (grey pluses) are from the
cloud particle size spectrometer (LCPS) aboard the Pioneer Venus probes
(Knollenberg & Hunten, 1980); (c) the H2SO4 chemical net production rate
profile (red) from Krasnopolsky (2012) and the derived H2SO4 net flux profile
(blue). Minus indicates a downward flux.

###  2.2 LTE for H2O and Non-LTE for H2SO4

The LTE condition states that the vapor should follow its saturation vapor
mixing ratio (SVMR, $q^{svp}_{i}$) in the cloud layer. The SVMR depends on the
local temperature and the cloud acidity (see Section 2.4). As shown in
Knollenberg & Hunten (1980), H2O condensation is sufficiently fast to keep its
curve following the SVMR and maintain an LTE condition such that:

$q^{g}_{2}=q^{svp}_{2}.$ (4)

But for H2SO4, LTE is not valid above 60 km, where the chemical production
rate is high (Knollenberg & Hunten, 1980). The vapor abundance is mainly
controlled by the vapor condensation and the chemical production, while the
vapor diffusion is negligible (Dai et al., 2022). To account for this non-LTE
effect, we introduce a vapor condensation rate that balances the net chemical
production of the H2SO4 vapor:

$P_{1}=C_{1}(q^{g}_{1}-q^{svp}_{1}),$ (5)

where $C_{1}$ is the condensation rate coefficient that critically depends on
the cloud particle size and their abundances (see Section 2.3). The term
$P_{1}/C_{1}=q^{g}_{1}-q^{svp}_{1}$ quantifies the supersaturation of the
H2SO4 vapor.

We introduce the H2O/H2SO4 molar ratio in the droplets
$m=q^{c}_{2}/q^{c}_{1}$. Substituting equations (3-5) to equations (1-2), we
can obtain:

$\displaystyle-\int^{\infty}_{z}P_{1}(z)dz=$ $\displaystyle-
K_{zz}n_{atm}\frac{d\left(q^{svp}_{1}+\frac{P_{1}}{C_{1}}\right)}{dz}+q^{c}_{1}n_{atm}v-$
(6) $\displaystyle\frac{1}{1+m}K_{zz}n_{atm}\frac{d[(1+m)q^{c}_{1}]}{dz},$

$\displaystyle\int^{\infty}_{z}P_{1}(z)dz=$ $\displaystyle-
K_{zz}n_{atm}\frac{dq^{svp}_{2}}{dz}+q^{c}_{2}n_{atm}v-$ (7)
$\displaystyle\frac{m}{1+m}K_{zz}n_{atm}\frac{d\left[\left(1+\frac{1}{m}\right)q^{c}_{2}\right]}{dz}.$

Note that $C_{1}$, $q^{svp}_{i}$, $v$, and $m$ are functions of $q^{c}_{1}$
and $q^{c}_{2}$ (see Section 2.3). Given $P_{1}$, $K_{zz}$ and $n_{atm}$, one
can solve for $q^{c}_{1}$ and $q^{c}_{2}$ in (6) and (7) without time
stepping.

The condensate settling and diffusional fluxes can be eliminated in (6) and
(7), which further simplifies the calculation. The governing equation of $m$
can be analytically obtained:

$\displaystyle
m=\frac{K_{zz}n_{atm}\frac{dq^{svp}_{2}}{dz}+\int^{\infty}_{z}P_{1}(z)dz}{K_{zz}n_{atm}\frac{d\left(q^{svp}_{1}+\frac{P_{1}}{C_{1}}\right)}{dz}-\int^{\infty}_{z}P_{1}(z)dz}.$
(8)

The main difference between our equation and equation (4) in KP94 is the
supersaturation term $P_{1}/C_{1}$ that is related to the H2SO4 condensation.
If the condensation is very fast, $C_{1}$ approaches infinity and the equation
(5) is reduced to that in the KP94 model. The detailed procedure to solve the
system will be given in Section 2.6.

###  2.3 Condensation Rate Coefficient, Cloud Particle Size and Settling
Velocity

The condensation rate coefficient of H2SO4 is adopted from Knollenberg &
Hunten (1980):

$\displaystyle C_{1}=\frac{2\pi
n_{p}D_{1}M_{1}n_{atm}D_{p}}{(M_{1}+M_{2}m)\left(1+\frac{2\lambda_{1}}{D_{p}\alpha}\right)},$
(9)

where $n_{p}$ is the cloud condensation nuclei (CCN, Fig. 1) number density
adopted from Gao et al. (2014), $M_{1}=98$ g mol-1 and $M_{2}=18$ g mol-1 are
the relative molecular weights, $D_{p}$ is the particle diameter. $\alpha=1$
is the accommodation coefficient (Seinfeld & Pandis, 2016, p. 500-502),
$D_{1}$ is the molecular diffusivity of H2SO4 (Hamill et al., 1977):

$\displaystyle D_{1}=\frac{\lambda_{1}}{3}\left(\frac{8kTN_{a}}{\pi
M_{1}}\right)^{\frac{1}{2}},$ (10)

where $k$ is the Boltzmann constant, $N_{a}$ is the Avogadro constant, $T$ is
the atmospheric temperature (Fig. 1) adopted from the observation of Pioneer
Venus at 45° latitude (Seiff et al., 1985), $\lambda_{1}$ is the mean free
path of the H2SO4 vapor (Hamill et al., 1977):

$\displaystyle\lambda_{1}=\left[\pi
n_{atm}\left(\frac{d_{1}+d_{atm}}{2}\right)^{2}\left(\frac{M_{atm}}{M_{1}+M_{atm}}\right)^{\frac{1}{2}}\right]^{-1},$
(11)

Where $d_{1}=4.4\times 10^{-8}$ cm (assumed in Hanson & Eisele, 2000) is the
molecular diameter of H2SO4, $d_{atm}=3.3\times 10^{-8}$ cm is the molecular
diameter of the atmospheric major component (carbon dioxide, CO2, Mehio et
al., 2014), $M_{atm}=44$ g mol-1 is the relative molecular weight of CO2.

The particle size is calculated given the condensate molar mixing ratios
$q^{c}_{1}$ and $q^{c}_{2}$:

$\displaystyle D_{p}=\left[\frac{6(M_{1}q^{c}_{1}+M_{2}q^{c}_{2})n_{atm}}{\pi
n_{p}\rho_{p}N_{a}}\right]^{\frac{1}{3}}.$ (12)

The particle sediments with the Stokes velocity:

$\displaystyle v=-\frac{2}{9}\frac{g\rho_{p}{D_{p}}^{2}}{4\eta}C_{c},$ (13)

where $g=870$ cm s-2 is the gravitational acceleration on Venus, $\eta$ is the
dynamic viscosity of CO2 gas that varies with temperature (Fenghour et al.,
1998), $C_{c}$ is the Cunningham slip correction factor following Seinfeld &
Pandis (2016, p. 371-372):

$\displaystyle
C_{c}=1+\frac{2\lambda_{1}}{D_{p}}\left(1.257+0.4e^{\frac{-1.1D_{p}}{2\lambda_{1}}}\right).$
(14)

###  2.4 Saturation Vapor Mixing Ratios

The SVMRs ($q^{svp}_{i}$) of both H2O and H2SO4 are calculated by:

$\displaystyle\ln P^{svp}_{i}$ $\displaystyle=\ln(q^{svp}_{i}kTn_{atm})$ (15)
$\displaystyle=\ln
P^{svp}_{i,p}+\frac{\mu_{i}-\mu^{0}_{i}}{R_{g}T}+\frac{4\sigma
M_{i}}{R_{g}T\rho_{p}D_{p}},$

where $P^{svp}_{i}$ is the saturation vapor pressure (SVP, in unit of bar) of
species $i$ over the sulfuric acid solution, $P^{svp}_{i,p}$ is the SVP (in
unit of bar) over pure species $i$, $\mu_{i}-\mu^{0}_{i}$ is the difference of
the chemical potential of species $i$ between the solution and its pure
condensate. This chemical potential difference varies with temperature and
acidity, and the data of it in this study is derived from Zeleznik (1991). The
last term comes from the Kelvin effect (Seinfeld & Pandis, 2016, p.419-423),
where $\sigma=71.11$ erg cm-2 is the surface tension of sulfuric acid (Myhre
et al., 1998, 80% sulfuric acid, 273 K), $\rho_{p}=1.8$ g cm-3 is the density
of the droplets. The $P^{svp}_{1,p}$ is from Kulmala & Laaksonen (1990):

$\displaystyle\ln P^{svp}_{1,p}=$ $\displaystyle 16.259-\frac{10156}{T_{0}}+$
(16) $\displaystyle
10156\left[-\frac{1}{T}+\frac{1}{T_{0}}+\frac{0.38}{T_{c}-T_{0}}\left(1+\ln\frac{T_{0}}{T}-\frac{T_{0}}{T}\right)\right],$

where $T_{0}=360$ K is a reference temperature, $T_{c}=905$ K is the critical
temperature. The $P^{svp}_{2,p}$ is from Tabazadeh et al. (1997):

$\displaystyle\ln
P^{svp}_{2,p}=e^{18.4524-\frac{3505.2}{T}-\frac{330918.6}{T^{2}}+\frac{12725068.3}{T^{3}}}.$
(17)

###  2.5 Cloud base Determination

Unlike the microphysical models, our simple model (and the KP94 model)
requires a known cloud base. We follow the definition of the cloud base
altitude in KP94 but use a different method to find it. The cloud base
altitude is where both the H2SO4 and H2O vapors start to condense together.
Below the cloud base, there is no condensate, and the H2SO4 and H2O vapors are
only controlled by the eddy diffusion and their chemical production/loss in
the upper clouds. Given their mixing ratios at the lower boundary well below
the cloud base (e.g., 40 km), the chemical-diffusion equilibrium vapor
profiles are shown in Fig. 2 ($q^{g\textquoteright}_{1}$ for H2SO4 and
$q^{g\textquoteright}_{2}$ for H2O). At the cloud base, both vapor mixing
ratios should be equal to their SVMRs that are controlled by the cloud acidity
(or $m$). Mathematically, one needs to find an altitude where a single value
$m$ exists to satisfy both following conditions:

$\displaystyle q^{g^{\prime}}_{1}(z)=q^{svp}_{1}(m,z),$ (18)

$\displaystyle q^{g^{\prime}}_{2}(z)=q^{svp}_{2}(m,z).$ (19)

Note that $q^{svp}_{i}$ also depends on $z$ because the temperature changes
with altitude. Graphically, one can solve for $m$ in equations (18) and (19)
respectively. The vertical profiles of the solution are shown in Fig. 2
($m\textquoteright_{1}$ for H2SO4 and $m\textquoteright_{2}$ for H2O). The
intersection between these two H2O/H2SO4 molar ratio profiles denotes the
altitude of the cloud base. One can show that the cloud base altitude obtained
in our method is fully consistent with that in KP94 but our method might be
more intuitively understandable.

Figure 2: The definition of the cloud base. The H2SO4 ($q^{g\textquoteright}_{1}$, red) and H2O ($q^{g\textquoteright}_{2}$, yellow) vapor profiles are only controlled by their chemical production/loss and diffusion. The $m\textquoteright_{1}$ profile (green) is a solution of equation (18) and the $m\textquoteright_{1}$ profile (blue) is a solution of equation (19). Their intersection (black circle) represents the cloud base altitude. Table 1: Model parameters and boundary conditions Parameter | Symbol | Nominal model
---|---|---
Atmospheric temperature | $T$ | Red line in Fig. 1(a)
Eddy diffusivity | $K_{zz}$ | Solid blue line in Fig. 1(a)
Atmospheric density | $n_{atm}$ | Red line in Fig. 1(b)
CCN number density | $n_{p}$ | Blue line in Fig. 1(b)
H2SO4 chemical net production rate | $P_{1}$ | Red line in Fig. 1(c)
H2SO4 net flux | $\Phi_{1}$ | Blue line in Fig. 1(c)
Molecular diffusivity of H2SO4 | $D_{1}$ | equation (10)
Mean free path of H2SO4 | $\lambda_{1}$ | equation (11)
Dynamic viscosity of CO2 | $\eta$ | Fenghour et al. (1998)
Relative molecular weight of H2SO4 | $M_{1}$ | 98 g mol-1
Relative molecular weight of H2O | $M_{2}$ | 18 g mol-1
Relative molecular weight of CO2 | $M_{atm}$ | 44 g mol-1
Droplet density | $\rho_{p}$ | 1.8 g cm-3
Gravitational acceleration | $g$ | 870 cm s-2
Accommodation coefficient | $\alpha$ | 1
Surface tension of sulfuric acid | $\sigma$ | 71.11 erg cm-2
Chemical potential | $\mu$ | Zeleznik (1991)
Molecular diameter of H2SO4 | $d_{1}$ | 4.4$\times$10-8 cm
Molecular diameter of CO2 | $d_{atm}$ | 3.3$\times$10-8 cm
Upper boundary (80 km) |  |
H2SO4 vapor | $dq^{g}_{1}/dz$ | 0
H2O vapor | $dq^{g}_{2}/dz$ | 0
H2SO4 condensate | $dq^{c}_{1}/dz$ | 0
H2O condensate | $dq^{c}_{2}/dz$ | 0
Lower boundary (40 km) |  |
H2SO4 vapor | $q^{g}_{1}$ | 3 ppm
H2O vapor | $q^{g}_{2}$ | 27 ppm
H2SO4 condensate | $q^{c}_{1}$ | 0 ppm
H2O condensate | $q^{c}_{2}$ | 0 ppm

###  2.6 Boundary condition and iterative procedure for solution

This model extends from 40 to 80 km. Both vapors and condensates are assumed
to be zero fluxes at the upper boundary. The H2SO4 and H2O vapor mixing ratios
at the lower boundary of the model domain (40 km) are fixed at 3 and 27 ppm,
respectively, agreeing with the observations of Venus Express (Marcq et al.,
2008; Oschlisniok et al., 2012, 2021) and ground-based spectroscopy (de Bergh
et al., 1991; de Bergh et al., 1995; Pollack et al., 1993; Arney et al.,
2014). We conclude the parameters and boundary conditions in Table 1.

The actual lower boundary for the condensates is the cloud base defined in
Section 2.5. Different from the D22 model where the droplets precipitate below
the cloud base and evaporate at lower altitudes, this model assumes no
condensate below the cloud base. In other words, the droplets evaporate so
rapidly within an infinitesimal distance that we neglect their downward
diffusion at the cloud base. This assumption is fine because the droplet mass
flux at the cloud base is primarily dominated by the gravitational settling of
large particles. However, it does cause a small difference between our model
results and that in D22, as shown in Section 3.

Our procedure to solve the system follows iterative steps below:

  1. 1.

Find out the cloud base altitude using the method in Section 2.5;

  2. 2.

Given $C_{1}$ (for the initial $C_{1}$, one can use infinity), solve for $m$
in equation (8) from the cloud base using the central difference scheme and
the scant method;

  3. 3.

For each altitude, calculate $q^{g}_{1}$ and $q^{g}_{2}$ (equations 4-5) based
on $m$ and $C_{1}$ in step (ii);

  4. 4.

Given $D_{p}$ (for the initial $D_{p}$, one can assume 1 µm), derive
sedimentation velocity $v$ using equation (13);

  5. 5.

Given the derived $m$ in step (ii), $q^{g}_{i}$ in step (iii), and $v$ in step
(iv), and assuming no particle diffusion at the cloud base, we can integrate
equations (6) and (7) from the cloud base to find out the profiles of
$q^{c}_{1}$ and $q^{c}_{2}$ respectively;

  6. 6.

Given $q^{c}_{1}$ and $q^{c}_{2}$, update $D_{p}$ using equation (12);

  7. 7.

Iterate steps (iv-vi) to get a converged $D_{p}$ profile;

  8. 8.

Given $D_{p}$ in step (vii), update $C_{1}$ using equation (9);

  9. 9.

Iterate steps (ii-viii) to find a converged solution of $m$, $q^{g}_{i}$ and
$q^{c}_{i}$.

The model reaches the steady state when the maximum variation of $m$ profile
is less than 0.001% at all altitudes. Using Python, this model is able to
reach the steady state within 10 seconds. The vertical resolution is 0.5 km
for the runs in this study, the same as D22 model. But the model behavior
remains the same in a higher resolution.

##  3 Results

Figure 3: The results of our nominal simulation and the comparisons to
observations and KP94 and D22 models. (a) the volume mixing ratios of H2SO4
vapor (black solid line), condensate (black dashed line), and SVMR (black
dotted line). The results from D22 and KP94 models are represented by red and
blue lines, respectively, with the same line styles. The H2SO4 vapor profile
above 60 km in KP94 is hardly retrieved from their linear-scale plots. Because
their vapor is assumed in LTE, here we calculate the profile using their
acidity and the thermodynamical data from Giauque et al. (1960). The
observations of the H2SO4 vapor mixing ratio are from Magellan Mission (brown,
Kolodner & Steffes, 1998), Venus Express (light green, Oschlisniok et al.
2012; pink, Oschlisniok et al. 2021), and Akatsuki Mission (gray, Imamura et
al. 2017); (b) the volume mixing ratios of H2O vapor, condensate, SVMR, as
well as comparisons to D22 and KP94 models. The observations are from Venus
Express (orange, Bertaux et al. 2007; light green, Fedorova et al. 2008) and
ground-based spectroscopy (cyan, de Bergh et al. 1991; dark cyan, de Bergh et
al. 1995; gray, Encrenaz et al. 2013; yellow, Encrenaz et al. 2020; brown,
Gurwell et al. 2007; pink, Meadows & Crisp 1996; sky blue, Sandor & Clancy
2005; (c) the distributions of the cloud acidity in our nominal simulation,
D22 and KP94. The estimated values are from Venus Express (light green, Arney
et al. 2014; pink, Barstow et al. 2012; yellow, Cottini et al. 2012; sky blue,
McGouldrick et al. 2021) and ground-based spectroscopy (orange, Hansen &
Hovenier 1974; dark cyan, Pollack et al. 1978); (d) the distributions of cloud
mass loading from our nominal simulation and D22. The observations are from
Pioneer Venus (orange, Knollenberg & Hunten 1980) and Venus Express (dark
cyan, Wilquet et al. 2009).

###  3.1 Nominal model results

We first compare our nominal model results with the observations and with the
simulations results in other models. For the model comparison, we mainly focus
on the KP94 model and the simple condensation model (D22). Our model adopts
the same condensation/evaporation physics in D22 and thus the two models are
very similar. The only difference, other than the numerical method, is that
this semi-analytical model pre-defined a cloud base property (i.e., altitude,
vapor mixing ratios, and the cloud acidity) using the method in Section 2.5,
while the cloud base is self-consistently simulated in the D22 condensation
model. In other words, our model imposed a discontinuity of the droplet flux
at the cloud base and neglected the cloud-base droplet diffusion. Because the
hydrogen flux should be conserved, the cloud mass loading and the
sedimentation flux at the cloud base in our model are a bit larger than those
in D22 (Fig. 3). The same argument also applies to the KP94 model. Moreover,
the KP94 model also assumes the LTE condition for the H2SO4 while our model
did not. We do not compare our results to microphysical models (e.g., Gao et
al., 2014; McGouldrick, 2017; Rossow, 1977) because the differences between
the condensation-only models and the full microphysical model have been
discussed in Dai et al. (2022) in detail.

The calculated H2SO4 vapor mixing ratio (Fig. 3) decreases from 5.5 ppm at the
cloud base to 1.7$\times$10-7 ppm at 60 km, then it increases to 1$\times$10-5
ppm at 66 km. The model results agree with the observations of Magellan
Mission (Kolodner & Steffes, 1998), Akatsuki Mission (Imamura et al., 2017),
Venus Express (Oschlisniok et al., 2012, 2021). Our model results agree pretty
well with a previous simulation in the condensation model D22 using time
stepping. Below 60 km, our model is also consistent with that in the KP94
model, implying that LTE is a good assumption below 60 km. Note that the KP94
used slightly different vapor mixing ratios at their lower boundary (47 km)
from our model.

The H2SO4 vapor does not follow LTE above 60 km. Our simulated H2SO4
supersaturation above 60 km is consistent with that in Dai et al. (2022). Fig.
4 illustrates the resultant H2SO4 condensation rate coefficient $C_{1}$, the
supersaturated H2SO4 mixing ratio $P_{1}/C_{1}$ ($=q^{g}_{1}-q^{svp}_{1}$),
and the SVMR of H2SO4 $q^{svp}_{1}$. The condensation rate coefficient
decreases exponentially with increasing altitude. This again confirms the
conclusion in Dai et al. (2022) that a slow condensation results in a large
vapor supersaturation in the upper cloud region. The supersaturated H2SO4
mixing ratio exceeds SVMR at 60 km, where the H2SO4 vapor pressure is
consequently about twice the saturation vapor pressure. Above 60 km, the large
chemical production of H2SO4 leads to an inversion of H2SO4 vapor profile.

Krasnopolsky (2011) predicted that the supersaturation of H2SO4 could occur
above 100 km by comparing the timescales of condensation and eddy mixing. But
they did not consider the H2SO4 chemical production which should be much
shorter. Here we give a simple estimate. The H2SO4 chemical net production
rate at 70 km is $P_{1}$(70 km)=1$\times$105 cm-3 s-1 (Krasnopolsky, 2012).
The H2SO4 vapor abundance at 70 km is $n_{1}$(70 km)=1$\times$107 cm-3,
adopted from Krasnopolsky (2011). Dai et al. (2022) proposed that chemical
production of H2SO4 vapor is much stronger than its diffusional divergence at
70 km. Then its chemical timescale is $n_{1}$(70 km)/$P_{1}$(70 km)=100 s at
70 km, which is much smaller than the condensation timescale of 800 s in
Krasnopolsky (2011). It indicates that the condensation is not efficient
enough to consume the quickly produced H2SO4 vapor, which leads to the
supersaturation. In the upper atmosphere, the lack of CCN could also lead to a
low condensation rate and subsequently induce the supersaturation
(McGouldrick, 2017).

Figure 4: The formation of the H2SO4 supersaturation, including the
distributions of the H2SO4 condensation rate coefficient $C_{1}$ (red), the
supersaturated H2SO4 mixing ratio $P_{1}/C_{1}$ (blue solid line), and the
SVMR of H2SO4 $q^{svp}_{1}$ (blue dotted line).

The derived H2O vapor mixing ratio (Fig. 3) decreases from 25 ppm at the cloud
base to 2.5 ppm at 70 km, agreeing with the observation of ground-based
spectroscopy (de Bergh et al., 1991; de Bergh et al., 1995; Gurwell et al.,
2007; Meadows & Crisp, 1996; Sandor & Clancy, 2005). It is about 1 ppm larger
than observations of Venus Express (Bertaux et al., 2007; Fedorova et al.,
2008) and ground-based spectroscopy (Encrenaz et al., 2013, 2020) above 70 km,
but the model results could be sensitive to the lower boundary condition (Dai
et al., 2022). Our H2O vapor mixing ratio is generally consistent with that in
D22 model within 30% ( 3 ppm versus 4 ppm). As mentioned earlier, our model
neglects part of the upward diffusion of the re-evaporated vapor below the
cloud base. The differences with Krasnopolsky (2015) are caused by different
thermodynamical properties in the SVMR calculations and different lower
boundary conditions. They adopted thermodynamical data from Giauque et al.
(1960), where the temperature variations on chemical potential were not
included in their study.

The vertical profile of the cloud mass loading (Fig. 3) agrees with
observations of LCPS (Knollenberg & Hunten, 1980) and estimates based on Venus
Express (Wilquet et al., 2009) except that it is underestimated by a factor of
two in 52-55 km and overestimated by a factor of two in 70-74 km. The
simulated mass loading is generally consistent with that in D22 model, except
at the cloud base. Our calculated cloud base altitude is 48 km, with a sharp
discontinuity in cloud mass, while that in the D22 model is smooth when the
cloud-base particle diffusion is included. Because of this different treatment
of the cloud base, the condensed H2SO4 is larger than their results in the
lower clouds. But note that both our model and the D22 model underestimated
the cloud mass by a factor of ten near the cloud base. This is because we did
not include the large Mode-3 particles in our study. Based on the observations
from Pioneer Venus, Toon et al. (1984) suspected that there could be large,
solid particles in the lower clouds that contribute to the mass loading.
Rimmer et al. (2021) proposed that the Mode-3 particles are made of hydroxide
salts such as NaOH and Ca(OH)2. Our study further confirms that H2SO4-H2O
system appears not sufficient to explain mass loading in the lower clouds.

The cloud acidity (Fig. 3) decreases from 98% at the cloud base to 78% at 62
km, then increases to 80% at 67 km. It is consistent with the estimate of
Barstow et al. (2012); McGouldrick et al. (2021); Pollack et al. (1978) but is
lower than Cottini et al. (2012) in the upper clouds and higher than Hansen &
Hovenier (1974) in the cloud deck. As mentioned above, our cloud acidity is
larger above 60 km because of less water vapor flux from below. The high
cloud-base acidity agrees with both the KP94 and D22 models. Our acidity in
60-66 km is smaller than KP94 model due to the different thermodynamical
treatments.

The cloud acidity tightly interacts with vapor vertical fluxes and the H2SO4
net flux (see equation 8). Fig. 5 shows the distributions of vapor diffusion
fluxes ($F^{g}_{1}$ for H2SO4 and $F^{g}_{2}$ for H2O), the H2SO4 net flux
($\Phi_{1}$), as well as the H2O/H2SO4 molar ratio in the droplets ($m$). The
H2O diffusion flux is relatively constant between 48-66 km, and potentially
decreases with altitude above 66 km. Its value remains on the same order of
magnitude as the H2SO4 net flux, indicating that eddy diffusion of H2O vapor
influences the cloud acidity over the cloud deck and the upper hazes. The
H2SO4 vapor diffusion flux dominates in the lower clouds but rapidly decreases
with altitude in the middle and upper clouds. Its value is over three orders
of magnitude smaller than H2SO4 net flux above 60 km. This indicates that
H2SO4 vapor diffusion largely regulates the cloud acidity in the lower clouds
but has much less effect above 60 km.

Figure 5: (a) the H2SO4 net flux $-\Phi_{1}$ (black, minus for log scale), the
H2SO4 vapor diffusion flux (red for $F^{g}_{1}$ and blue for $-F^{g}_{1}$),
and the H2O vapor diffusion flux ($F^{g}_{2}$, yellow); (b) the derived
droplet H2O/H2SO4 molar ratio $m$. Figure 6: The mixing ratios of H2SO4 vapor
(the first column), H2O vapor (the second column), condensed H2SO4 (the third
column), and condensed H2O (the fourth column) as functions of vapor mixing
ratios at the lower boundary. We use 70 km (the first row), 66 km (the second
row), and 50 km (the third row) to represent the cloud top, the H2SO4
production peak altitude, and the lower clouds, respectively. Figure 7: The
cloud acidity (the first column), cloud mass loading (the second column),
particle diameter (the third column), and H2SO4 supersaturated mixing ratio
(the fourth column) as functions of vapor mixing ratios at the lower boundary.
The rows are the same as those in Fig. 6.

###  3.2 Exploration of a large parameter space

After the validation in Section 3.1, we conclude that the semi-analytical
model agrees with the D22 model pretty well. Also, the new model is much
faster and allows us to explore the system behavior in a large parameter
space. Here we perform the survey for different H2SO4 and H2O vapor mixing
ratios at the lower boundary. We vary H2SO4 from 3 to 20 ppm, H2O from 27 to
37 ppm while keeping the same chemical production rate of H2SO4. The simulated
behaviors of vapors and condensates are shown in Fig. 6, the cloud acidity,
cloud mass loading, particle diameter, and H2SO4 supersaturation are shown in
Fig. 7. We use 70 km, 66 km, and 50 km to represent the cloud top, the H2SO4
production peak altitude, and the lower clouds, respectively. The difference
between the cloud top and the H2SO4 production peak altitude is small. But the
responses to the boundaries are quite different between the upper clouds and
the lower clouds.

In the upper cloud region, both H2O vapor and condensed H2O mixing ratios
increase with the increase of the H2O vapor at the lower boundary as more
water vapor is diffused to the upper atmosphere. However, while the H2SO4
vapor mixing ratio increases, the condensed H2SO4 decreases as the H2O vapor
increases at the boundary. Although the variation is small, this behavior
suggests an interesting response to the system. As the condensed water
increases, the cloud particle size increases but the cloud acidity decreases.
The former facilitates the condensation rate of H2SO4, but the latter
decreases it (see equation 9). That the condensed H2SO4 decreases implies that
the acidity effect is more important. A smaller H2SO4 condensation rate leads
to more H2SO4 vapor and less condensed H2SO4.

Opposite to H2O, the upper-cloud H2SO4 vapor mixing ratio decreases with the
increasing H2SO4 mixing ratio at the lower boundary, although the change is
very small. The H2SO4 in the upper clouds is regulated by condensation and
chemistry rather than diffusion. The increased H2SO4 at the lower boundary
increases the lower-cloud condensed H2SO4 and cloud acidity, which increases
the absorption of the H2O vapor and weakens its upward transport. As mentioned
above, a smaller H2O flux would enhance the H2SO4 condensation in the upper
clouds. As a result, the mixing ratio of H2SO4 vapor decreases, and the
condensed H2SO4 increases.

In the lower cloud region, the H2SO4 vapor mixing ratio, cloud mass loading,
and particle diameter show an opposite behavior to the upper clouds, while
those of the H2O vapor and condensate are similar. In the lower clouds, the
H2SO4 vapor mixing ratio is controlled by its SVMR rather than
supersaturation. As the H2O vapor mixing ratio at the lower boundary
increases, the lower-cloud acidity decreases, decreasing the H2SO4 SVMR. An
increase of H2SO4 vapor mixing ratio at the lower boundary enhances the H2SO4
condensation and the mass loading in the lower clouds. We note that these
variations are significant when the boundary H2SO4 vapor mixing ratio is
small, as the result of the cloud base altitude increases with decreasing
boundary H2SO4 vapor mixing ratio (see Fig. 2).

Next, we explore the system response to varying eddy diffusivities. The
diffusivity in the middle and upper atmosphere is constrained by observations
(Lane & Opstbaum, 1983; Mahieux et al., 2021; Woo & Ishimaru, 1981), but it is
largely unknown in the lower atmosphere. Krasnopolsky & Pollack (1994)
inferred smaller diffusivity below the clouds. Here we test different
diffusivities in the lower atmosphere. We increase (or decrease) the standard
diffusivity profile below 54 km by varying the diffusivity at 40 km to 2000,
4000, 6000, and 8000 cm2 s-1, respectively (see Fig. 1a). With these
variations, the standard diffusivity at 48 km (near the cloud base) is
increased by -13.4%, 16.5%, 38.6%, and 56.8%, respectively. We name these
cases with their diffusivity at 40 km ($K_{zz40}$). The results are shown in
Fig. 8.

Figure 8: The simulations with different eddy diffusivities in the lower
atmosphere. We increase (or decrease) the standard diffusivity profile below
54 km by varying the diffusivity at 40 km to 2000, 4000, 6000, and 8000 cm2
s-1, respectively (see Fig. 1a). With these variations, the standard
diffusivity at 48 km (near the cloud base) is increased by -13.4%, 16.5%,
38.6%, and 56.8%, respectively. We name these cases with their diffusivities
at 40 km ($K_{zz40}$). The line styles, observations and estimates are the
same as those in Fig. 3

With the increase of the lower atmospheric diffusivity, the cloud base is
lifted. It increases the lower-cloud mass loading by strengthening the re-
condensation of the H2SO4 vapor that evaporated from cloud-base droplets (Dai
et al., 2022). But it has a minor effect on the mass loading of the middle and
upper clouds because the upper cloud mass mainly originates from the chemical
production of H2SO4 rather than the transport from below. Moreover, with the
strengthened vertical transport of H2O vapor from the lower atmosphere, its
mixing ratio increases throughout the cloud deck and the upper hazes. This is
consistent with Krasnopolsky (2015). Besides, the H2O condensation is
strengthened with increasing H2O vapor mixing ratio. It also decreases the
cloud acidity, especially in the upper clouds where the condensation is
efficient.

Lastly, we explore the effect of different CCN profiles. The CCN profile is
prescribed in our model since nucleation is not included, which could
influence the condensation efficiency and impact the cloud formation. Our
nominal CCN profile differs from observations by less than four times (see
Fig. 1). But here we scale the CCN number density profile by up to three
orders of magnitude to see the system response. The results are shown in Fig.
9. The increase of the CCN increases the mixing ratios of both condensates by
increasing the condensation rate. Thus, the supersaturation of H2SO4
decreases, agreeing with Dai et al. (2022). Surprisingly, even if the CCN
increases by three orders of magnitude, it is not able to entirely remove the
supersaturation of H2SO4 in the upper clouds. Besides, few variations are seen
in the cloud acidity and H2O vapor because the condensation is also enhanced
as the CCN number density increases.

Figure 9: The simulations with different CCN number density profiles. The
standard CCN number density profile is scaled by -3 (blue), -2 (cyan), -1
(green), 1 (golden), 2 (light brown), and 3 (red) orders of magnitude,
respectively. The nominal results (black) are represented for comparison. The
line styles, observations and estimates are the same as those in Fig. 3.

##  4 Conclusions

This study provided a fast, semi-analytical gas-cloud model on Venus, based on
the conservation of the hydrogen flux and the interaction between cloud
acidity and vapor distributions for H2SO4 and H2O. The physics in this model
is based on the condensation cloud model in Dai et al. (2022). The calculation
of the vapors is built on the framework of Krasnopolsky & Pollack (1994) and
Krasnopolsky (2015) but our model only assumes LTE for H2O but not for H2SO4.
Moreover, we include the feedback of cloud condensation to the vapor
distributions, allowing us to explore the supersaturation of the H2SO4 vapor
above 60 km. Similar to the KP94 model, our model can be directly solved using
a semi-analytical method without time stepping.

Our simulation agrees with the observed vertical distributions of H2SO4 and
H2O vapors, the cloud mass loading, and the estimate of cloud acidity. The
results are consistent with Dai et al. (2022) including a highly
supersaturated H2SO4 layer above 60 km, different from that in the previous
LTE model in Krasnopolsky & Pollack (1994). We confirm the conclusion in Dai
et al. (2022) that a slow condensation results in such a large vapor
supersaturation in the upper cloud region. Based on the hydrogen flux
conservation, we found that the vertical transport of H2O vapor affects the
cloud acidity throughout the clouds but that of H2SO4 vapor is only important
in the lower clouds.

The fast model allows us to explore the behavior of the sulfuric acid gas-
cloud system in a large parameter space. We found that the increase of H2O
vapor mixing ratio at the lower boundary decreases the H2SO4 condensation rate
in the upper clouds, suggesting that the effect of acidity is more important
than particle size on the upper-cloud H2SO4 condensation rate. The increase of
H2SO4 mixing ratio at the lower boundary enhances its condensation in the
lower clouds, increasing the lower-cloud mass loading. But it increases the
absorption of H2O vapor in the lower clouds and decreases its upward
transport, leading to a similar behavior in the upper clouds to decreasing the
H2O at the lower boundary. Also, increasing the H2SO4 at the lower boundary
lifts the cloud base altitude. Increasing lower atmospheric diffusivity could
also lift the cloud base and increases the lower-cloud mass loading and the
upward transport of H2O vapor, but it has few effects on H2SO4 vapor.
Increasing CCN number density by three orders of magnitude significantly
decreases the H2SO4 supersaturation in the upper clouds but is not able to
entirely remove it, further confirming that the H2SO4 supersaturation is
robust.

On Venus, clouds provide a vital link between the upper atmosphere and the
lower atmosphere. Our study shows that the observations of the cloud deck
properties can be used to constrain the lower atmospheric conditions, and vice
versa. Clouds are also closely connected with the atmospheric chemistry in
both gas and liquids (e.g., Rimmer et al., 2021), which need to be
investigated in further detail. On extra-solar terrestrial planets, it was
proposed that bright sulfuric acid clouds can be used to distinguish planets
with and without the surface ocean (e.g., Loftus et al., 2019). Since both the
complicated microphysics and the time integration are avoided in our semi-
analytical model, the model is 100 times faster than the previous condensation
model (Dai et al., 2022) and 1000 times faster than full microphysics models
(e.g., Gao et al., 2014; McGouldrick, 2017). It is fast enough to be coupled
with climate models and chemistry models to understand the cloudy atmospheres
of Venus and Venus-like extra-solar planets.

## Acknowledgements

This work is supported by the Chinese Academy of Sciences B-type Strategic
Priority Program No. XDB41000000. L. X. Zhang acknowledges support from the
National Science Foundation through grant AST1740921 and J. Cui acknowledges
support from the National Natural Science Foundation of China through grant
42030201.

## Data Availability

The data in this study is publicly available online at
https://doi.org/10.6084/m9.figshare.19590562.v1 (Dai, 2022).

## References

  * Arney et al. (2014) Arney G., Meadows V., Crisp D., Schmidt S. J., Bailey J., Robinson T., 2014, Journal of Geophysical Research (Planets), 119, 1860
  * Barstow et al. (2012) Barstow J. K., et al., 2012, Icarus, 217, 542
  * Bertaux et al. (2007) Bertaux J.-L., et al., 2007, Nature, 450, 646
  * Bierson & Zhang (2020) Bierson C. J., Zhang X., 2020, Journal of Geophysical Research (Planets), 125, e06159
  * Cottini et al. (2012) Cottini V., Ignatiev N. I., Piccioni G., Drossart P., Grassi D., Markiewicz W. J., 2012, Icarus, 217, 561
  * Dai (2022) Dai L., 2022, 10.6084/m9.figshare.19590562.v1
  * Dai et al. (2022) Dai L., Zhang X., Shao W. D., Bierson C. J., Cui J., 2022, Journal of Geophysical Research (Planets), 127, e07060
  * de Bergh et al. (1991) de Bergh C., Bezard B., Owen T., Crisp D., Maillard J. P., Lutz B. L., 1991, Science, 251, 547
  * de Bergh et al. (1995) de Bergh C., Bezard B., Crisp D., Maillard J. P., Owen T., Pollack J., Grinspoon D., 1995, Advances in Space Research, 15, 79
  * Encrenaz et al. (2013) Encrenaz T., et al., 2013, Astronomy & Astrophysics, 559, A65
  * Encrenaz et al. (2020) Encrenaz T., et al., 2020, Astronomy & Astrophysics, 639, A69
  * Esposito et al. (1983) Esposito L., Knollenberg R., Marov M. Y., Toon O., Turco R., 1983, The clouds and hazes of Venus, Venus DM Hunten, L. Colin, TM Donahue, VI Moroz, 484–564
  * Fedorova et al. (2008) Fedorova A., et al., 2008, Journal of Geophysical Research (Planets), 113, E00B22
  * Fenghour et al. (1998) Fenghour A., Wakeham W. A., Vesovic V., 1998, Journal of Physical and Chemical Reference Data, 27, 31
  * Gao et al. (2014) Gao P., Zhang X., Crisp D., Bardeen C. G., Yung Y. L., 2014, Icarus, 231, 83
  * Giauque et al. (1960) Giauque W., Hornung E., Kunzler J., Rubin T., 1960, Journal of the American Chemical Society, 82, 62
  * Gurwell et al. (2007) Gurwell M. A., Melnick G. J., Tolls V., Bergin E. A., Patten B. M., 2007, Icarus, 188, 288
  * Hamill et al. (1977) Hamill P., Toon O. B., Kiang C. S., 1977, Journal of Atmospheric Sciences, 34, 1104
  * Hansen & Hovenier (1974) Hansen J. E., Hovenier J. W., 1974, Journal of Atmospheric Sciences, 31, 1137
  * Hanson & Eisele (2000) Hanson D. R., Eisele F., 2000, Journal of Physical Chemistry A, 104, 1715
  * Imamura & Hashimoto (2001) Imamura T., Hashimoto G. L., 2001, Journal of the Atmospheric Sciences, 58, 3597
  * Imamura et al. (2017) Imamura T., et al., 2017, Earth, Planets and Space, 69, 137
  * James et al. (1997) James E. P., Toon O. B., Schubert G., 1997, Icarus, 129, 147
  * Knollenberg & Hunten (1980) Knollenberg R. G., Hunten D. M., 1980, Journal of Geophysical Research, 85, 8039
  * Kolodner & Steffes (1998) Kolodner M. A., Steffes P. G., 1998, Icarus, 132, 151
  * Krasnopolsky (1985) Krasnopolsky V. A., 1985, Planetary and Space Science, 33, 109
  * Krasnopolsky (2011) Krasnopolsky V. A., 2011, Icarus, 215, 197
  * Krasnopolsky (2012) Krasnopolsky V. A., 2012, Icarus, 218, 230
  * Krasnopolsky (2015) Krasnopolsky V. A., 2015, Icarus, 252, 327
  * Krasnopolsky & Parshev (1981) Krasnopolsky V. A., Parshev V. A., 1981, Nature, 292, 610
  * Krasnopolsky & Pollack (1994) Krasnopolsky V. A., Pollack J. B., 1994, Icarus, 109, 58
  * Kulmala & Laaksonen (1990) Kulmala M., Laaksonen A., 1990, The Journal of Chemical Physics, 93, 696
  * Lane & Opstbaum (1983) Lane W. A., Opstbaum R., 1983, Icarus, 54, 48
  * Loftus et al. (2019) Loftus K., Wordsworth R. D., Morley C. V., 2019, The Astrophysical Journal, 887, 231
  * Mahieux et al. (2021) Mahieux A., et al., 2021, Icarus, 361, 114388
  * Marcq et al. (2008) Marcq E., Bézard B., Drossart P., Piccioni G., Reess J. M., Henry F., 2008, Journal of Geophysical Research (Planets), 113, E00B07
  * McGouldrick (2017) McGouldrick K., 2017, Earth, Planets and Space, 69, 161
  * McGouldrick & Toon (2007) McGouldrick K., Toon O. B., 2007, Icarus, 191, 1
  * McGouldrick et al. (2021) McGouldrick K., Peralta J., Barstow J. K., Tsang C. C. C., 2021, The Planetary Science Journal, 2, 153
  * Meadows & Crisp (1996) Meadows V. S., Crisp D., 1996, Journal of Geophysical Research, 101, 4595
  * Mehio et al. (2014) Mehio N., Dai S., Jiang D.-e., 2014, Journal of Physical Chemistry A, 118, 1150
  * Mills (1998) Mills F. P., 1998, PhD thesis, California Institute of Technology
  * Myhre et al. (1998) Myhre C. E., Nielsen C. J., Saastad O. W., 1998, Journal of Chemical & Engineering Data, 43, 617
  * Oschlisniok et al. (2012) Oschlisniok J., Häusler B., Pätzold M., Tyler G. L., Bird M. K., Tellmann S., Remus S., Andert T., 2012, Icarus, 221, 940
  * Oschlisniok et al. (2021) Oschlisniok J., Häusler B., Pätzold M., Tellmann S., Bird M. K., Peter K., Andert T. P., 2021, Icarus, 362, 114405
  * Parkinson et al. (2015) Parkinson C. D., et al., 2015, Planetary and Space Science, 113, 205
  * Pollack et al. (1978) Pollack J. B., Strecker D. W., Witteborn F. C., Erickson E. F., Baldwin B. J., 1978, Icarus, 34, 28
  * Pollack et al. (1993) Pollack J. B., et al., 1993, Icarus, 103, 1
  * Rimmer et al. (2021) Rimmer P. B., Jordan S., Constantinou T., Woitke P., Shorttle O., Hobbs R., Paschodimas A., 2021, The Planetary Science Journal, 2, 133
  * Rossow (1977) Rossow W. B., 1977, Journal of Atmospheric Sciences, 34, 417
  * Rossow & Gierasch (1977) Rossow W. B., Gierasch P. J., 1977, Journal of Atmospheric Sciences, 34, 405
  * Sandor & Clancy (2005) Sandor B. J., Clancy R. T., 2005, Icarus, 177, 129
  * Seiff et al. (1985) Seiff A., et al., 1985, Advances in Space Research, 5, 3
  * Seinfeld & Pandis (2016) Seinfeld J. H., Pandis S. N., 2016, Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. John Wiley & Sons
  * Shao et al. (2020) Shao W. D., Zhang X., Bierson C. J., Encrenaz T., 2020, Journal of Geophysical Research (Planets), 125, e06195
  * Tabazadeh et al. (1997) Tabazadeh A., Toon O. B., Clegg S. L., Hamill P., 1997, Geophysical Research Letters, 24, 1931
  * Toon et al. (1984) Toon O. B., Ragent B., Colburn D., Blamont J., Cot C., 1984, Icarus, 57, 143
  * Wilquet et al. (2009) Wilquet V., et al., 2009, Journal of Geophysical Research (Planets), 114, E00B42
  * Winick & Stewart (1980) Winick J. R., Stewart A. I. F., 1980, Journal of Geophysical Research, 85, 7849
  * Woo & Ishimaru (1981) Woo R., Ishimaru A., 1981, Nature, 289, 383
  * Yung & Demore (1982) Yung Y. L., Demore W. B., 1982, Icarus, 51, 199
  * Zeleznik (1991) Zeleznik F. J., 1991, Journal of Physical and Chemical Reference Data, 20, 1157
  * Zhang et al. (2012) Zhang X., Liang M. C., Mills F. P., Belyaev D. A., Yung Y. L., 2012, Icarus, 217, 714