Two saturated states of the vertical shear instability in protoplanetary disks with vertically varying cooling times

Yuya Fukuhara1
Satoshi Okuzumi1
Tomohiro Ono12

1Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan
2School of Natural Sciences, Institute for Advanced Study, Princeton, NJ 08544, USA

protoplanetary disks — hydrodynamics — turbulence — methods: numerical

Turbulence in protoplanetary disks plays an important role in dust evolution and planetesimal formation.
The vertical shear instability (VSI) is one of the candidate hydrodynamic mechanisms that can generate turbulence in the outer disk regions.
The VSI requires rapid gas cooling in addition to vertical shear.
A linear stability analysis suggests that the VSI may not operate around the midplane where gas cooling is inefficient.
In this study, we investigate the nonlinear outcome of the VSI in disks with a linearly VSI-stable midplane region.
We perform two-dimensional global hydrodynamical simulations of an axisymmetric disk with vertically varying cooling times.
The vertical cooling time profile determines the thicknesses of the linearly VSI-stable midplane layer and unstable layers above and below the midplane.
We find that the thickness of the midplane stable layer determines the vertical structure of VSI-driven turbulence in the nonlinear saturated state.
We identify two types of final saturated state: (1) T states characterized by vertical turbulent motion penetrating into the VSI-stable midplane layer and (2) pT states characterized by turbulent motion confined in the unstable layers.
The pT states are realized when the midplane VSI-stable layer is thicker than two gas scale heights.
We also find that the VSI-driven turbulence is largely suppressed at all heights when the VSI-unstable region lying above and below the midplane is thinner than two gas scale heights.
We present empirical formulas that predict the strength of VSI-driven turbulence as a function of the thicknesses of the unstable and stable layers.
These formulas will be useful for studying how VSI-driven turbulence and dust grains controlling the disk cooling efficiency evolve simultaneously.

§ INTRODUCTION

The initial stage of planet formation is the formation of kilometer-sized planetesimals from micron-sized dust grains in protoplanetary disks (see  and  for recent reviews).
Micrometer-sized grains initially grow into larger aggregates through coagulation.
Sufficiently large aggregates settle onto the disk midplane and form a dense dust layer <cit.>.
Large aggregates can also concentrate radially toward local maxima of the disk gas pressure (e.g., ).
These dust concentration processes may lead to planetesimal formation through the streaming instability <cit.> and gravitational instabilities <cit.>.
If the aggregates are sticky enough, they may form planetesimals directly by successive coagulation (e.g. ).

In these planetesimal formation processes, gas disk turbulence can play many important roles, including both positive and negative ones.
Large long-lived vortices created by turbulence may lead to dust concentration and subsequent planetesimal formation through gravitational collapse (e.g., ).
On the other hand, strong turbulence induces a large impact velocity between dust particles (e.g., ) that may inhibit planetesimal formation through coagulation (e.g., ).
Turbulence can also diffuse the midplane dust layer and prevent the onset of gravitational and streaming instabilities (e.g., ).
Therefore, the strength of turbulence would determine how planetesimals form.

Recent radio interferometric observations with the Atacama Large Millimeter–submillimeter Array (ALMA) and other interferometers have provided important constraints on turbulent intensity in the outer region of protoplanetary disks.
Molecular line emission observations point to weak or no turbulence in the upper layers <cit.> and regions within one scale height above and below the midplane <cit.> in the outer disk region (for a recent review, see ).
Strong nonthermal gas motion is detected in a disk around DM Tau <cit.>.
Observations of dust rings and gaps in massive and large disks (e.g., ) can also be used to indirectly constrain the turbulence intensity because the morphology of the rings and gaps can be affected by turbulent dust diffusion.
The well-defined morphology of the dust gaps in the disks around HL Tau <cit.> and Oph163131 <cit.> suggests a significantly low level of turbulence.
On the other hand, one of the two major dust rings in the HD 163296 disk shows a high degree of vertical dust diffusion <cit.>, potentially implying strong turbulent diffusion in that location.
Taken together, the observations so far suggest that the outer disk regions are mostly laminar but can become strongly turbulent under some circumstances.

To better understand when and where strong turbulence emerges, it is essential to study what mechanisms drive disk turbulence.
Regarding the outer disk region, it was previously believed that the magnetorotational instability (MRI; ) is the most viable mechanism for driving turbulence (e.g., ).
However, recent theoretical studies have shown that the viability of the MRI in the outer disk is limited by ambipolar diffusion <cit.>.
Instead, recent studies point out the importance of purely hydrodynamic instabilities (for reviews, see  and ), in particular the vertical shear instability <cit.>.
This is a hydrodyamical instability similar to the Goldreich–Schubert–Fricke instability in differentially rotating stars <cit.>.
This instability requires a vertical gradient in the gas orbital velocity and rapid cooling of disk gas <cit.>.
Turbulence driven by the VSI has a predominant vertical motion (e.g., ) that can prevent dust settling  <cit.>.
Because the VSI requires rapid gas cooling, it is most active in outer disk regions with low optical depths <cit.>.
In the outer regions, the VSI can take over the MRI <cit.>.

Despite its importance, it is yet to be fully understand in what conditions the VSI produces turbulence around the midplane, where dust evolution and planetesimal formation mainly take place.
According to local linear analysis, the VSI tends to be suppressed in a midplane region with a high optical depth and hence a low cooling rate <cit.>.
However, <cit.> recently showed that such an optically thick midplane region can become turbulent if the VSI operates above the region.
They performed hydrodynamical simulations using a vertically varying cooling rate profile and found that a vertical gas motion generated by the VSI crosses the midplane in the final saturated state.
Yet, it remains unclear how strongly their results depend on the assumed vertical profile of the gas cooling rate.
In general, the disk cooling profile is determined by the size and spatial distribution of dust grains that dominate the disk opacity <cit.>.
Therefore, the cooling profile can change significantly as the grains  grow and/or settle toward the midplane <cit.>.
Very recently, <cit.> also showed that depletion of sub-micrometer-sized dust grains due to coagulation increases the gas cooling time and makes the disk VSI-stable.
These aspects are particularly important in the context of planetesimal formation, where dust evolution is inevitable.

In this paper, we investigate systematically how the saturated state of VSI-driven turbulence depends on the vertical profile of the disk cooling rate.
We perform global two-dimensional (2D) hydrodynamical simulations of an axisymmetric protoplanetary disk with a parameterized vertical profile of the disk gas cooling timescale.
We show that, depending on the thickness of the VSI-stable midplane region, a vertical gas motion generated by the VSI above the midplane either does or does not penetrate into the midplane.
We also provide empirical formulas that predict the saturated level of VSI-driven turbulence from a given vertical profile of the gas cooling timescale.
Such empirical relations will be useful to study how the saturated state of VSI-driven turbulence varies with long-term evolution of the background disk.

This paper is organized as follows.
In section <ref>, we describe our numerical hydrodynamic simulation setup, our disk cooling (thermal relaxation) model, and a simulation analysis method.
We then present the main results in section <ref> and discuss the implications of our study in section <ref>.
Section <ref> presents a summary.

§ NUMERICAL METHOD

We perform global 2D hydrodynamical simulations of an axisymmetric protoplanetary disk in spherical polar coordinates $(r,~\theta,~\phi)$.
We also use the cylindrical radius $R = r\sin{\theta}$ and vertical coordinate $z = r\cos{\theta}$.
The solved continuity, motion, and energy equations of hydrodynamics are

\begin{equation}\label{eq:continuity}
\frac{\pd \rho_{\rm g}}{\pd t}+ \nabla \cdot \paren{\rho_{\rm g}\bm v} = 0,
\end{equation}

\begin{equation}\label{eq:EOM}
\frac{\pd \rho_{\rm g} {\bm v}}{\pd t} + \nabla\cdot\left(\rho_{\rm g}{\bm v}{\bm v} \right) = -\nabla P-\rho_{\rm g}\nabla \Phi,
\end{equation}

\begin{equation}\label{eq:energy_eq}
\frac{\pd E}{\pd t}+\nabla \cdot \left[{\bm v}\paren{E+P} \right] = -\rho_{\rm g}{\bm v}\cdot\nabla \Phi -\frac{U\paren{\rho_{\rm g},~T}-U\paren{\rho_{\rm g},~T_0}}{t_{\rm relax}},
\end{equation}

where $\rho_{\rm g}$ is the gas density, $\bm v$ is the gas velocity, $P$ is the gas pressure, $\Phi$ is the gravitational potential, $E$ is the total energy per unit volume, $U$ is the internal energy per unit volume, $T$ is the temperature, and $T_0$ is the temperature before the cooling.
The gas velocity has three components for the radial, meridional, and azimuthal in spherical polar coordinates ${\bm v} = (v_r,~v_\theta,~v_\phi)$.
The gravitational potential is given by $\Phi = -GM/R$, where $G$ is the gravitational constant and $M$ is the mass of the central star.
The total energy per unit volume $E$ and internal energy per unit volume $U$ are

\begin{equation}
E = U\paren{\rho_{\rm g},~T} + \frac{1}{2}\rho_{\rm g} \mbox{\bm $v$}^2,
\end{equation}

\begin{equation}
U\paren{\rho_{\rm g},~T} = \frac{P\paren{\rho_{\rm g},~T}}{\gamma-1},
\end{equation}

respectively, where $\gamma$ is the heat capacity ratio and is taken to be 1.4.
In equation (<ref>), $t_{\rm relax}$ is the thermal relaxation (cooling) timescale.
In this study, the thermal relaxation timescale is an important parameter that determines the strength of VSI-driven turbulence and is set with a spatial distribution.
We describe more details of $t_{\rm relax}$ model in section <ref>.

To solve the hydrodynamical equations (equations (<ref>)–(<ref>)), we use the Athena++ code <cit.> with the combination of the Harten–Lax–Van Leer (HLLC) approximate Riemann solver <cit.>, the reconstruction scheme of a second-order piecewise linear method <cit.>, and second-order Runge–Kutta time integrator.
The Courant–Friedrichs–Lewy (CFL) number is set to 0.3.

§.§ Simulation setup

Parameters choices of global 2D hydrodynamical simulations.

Parameter     Symbol     Value

Reference gas density     $\rho_0$     1.0

Reference radius     $R_0$     1.0

Radial power–law index for the gas density     $p$     $-2.25$

Reference sound speed     $c_0$     0.05

Radial power–law index for the temperature     $q$     $-0.5$

Reference gas scale height     $H_0$     0.05

We assume that the disk is initially in vertical hydrostatic equilibrium and give the initial gas density profile as

\begin{equation}\label{eq:initial_gas_density}
\rho_{\rm g} = \rho_0 \paren{\frac{R}{R_0}}^p \exp{\left(-\frac{z^2}{2H_{\rm g}^2} \right)},
\end{equation}

where $\rho_0$ is the reference gas density, $R_0$ is the reference radius, and $p$ is the radial power-law index for the gas density.
The sound speed is given by

\begin{equation}
c_{\rm s} = c_0\paren{\frac{R}{R_0}}^{q/2},
\end{equation}

where $c_0$ is the reference sound speed and $q$ is the radial power-law index of the temperature.
The disk scale height is

\begin{equation}
H_{\rm g} = \frac{c_{\rm s}}{\Omega_{\rm K}} = H_0\paren{\frac{R}{R_0}}^{(q+3)/2},
\end{equation}

where $\Omega_{\rm K} = \sqrt{GM/R^3}$ is the Keplerian frequency and $H_0 = c_0/\Omega_{\rm K}$ is the reference gas scale height.
Our parameter choices are summarized in table <ref>.

The initial velocities are set to $v_r = v_\theta = 0$ and $v_\phi = R\Omega\paren{R,~z}$, where $\Omega\paren{R,~z}$ is the gas angular velocity. For $H_{\rm g} \ll R$, which is the case for the disk model in our simulations, $\Omega\paren{R,~z}$ can be approximated as <cit.>

\begin{equation}\label{eq:gas_angular_velocity}
\Omega\paren{R,~z} = \Omega_{\rm K} \left[ 1+\frac{1}{2}\paren{\frac{H_{\rm g}}{R}}^2 \paren{p+q+\frac{q}{2}\frac{z^2}{H_{\rm g}^2}} \right].
\end{equation}

From equation (<ref>), the vertical shear of the gas rotation velocity $\pd(R\Omega)/\pd z$ is given by

\begin{equation}\label{eq:vertical_shear}
\frac{\pd (R\Omega)}{\pd z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K}.
\end{equation}

This physical quantity is the driving force behind the dynamics of the VSI and characterizes the strength of VSI-driven turbulence.
We also add small cellwise random velocities with an amplitude of $10^{-5} = 2 \times 10^{-4}c_0$ to the initial velocity field.

At all computational boundaries, we fix both density and pressure to the initial values.
For the velocity components normal to the boundaries, we apply the outflow boundary conditions preventing inflow at the inner and outer radial boundaries and reflecting boundary conditions at the upper and lower meridional boundaries.

Grid cells have logarithmically and linearly uniform spacings in the radial and meridional directions, respectively.
The radial and meridional domains cover $0.5 \leq r \leq 2.5$ and $\pi/2 - 5H_0/R_0 \leq \theta \leq \pi/2 + 5H_0/R_0$.
Because a resolution of 100 cells or more per scale height is necessary to resolve VSI-driven turbulence <cit.>, we adopt a resolution of 128 cells per scale height in both radial and vertical directions.
Therefore, our simulations use 4160$\times$1280 grid cells[The highest resolution simulation in <cit.> has 4096 cells in the radial domain extending from 0.5 to 5.0 and thereby has $\sim 100$ cells per gas scale height in the radial direction.

We adopt the code units $M = G = R_0 = 1$.
In this unit system, the orbital period at $R=R_0$ is $P_{\rm in} = 2\pi$.

§.§ Cooling model

We parameterize the thermal relaxation timescale $t_{\rm relax}$ in terms of the dimensionless cooling (thermal relaxation) time

\begin{equation}\label{eq:beta}
\beta = t_{\rm relax}\Omega_{\rm K}.
\end{equation}

We use the terms “thermal relaxation" and “cooling" interchangeably throughout this study.

The cooling time varies vertically depending on the optical depth from the surface <cit.>.
In a region that is optically thin to its own thermal emission, the cooling time is primarily determined by the collision time for gas molecules and grains, and therefore scales inversely with the density and cross section of dust grains.
In an optically thick region, the diffusion approximation applies, and the cooling time scales linearly with the opacity.
In both cases, the exact vertical dependence of the cooling time would depend on the size and vertical distribution of the grains.
To study the $\beta$ dependence of our simulations systematically, we parameterize $\beta$ as

\begin{equation}\label{eq:beta_model}
\beta (R,z) = \left[\beta_0\exp{\paren{\frac{z^2}{a^2H_{\rm g}^2}}} + \beta_1\exp{\paren{-\frac{z^2}{b^2H_{\rm g}^2}}} \right]\left(\frac{R}{R_0}\right)^{(q+1)/2},
\end{equation}

where $\beta_0$, $a$, $\beta_1$, and $b$ are dimensionless numbers that characterize the $\beta$ profile.
The first and second terms in the square bracket mimic the increase and decrease of $\beta$ with $|z|$ in the optically thin (high-$|z|$) and thick (low-$|z|$) regions, respectively, with $aH_{\rm g}$ and $bH_{\rm g}$ representing the vertical length scales over which $\beta$ varies.
The radial dependence of $\beta$ is chosen to be the same as one of $H_{\rm g}/R$ to investigate the dependence of the radial constant thicknesses of VSI-unstable and stable layers determined by the critical cooling time of the VSI (see also section <ref>).
In reality, $\beta$ can have a different radial dependence from this (e.g., ).
In the simulations, we take $\beta$ not to exceed $10$ to avoid excessively large values of $\beta$ at high-$|z|$.

In this study, we treat $a$ and $b$ as free parameters.
In reality, these values can vary with dust growth and settling.
The value of $a$ can increase or decrease depending on dust growth and settling.
This is because the optically thick area around the midplane widens for a certain amount of dust concentration due to the presence of regions that change from a high to low optical depth, and narrows beyond that point.
On the other hand, dust growth and settling cause dust depletion and increase cooling time at high altitudes <cit.>, which corresponds to a decrease in $b$.

Vertical $\beta$ profile for runs of $(a,~\beta_1) = (2.10,~0)$ and $(0.75,~0)$ (upper panels) and $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ and $(2.1,~0.5,~0.9)$ (lower panels) at $R=1.0$. The gray regions are marked as “unstable layer" and other white regions are marked as “stable layer", which are determined by applying the linear criterion (equation (<ref>)). The dotted lines represent $\beta = \beta_{\rm gc}$. We refer to the thicknesses of the unstable layer (all panels) and stable midplane layer (only lower panels) as $\Delta L_{\rm u}$ and $\Delta L_{\rm s}$, respectively.

We perform 46 simulations with different values of $a$, $\beta_1$, and $b$ with $\beta_0 = 2\times 10^{-3}$ as summarized in table <ref> of appendix <ref>.
We also perform one simulation with the locally isothermal equation of state.
We take $\beta$ to be constant in time.
Figure <ref> illustrates the vertical profile of $\beta$ for four runs.
The $\beta$ profiles for all our runs are shown in figure <ref> in appendix <ref>.

§.§ Linearly unstable and stable layers

The prescribed vertical cooling rate profile determines where the linear growth of the VSI occurs.
Following <cit.> and  <cit.> expected linear VSI growth in regions that fulfill the linear instability criterion, we use

\begin{equation}\label{eq:global_criterion}
\beta \lesssim \beta_{\rm gc},
\end{equation}

\begin{equation}\label{eq:global_cooling_time}
\beta_{\rm gc} = \frac{|q|}{\gamma-1}\frac{H_{\rm g}}{R}
\end{equation}

is the dimensionless critical cooling time <cit.>.
In this study, we refer to such regions the linearly unstable layers[<cit.> called these the VSI zones.].
<cit.> originally proposed equation (<ref>) as a criterion for the vertically global instability (see also appendix <ref>).
However, we find that this criterion well predicts where the VSI starts to grow in our simulations with vertically varying cooling times (see section <ref>).

All our simulations have linearly unstable layers at some heights (see  figure <ref> for illustrative examples).
Runs with $\beta_1 = 0$ have a cooling time monotonically decreasing toward the midplane, yielding a single unstable layer at $|z| < z_{\rm u}$, where $z_{\rm u}$ is the height of the unstable layer's upper bundary.
Runs with $\beta_1 > 0$ have two linearly unstable layers sandwiching a midplane region where equation (<ref>) breaks down.
We call this midplane region the linearly stable layer and denote the height of its boundary by $z_{\rm s}$.
When the linearly stable layer is absent, we set $z_{\rm s} = 0$.

We define the thicknesses of the linearly stable and unstable layers as $\Delta L_{\rm s} = 2z_{\rm s}$ and $\Delta L_{\rm u} = 2z_{\rm u} - \Delta L_{\rm s}$. When $z_{\rm s} > 0$, $\Delta L_{\rm u}$ accounts for the thicknesses of the two separated unstable layers lying at $z < 0$ and $z > 0$ (see figure <ref>).
In general, $\Delta L_{\rm u}$ decreases with decreasing $a$, and $\Delta L_{\rm s}$ increases with increasing $b$.
The values of $\Delta L_{\rm u}$ and $\Delta L_{\rm s}$ for all runs are summarized in table <ref> of appendix <ref>.
As we show in section <ref>, $\Delta L_{\rm u}$ and $\Delta L_{\rm s}$ are key quantities that dictate the saturated state of VSI-driven turbulence.

§.§ Turbulence diagnostics

We quantify the strength of VSI-driven turbulence using the time average of the squared vertical velocity $\langle v_z^2\rangle$, where $v_z = v_r\cos{\theta}-v_\theta\sin{\theta}$ is the vertical velocity.
The bracket $\langle\cdots\rangle$ denotes the time average.
In our simulations, the time averaging is performed after the system relaxes into a quasi-steady state.
The mean squared vertical velocity is related to the vertical diffusion coefficient for gas and small dust particles <cit.>.
We discuss this in more detail in section <ref>.

Turbulence also transports the disk's radial angular momentum.
The efficiency of angular momentum transport is measured by the Reynolds stress $\langle\rho_{\rm g}\delta v_r\delta v_\phi\rangle$, where $\delta v_r = v_r -\langle v_r \rangle$ and $\delta v_\phi = v_\phi -\langle v_\phi \rangle$ are the dispersions of the radial and azimuthal velocity, respectively.

In this study, we use the dimensionless Reynolds stress $\alpha_{r\phi}$ defined by

\begin{equation}
\alpha_{r\phi} = \frac{\langle\rho_{\rm g}\delta v_r\delta v_\phi\rangle}{\langle P \rangle}.
\end{equation}

We compute the spatial distribution of $\alpha_{r\phi}$ from our simulation results to estimate any possible radial angular momentum transport caused by the VSI-driven turbulent motions.
Furthermore, to see a global effect of turbulence, we calculate the vertical average of $v_z^2$ and $\alpha_{r\phi}$ from $z = -z_{\rm u}$ to $z = z_{\rm u}$, which is denoted by overbars.

The time that it takes for the system to reach a quasi-steady state differs from one run to another (see figure <ref>), which may depend on the unstable layer thickness.
Therefore, we stop a run at 400 orbits if the quasi-steady state has already been reached by 250 orbits; otherwise, we continue the run until 800 orbits.
The stopping times for all runs are summarized in table <ref> of appendix <ref>.
For simulations with shorter and longer runtimes, time averaging is performed over periods of 250–400 and 650–800 orbits, respectively.

§ RESULTS

In this section, we present our simulation results to study how the vertical profile of the cooling time affects VSI-driven turbulence in protoplanetary disks.
We begin by defining two saturated states of VSI-driven turbulence and then analyze how the saturated state depends on the thicknesses of the unstable and stable layers in section <ref>.
We construct empirical formulas of $\langle v_z^2\rangle|_{\rm mid}$, $\langle \overline{v_z^2}\rangle$, and $\overline{\alpha_{r\phi}}$ as a function of the thicknesses of the unstable and midplane stable layers, i.e., $\Delta L_{\rm s}$ and $\Delta L_{\rm u}$, in section <ref>.

§.§ Two saturated states of VSI-driven turbulence

Vertical velocity $v_z/c_{\rm s}$ for runs with $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ (left panel) and $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (right panel), which relaxes into the T state and pT state, respectively, as a function of $R$ and $z$ at the end of the simulations. The dashed lines indicate the unstable layer's upper boundaries at $z = \pm z_{\rm u}$ (uppermost and lowermost lines) and stable midplane layer's upper boundaries at $z = \pm z_{\rm s}$ (two middle lines).

We find that the width of the linearly VSI-stable layer at the midplane determines the vertical structure of VSI-driven turbulence in a steady state.
Specifically, we identify two possible saturated states of turbulence.
In the first class of saturated states, which we call the T (turbulent) states, the vertical gas motion generated in the linearly VSI-unstable layers penetrates into the linearly VSI-stable midplane.
This state was already seen in the simulations by <cit.>.
In the second class, which we call the pT (partially turbulent) states, the vertical gas motion is well confined in the VSI-unstable layers, leaving the VSI-stable midplane layer only weakly turbulent.

Figure <ref> illustrates the two saturated states.
Here, we present the two-dimensional maps of the vertical velocity $v_z$, normalized by the sound speed $c_{\rm s}$, at the end of the simulations for $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ and $(2.1,~0.5,~0.9)$.
The run with $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ relaxes into the T state with a vertically uniform gas motion through $|z| < z_{\rm u} = 3.9 H_{\rm g}$.
This state is not expected from the linear stability analysis, which predicts that the midplane region of $|z| < z_{\rm s} = 0.9H_{\rm g}$ is linearly VSI-stable.
The run with $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ relaxes into the pT state with little vertical gas motion inside the linearly VSI-stable layers ($z_{\rm s} = 1.3H_{\rm g}$ for this case).

Our simulations do not clearly reproduce a finer spatial profile of vertical velocity caused by the secondary parametric instability, which can be related to the nonlinear saturation process of the VSI.
This is because our simulation resolution is insufficient to resolve the parametric instability that may require $\sim 300$ cells per gas scale height in the radial direction <cit.>, while it is sufficient to resolve the dominant VSI modes <cit.>.

Vertical velocity $v_z/c_{\rm s}$ for runs with $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ (T state; left panel) and $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (pT state; right panel) as a function of time and $z/H_{\rm g}$ at $R=1.0$. The dashed lines indicate $z = \pm z_{\rm u}$ (top and bottom lines) and $z = \pm z_{\rm s}$ (two middle lines).

To see how the final saturated states are reached, we plot in figure <ref> the vertical profiles of $v_z$ at $R=1.0$ for the two cases displayed in figure <ref>.
In both cases, the vertical gas motion starts to develop near the unstable layer's upper boundary at $|z|= z_{\rm u}$, where the vertical shear is strong.
For $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ (T state), the vertical flows developed in the VSI-unstable layers overshoot the stable midplane layer and eventually form a unified flow (see also <cit.>).
For $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (pT state), one can see that the gas vertical motion is well confined in the unstable layers.
The vertical profile of $\langle v_z^2\rangle$ is steady over the $t = 400$–$800$ orbits, suggesting that our simulation captures the final saturated state.
Figures <ref> in appendix <ref> presents the time evolution of turbulence diagnostics for the two runs presented here.

Time and vertical mean squared vertical velocity, $\langle \overline{v_z^2}\rangle$ vs. time mean squared vertical velocity at the midplane, $\langle v_z^2 \rangle|_{\rm mid}$, at $R=1.0$ for all runs presented in this study. The symbols correspond to the T states (circles) and pT states (crosses). The dashed lines show $\langle v_z^2 \rangle|_{\rm mid} = C_1 \times\langle \overline{v_z^2}\rangle$ with $C_1 = 1$, $0.1$, and $0.01$.

We define the T and pT states more quantitatively by using the time-averaged $v_z^2$ at the midplane, $\langle v_z^2 \rangle|_{\rm mid}$, and $v_z^2$ averaged both in time and in the vertical direction, $\langle \overline{v_z^2}\rangle$.
The ratio between the two quantities reflects how strongly the vertical gas flow penetrates into the midplane region.
Simulations exhibiting strongly overshooting vertical gas flows (T states) yield $\langle v_z^2 \rangle|_{\rm mid} \approx \langle\overline{v_z^2}\rangle$, whereas those with a less turbulent midplane region than in the unstable layers (pT states) yield $ \langle v_z^2 \rangle|_{\rm mid} \ll \langle\overline{v_z^2}\rangle$.
Figure <ref> shows $\langle v_z^2 \rangle|_{\rm mid}$ versus $\langle\overline{v_z^2} \rangle$ of all runs presented in this study.
A majority of our simulations result in either $\langle v_z^2 \rangle|_{\rm mid} \approx \langle \overline{v_z^2}\rangle$ or $\langle v_z^2 \rangle|_{\rm mid} \approx 0.01$–$0.1\langle \overline{v_z^2}\rangle$.
In the following, we refer to the T and pT states as saturated states with $\langle v_z^2 \rangle|_{\rm mid} > 0.1\langle \overline{v_z^2}\rangle$ and $\langle v_z^2 \rangle|_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$, respectively.
The states of all our simulation runs are summarized in table <ref> of appendix <ref>.

Time mean squared vertical velocity at the midplane $\langle v_z^2 \rangle|_{\rm mid}$ (upper panel) and its vertical average $\langle \overline{v_z^2}\rangle$ (lower panel) at $R=1.0$ from all simulation runs, mapped in the $\Delta L_{\rm u}$–$\Delta L_{\rm s}$ plane. The dashed line indicates the simulations relaxing to pT states ($\langle v_z^2 \rangle|_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$).

From the examples shown in figures <ref> and <ref>, one can expect that the thickness of the linearly VSI-stable layer, $\Delta L_{\rm s}$, determines the final saturated state.
To test this hypothesis, we map in figure <ref> the values of the turbulence diagnostics $\langle v_z^2 \rangle|_{\rm mid}$ and $\langle \overline{v_z^2}\rangle$ from all our simulations against $\Delta L_{\rm u}$ and $\Delta L_{\rm s}$.
We find that pT states ($\langle v_z^2 \rangle|_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$) are realized when the linearly stable midplane layer is as wide as $\Delta L_{\rm s} \gtrsim 2H_{\rm g}$.

Another important finding from figure <ref> is that the thickness of the linearly unstable layer, $\Delta L_{\rm u}$, determines the vertically averaged saturation level $\langle \overline{v_z^2}\rangle$.
Turbulence is largely suppressed at all heights in the cases of $\Delta L_{\rm u} \lesssim 2H_{\rm g}$.
One can see that $\langle \overline{v_z^2}\rangle$ as well as $\langle v_z^2 \rangle|_{\rm mid}$ decreases sharply from $\sim 2\times 10^{-3}c_0^2$ to $\ll 10^{-3}c_0^2$ as $\Delta L_{\rm u}$ falls below $2H_{\rm g}$.

Time and vertical mean squared vertical velocity, $\langle \overline{v_z^2}\rangle$, vs. vertical mean Reynolds stress, $\overline{\alpha_{r\phi}}$, at $R=1.0$ for all runs. The symbols correspond to the T states (circle) and pT states (crosses). The dashed lines show $\overline{\alpha_{r\phi}} = C_2\times \langle \overline{v_z^2}\rangle/c_0^2$ with $C_2 = 1$, $0.1$, and $0.01$.

For both the T and pT states, the vertical averaged Reynolds stress $\overline{\alpha_{r\phi}}$ is tightly correlated with $\langle \overline{v_z^2}\rangle$.
Figure <ref> plots $\langle \overline{v_z^2}\rangle$ versus $\overline{\alpha_{r\phi}}$ at $R=1.0$ from all runs presented in this study.
This figure shows that most simulations produce $\overline{\alpha_{r\phi}} = 0.01$–$0.1\langle \overline{v_z^2}\rangle/c_0^2$.

Our simulation results may change if we perform a three-dimensional (3D) simulation.
However, the 3D effecs may be minor because <cit.> already showed that the average vertical velocity in a 3D simulation is comparable to that in 2D simulations.

§.§ Empirical formulas for turbulent quantities

The results presented in the previous subsection show that important turbulent quantities $\langle v_z^2 \rangle|_{\rm mid}$, $\langle \overline{v_z^2}\rangle$, and $\overline{\alpha_{r\phi}}$ are all tightly correlated with $\Delta L_{\rm s}$ and $\Delta L_{\rm u}$.
This suggests that one can predict these quantities for general cases as a function of $\Delta L_{\rm s}$ and $\Delta L_{\rm u}$ without having to perform further simulations.
Here, we construct such formulas based on our simulation results.

Because there are two types of saturated states, we consider a fitting function for $\langle v_z^2 \rangle|_{\rm mid}$ of the form

\begin{equation}\label{eq:fitting_vz2}
\left.\frac{\langle v_z^2 \rangle}{c_{\rm s}^2}\right|_{\rm mid} = f_{\rm T}(\Delta L_{\rm u},~\Delta L_{\rm s}) + f_{\rm pT}(\Delta L_{\rm u},~\Delta L_{\rm s}),
\end{equation}

where $f_{\rm T}$ and $f_{\rm pT}$ represent $\langle v_z^2 \rangle|_{\rm mid}$ for the T and pT states, respectively.
For the T states, figure <ref> shows that $\langle v_z^2 \rangle|_{\rm mid}$ is approximately constant except at $\Delta L_{\rm s} \sim 2H_{\rm g}$, which is the boundary between the T and pT states, and at $\Delta L_{\rm u} \sim 2H_{\rm g}$, where $\langle v_z^2 \rangle|_{\rm mid}$ sharply drops.
We reproduce these features with a simple function

\begin{eqnarray}\label{eq:fitting_vz2_fT}
f_{\rm T} = 2\times 10^{-3} \varsigma_{3.5}\left(x_1\right) \varsigma_{25}\left(x_2\right),
\end{eqnarray}

where $\varsigma_{c}(x)$ is the sigmoid function defined by

\begin{equation}
\varsigma_c(x) = \frac{1}{1+\exp{(-cx)}} = \frac{\tanh{(cx/2)}+1}{2}.
\end{equation}

The sigmoid function approaches unity and zero in the limits of $x \gg 0$ and $x \ll 0$, respectively.
Therefore, we adopt this function to reproduce the characteristics of the turbulence diagnostics that are almost constant within the region of $\Delta L_{\rm u} > 2H_{\rm g}$ and $\Delta L_{\rm s} < 2H_{\rm g}$ and change sharply and continuously at some boundaries.
The arguments of the first and second sigmoid functions in $f_{\rm T}$ express the boundaries that $\langle v_z^2 \rangle|_{\rm mid}$ decreases sharply around $\Delta L_{\rm u} \sim 2H_{\rm g}$ and $\Delta L_{\rm s} \sim 2H_{\rm g}$, respectively.
We set $x_1 = 3.5\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s}-4.8$ and $x_2 = 0.07\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s}+1.8$, where $\Delta \tilde{L}_{\rm s} = \Delta L_{\rm s}/H_{\rm g}$ and $\Delta \tilde{L}_{\rm u} = \Delta L_{\rm u}/H_{\rm g}$.
For $f_{\rm pT}$, we use

\begin{equation}\label{eq:fitting_vz2_fpT}
f_{\rm pT} = 2\times 10^{-5} \varsigma_{7}\left(x_3\right),
\end{equation}

where $x_3 = \ln(\max\{\Delta \tilde{L}_{\rm u}-2.5,~0\})-\Delta \tilde{L}_{\rm s}+1.8$ that expresses the curve boundary of the drop in the pT states toward large $\Delta L_{\rm s}$.
We use the natural logarithmic function in $x_3$ to represent this curve boundary.

Upper panel: comparison of $\langle v_z^2 \rangle|_{\rm mid}$ from simulations (points) and the empirical formula (equations (<ref>)–(<ref>); background) on the $\Delta L_{\rm u}$–$\Delta L_{\rm s}$ plane. The triangles show the simulation results of <cit.>. The dashed lines are contours of $\langle v_z^2 \rangle|_{\rm mid} = 10^{-3}$, $10^{-4}$, $10^{-5}$, and $10^{-6}$ from the formula. The dotted lines show $3.5\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s}-4.8 = 0$, $0.07\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s}+1.8 = 0$, and $\ln(\max\{\Delta \tilde{L}_{\rm u}-2.5,~0\})-\Delta \tilde{L}_{\rm s}+1.8 = 0$. Lower panel: same as the upper panel, but comparing $\langle \overline{v_z^2}\rangle$ from the simulations and from equations (<ref>)–(<ref>).

Comparison of $\langle v_z^2\rangle_{\rm mid}$ from simulations  and empirical formula (equations (<ref>)–(<ref>); dashed lines) as a function of $\Delta L_{\rm s}$ The circles, crosses, and triangles are from simulations with $\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s} = 10.0$, $5.6$, and $3.7$, respectively.

The upper panel of figure <ref> compares the values of $\langle v_z^2 \rangle|_{\rm mid}$ from our simulations with those from the empirical formula.
This formula reproduces the simulation results of $\langle v_z^2 \rangle|_{\rm mid}/c_0^2 \sim 10^{-3}$ (T state) and $\sim 10^{-5}$ (pT state).
The sigmoid functions in $f_{\rm T}$ and $f_{\rm pT}$ are also useful to represent the sharp decrease in $\Delta L_{\rm u} \lesssim 2H_{\rm g}$ and $\Delta L_{\rm s} \sim 2H_{\rm g}$.
Furthermore, this formula also represents the simulation results of <cit.>.
To quantify the errors between the simulation results and formula, we show in figure <ref> the values of $\langle v_z^2 \rangle|_{\rm mid}$ for some simulations and their corresponding values for the empirical formula as a function of $\Delta L_{\rm s}$.
This figure indicates that the formula is accurate to less than an order of magnitude in $\langle v_z^2 \rangle|_{\rm mid}$ at $\Delta L_{\rm s} \lesssim 2 H_{\rm g}$.
Furthermore, the formula replicates the sharp drop at $\Delta L_{\rm s} > 2H_{\rm g}$ and $\langle v_z^2 \rangle|_{\rm mid}/c_0^2 \sim 10^{-5}$ for the pT state's simulations of $\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s} = 10.0$.

To verify the accuracy of the formula, we calculate the root-mean-squared error of the formula with respect to $\log_{10}\langle v_z^2 \rangle|_{\rm mid}$.
For simulations with $\langle v_z^2 \rangle|_{\rm mid} > 10^{-5}$, the error is 0.48 dex, which means that the formula has an accuracy of less than an order of magnitude.
The error increases to 0.95 dex if we include all simulations.

Similarly, we propose an empirical formula for $\langle \overline{v_z^2}\rangle$ given by

\begin{equation}\label{eq:fitting_vz2_int}
\frac{\langle \overline{v_z^2}\rangle}{c_{\rm s}^2} = g_{\rm T}(\Delta L_{\rm u},~\Delta L_{\rm s}) + g_{\rm pT}(\Delta L_{\rm u},~\Delta L_{\rm s}),
\end{equation}

where $g_{\rm T}$ and $g_{\rm pT}$ represent $\langle \overline{v_z^2}\rangle$ in the T and pT states, respectively.
Because the difference between $\langle v_z^2 \rangle|_{\rm mid}$ and $\langle \overline{v_z^2}\rangle$ appears in only the pT state (see figure <ref>), we fit $\langle \overline{v_z^2}\rangle$ by varying only the coefficients of $f_{\rm T}$ and $f_{\rm pT}$.
Therefore, we determine $g_{\rm T}$ and $g_{\rm pT}$ as

\begin{equation}\label{eq:fitting_vz2_int_gT}
g_{\rm T}(\Delta L_{\rm u},~\Delta L_{\rm s}) = 0.5\cdot f_{\rm T}(\Delta L_{\rm u},~\Delta L_{\rm s}),
\end{equation}

\begin{equation}\label{eq:fitting_vz2_int_gpT}
g_{\rm pT}(\Delta L_{\rm u},~\Delta L_{\rm s}) = 50\cdot f_{\rm pT}(\Delta L_{\rm u},~\Delta L_{\rm s}).
\end{equation}

The lower panel of figure <ref> compares $\langle \overline{v_z^2}\rangle$ from the formula with those from the simulations.
The formula for $\langle \overline{v_z^2}\rangle$ reproduces the high turbulence level ($\langle \overline{v_z^2}\rangle \sim 10^{-3}c_0^2$) and a sharp decrease at $\Delta L_{\rm u} \lesssim 2H_{\rm g}$.

Ratio of $\overline{\alpha_{r\phi}}$ and $\langle \overline{v_z^2}\rangle/c_0^2$, $C_2$, as a function of $\Delta L_{\rm u}$. The points and crosses plot the simulation results for all runs presented in this study and for <cit.>, respectively. The dashed line shows the fitting function in equation (<ref>).

Comparison of $\overline{\alpha_{r\phi}}$ from simulation results (circles, crosses, and triangles) and empirical formula (equation (<ref>); dashed lines) as a function of $\Delta L_{\rm s}$. The circles, crosses, and triangles are from simulations with $\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s} = 10.0$, $5.6$, and $3.7$, respectively.

The tight correlation between $\overline{\alpha_{r\phi}}$ and $\langle \overline{v_z^2}\rangle$ (figure <ref>) motivates us to consider a formula for $\overline{\alpha_{r\phi}}$ of the form

\begin{eqnarray}\label{eq:fitting_alpha}
\overline{\alpha_{r\phi}} &=&  C_2\left(\Delta L_{\rm u}\right)\frac{\langle \overline{v_z^2}\rangle}{c_0^2} \nonumber \\
&=& C_2\left(\Delta L_{\rm u}\right)\left[g_{\rm T}(\Delta L_{\rm u},~\Delta L_{\rm s}) + g_{\rm pT}(\Delta L_{\rm u},~\Delta L_{\rm s}) \right],
\end{eqnarray}

where $C_2\left(\Delta L_{\rm u}\right)$ represents the ratio of $\overline{\alpha_{r\phi}}$ to $\langle \overline{v_z^2}\rangle/c_0^2$.
Based on the simulation results shown in figure <ref>, we propose

\begin{equation}\label{eq:C2_fit}
C_2\left(\Delta L_{\rm u}\right) = 0.01 + 0.09\varsigma_1 \left( \Delta \tilde{L}_{\rm u} - 3.5 \right).
\end{equation}

This function smoothly decreases from $0.01$ to $0.1$ as $\Delta L_{\rm u}$ decreases.
Figure <ref> shows that equation (<ref>) reproduces the ratio $\overline{\alpha_{r\phi}}/(\langle \overline{v_z^2}\rangle/c_0^2)$ from the our simulations as well as simulations of <cit.> to within one order of magnitude.
The root-mean-squared error of the formula for $\log_{\rm 10}[\overline{\alpha_{r\phi}}/(\langle \overline{v_z^2}\rangle/c_0^2)]$ is 0.24 dex.
If we limit the comparison to cases with $\overline{\alpha_{r\phi}} \gtrsim 10^{-5}$, the formula is accurate to within a factor of a few.
This can be seen in figure <ref>, where we plot the values of $\overline{\alpha_{r\phi}}$ from the simulations and the formula as a function of $\Delta L_{\rm s}$.

§ DISCUSSION

§.§ Estimating dust vertical diffusion coefficient

The predominantly vertical gas motion in VSI-driven turbulence causes strong vertical dust diffusion <cit.>.
Our simulations show that the gas velocity dispersion at the midplane varies with the thicknesses of the VSI-stable and unstable layers.
Qualifying how the dust diffusion coefficient varies will be useful for studying planetesimal formation (e.g., ) and testing theory with millimeter observations of dust rings and gaps (e.g., ).
A direct measurement of the dust vertical diffusion coefficient requires a calculation of dust grains' motion in hydrodynamical simulations, which is beyond the scope of this study.
Here, we indrectly estimate the dust vertical diffusion coefficient using the gas vertical velocity dispersion measured in our simulations.

Formally, the dust vertical diffusion coefficient is defined as

\begin{equation}
D_z = \frac{1}{2}\frac{d\langle z^2 \rangle_{\rm p}}{dt},
\end{equation}

where $\langle z^2 \rangle_{\rm p}$ is the ensemble average of $z^2$ for dust grains.
For grains whose stopping time is shorter than the orbital period and the correlation time $\tau_{\rm corr}$ of turbulence, $D_z$ reduces to the gas vertical diffusion coefficient, which can be estimated as <cit.>

\begin{equation}
D_z \sim \langle v_z^2 \rangle \tau_{\rm corr}.
\end{equation}

It is useful to normalize $D_z$ as

\begin{equation}\label{eq:alpha_z}
\alpha_z \equiv \frac{D_z}{c_{\rm s} H_{\rm g}} \sim   \frac{\langle v_z^2 \rangle}{c_{\rm s}^2}\cdot \tau_{\rm corr}\Omega_{\rm K},
\end{equation}

where we have used $H_{\rm g} = c_{\rm s}/\Omega_{\rm K}$.
<cit.> and <cit.> independently estimate $\tau_{\rm corr}$ and report  $\tau_{\rm corr}\Omega_{\rm K} \sim 0.2$ and $\sim 20$, respectively[<cit.> report a dimensionless vertical diffusion coefficient and a mean squared vertical velocity of $\alpha_z = 5.4 \times 10^{-3}$ and $\langle v_z^2 \rangle^{1/2} = 0.0166c_{\rm s}$, respectively. These values and the relation in equation (<ref>) lead to $\tau_{\rm corr}\Omega_{\rm K} \sim 20$.].

We are particularly interested in the value of $\alpha_z$ at the midplane, $\alpha_{z,\rm mid}$, because the thickness of the dust sedimentary layer at the midplane scales as $\alpha_{z,\rm mid}^{-1/2}$ <cit.>.
Our simulations show $\langle v_z^2 \rangle|_{\rm mid}/c_{\rm s}^2 \sim 2\times 10^{-3}$ (section <ref>) for fully developed VSI-driven turbulence in T states with $\Delta  L_{\rm u} \gtrsim 2H_{\rm g}$.
For this case, equation (<ref>) predicts a dimensionless vertical diffusion coefficient of $\alpha_{z, \rm mid} \sim 4 \times 10^{-4\dots-2}$, with the uncertainty originating from that of $\tau_{\rm corr}$.
This predicted value is higher than the dimensionless Reynolds stress in VSI-driven turbulence ($\alpha_{r\phi}\sim 2\times 10^{-4}$), reflecting the predominantly vertical motion of VSI-driven turbulence.
For the cases of $\Delta L_{\rm u} < 2H_{\rm g}$ and $\Delta L_{\rm s} > 2H_{\rm g}$, equation (<ref>) predicts much smaller diffusion coefficients of $\alpha_{z, \rm mid} \ll 10^{-4}$ and $\alpha_{z, \rm mid} \sim 4 \times 10^{-6\dots-4}$, respectively.
Implications of the suppressed VSI-driven turbulence for dust evolution and disk observations are discussed in the following section.

§.§ Implications for dust evolution and observations of protoplanetary disks

The suppression of VSI-driven turbulence at the midplane in the cases of $\Delta L_{\rm u} < 2H_{\rm g}$ or $\Delta L_{\rm s} > 2H_{\rm g}$ has important implications for dust settling and planetesimal formation in the outer regions of protoplanetary disks.
Weak turbulence yields low relative velocities of dust particles, which is preferred for dust growth through coagulation without collisional fragmentation and erosion (e.g., ).
Furthermore, the weak turbulent diffusion ($\alpha_z \lesssim 10^{-4}$) in suppressed turbulence would promote planetesimal formation through the streaming and gravitational instabilities (e.g., ).
Therefore, disk regions with $\Delta L_{\rm u} < 2H_{\rm g}$ or $\Delta L_{\rm s} > 2H_{\rm g}$ would be preferential sites for planetesimal formation.

The suppression of VSI-driven turbulence may also explain the low level of vertical dust diffusion inferred from ALMA observations of some protoplanetary disks.
The dust rings around HL Tau and Oph 163131 exhibit well-separated morphology in millimeter images, indicating that the large dust particles in the rings have settled onto the midplane.
Assuming that millimeter-sized particles dominate the millimeter emission, these observations point to small vertical diffusion coefficients of $\alpha_z \sim$ a few $10^{-4}$ for HL Tau <cit.> and of $\alpha \lesssim 10^{-5}$ for Oph 163131 <cit.>.
These estimated values of $\alpha_z$ are consistent with suppressed VSI-driven turbulence in the cases of $\Delta L_{\rm u} < 2H_{\rm g}$ or $\Delta L_{\rm s} > 2H_{\rm g}$ (see section <ref>).
Therefore, we hypothesize that VSI-driven turbulence is indeed suppressed in the outer regions of these disks.
Testing this hypothesis requires detailed modeling of these disks' cooling structure.

§.§ Need for self-consistent modeling of dust and VSI evolution

Because dust particles control disk cooling, whether the condition $\Delta L_{\rm u} < 2H_{\rm g}$ or $\Delta L_{\rm s} > 2H_{\rm g}$ for the suppression of VSI-driven turbulence is realized would depend on the size and amount of the particles.
Dust growth and settling lead to a VSI-unstable region that has smaller $\Delta L_{\rm u}$, i.e., that is more confined around the midplane <cit.>.
Depletion of small grains that dominate the gas disk cooling can also increase the cooling time <cit.> and thereby make the VSI-unstable region smaller.
These effects may result in suppression of VSI-driven turbulence.
On the other hand, an increase in the dust surface density leads to a wider optically thick region around the midplane.
This may make the region around the midplane cooling less efficient and consequently suppress VSI-driven turbulence at the midplane.
Assessing whether $\Delta L_{\rm u} < 2H_{\rm g}$ or $\Delta L_{\rm s} > 2H_{\rm g}$ can be realized under realistic conditions requires self-consistent modeling of dust growth, dust surface density evolution,  and disk cooling that takes into account gas–dust thermal coupling.

Moreover, because a change in the saturated state of VSI-driven turbulence would also affect dust evolution, the evolution of dust and VSI should be coupled.
For instance, suppression of VSI-driven turbulence by dust growth, if it really occurs, would reduce the collision velocity between the dust particles and thus promote their further growth.
A similar positive feedback can also be expected for suppression of VSI-driven turbulence by dust settling.
These positive feedback effects can be important for understanding planet formation and turbulence in outer disk regions.

To quantify these effects, the empirical formulas presented in section <ref> will be useful.
These formulas represent the correlation between the thicknesses of the VSI-unstable and stable layers and the VSI-driven turbulence intensity.
The size of dust particles as well as the spatial distribution controls gas cooling and determines the thicknesses of the VSI-unstable and stable layers.
Our formulas can be used to predict how the saturated level of VSI-driven turbulence would evolve with the long-term evolution of dust.
This will be done in our future work.

§.§ Limitation of our simulations

Our simulations are subject to two important limitations that should be addressed in future work.
First, the cooling rates adopted in our simulations are vertically varying but constant in time.
In reality, because VSI-driven turbulence alters the dust distribution (e.g., ), the cooling time distribution determined by dust can change with turbulence.
Therefore, turbulence and cooling time would evolve simultaneously until turbulence and dust profile are in equilibrium states.
If the dust evolution timescale is longer than the timescale for VSI-driven turbulence saturation, we can investigate the co-evolution of dust and turbulence using empirical formulas presented in section <ref>.
This is because it is reasonable to assume that the size and spatial distribution of dust grains will not change while the VSI develops turbulence.
On the other hand, if this assumption breaks down, the stability of this system should be studied in hydrodynamical simulations that include the thermal coupling between gas and dust.
We plan to address these open issues in future work.

Second, our simulations neglect the effects of dust and magnetic fields on gas disk dynamics.
Dust would increase the effective buoyancy frequency of the gas <cit.> and weaken the VSI <cit.>.
Magnetic fields threading the disk may also suppress the VSI either directly through magnetic tension or indirectly through MRI turbulence <cit.>.
The roles of magnetic fields in VSI suppression can be positive or negative depending on non-ideal magnetohydrodynamical effects (ambipolar diffusion, Ohmic resistivity, and Hall effect;  ).
We plan to quantify these effects using simulations including dust feedback and magnetic field in the future.

§ SUMMARY

In this study, we have investigated how the saturated state of VSI-driven turbulence depends on the vertical profile of the disk cooling rate.
We have performed global two-dimensional hydrodynamical simulations of an axisymmetric protoplanetary disk with vertically varying cooling times.
Our key findings are summarized as follows.

* The thickness of the linearly VSI-stable layer at the midplane determines the vertical structure of VSI-driven turbulence in a steady state (figures <ref> and <ref>). We have identified two final saturated states of turbulence. In the first state, the vertical gas motion generated in the linearly VSI-unstable layers penetrates the VSI-stable midplane layer (T states). In the second state, the vertical gas motion is well confined in the unstable layers (pT states), leaving the stable midplane layer only weakly turbulent. Using the time averaged squared vertical velocity at the midplane $\langle v_z^2 \rangle|_{\rm mid}$ and $v_z^2$ averaged both in time and in the vertical direction $\langle \overline{v_z^2}\rangle$, we refer to the T and pT states with $\langle v_z^2 \rangle|_{\rm mid} > 0.1\langle \overline{v_z^2}\rangle$ and $\langle v_z^2 \rangle|_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$, respectively (figure <ref>).

* The pT states are realized when the thickness of the VSI-stable midplane layer $\Delta L_{\rm s}$ is larger than two gas scale heights (figure <ref>). When the thickness of the VSI-unstable layer $\Delta L_{\rm u}$ is thinner than $2H_{\rm g}$, VSI-driven turbulence is also largely suppressed at all heights. The turbulence diagnostic value $\langle \overline{v_z^2}\rangle$ as well as $\langle v_z^2 \rangle|_{\rm mid}$ decreases sharply from $\sim 2\times 10^{-3}c_0^2$ to $\ll 10^{-3}c_0^2$, where $c_0$ is the sound speed, as $\Delta L_{\rm u}$ falls below $2H_{\rm g}$ (figure <ref>).

* For both the T and pT states, the vertical averaged Reynolds stress $\overline{\alpha_{r\phi}}$ and $\langle \overline{v_z^2}\rangle$ are connected by $\overline{\alpha_{r\phi}} = 0.01$–$0.1\langle \overline{v_z^2}\rangle/c_0^2$ (figure <ref>).

* We propose empirical formulas for the turbulence diagnostics $\langle v_z^2 \rangle|_{\rm mid}$, $\langle \overline{v_z^2}\rangle$, and $\overline{\alpha_{r\phi}}$ as a function of $\Delta L_{\rm u}$ and $\Delta L_{\rm s}$ (equations (<ref>), (<ref>), and (<ref>), respectively). These formulas reproduce the tightly correlation of turbulence diagnostics with $\Delta L_{\rm u}$ and $\Delta L_{\rm s}$ (figure <ref>) and have an accuracy of less than an order of magnitude in strong turbulence (figures <ref> and <ref>). These formulas will be useful for predicting how the states of VSI-driven turbulence vary with the long-term evolution of dust.

Our results suggest that the suppression of VSI-driven turbulence at the midplane in the cases of $\Delta L_{\rm u} < 2H_{\rm g}$ or $\Delta L_{\rm s} > 2H_{\rm g}$ can lead to the strongly vertical settling of dust particles.
This effect may promote planetesimal formation through dust coagulation and the gravitational and streaming instabilities in outer disk regions.
This effect may also explain the low level of vertical dust diffusion observed in the dust rings of some protoplanetary disks.

For simplicity, the present study has modeled the vertical cooling rate profile with a parameterized analytic function.
In reality, gas cooling is regulated by dust particles, and therefore the disks' cooling structure should depend on the particles' size and spatial distribution.
This implies that the evolution of dust and VSI can be coupled because VSI-driven turbulence can affect dust evolution.
For instance, dust growth would make the VSI-unstable region confined around the midplane and consequently suppress VSI-driven turbulence.
This turbulence suppression by dust growth may reduce the collision velocity between dust grains and thus promote their further growth.
To quantify this positive feedback effect, the empirical formulas presented in this study may be useful.
Verification of planetesimal formation by this positive feedback requires studying how the VSI and dust co-evolve in the future.

We thank Mario Flock for a discussion about the correlation time of VSI-driven turbulence.
We also thank the anonymous referee for helpful comments.
This work was supported by JSPS KAKENHI Grant Numbers JP18H05438, JP20H00182, JP20H01948, JP20J01376, and JP22J22593.
Numerical computations were carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

§ PARAMETER CHOICES AND FINAL SATURATED STATES FOR ALL RUNS

Parameter choices of our simulations with $\beta$ cooling models and final saturated states for all runs.
= 1.3mm

Run Name     $a$     $\beta_1$     $b$     $\Delta L_{\rm u}/H_{\rm g}$     $\Delta L_{\rm s}/H_{\rm g}$     Time [$P_{\rm in}$]     T/pT

isothermal     -     -     -     -     -     400     -

a270-b000     2.7     0.0     -     10.0     0.0     400     T

a270-b015     2.7     0.5     0.15     9.6     0.4     400     T

a270-b030     2.7     0.5     0.3     9.1     0.9     400     T

a270-b050     2.7     0.5     0.5     8.6     1.5     400     T

a270-b060     2.7     0.5     0.6     8.3     1.8     400     T

a270-b070     2.7     0.5     0.7     8.0     2.0     400     T

a270-b080     2.7     0.5     0.8     7.7     2.3     400     T

a270-b090     2.7     0.5     0.9     7.4     2.6     800     pT

a270-b100     2.7     0.5     1.0     7.1     2.9     800     pT

a270-b110     2.7     0.5     1.1     6.8     3.2     800     pT

a270-b120     2.7     0.5     1.2     6.5     3.5     800     pT

a270-b130     2.7     0.5     1.3     6.2     3.8     800     pT

a210-b000     2.1     0.0     -     7.8     0.0     400     T

a210-b015     2.1     0.5     0.15     7.3     0.4     400     T

a210-b030     2.1     0.5     0.3     6.9     0.9     400     T

a210-b050     2.1     0.5     0.5     6.3     1.5     400     T

a210-b060     2.1     0.5     0.6     6.0     1.8     400     T

a210-b070     2.1     0.5     0.7     5.8     2.1     400     T

a210-b080     2.1     0.5     0.8     5.5     2.3     800     T

a210-b090     2.1     0.5     0.9     5.2     2.6     800     pT

a210-b100     2.1     0.5     1.0     4.9     2.9     800     pT

a210-b110     2.1     0.5     1.1     4.6     3.2     800     pT

a150-b000     1.5     0.0     -     5.6     0.0     400     T

a150-b015     1.5     0.5     0.15     5.1     0.4     400     T

a150-b030     1.5     0.5     0.3     4.7     0.9     400     T

a150-b050     1.5     0.5     0.5     4.1     1.5     400     T

a150-b060     1.5     0.5     0.6     3.8     1.8     400     T

a150-b070     1.5     0.5     0.7     3.5     2.1     800     T

a150-b080     1.5     0.5     0.8     3.2     2.4     800     T

a150-b090     1.5     0.5     0.9     2.9     2.7     800     T

a100-b000     1.0     0.0     -     3.7     0.0     400     T

a100-b015     1.0     0.5     0.15     3.3     0.4     400     T

a100-b030     1.0     0.5     0.3     2.8     0.9     400     T

a100-b050     1.0     0.5     0.5     2.3     1.5     800     T

a100-b060     1.0     0.5     0.6     1.9     1.8     800     T

a100-b070     1.0     0.5     0.7     1.6     2.1     800     T

a100-b080     1.0     0.5     0.8     1.3     2.4     800     T

a075-b000     0.75     0.0     -     2.9     0.0     800     T

a075-b015     0.75     0.5     0.15     2.3     0.4     800     T

a075-b030     0.75     0.5     0.3     1.9     0.9     800     T

a075-b040     0.75     0.5     0.4     1.6     1.2     800     T

a075-b050     0.75     0.5     0.5     1.3     1.5     800     T

a050-b000     0.5     0.0     -     1.8     0.0     800     T

a050-b015     0.5     0.5     0.15     1.4     0.4     800     T

a050-b030     0.5     0.5     0.3     1.0     0.9     800     T

a025-b000     0.25     0.0     -     0.9     0.0     800     T

Vertical profiles of dimensionless cooling time $\beta$ as a function $z/H_{\rm g}$ at $R=1.0$ for all runs presented in this study.

We perform 46 hydrodynamical simulations with different values of $a$, $\beta_1$, and $b$ that determine the vertical profile of the cooling time $\beta$ (equation (<ref>)) and thicknesses of the unstable and midplane stable layers (equation (<ref>)).
The values of $a$, $\beta_1$, $b$, $\Delta L_{\rm u}$, and $\Delta L_{\rm s}$ for all runs presented in this study are summarized in table <ref>.
Figure <ref> illustrates the vertical $\beta$ profiles for all runs.

We stop a run at a different time for each run because the time, until the system reaches in quasi-steady state, differs from one run to another.
The runtimes and final saturated states of all our simulations are also summarized in table <ref>.

§ EXAMPLES OF TURBULENCE TIME EVOLUTION AND VERTICAL STRUCTURE FOR RUNS WITH T AND PT STATES

Time evolution of the squared vertical velocity at the midplane $v_z^2|_{\rm mid}$ and $v_z^2$ averaged in vertical direction $\overline{v_z^2}$ for runs with $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ (T state; upper panel) and $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (pT state; lower panel) at $R=1.0$.

In figure <ref>, we display the two-dimensional maps and vertical profile's time evolution of vertical velocity for runs with $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ and $(2.1,~0.5,~0.9)$, which reach the T state and pT state, respectively.
These two runs also differ in their time evolution of $v_z^2|_{\rm mid}$ and $\overline{v_z^2}$.
Figure <ref> shows time evolution of $v_z^2|_{\rm mid}$ and $\overline{v_z^2}$ for the two runs presented in section <ref>.
For $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ (T state), the vertical flows develop at $\lesssim 50 P_{\rm in}$, corresponding to the growth rate $\Gamma_{\rm VSI}$ being $\sim 10^{-2}\Omega_{\rm K}$.
On the other hand, for $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (pT state), the vertical gas motion grows at $\lesssim 300 P_{\rm in}$ with $\Gamma_{\rm VSI} \sim 10^{-3} \Omega_{\rm K}$.
The values of $v_z^2|_{\rm mid}$ and $\overline{v_z^2}$ are steady over $t = 250$–$400$ orbits and $t = 650$–$800$ orbits for runs with the T state and pT state, respectively, which are consistent with the vertical velocity's time evolution shown in figure <ref>.

Vertical profile of time mean squared vertical velocity $\langle v_z^2 \rangle$ (upper panel) and Reynolds stress $\alpha_{r\phi}$ (lower panel) for runs with the isothermal state, $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ (T state), and $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (pT state) at $R=1.0$. The dashed dotted lines represent the height of the unstable layer's upper boundaries $z = \pm z_{\rm u}$, which is the same for these two runs. The dotted lines represent the height of the midplane stable layer's uuper boundaries $z = \pm z_{\rm s}$ for each run.

These two runs reach different quasi-steady states, resulting in different turbulence profiles (see figures <ref> and <ref>).
To evaluate this quantitatively, we plot in figure <ref> vertical profiles of $\langle v_z^2 \rangle$ and $\alpha_{r\phi}$ for these two cases.
In both cases, $\langle v_z^2 \rangle$ and $\alpha_{r\phi}$ have a vertical profile with peaks near the unstable layer's upper boundary at $z = \pm z_{\rm u}$.
The run with $(a,~\beta_1,~b) = (2.1,~0.5,~0.6)$ (T state) has almost the same turbulent structure as the run of the isothermal state, which is the fundamental ideal state for the VSI.
In contrast, for $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (pT state), $\langle v_z^2 \rangle$ as well as $\alpha_{r\phi}$ decreases sharply near the stable midplane layer's upper boundary at $z = \pm z_{\rm s}$.
This vertical distribution of $\langle v_z^2 \rangle$ yields the difference of more than an order of magnitude between $\langle v_z^2 \rangle|_{\rm mid}$ and $\langle \overline{v_z^2} \rangle$ (see figure <ref>).
The increase in $\langle v_z^2 \rangle$ around the midplane may be related to the linear VSI local criterion in appendix <ref>.
Additionally, figure <ref> indicates that $\alpha_{r\phi}$ is an order-of-magnitude smaller than $\langle v_z^2 \rangle$ for both the T and pT states.
This is consistent with the correlation between $\overline{\alpha_{r\phi}}$ and $\langle \overline{v_z^2}\rangle$ shown in figure <ref>.

§ LOCAL CRITERION OF LINEAR VSI

So far we have used the global criterion of the linear VSI in equation (<ref>) to determine the thicknesses of the unstable and stable layers.
However, the slight rise of $\langle v_z^2 \rangle$ near the midplane of the run with $(a,~\beta_1,~b) = (2.1,~0.5,~0.9)$ (pT state) shown in the upper panel of figure <ref> may be due to an unstable midplane layer with a slightly thin determined by a local criterion of the VSI.

The onset of the VSI requires reducing buoyancy by rapid disk cooling.
When this effect is compared with the vertical shear that drives the VSI, the local criterion is given by <cit.>

\begin{equation}\label{eq:local_criterion}
\beta \leq \beta_{\rm lc} =  \frac{|\pd_z(R\Omega)|}{N_z^2}\Omega_{\rm K}.
\end{equation}

Here, $\beta_{\rm lc}$ is the local critical dimensionless cooling time and $N_z$ is the Brunt–Väisälä frequency defined by

\begin{equation}
N_z^2 \equiv -\frac{1}{\rho_{\rm g} C_P}\cdot\frac{\partial P}{\partial z}\cdot\frac{\partial s}{\partial z},
\end{equation}

where $C_P$ and $s$ are the specific heat at constant pressure and the specific entropy, respectively.
The specific entropy is given by $s = C_V\ln{\paren{P/\rho_{\rm g}^\gamma}}$, where $C_V$ is the specific heat at constant volume and $\gamma$ is the heat capacity ratio.
The local criterion can provide the global criterion in equation (<ref>) at $z = \gamma H_{\rm g}/2$ <cit.>.

The presence of the thin unstable midplane layer can be confirmed by looking at the vertical dependence of the local criterion.
Using the ideal gas law $P = k_{\rm B}\rho_{\rm g}T/m_{\rm g}$, where $k_{\rm B}$ is the Boltzmann constant and $m_{\rm g}$ is the mean molecular mass of the gas, and assuming vertically isothermal and vertical hydrostatic equilibrium (equation (<ref>)), we get $N_z^2 \propto (\pd_z \ln\rho_{\rm g})^2 \propto z^2$.
With this assumption and $\pd_z (R\Omega) \propto z$ (see equation (<ref>)), the local dimensionless critical cooling time can be analytically performed, resulting in $\beta_{\rm lc} \propto |z|^{-1}$.
Therefore, in regions where $z$ is sufficiently small (but not zero), $\beta_{\rm lc}$ increases, and thus the local criterion in equation (<ref>) is satisfied for any finite $\beta$.
As a result, the thin VSI-active layer determined by the local criterion always exists near the midplane[For constant $\beta = 1$ in vertical direction, a thickness of the midplane unstable layer is $\sim 0.1H_{\rm g}$.], although the VSI does not operate the midplane ($z = 0$) where the vertical shear vanishes.

This thin unstable layer near the midplane may be weakly turbulent by the VSI.
This may produce the slight vertical shading visible near the midplane in the right panel of figure <ref> and the slight increase of $\langle v_z^2 \rangle$ near the midplane in the upper panel of figure <ref>.
Therefore, to discuss exactly where the VSI operates turbulence, it may be necessary to apply the local criterion as well as the global criterion.


# Towards Autoformalization of Mathematics and Code Correctness: Experiments
with Elementary Proofs

Garett Cunningham
School of EECS
Ohio University
Athens, OH 45701
<EMAIL_ADDRESS>

&Razvan C. Bunescu
Department of Computer Science
UNC Charlotte
Charlotte, NC 28223
<EMAIL_ADDRESS>
&David Juedes
School of EECS
Ohio University
Athens, OH 45701
<EMAIL_ADDRESS>

###### Abstract

The ever-growing complexity of mathematical proofs makes their manual
verification by mathematicians very cognitively demanding. Autoformalization
seeks to address this by translating proofs written in natural language into a
formal representation that is computer-verifiable via interactive theorem
provers. In this paper, we introduce a semantic parsing approach, based on the
Universal Transformer architecture, that translates elementary mathematical
proofs into an equivalent formalization in the language of the Coq interactive
theorem prover. The same architecture is also trained to translate simple
imperative code decorated with Hoare triples into formally verifiable proofs
of correctness in Coq. Experiments on a limited domain of artificial and
human-written proofs show that the models generalize well to intermediate
lengths not seen during training and variations in natural language.

##  1 Introduction

To the uninitiated, the notion of mathematical proof represents simply an
argument written by people to convince others of mathematical truth. However,
in a real sense, mathematical proof must have formal underpinnings that go
beyond the written argument. Arguments that lack such underpinnings might have
fatal errors or even logical inconsistencies (see, for example, Russell’s
Paradox (Irvine and Deutsch, 2021)). Nevertheless, mathematical arguments
written in natural language are the norm and they have great value.

In Tymoczko (1979)’s well-known paper that discusses a somewhat controversial
(at the time) proof of the Four Color Theorem (Appel and Haken, 1977; Appel et
al., 1977), he explores “what is a mathematical proof?” He posits that all
mathematical proofs must be (i) convincing, (ii) surveyable, and (iii)
formalizable. The first two points are for the reader—proofs must be
convincing to and comprehensible by mathematicians. For the third point, he
notes that, “Most mathematicians and philosophers believe that any acceptable
proof can be formalized. We can always find an appropriate formal language and
theory in which the informal proof can be embedded and ‘filled out’ into a
rigorous formal proof.” For most mathematicians, this third part is crucial
for ensuring that subtle, but fatal, errors in logic do not exist in
mathematical proof.

Great progress has been made since the 1970’s in fully formalizing significant
mathematical results. For instance, the Feit-Thompson Theorem (Gonthier et
al., 2013; Gonthier, 2013) and the Four Color Theorem (Gonthier, 2008) have
been formally verified using the proof assistant Coq Bertot and Castéran
(2013), and the Kepler Conjecture (Hales, 2005; Hales et al., 2017) has been
formally verified using the proof assistants Isabelle and HOL Light (Nipkow et
al., 2002). Moreover, proof assistants have demonstrated immense utility for
software verification, such as the full certification of a C compiler (Leroy,
2009). Proofs demonstrating the correct behavior of code share a similar
structure to proofs in pure mathematics, where systems like Hoare logic
replace standard first-order logic. Thus, Tymoczko’s criteria for mathematical
proof can be extended to the verification of programs. For many experts, LaTeX
provides an excellent tool for satisfying the first two criteria. In addition,
carefully written LaTeX (Higham, 2020) provides a rich structure for
establishing the third criterion.

The vast majority of modern mathematics is expressed using natural language
(NL), with the overwhelming majority typeset in LaTeX. Fully formalizing
mathematics using proof assistants is still a difficult and time consuming
task. This paper takes some preliminary steps toward bridging this gap by
exploring how modern machine learning techniques can be used to convert
carefully written LaTeX into equivalent, and formally verified mathematics in
Coq, a process referred to as _autoformalization_ in the literature (Szegedy,
2020).

Wang et al. (2018, 2020) explored the similar task of translating mathematical
statements from LaTeX into Mizar, using LSTM-based models with attention. To
generate aligned LaTeX-Mizar pairs, they use a tool Bancerek (2006) that
translates top-level Mizar statements into artificial LaTeX sentences, a task
that is facilitated by the fact that Mizar is human readable and similar in
length with the corresponding LaTeX version. Carman (2021) evaluated the
competency of LSTMs toward formalizing a restricted set of artificially
generated theorems about simple arithmetic expressions, reporting reasonable
success over expression lengths seen during training. More recently, Wu et al.
(2022) evaluated Codex and PaLM on a significantly more limited, but human-
written set of theorems in algebra and number theory.

In contrast to prior work, we address the autoformalization of both theorems
and their proofs, and extend the scope to proofs of code correctness. We use a
number of manually written mathematical statements to abstract a complex
grammar that is then used to generate a dataset of substantially longer and
more diverse mathematical theorems and proofs. We develop an architecture
based on the Universal Transformer (Dehghani et al., 2018) and adapt a copying
mechanism (Gu et al., 2016) to handle arbitrary numbers and variable names at
test time. The models are evaluated extensively on their ability to
systematically generalize to statement lengths not seen during training, for
which we report sequence-level accuracy as well as a _semantic-level_ accuracy
calculated by combining sequence-level accuracy for the theorem and running
Coq to determine if the generated proof is correct. Code and data are made
available at https://github.com/gc974517/autoformalization.

##  2 Dataset of Theorems and Proofs

We create two independent datasets of mathematical statements that overall
correspond to four classes of theorems and proofs: the first dataset contains
three classes of arithmetic statements (even-odd, composites, and powers),
described in detail in Section 2.1, and the second dataset containing
statements about code correctness via Hoare logic (poly), described in detail
in Section 2.2. In each example, the input theorem-proof pair is given in
LaTeX, whereas the formalized output is represented in Coq. This work focuses
on the proof assistant Coq Bertot and Castéran (2013) because (a) there is a
rich set of mathematical libraries that have been developed for it, (b) it has
been used successfully to reason about significant computation artifacts, such
as the ComperCert C compiler Leroy (2009)), and (c) it benefits from a rich
set of training material for the proof assistant related to software
verification Pierce et al. (2010).

Each class of examples demonstrates features necessary for the successful
autoformalization of mathematical theorems and proofs. For example, powers and
composites examples may define useful terminology to make the theorems
shorter, e.g. proving that 4 is a square, or conversely they may state
theorems directly without any preliminary definitions, e.g. proving $\exists
n.\,n^{2}=4$. As shown in Figures 2.1 and 2.1, this corresponds in Coq to
aliasing propositions using the Definition keyword. Additionally, the examples
in the dataset provide a stress test of the copying mechanism described in
Section 3.1, testing its ability to learn the correct order and number of
terms to include in mathematical expressions, as well as their placement in
theorems and proofs, in a way that generalizes to arbitrary tokens in
mathematical language.

For each of the four classes of theorems and proofs, we manually created a few
examples ourselves in order to guide the construction of a complex grammar
that is then used to generate a dataset of substantially longer and more
diverse mathematical theorems and proofs. Each dataset is generated using its
corresponding grammar in an identical way. First, a random seed is sampled
that controls the overall structure of the theorem, proof, and definition, if
any. Then, the skeleton structure of the proof is completed with phrases that
are sampled from a separate context-free grammar. The coarse control of the
skeleton structure allows the construction of examples with interesting
features like sublemmas, forward or backward proof direction, coreference, or
additional conditions for the theorem, among others.

LaTeX Input Sequence
---

######  Theorem.

$28M+308$ is even.

###### Proof.

We know the summation between even numbers in $\mathbb{N}$ will be an even
number. Observe that $308$ is known to be even. Additionally, note that the
pair $M\times 28$ is trivially even. This is true because the coefficient $28$
is even. ∎
Coq Output Sequence
⬇ Require Import Arith. Theorem M28_308: forall M : nat, Nat.even(28 * M +
308) = true. Proof. intros. repeat rewrite Nat.even_add. assert (H1: Nat.even
308 = true). { auto. } assert (H2: Nat.even 28 = true). { auto. } assert (H3:
Nat.even (28 * M) = true). { rewrite Nat.even_mul. rewrite H2. auto. } rewrite
H1. rewrite H3. auto. Qed.

Figure 1: Generated example from the even-odd set.

LaTeX Input Sequence
---

######  Theorem.

$450+a\cdot 192+j\cdot 462$ is guaranteed to be even for any natural terms
$j$, and $a$.

###### Proof.

It can be justified that $192\cdot a+j\cdot 462$ is trivially even. Note that
$192a$ is an even number in $\mathbb{N}$ because multiplying between an even
integer with an arbitrary number in $\mathbb{N}$ is guaranteed to be even.
Likewise, $462j$ is trivially an even number in $\mathbb{N}$. The claim is
proven as a consequence of the fact that the sum of even numbers with an even
number will be in itself an even number. Therefore, our theorem holds. ∎
Coq Output Sequence
⬇ Require Import Arith. Theorem a450_192j_450_even: forall j a : nat, Nat.even
(192 * a + 462 * j + 450) = true. Proof. intros. rewrite Nat.even_add. assert
(H1: Nat.even (192 * a) = true). { rewrite Nat.even_mul. auto. } assert (H2:
Nat.even (462 * j) = true). { rewrite Nat.even_mul. auto. } assert (H3:
Nat.even (192 * a + 462 * j) = true). { repeat rewrite Nat.even_add. rewrite
H1. rewrite H2. auto. } rewrite H3. auto. Qed.

Figure 2: Instance of sublemma use in the even-odd dataset. The proof that the
sum of non-constant terms is even (assertion H3) is given before proving the
theorem.

Many of the difficulties in formalizing mathematical statements from NL into
Coq stem from the wide variability in the level of detail of mathematical
proofs, and the frequent mismatch between what is considered an acceptable
inference step in NL proofs vs. an inference step in Coq. Furthermore, there
may be multiple Coq proofs for any given theorem, at different levels of
granularity. We address this ambiguity by requiring the structure of the Coq
proof to match the overall structure of the NL proof. This is achieved by
quasi-synchronously generating the LaTeX and Coq versions of mathematical
statements, while still allowing for some simple re-orderings in order to
improve generalization performance, e.g. swapping arguments of commutative
operations.

In total, the grammar-based method for generating examples can theoretically
produce over 283 million unique arithmetic examples and over 491,000 unique
code examples, before considering variations in phrasing by sampling from the
context-free grammar.

###  2.1 Arithmetic Statements

We generated three classes of mathematical statements, i.e. theorem-proof
pairs:

  * •

even-odd: an expression is even or odd.

  * •

composites: a number is composite.

  * •

powers: a number is an integer power of $n$.

even-odd examples contain arithmetic expressions of $n$ variables with even
coefficients that are summed with a constant term, meaning that the parity of
this constant determines the parity of the whole expression. Proofs make use
of this fact with varying rigor based on our manually designed grammar, an
example of which is shown by Figure 2. The Coq program is generated
concurrently with the paired LaTeX example. The example shown in Figure 2
illustrates the generation and use of prior facts to prove an implicit
sublemma, in both the natural language and matching Coq version.

Examples of theorems and proofs for powers and composites share a similar
structure in both their LaTeX and Coq forms, as shown in Figures 2.1 and 2.1,
respectively. The theorems assert the existence of a natural number such that
a defining property holds and their proofs are constructive, with the
distinction that examples for composites prove factorization into $n$ factors.

LaTeX Input Sequence
---

######  Definition.

We define that $w\in\mathbb{N}$ is a _composite_ natural number if taking some
$R$, $Q\in\mathbb{N}$ we have $Q$, $R\geq 2$ and $Q\times R=w$.

######  Theorem.

$35$ is a _composite_ whole number.

###### Proof.

Remember that a composite natural number is the multiplication between $Q$ and
$R$ such that $Q$ and $R\geq 2$. Allow $R=7$, $Q=5$. We justify the result is
valid as $35=Q\cdot R$. ∎
Coq Output Sequence
⬇ Require Import Lia. Definition composite (w : nat) := exists R Q : nat, (Q
>= 2) /\ (R >= 2) /\ (Q * R = w). Theorem w_composite: composite 35. Proof.
unfold composite. exists 7. exists 5. lia. Qed.

Figure 3: Generated composites example.

LaTeX Input Sequence
---

######  Definition.

We say $o$ is a _square_ whenever there exists some whole number $Z$ such that
$Z\geq 2$ and $o=Z^{2}$.

######  Theorem.

$o=64$ is a _square_.

###### Proof.

Let $Z=8$. Observe that $64=8^{2}$. Also notice $Z=8$ is more than or equal to
2. This yields 64 is a square whole number. ∎
Coq Output Sequence
⬇ Require Import Lia. Definition square (o : nat) := exists Z : nat, (Z >= 2)
/\ (o = Z^2). Theorem square_64: square 64. Proof. unfold square. exists 8.
assert (H1: 8 >= 2). { lia. } repeat split. apply H1. Qed.

Figure 4: Generated example from the powers set.

For both training and testing, we generate 5,000 even-odd, 5,000 composites,
and 2,000 powers examples. We train on values of $n\in\\{2,3,5,7,9\\}$ and
test on values $n\in\\{2,3,\ldots,12\\}$, where $n$ represents the number of
variables in the arithmetic expression, the number of factors, or the power,
respectively. This is done in order to evaluate the model’s ability to
generalize to unseen arithmetic expression lengths and numbers of factors.

LaTeX Input Sequence |  | Coq Output Sequence
---|---|---

Figure 5: Generated poly example: [Left] the Hoare logic proof; [Right] the
code correctness proof in Coq.

####  2.1.1 Handwritten Examples

We also created a small collection of 45 human-written LaTeX theorem-proof
pairs to evaluate performance on examples outside of our manually generated
grammar. These are distinct from the original manually written examples that
were used to guide the development of the generative grammar. There are 15
examples for each type of proof from the arithmetic set, using the same
vocabulary with a number of unseen grammatical structures.

###  2.2 Code Correctness Statements

We create a dataset of correctness proofs about short programs written in the
imperative programming language _Imp_ (Pierce et al., 2018), which we call
poly. The programs represent various algorithms for evaluating a polynomial,
and their proofs of correctness verify that the programs correctly model the
polynomial as a mathematical function. Proofs are conducted as either fully
decorated programs or as sequences of Hoare triples with natural language
justifying steps in between. An example is shown in Figure 2.1.

For both training and testing data, we generate 5,000 examples. We train on
programs containing 2, 3, 5, 7, 9, and 11 lines, then test on programs
containing from 2 up to 14 lines to evaluate the model’s ability to generalize
to novel program lengths.

##  3 Semantic Parsing Architecture

To formalize LaTeX statements into Coq, we developed an encoder-decoder
architecture based on the Universal Transformer Dehghani et al. (2018).
Similar to Csordás et al. (2021), we do so by adding recursive passes into the
encoder and decoder of a base Transformer (Vaswani et al., 2017), thus making
the model analogous to a Universal Transformer without adaptive computation
time (ACT). Further, we introduce a copying mechanism and support for out-of-
vocabulary mathematical terms.

###  3.1 Copying Mechanism

Mathematical language contains features uncommon or non-existent in natural
language, such as numbers, variables, and carefully defined terminology. In
order to address the use of general mathematical jargon, these tokens are
replaced in the LaTeX input with generic forms denoting their usage, such as
<var1> up to <varN> for variables, which effectively ensures generalization to
variable renaming Ferreira et al. (2022), <nat1> up to <natN> for numbers, or
<def> for definitions, coupled with the use of a copying mechanism adapted
from Gu et al. (2016). Note that a different generic token is introduced for
each unique numerical constant or variable literal in the theorem and its
proof, and the corresponding generic token is used in the Coq version. For
example, considering the $\langle$LaTeX, Coq$\rangle$ pair in Figure 2.1,
<nat1>, <nat2>, <nat3>, and <nat4> would be used to replace the constants 2,
35, 7, and 5 respectively, everywhere in the LaTeX and Coq statements.
Similarly, <var1>, <var2>, and <var3> were used to replace variable literals
$w$, $R$, and $Q$. This is in contrast to using just two generic tokens <nat>
and <var> everywhere, which would make all numbers coreferent and all
variables coreferent. Preliminary experiments validated the utility of
encoding these distinctions while maintaining the correct coreference in both
LaTeX and Coq statements.

Overall, by using generic tokens for numbers, variables, and definitions, only
a limited set of embeddings need to be trained and the model is forced to
utilize contextual information in order to appropriately copy tokens into the
Coq output. In this way, the model has the ability to generalize to unseen
numbers or variable and definition names.

The original CopyNet (Gu et al., 2016) used an encoder-decoder architecture
with a copying mechanism to calculate the probabilities of generating in-
vocabulary tokens vs. copying tokens from the input sequence to the output.
Our autoformalization task guarantees mutual exclusivity between generating
(g) and copying (c) tokens, which allows using a simplified formula for
calculating the probability of producing a token $y_{t}$ at time step $t$.
Letting $\mathcal{V}_{c}$ denote the Coq vocabulary, $X$ denote the input
sequence of LaTeX tokens, and $\mathcal{X}$ denote the collection of unique
tokens in $X$, we calculate the probability of producing $y_{t}$ as:

$\displaystyle p(y_{t})$
$\displaystyle=\begin{cases}p(y_{t},\text{g})=\displaystyle\frac{1}{Z_{t}}e^{\psi_{g}(y_{t})},\hfill
y_{t}\in\mathcal{V}_{c}\\\
p(y_{t},\text{c})=\displaystyle\frac{1}{Z_{t}}\displaystyle\sum_{x_{j}\in
X:x_{j}=y_{t}}\\!\\!\\!e^{\psi_{c}(x_{j})},\hfill
y_{t}\in\mathcal{X}\end{cases}$

where
$Z_{t}=\displaystyle\sum_{y_{t}\in\mathcal{V}_{c}}e^{\psi_{g}(y_{t})}+\sum_{x_{j}\in
X}e^{\psi_{c}(x_{j})}$. The scoring functions are given by
$\psi_{g}(y_{t})=\mathbf{v}_{y_{t}}^{\top}\mathbf{W}_{o}\mathbf{s}_{t}$ and
$\psi_{c}(x_{j})=\tanh{\left(\mathbf{h}_{j}^{\top}\mathbf{W}_{c}\right)}\mathbf{s}_{t}$,
where $\mathbf{v}_{y_{t}}$ is a one-hot encoding of $y_{t}$, $\mathbf{h}_{j}$
is the hidden encoder state for the input token $x_{j}$, $\mathbf{s}_{t}$ is
the decoder state at step $t$, and $\mathbf{W}_{o}$ and $\mathbf{W}_{c}$ are
learnable parameters.

###  3.2 Encoder-Decoder Architecture

We diverge from the standard Transformer architecture in a few crucial ways:

  * •

Probabilities are calculated via $p(y_{t})$ above.

  * •

Absolute positional encodings are removed.

  * •

Self-attention uses relative positional representations as in Shaw et al.
(2018).

  * •

Stacks of $N$ encoder/decoder blocks have $T$ recurrent passes.

All other aspects of the model remain unchanged from the original Transformer.
We emphasize relative positional information over absolute in our model
architecture. Preliminary evaluations on the even-odd dataset showed that
Transformer models that use absolute positional encodings obtain 0% sequence-
level accuracy on expression lengths that are not seen at training time.
Removing reliance on absolute position resolves this type of systematic
generalization. The use of relative positional encodings for the Transformer-
based models was thus essential for achieving stronger systematic
generalization, which also agrees with the findings of Csordás et al. (2021)
on other NLP tasks.

##  4 Experimental Evaluations

To evaluate the performance of trained models, we ran two primary experiments:
first on the collection of arithmetic examples, then on the collection of code
correctness examples. All models are evaluated in terms of sequence-level
accuracy, where an example is considered correctly processed only if the
generated Coq sequence for both the theorem and its proof perfectly matches
token by token the ground truth sequence. We also report semantic-level
accuracy, for which the generated Coq theorem needs to attains a perfect
sequence-level accuracy and the Coq engine verifies that the generated Coq
proof truly proves the generated Coq theorm, regardless of whether it matches
the ground truth version of the proof. This emphasizes that the model was able
to capture the general meaning of the natural language proof by correctly
translating the theorem and successfully proving it using the natural language
version as a guide.

All experiments were performed on one NVIDIA RTX-A6000 GPU with 48GB of
memory.

###  4.1 Arithmetic Statements

We evaluate a Transformer model on the full data combining even-odd +
composites + powers and using both the theorem and its proof in each sequence.
We tune a model with embedding and state sizes of 32, a feed forward width of
256, 4 encoder and decoder blocks with 4 recurrent passes, 4 attention heads,
and a clipping value of 2 for self-attention. We trained this model over
minibatches of size 20, optimized with Adam using $\beta_{1}=0.9$,
$\beta_{2}=0.98$, $\varepsilon=1e-9$, and an initial learning rate of 0.001,
annealed by a factor of $1/\sqrt{10}$ based on training loss plateaus with a
patience of 5 epochs.

The results in Table 1 show that the model generalizes well to the
intermediate lengths of $\\{4,6,8\\}$, with a small number of correctly
translated examples longer than the maximum of 9 used in training. Otherwise,
the model fails to generalize to longer unseen lengths, which is not
surprising, given that Transformer models are known to fail dramatically at
systematic generalization on longer inputs for various NLP tasks Csordás et
al. (2021), or to incur substantial decrease in accuracy for longer symbolic
integration problems Welleck et al. (2022). Switching to semantic-level
evaluation leads to a significant increase in accuracy for composites, with a
more modest increase for even-odd.

| even-odd | composites | poly
---|---|---|---
$n$ | Seq | Sem | Seq | Sem | Both
2 | 99.6 | 99.8 | 76.7 | 97.6 | 100.0
3 | 99.4 | 99.6 | 64.6 | 94.2 | 100.0
4 | 99.4 | 99.4 | 56.1 | 93.9 | 82.1
5 | 99.2 | 99.6 | 54.9 | 94.4 | 99.2
6 | 98.8 | 98.8 | 57.1 | 94.3 | 45.1
7 | 99.1 | 99.5 | 58.5 | 93.4 | 96.5
8 | 93.8 | 94.0 | 53.5 | 88.3 | 15.7
9 | 98.6 | 98.6 | 53.7 | 93.7 | 98.2
10 | 7.0 | 7.0 | 1.2 | 1.6 | 35.6
11 | 0.0 | 0.0 | 0.0 | 0.0 | 93.5
12+ | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
powers | Seq = 100.0 | Sem = 100
Table 1: Sequence-level (Seq) and semantic-level (Sem) accuracy (%) on test
examples, split by expression length, with the exception of powers.

###  4.2 Code Correctness Statements

We extend our scope to include data representing proofs of program correctness
using the language of Hoare logic. We train a separate model with the same
embedding and state sizes, feed forward width, and learning rates as in
Section 4.1. Depth is increased to 8 encoder and decoder blocks with 8
recurrent passes, 8 attention heads, and a clipping value of 8. The model is
trained over minibatches of size 1 with Adam, with a patience of 3 epochs.

The poly results shown in Table 1 demonstrate that the model is able to
generalize to program line counts of $\\{4,6,8,10\\}$ unseen during training
with diminishing returns as the program length grows, eventually failing to
generalize for lengths longer than the maximum seen in training. We observe
that increasing the depth of the model significantly improved generalization.
A model with identical hyperparameters to the arithmetic experiment yielded
less then half the sequence-level accuracy for intermediate program lengths.
Therefore, further increasing the depth of the model could push performance
closer to optimal generalization to intermediate lengths at the cost of
significantly more computing resources. Additionally, poly examples are far
less prone to non-fatal token swapping errors. We observe that semantic-level
accuracy is identical to sequence-level, as all copying errors compromised the
validity of the proof. Therefore, accuracies are shown as one column (Both).

###  4.3 Handwritten Examples

We also evaluate the semantic-level accuracy of the trained models on the
collection of 45 human-written LaTeX theorem-proof pairs. This is done by
manually verifying that the generated Coq theorem corresponds to the LaTeX
version and that the subsequent proof is correct according to the Coq
interpreter. The fully trained model achieved 53.3% for both even-odd and
composites, and 73.3% for powers.

Mistakes in almost all cases are confined to the mishandling of out-of-
vocabulary tokens, such as mis-copying a variable within a definition or the
omission of an assertion in the proof tied to a term. The model otherwise
generated syntactically sound Coq code. Mistakes strongly correlate with
examples that deviate significantly from the grammatical structure of the
artificial data. Thus, pre-trained language models as evaluated by Wu et al.
(2022) or pre-training new models on mathematical corpora like MATH (Hendrycks
et al., 2021) may serve to alleviate the problems caused by the scarcity of
aligned natural and formal mathematics data.

##  5 Concluding Remarks

As we have seen, it is feasible to train machine learning models to perform
autoformalization over very restricted domains of math and code correctness
proofs. These models show capability to systematically generalize to new
expression lengths and program sizes. Moreover, these models were able to
translate previously unseen hand written natural language examples, albeit
with lower accuracy. We are hopeful that this approach can be applied to
autoformalization of a larger segment of mathematics and code verification.

As mentioned by Szegedy (2020), "Autoformalization is not just a challenge:
successful autoformalization would represent a breakthrough for general AI
with significant implications in various domains." We see an especially
significant impact in education, where integration of autoformalization into
proof assistants for introductory mathematics and software verification
courses would enable the detection of missing steps or misconceptions in
students’ proofs.

## References

  * Appel and Haken (1977) K. Appel and W. Haken. 1977.  Every planar map is four colorable, part I: discharging. _Illinois Journal of Mathematics_ , 21(3):429 – 490.
  * Appel et al. (1977) K. Appel, W. Haken, and J. Koch. 1977.  Every planar map is four colorable. II: Reducibility.  _Ill. J. Math._ , 21:491–567.
  * Bancerek (2006) Grzegorz Bancerek. 2006.  Automatic translation in formalized mathematics.  _Mechanized Mathematics and Its Applications_ , 5(2):19–31.
  * Bertot and Castéran (2013) Yves Bertot and Pierre Castéran. 2013.  _Interactive theorem proving and program development: Coq’Art: the calculus of inductive constructions_.  Springer Science & Business Media.
  * Carman (2021) Benjamin Andrew Carman. 2021.  _Translating LaTeX to Coq: A Recurrent Neural Network Approach to Formalizing Natural Language Proofs_.  Ph.D. thesis, Ohio University.
  * Csordás et al. (2021) Róbert Csordás, Kazuki Irie, and Juergen Schmidhuber. 2021.  The devil is in the detail: Simple tricks improve systematic generalization of transformers.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 619–634.
  * Dehghani et al. (2018) Mostafa Dehghani, Stephan Gouws, Oriol Vinyals, Jakob Uszkoreit, and Łukasz Kaiser. 2018.  Universal transformers.  _arXiv preprint arXiv:1807.03819_.
  * Ferreira et al. (2022) Deborah Ferreira, Mokanarangan Thayaparan, Marco Valentino, Julia Rozanova, and Andre Freitas. 2022.  To be or not to be an Integer? Encoding Variables for Mathematical Text.  In _Findings of the Association for Computational Linguistics: ACL 2022_ , pages 938–948, Dublin, Ireland. Association for Computational Linguistics.
  * Gonthier (2013) G. Gonthier. 2013.  Engineering mathematics: The odd order theorem proof.  In _Conference Record of the Annual ACM Symposium on Principles of Programming Languages_ , pages 1–2.  Cited By 12.
  * Gonthier (2008) Georges Gonthier. 2008.  Formal proof – the four-color theorem.  _Notices of the AMS_ , 55(11):1382–1393.
  * Gonthier et al. (2013) Georges Gonthier, Andrea Asperti, Jeremy Avigad, Yves Bertot, Cyril Cohen, François Garillot, Stéphane Le Roux, Assia Mahboubi, Russell O’Connor, Sidi Ould Biha, Ioana Pasca, Laurence Rideau, Alexey Solovyev, Enrico Tassi, and Laurent Théry. 2013.  A machine-checked proof of the odd order theorem.  In _Interactive Theorem Proving_ , pages 163–179, Berlin, Heidelberg. Springer Berlin Heidelberg.
  * Gu et al. (2016) Jiatao Gu, Zhengdong Lu, Hang Li, and Victor OK Li. 2016.  Incorporating copying mechanism in sequence-to-sequence learning.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1631–1640.
  * Hales et al. (2017) Thomas Hales, Mark Adams, Gertrud Bauer, Tat Dat Dang, John Harrison, Le Truong Hoang, Cezary Kaliszyk, Victor Magron, Sean McLaughlin, Tat Thang Nguyen, and et al. 2017.  A formal proof of the Kepler conjecture.  _Forum of Mathematics, Pi_ , 5:e2.
  * Hales (2005) Thomas C. Hales. 2005.  A proof of the Kepler conjecture.  _Ann. Math. (2)_ , 162(3):1065–1185.
  * Hendrycks et al. (2021) Dan Hendrycks, Collin Burns, Saurav Kadavath, Akul Arora, Steven Basart, Eric Tang, Dawn Song, and Jacob Steinhardt. 2021.  Measuring mathematical problem solving with the math dataset.  _NeurIPS_.
  * Higham (2020) Nicholas J. Higham. 2020.  _Handbook of Writing for the Mathematical Sciences_, third edition.  Society for Industrial and Applied Mathematics, Philadelphia, PA, USA.
  * Irvine and Deutsch (2021) Andrew David Irvine and Harry Deutsch. 2021.  Russell’s Paradox.  In Edward N. Zalta, editor, _The Stanford Encyclopedia of Philosophy_ , Spring 2021 edition. Metaphysics Research Lab, Stanford University.
  * Leroy (2009) Xavier Leroy. 2009.  Formal verification of a realistic compiler.  _Communications of the ACM_ , 52(7):107–115.
  * Nipkow et al. (2002) Tobias Nipkow, Lawrence C Paulson, and Markus Wenzel. 2002.  _Isabelle/HOL: a proof assistant for higher-order logic_, volume 2283.  Springer Science & Business Media.
  * Pierce et al. (2010) Benjamin C Pierce, Chris Casinghino, Marco Gaboardi, Michael Greenberg, Cătălin Hriţcu, Vilhelm Sjöberg, and Brent Yorgey. 2010\.  Software foundations.  _Webpage: http://www. cis. upenn. edu/bcpierce/sf/current/index. html_.
  * Pierce et al. (2018) Benjamin C. Pierce, Arthur Azevedo de Amorim, Chris Casinghino, Marco Gaboardi, Michael Greenberg, Cǎtǎlin Hriţcu, Vilhelm Sjöberg, Andrew Tolmach, and Brent Yorgey. 2018.  _Programming Language Foundations_.  Software Foundations series, volume 2. Electronic textbook.  Version 5.5. http://www.cis.upenn.edu/~bcpierce/sf.
  * Shaw et al. (2018) Peter Shaw, Jakob Uszkoreit, and Ashish Vaswani. 2018.  Self-attention with relative position representations.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 2 (Short Papers)_ , pages 464–468, New Orleans, Louisiana. Association for Computational Linguistics.
  * Szegedy (2020) Christian Szegedy. 2020.  A promising path towards autoformalization and general artificial intelligence.  In _Intelligent Computer Mathematics. CICM 2020. Lecture Notes in Computer Science_ , volume 12236, pages 3–20. Springer.
  * Tymoczko (1979) Thomas Tymoczko. 1979.  The four-color problem and its philosophical significance.  _Journal of Philosophy_ , 76(2):57–83.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  _Advances in Neural Information Processing Systems_ , 30:5998–6008.
  * Wang et al. (2020) Qingxiang Wang, Chad Brown, Cezary Kaliszyk, and Josef Urban. 2020.  Exploration of neural machine translation in autoformalization of mathematics in Mizar.  In _Proceedings of the 9th ACM SIGPLAN International Conference on Certified Programs and Proofs_ , CPP 2020, page 85–98, New York, NY, USA. Association for Computing Machinery.
  * Wang et al. (2018) Qingxiang Wang, Cezary Kaliszyk, and Josef Urban. 2018.  First experiments with neural translation of informal to formal mathematics.  In _Intelligent Computer Mathematics_ , pages 255–270, Cham. Springer International Publishing.
  * Welleck et al. (2022) Sean Welleck, Peter West, Jize Cao, and Yejin Choi. 2022.  Symbolic Brittleness in Sequence Models: On Systematic Generalization in Symbolic Mathematics.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 36(8):8629–8637.  Number: 8.
  * Wu et al. (2022) Yuhuai Wu, Albert Q Jiang, Wenda Li, Markus N Rabe, Charles Staats, Mateja Jamnik, and Christian Szegedy. 2022.  Autoformalization with large language models.  _arXiv preprint arXiv:2205.12615_.