{'fraction_non_alphanumeric': 0.03712856121719595, 'fraction_numerical': 0.017808638823649546, 'mean_word_length': 4.444901590125653, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Prospects for Hurricane-like Vortices in Protoplanetary Disks

When ice on the surface of dust grains in protoplanetary disk sublimates, it adds its latent heat of water sublimation to the surrounding flow. Drawing on the analogy provided by tropical cyclones on Earth, we investigate if this energy source is sufficient to sustain or magnify anticyclonic disk vortices that would otherwise fall victim to viscous dissipation. An analytical treatment, supported by exploratory two-dimensional simulations, suggests that even modestly under-saturated flows can extend the lifetime of vortices, potentially to a degree sufficient to aid particle trapping and planetesimal formation. We expect the best conditions for this mechanism to occur will be found near the disk's water ice line if turbulent motions displace gas parcels out of thermodynamic equilibrium with the dust mid-plane.

In the consensus theoretical framework, a key location in a protoplanetary disk is the water ice line, which marks the innermost radius at which water vapor and ice can coexist in equilibrium. Depending on the partial pressure of water vapor, the ice line lies at a temperature of ∼ 150 K (e.g., Podolak & Zucker 2004;Lecar et al. 2006;Martin & Livio 2012). At the ice line, if the system is perturbed from phase equilibrium, ice will sublimate to vapor, or conversely, water vapor will deposit onto the silicate surfaces of dust, depending on whether the flow is under-or super-saturated. In the process, the flow is provided with, or deprived of enthalpy given by the latent heat of sublimation of water.
Our goal with this paper is to investigate the possibility that such a thermodynamic disequilibrium between the gas disk and the particle mid-plane can maintain, enhance, and spur the formation of vortices. We draw on the analogy to the well-understood process by which the thermodynamic disequilibrium between the ocean and the atmosphere provides the energy source for tropical cyclones on Earth (Kleinschmidt 1951;Emanuel 1991). There, strong surface winds roil the ocean surface and in the process pick up water which increases the moisture content of the previously under-saturated near-surface air. This exchange transfers enthalpy in form of the water's latent heat of evaporation from the ocean to the atmosphere, and the enthalpy is made accessible as potential energy for thermodynamic transformation. In particular, once the flow has acquired enough H 2 O to achieve saturation, condensation and cloud formation converts latent heat to sensible heat which in turn drives vigorous upwards convection in the eye-wall of the hurricane. This, in turn, decreases surface-level air pressure and accelerates surface-level air into vortical motion via Coriolis deflection around the storm's pressure minimum. The result is strong azimuthal circulation due to conservation of angular momentum, and likewise a transverse circulation due the mass conservation, in which the surface level inflow is balanced by a radial outflow at the tropopause. At this high-altitude (h ∼ 10 − 15 km) level, the atmosphere becomes optically thin in the infrared and the flow can radiate excess heat. Positive feedback is achieved due to faster wind speeds increasing the rate of enthalpy transfer from ocean to atmosphere, thus intensifying convection, ultimately leading to a mature hurricane, which is maintained as long as the hurricane is over the ocean and can draw in sufficiently under-saturated air. In steadystate, the frictional losses, namely dissipation at the ocean-atmosphere interface and radiative losses at the tropopause, are in balance with the latent heat flux.
A similar disequilibrium between ice-coated dust grains in a mid-plane layer and surrounding gas flow in protoplanetary disks will be associated with a source of enthalpy that scales with the level of under-saturation. Our goal is to carry out an initial investigation of the prospects for this mechanism to operate in protoplanetary disks. We model hurricane-like heat engines in two-dimensional hydrodynamic simulations by parameterizing both the thermodynamic disequilibrium and the mixing process.
The paper is structured as follows. Sect. 2 outlines the physical model of latent heat flux driven vortices, applied to protoplanetary disks. In Sect. 3 we present our numerical implementation, while the obtained results are in Sect. 4. Finally, we discuss our work in Sect. 5

VORTEX MODEL
The foregoing discussion suggests that the structure and the thermodynamical properties of protostellar disks allow one to draw strong physical analogies with terrestrial hurricanes. In this section we elaborate the specific details of a model that is based on this analogy.

Latent heat flux
When a vortex near the ice line of a protoplanetary disk mixes with icy dust grains from the mid-plane, any under-saturation will lead to sublimation which increases the flow's water vapor content, and thus is associated with a net enthalpy transport from dust at the mid-plane to the gas flow that is proportional to the specific latent heat of water sublimation L s . In analogy to Hurricanes, we characterize the bulk latent heat flux LHF (units erg/s) over a surface dA as (compare to e.g., Liu et al. 2011) where |v| is the flow's wind speed, and ρ is the total gas density. C E is the dimensionless turbulent exchange coefficient which describes how efficient humidity is transferred to the flow. Note that the turbulent exchange coefficient itself is believed to depend on flow velocity (e.g., Cronin et al. 2019). We discuss this further in Sect. 5.2. ∆q is the under-saturation of the flow given by where q is the vapor mixing ratio (also known as specific humidity) of the under-saturated flow that mixes with the dust-layer. It is defined as the ratio of vapor density (absolute humidity) ρ v over total gas density ρ, i.e.
where ρ d is the density of "dry" gas. In Eq.
(2), q * is the vapor mixing ratio at saturation.

Phase equilibrium and saturation
Phase equilibrium between water ice and vapor is characterized by the relative humidity also known as saturation ratio, being equal to unity, i.e. RH = 1. Here, P v is the pressure of the water vapor phase, for an ideal gas given by P v = ρ v R v T . It is thus a partial pressure of the total gas pressure P = ρR g T , where R v = R/M H2O = 4.615 · 10 6 ergs/g/K, and R g = R/M g = k B /(µm p ) = 3.568 · 10 7 ergs/g/K are specific gas constants for water vapor and the total gas respectively, with mean molecular weight of the gas µ = 2.33 (assuming solar composition), molar mass of water M H2O = 18.015 g/mol, proton mass m p , universal gas constant R, and Boltzmann constant k B . P v, * is the saturation vapor pressure. We can deduce P v, * from the curve of vapor-ice coexistence for water, the slope of which is given by the Clausius-Clapeyron relation (Clausius 1850) dP where ∆v is the change in specific volume associated with the phase transition, and L s is the latent heat of water sublimation. Following Ros et al. (2019), we use  Figure 1. Visualization of key properties of a steady-state disk between 1 AU and 5 AU. Our fiducial model places the water ice line at ∼ 3.6 AU (blue line). Shown are as a function of distance to star R, from left to right, top row: temperature profile T and surface density profile Σ based on the MMSN (Hayashi 1981), disk aspect ratio H/R, mid-plane gas density ρ0 and (total) water density using a fiducial H2O abundance of 1.2 · 10 −4 (Furuya et al. 2013). Bottom row: total pressure P and partial vapor pressure Pv which is the minimum of saturation vapor pressure Pv, * and maximum partial pressure Pv,max = ρH 2 ORvT , water vapor and ice density, relative humidity RH (compare to Ros et al. 2019), and lastly vapor mixing ratio at saturation q * and actual mixing ratio q. We argue that an inwardly transported fluid parcel would be under-saturated with relative to its new environment.
This implies for the saturation vapor mixing ratio, again using an ideal equation of state, It is important to highlight that P v, * and q * are not physical properties of the system state, but only quantify partial pressure P and mixing ratio q respectively necessary to achieve phase equilibrium.
The system can only establish equilibrium, however, if there is enough water available to satisfy P v = P v, * . Thus, the (water) iceline can be defined as the furthestin location where the saturation criterion can be met. Outside the iceline, excess vapor will be deposited onto dust grains, such that P v = P v, * ⇔ RH = 1 to maintain phase equilibirum (Ros et al. 2019). Inside the ice line, all the available H 2 O is in the vapor phase, yet the resulting partial pressure is less than the required P v, * , resulting in a phase disequilibrium, in particular in an under-saturated flow.
The iceline location depends on temperature and density profiles and the water abundance. We proceed with a fiducial disk model where temperature scales as and the gas surface density scales as with distance to star R. For the MMSN (Hayashi 1981), R 0 = 1 AU, T 0 = 280 K, β T = 1/2, Σ 0 = 1000 g/cm 2 , and β Σ = 1. Vertical hydrostatic equilibrium implies that the mid-plane gas density is ρ 0 = Σ/( √ 2πH). H = c s /Ω is the gas pressure scale height given by the ratio of local speed of sound c s = γk B T /(µm p ) and Keplerian angular frequency Ω = GM /R 3 . Here, we introduced adiabatic index γ = 5/3, gravitational constant G and solar mass M . With these power law scalings, we can write the vapor mixing ratio at saturation q * as a function of R as For the water abundance, we assume a fiducial value of n H2O /n H = 1.2 · 10 −4 (Furuya et al. 2013), which implies for the total (vapor + ice) water density ρ H2O = (n H2O /n H )(M H2O /M H )ρ 0 . The vapor pressure is therefore with the iceline located where the two arguments of the minimum function are equal. The water vapor and ice densities are then just An alternate approach to calculating ice and vapor densities is to assume an ice fraction in dust grains (of order 0.5 in models by Min et al. 2011), and a dust-to-gas ratio (of order 10 −2 in the ISM, Savage & Jenkins 1972). This leads to a water abundance of n H2O /n H ∼ 2 · 10 −4 that is consistent with our chosen value. Fig. 1 shows the radial variations of the properties introduced in this section. For the utilized MMSN disk model, the iceline is located at R 0 ∼ 3.6 AU, which is somewhat further out than what is obtained in more sophisticated considerations of the temperature profile (e.g, Kennedy & Kenyon 2008;Martin & Livio 2012).
We note, that the microphysics that govern H 2 O deposition and sublimation in disks are not well constrained (for a review see e.g., Henning & Semenov 2013). We assume that the associated timescales are short compared to gas transport. In addition, our model neglects chemical reactions. For example, UV photons can photodissociatie H 2 O ice (Gerakines et al. 1996), in particular if dust grains are lifted to the disk surface via turbulent mixing (Furuya et al. 2013). Likewise, water formation through, e.g. hot neutral chemistry (e.g., Du & Bergin 2014) is not considered. As a result, in our model, the total water density can only change via radial advection.

Phase disequilibria due to radial mixing
Indeed, any such radial transport of gas parcel, for example epicyclic oscillations or turbulent motions, introduces thermodynamic disequilibrium. In particular, if a fluid parcel at R 1 in phase equilibrium with vapor mixing ratio q 1 = q * (R 1 ) were to displaced radially inward to radius R 2 < R 1 , it would be under-saturated relative to its new environment where the ambient saturation mixing ratio is q * (R 2 ) > q 1 . The parcel would move towards phase equilibrium by sublimating water ice, which would increase the flow's enthalpy in form of latent heat of water sublimation to the flow according to Eq. (1). Since the amount of available ice -a requirement for any latent heat transfer -increases outwards, but the gradient in saturation vapor pressure increases inwards (see Fig. 1), the location just outside the water iceline constitutes a sweet spot where the achieved latent heat flux is maximal. For example for our fiducial model, a radial dispacement of one gas scale height leads to ∆q ∼ 10 −4 just outside the ice line, but ∆q ∼ 10 −14 at 10 AU which is negligible.
In our analysis, we remain largely agnostic to the origin of radial displacements causing disequilibria, and treat ∆q as a free parameter and will investigate values of 10 −6 ≤ ∆q ≤ 10 −4 .

Vortex energetics
The latent heat flux increases the water vapor mixing ratio and thus the total enthalpy of the flow where C P = dh/dT is the specific heat capacity at constant pressure, and we assumed the ice content to be neglible. The first law of thermodynamics for a canonical ensemble dh = T ds + 1/ρdP , with specific entropy s, becomes Using Bernoulli's equation, the full energy equation is given by (Emanuel 1991) The first term represents the total internal energy of the flow including from left to right the kinetic energy, the gravitational potential energy, the sensible heat, and the latent heat stored in water vapor. The second term is the heat input that must be provided by the ice nuclei, and the third term is any frictional dissipation of energy along the flow path dl. This may entail shedding of Rossby waves, and any turbulent dissipation induced by (for example) Keplerian shear. We constrain our first-step investigation to optically thick disks, which are primarily heated through stellar irradiation rather than accretion heating (Armitage 2015). As such, the disk vortices that we model do not experience radiative losses.
In a steady-state, the flow's internal energy is constant, and integration of Eq. (15) implies that along any stream line, heating is balanced by friction The heating is provided by the enthalpy flux associated with the change in vapor mixing ratio as the flow picks up water vapor (compare with Emanuel 1991). Using the first law in Eq. (14), the heating can be written Here, we have assumed that while the flow itself is adiabatic, its mixing with the dust layer proceeds isothermally. In other words, the disequilibrium is characterized by phase imbalance alone, not by a temperature gradient between ice grains and gas flow. Eq. (16) implies that along a streamline and thus, if the streamline is closed, In a steady-state, the energy dissipation rate of the vortex is exactly balanced by the latent heat flux, i.e.
where A is the total area of contact between vortex flow and ice layer.
In order for this mechanism to operate, the acquired latent heat of sublimation -an internal energy -must be converted into work. We model this "piston" in an ad-hoc manner by prescribing an increase in the flow momentum, such that ∂v 2 /∂t ∝ LHF. This is in direct analogy to terrestrial hurricanes, where latent heatenhanced convection locally increases the atmospheric scale height, thereby lowering the ambient pressure and advecting in more air. We draw further attention to Sect. 5.1, where we discuss a proposed three-dimensional structure of this mechanism in protoplanetary disks that includes convection.

Summary of model assumptions
It is worth reiterating the assumptions that inform Eq. (20).
• The enthalpy transfer from the dust mid-plane to gas flow is driven by the thermodynamic disequilibrium associated with an under-saturation of the gas given by ∆q. The efficiency of the latent heat transfer characterized by constant turbulent exchange coefficient C E (see Sect. 5.2).
• The thermodynamic disequilibrium, specifically the under-saturation ∆q can be maintained. For closed streamlines, this requires the transport of material to regions where deposition and desaturation occurs, i.e. to greater R where q * is smaller. For open systems, this can be achieved if the vortex can draw-in new under-saturated gas (compare to terrestrial hurricane in Emanuel 1991).
• The dust-midplane acts as a heat bath that keeps the flow at constant temperature as it mixes with the flow. The mixing process is therefore isothermal.
• Along a closed streamline, d ln P = 0, i.e. any expansion is followed by compression of equal magnitude. A special case satisfying this assumption is a flow, in which heat input is exclusively associated with an increase in internal energy and not with any work (or temperature), i.e. an isothermal and isobaric phase transition. The increase in internal energy is associated with an increase in flow speed to achieve positive feedback.

NUMERICAL SETUP
We employ Athena++ 1 (Stone et al. 2020) which solves hydrodynamic equations using a high-order Godunov scheme. Hereby, the hydrodynamic equations are solved on a Eulerian grid, and in the shearing sheet reference frame (Stone & Gardiner 2010). The x, and y coordinates correspond to radial and azimuthal directions respectively. We utilize an adiabatic equation of state. As a first step, we perform two-dimensional shearing sheet simulations with a fiducial N x ×N y = 256×256 mesh. In Sect. 4.1 we extend the calculations to explore variations in the grid resolution. The physical dimensions are L x = L y = 2H.
The full system of equation reads where Terms containing the Keplerian angular frequency Ω are associated with the shearing sheet and include the Coriolis and centrifugal forces, as well as the linearized background shear −3/2Ωx. In the momentum equation, f damp is a damping coefficient associated with soak zones, and f heat is the driving coefficient corresponding to the heat-work conversion outlined in the previous section. We construct f heat to distribute the momentum provided by the enthalpy of sublimation such that v x /v y = ∆v x /∆v y . Thus in each zone, there is an increase in specific momentum given by where LH = LHF · dt = L s C E ρ|v|∆qdAdt is the specific latent heat acquired in a given time step dt within area dA. In each zone, we set dA = L x /N x · L y /N y = 6.1 · 10 −5 H 2 as the domain of contact with the dust midplane. Deviations of the boundary layer area from this simple shape are absorbed in the turbulent exchange coefficient, C E .

Scaling Relations
In code units, we set Ω −1 , c s , the unperturbed background density ρ 0 , and the unperturbed temperature T 0 to unity, which implies that the gas scale height, H, the unperturbed pressure, P 0 , and the reference volume, V = H 3 are also unity in code units. Using Eq. (8), and µ = 2.33 the sound speed can be expressed in cgs units with c s = 9.96 · 10 4 cm s The energy per unit mass corresponding to a value of unity in code units is given by c 2 s /2 = 4.96 · 10 9 (R/1AU ) −1/2 erg/g. We can therefore express the latent heat of water sublimation, L s = 3 · 10 10 erg/g (Datt 2011) in terms of code units as

Vortex initial conditions
Adams & Watkins (1995) provide circular symmetric vortex solutions for a vortex in geostrophic balance in the limit of zero shear. Similarly to Seligman & Laughlin (2017), we seed our simulation with an extended vortex which is characterized by a vortex strength (circulation) of Γ (units cm 2 /s) and a characteristic radius R c , within which the vortensity (vorticity over surface density) is set constant. If is the radial coordinate, and k R = Ω/c s = 1/H the Rossby wavenumber, the velocity profile is given by (Adams & Watkins 1995) where I α and K α are modified Bessel functions of the first and second kind respectively. The corresponding density profile that in the absence of shear would establish geostrophic balance is with stream function and Density and velocity profiles for Γ = 0.1, and R c = 0.1H are shown in Fig. 2. Azimuthal velocities are maximal at = R c , which is also the location of the steepest density gradient.
The solution of Adams & Watkins (1995) is stable in the absence of shear. Simulations by Godon & Livio (2000); Seligman & Laughlin (2017) indicate that in the presence of background shear, the vortex adjusts, sheds Rossby waves, and is forced into an oval shape. We discuss this readjustment further in Sect. 4.1. Work by Lyra (2021) has recently presented disk vortex solutions in elliptical coordinates, but an analytical solution that incorporates an equilibrium vortex into a background Keplerian shear has not yet been found. Indeed, the elliptical instability (Lesur & Papaloizou 2009) caused by resonance of vortex turnover time and epicyclic frequency is expected to trigger for all two-dimensional vortices in disks, suggesting a physical equilibrium solution may not exist.

Boundary conditions and soak zones
We adopt boundary conditions that are periodic in y and shear-periodic in x. In the absence of a latent heat flux, we expect seeded vortices to dissipate over time due to shedding of Rossby waves and friction with the shear flow generated by numerical viscosity inherent to the computational method. (We do not include an explicit viscous dissipation.). In order to prevent Rossby waves from propagating across the shear-periodic boundaries, reappearing in the domain, and interfering with the vortex, we dedicate the 5% of the grid cells closest to the radial simulation boundaries as soak zones. These zones act to remove momentum from the simulation, characterized by f damp , thereby improving the isolation of the vortex.

Diagnostic measures
As a measures of vortex dissipation, we take the maximum value of the flow's vorticity as well as total kinetic energy E kin . We expect our vorticities to decay exponentially (compare to e.g., Godon & Livio 1999). Hence, we define the exponential decay time τ via where ω 0 is the vortex's initial maximum vorticity and ω eq some equilibrium vorticity.

Vortex adjustment and dissipation
Before including the heat engine implied by Eq.
(1), we test our numerical setup by evolving the vortex initial condition (with Γ = 0.1, R c = 0.1H) for a variety of different resolutions for 100 orbits and measure the ensuing vortex dissipation rate. Evolution of vorticity and total kinetic energy as a function of grid resolution is shown in Fig. 3. The lowest resolution simulation cannot maintain the vortex. Geostrophic balance is never established and the vortex effectively vanishes within 10 orbits. (Note that there is some vorticity even in the absence of the vortex due to the shear flow.) The other simulations maintain a vortex over the course of 100 orbits. All vortices lose kinetic energy and vorticity, the rate of which decreases smoothly with increasing resolution. The right panel of Fig. 3 includes two snapshots the vorticity of the highest-resolution simulation 512 × 512. Dissipation is manifested in drag (decrease in maximum vorticity) and diffusion (increase in radial extent).
The resolution also impacts the adjustment time required for the vortex to adapt to the Keplerian shear flow t adj . We estimate the relative energy dissipation rate viaĖ where t max is the vortex dissipation time of 10 orbits for the 64 × 64 simulation, and the simulation runtime of 100 orbits for the others. The resulting dissipation rates are also shown in the left panel of Fig. 3. We fitted an exponential decay of formĖ/E = m exp(−tN ) + c. The dissipation rates are converging to a value of orderĖ/E ∼ 10 −4 Ω. Given the dissipation rates in Fig. 3, Eq. (20) asserts an order of magnitude for the required under-saturation ∆q req to maintain the flow, i.e. ∆q req ∼Ė/(ρ|v|AC E L s ). Our seeded vortex has a total kinetic energy of order E ∼ 10 −4 in code units, implyingĖ ∼ 10 −8 for the fiducial simulation (see Fig. 3). For C E = 0.1, ρ ∼ 1, a vortex area of A ∼ R 2 c = 10 −2 H 2 , and velocities of order |v| ∼ 10 −2 c s therein, the required under-saturation that is to be sustained is of order ∆q ∼ 10 −4 , which is modest. Fig. 3 suggests that simulations with resolution less than 512 × 512 are not converged, and the observed vortex decay is mainly produced by numerical dissipation. This numerically-induced effective friction is therefore present in our fiducial 256 × 256-resolution simulations, and we proceed with the ansatz that the numerical dissipation in the calculations would be mirrored by an effective physical viscosity within a real protostellar disk.
Indeed, if the dissipation of a real disk vortex is driven by diffusion at the interface between the vortex flow and the Keplerian shear flow, then we expect the dissipation time, t d , to relate to the diffusion time. For a vortex of size, s v , and (numerical) kinematic viscosity ν, this is given by where we have expressed ν in terms of the commonly used α-parameterization (Shakura & Sunyaev 1973) Using s v = R c = 0.1H, Eq. (35) implies an effective turbulence parameter of α ≈ 4 · 10 −4 , which is in line with the dissipation required to generate observed time scales for protostellar disk evolution (e.g., Flock et al. 2017). It is also consistent with the results of Godon & Livio (1999), who found that the exponential decay time, τ , depends on the viscosity parameter α. For α = 10 −4 , the anticyclonic vortices in Godon & Livio (1999) propagate for τ ∼ 50 orbits, a duration comparable to that observed in our fiducial simulation with 256 × 256 resolution.    Fig. 3, but for a fixed resolution of 256 × 256. Each color corresponds to a different under-saturation ∆q which drives the heat engine that is turned on after 5 orbits. The gray shaded region marks this period in which the vortex adjusts to the shear flow. The right panel includes two snapshots of the vorticity of the simulation with the strongest under-saturation. Note that the vorticity colorbars have a different scaling to those in the right panel of Fig. 3, reflecting the now much stronger vortices.

Vortex preservation and strengthening
Specifically we examine the influence of water vapor under-saturations in the range 10 −6 ≤ ∆q ≤ 10 −4 . Simulation run-times are 20 orbits. For simulations with ∆q 10 −5 , the vortex is maintained and both vorticity and vortex kinetic energy increase with time. The efficiency with which the vortex strengthens is somewhat less than what our order of magnitude estimate in Sect. 4.1 predicted. While vortices in the simulations with ∆q 10 −5 dissipate (albeit at a slower rate than in the reference simulation without heat engine) the vortices seeded into simulations with ∆q 10 −5 increase both their vorticity and their size. An example of this growth can be seen for ∆q = 10 −4 in the right panel of Fig. 4. Increasing vorticity and vortex velocities also lead to a relative increase in viscous dissipation of vortical energy, which is why neither kinetic energy nor vorticity increase linear in ∆q.
We deliberately chose under-saturations ∆q to keep velocities sub-sonic within the run time of 20 orbits. This is in contrast to Les & Lin (2015) who investigated vortices in the context of disk gaps. Their vortices are formed and maintained by the the Rossby Wave instability (Lovelace et al. 1999) and exhibit shock-inducing velocities, which generate dissipation. Our vortices would undergo a similar fate if they received energy input sufficient to drive super-sonic velocities. For example, the vortex in the simulation with ∆q = 10 −4 (yellow line in Fig. 4) increases its peak velocities at a rate of 7.5 · 10 −3 c s /orbit. If maintained, such boosting would lead to super-sonic velocities at a run time of about 120 orbits.

Vortex formation via latent heat flux
Anticyclonic vortices are believed to form through a number of disc processes. Suggested mechanisms have included the Rossby wave instability (Lovelace et al. 1999) and the Vertical Shear Instability (VSI, e.g., Nelson et al. 2013;Manger & Klahr 2018). The latter phenomenon draws its energy from the vertical gradient in orbital frequency caused by radial variations in temperature and entropy. In this section, we explore the potential of the latent heat flux to generate vortices by tapping the thermodynamic disequilibrium between the gas flow and the icy layer in the presence of a background relative drift between the two layers. Note that because the latent heat flux is driven by a phase disequilibrium, and not by a temperature gradient, it is distinct from the VSI or the thermal instability investigated by Owen (2020).
So far, we have treated the dust layer as a distinct laminar surface. In a more realistic disk, the coupling of dust and gas (characterized by the stopping time t s , a proxy for the particle size) leads to an equilibrium drift of particles both azimuthally and radially, even if v = 0 in the Keplerian frame. For a pressure gradient and local dust-to-gas ratio, µ, the equilibrium relative streaming velocity w in the Keplerian frame is given by (Nakagawa et al. 1986;Squire & Hopkins 2018) where v K = RΩ and St = t s Ω is the dimensionless stopping time or Stokes number. Typically, St ∼ 10 −3 − 10 −2 1 (see e.g., Testi et al. 2014), and our subsequent discussion assumes that the stopping time is in this regime. Vertical stratification of the particle mid-plane leads to height-dependent dust-to-gas ratio, µ, and an attendant relative velocity w that decreases with height: The dust-rich mid-plane is less coupled to the gas flow and has a near-Keplerian orbital velocity in equilibrium. By contrast, the upper layers are more coupled to the gas flow and orbit close to the sub-Keplerian gas velocity of (1 − η)v K .
This intrinsic vertical shear has long been thought to drive Kelvin Helmholtz Instability (KHIs, Weidenschilling 1980), which mix the dust layer with the gas flow, and in the process prevent settling of the dust to a razor-thin layer. We envision particle mixing in the vortex to proceed in a similar manner, as particles in the dust-rich mid-plane likewise tend to collectively decouple from the gas flow even it deviates from the equilibrium sub-Keplerian stream lines. This decoupling becomes particularly prevalent for dust-rich disk regions with super-solar metallicities. There, the KHI becomes increasingly ineffective in turning over the nowmassive particle stream, leading to a development of a high-density mid-plane cusp (Sekiya 1998;Youdin & Shu 2002;Gómez & Ostriker 2005). We discuss mixing in a disk subject to KHI further in Sect. 5.2.
The underlying equilibrium stream in Eq. (38) imposes asymmetry onto an anticyclonic vortex supported by the Keplerian shear. The relative velocity between particles and vortical flow and thus the latent heat flux in Eq. (1) is decreased (increased) on the side towards (away from) the star and trailing (leading) the vortex's orbit, depending on if the Nakagawa drift is parallel (or anti-parallel) to direction of the vortex stream lines.
For the vortices simulated here, peak velocities are around 10% of the sound speed. For fiducial values of St = 0.1, µ = 0.02, η = 0.005, and v K = 0.05 c s , the radial and azimuthal components of the background drift, w are of order 1% and 0.1% of the sound speed respectively. These deviations are an order of magnitude smaller than vortex flow velocities, and are thus of minor importance for the overall structure of the flow.
For weaker vortices, however, the background drift provides a significant contribution to the flow speed that drives mixing. In fact, even without a stable preexisting vortex, the background drift alone may be able to bootstrap a vortex flow, if the gas flow is sufficiently under-saturated. In order to investigate this we seeded a simulation with a random vortex field consisting of 100 vortices of with circulation 0.01 ≤ Γ ≤ 0.2 and size 0.01 ≤ R c /H ≤ 0.1, uniform-randomly distributed. We test this initial condition for setups with background drift both enabled and disabled, and compare to a reference simulation without latent heat flux. The background drift is implemented only insofar that the relative velocity in Eq. (38) is added to the relative velocity in Eq.(1). Hereby, we choose St = 0.1 and µ = 0.02 which correspond to marginally coupled dust particles and a typical dust-to-gas ratio (e.g., Birnstiel et al. 2012). We test two under-saturations of ∆q = 10 −5  Figure 5. Evolution of the kinetic energy for random vortex initial condition for five different setups. Except for the run where LHF is turned off (red line), the adjustment time is set to one orbit (gray shaded region). Snapshots throughout the evolution of the absolute vorticity are shown for the simulation with ∆q = 10 −4 and drift included. and ∆q = 10 −4 , both of which were enough to maintain the vortex in Fig. 4. For the two simulations that include the latent heat flux, the flux is enabled after a vortex adjustment time of one orbit. Simulations were run for a total duration of 20 orbits. Fig. 5 shows the evolution of the normalized kinetic energy. In the reference simulation where LHF is disabled, the initial condition dissipates and no lasting vortices survive (or form). Simulations with latent heat flux and drift (solid lines in Fig. 5) are able to develop vortices from the initial condition. Fig. 5 includes four snapshots of the vorticity throughout the evolution of the ∆q = 10 −4 simulation. At 1 orbit, when the latent heat flux is turned on, the initial vortex field has been sheared out by the Keplerian flow and a few anticyclonic vortices can be identified. Cyclonic vortices have been destroyed by the shear flow. The surviving vortices are subsequently strengthened by the latent heat flux, interact with each other, and ultimately merge, forming one large vortex spanning about a scale height in radial direction at 20 orbits. The maximum vorticity increases about a factor of three. The simulations with latent heat flux turned on follow the precedent set in Fig. 4 where higher under-saturations, ∆q, allow for stronger vortices. This is also seen in simulations without drift (dashed lines in Fig. 5). Only the simulation with the greater latent heat flux is able to produce stable vortices, while the simulation with ∆q = 10 −5 dissipates in a manner similar to that observed in the reference simulation.

DISCUSSION
Our analytical considerations and a variety of rudimentary numerical simulations point toward the possibility that hurricane-like vortices may exist in protoplanetary disks. In Sect. 4.1, we measured dissipation rates of the vortex initial condition in Adams & Watkins (1995) for a number of resolutions. In Sect. 4.2, we quantified the thermal disequilibrium between undersaturated gas flow and ice-coated dust required to maintain, or strengthen a dissipating vortex. As predicted in Eq.(20), the required disequilibrium and the associated latent heat flux depends on the dissipation rate of the vortex. In Sect. 4.3, we investigated the potential of the latent heat flux to form strong vortices from a random vortex initial condition.
To better understand the applicability of our work to real protoplanetary disks, it is thus worth discussing applicability of our model to real disks and the feasibility of latent heat fluxes required to balance expected vortex dissipation rates. We also discuss avenues for a mixing model, and the impact on dust trapping and planetesimal formation.

Vortex streamlines in two and three dimensions
We restricted our first-step numerical approach to two dimensions, yet hurricanes on Earth are characterized by three dimensional streamlines. In addition to the strong  (Emanuel 1991). y = 0 is the center of the vortex. The bottom part shows the run of pressure and circular velocities as seen in Fig. 6. vortical flow, there is a transverse circulation. There, near-surface air streams radially inwards and is subject to the latent heat flux. The conversion of latent to sensible heat via rainfall drives vigorous upwards convection in the eyewall, which decreases the central pressure thus drawing in more ambient air. This strengthens the winds and thus also the latent heat flux, thus leading to the growth. In the eye wall, the warm air expanses upwards and cools until the temperature of the convective flow transitions smoothly into the ambient temperature profile. In order to mechanically maintain the flow, and allow under-saturated air into the system, there is a radial outflow at the tropopause, where excess heat is radiated away. With this mechanism in mind, it is worthwhile considering plausible forms for three-dimensional streamlines that create analogous heat engines in protoplanetary disks.
The key difference between the Earth case and the protoplanetary disk case is the strong Keplerian shear in the latter, which rapidly tears radially extended flow structures into azimuthally elongated bands, and also rapidly dissipates retrograde (cyclonic) vortices (see e.g., Godon & Livio 1999). By contrast, prograde (anticyclonic) vortices are supported by the shear flow and dissipate significantly more slowly in the absence of an energy source. Moreover, there are a number of physical mechanisms that have long been posited to produce large-scale vortices (see Sect. 4.3). The novelty of our approach derives from our focus on the applicability of the latent-heat flux driven heat engine to anticyclonic vortices.
An anticyclonic vortex in geostrophic balance is a region of increased pressure relative to the ambient gas pressure in the disk. Fig. 6 shows a pressure map of the final snapshot of the ∆q = 10 −4 simulation in Fig. 4 with velociy flow field plotted on top. While pressure in the vortex center is indeed higher than the ambient pressure (required to maintain geostrophic balance), the pressure maximum occurs at around (x, y) ∼ (0, 0.5H), where circular velocities are high and the shear contribution is zero. It is here where we draw the analogy to a Hurricane on Earth and expect convection away from the mid-plane to mechanically maintain an additional transverse circulation.
We show a schematic of our proposed transverse circulation in Fig. 7. In analogy to the hurricane on Earth, the vortex flow picks up enthalpy of water sublimation by mixing with icy dust grains as gas parcels spiral towards the local minimum where azimuthal velocities are expected to be maximum (A). Note that we do not expect the inflow to proceed axisymmetric, but mainly azimuthally displaced from the vortex center as this where the pressure gradient ∂P/∂ is maximal (see Fig. 6). In order for the additional enthalpy to drive vertical convection (B) and maintain the pressure gradient, it must be energetically favorable for the additional heat to be carried dynamically rather than radiatively. For convective instability (e.g., Lin & Papaloizou 1980) where the radiative temperature gradient is given by (e.g., Kippenhahn & Weigert 1990) where a = 7.57 · 10 −15 erg cm −3 K −4 is the radiation constant, κ = κ(T ) the opacity, and l/m the energy flux per unit mass. For convection, it must satisfy where we have inserted T = 150 K, P = 10 erg/cm 3 (see Fig. 1), and a fiducial mean opacity of order 1 cm 2 /g corresponding to an optically thick mid-plane (e.g. Henning & Stognienko 1996). The latent heat flux per unit mass is of order LHF/m ∼ L s C E |v|∆q/L, where L is a characteristic length scale over which the latent heat is distributed. For C E = 0.1, |v|/c s = 0.1, ∆q = 10 −4 , L/H = 0.1, and at R = 3.6 AU, our fiducial model yields LHF/m ∼ 0.1 erg/g/s. Given these rough estimates, this energy flux density would satisfy Eq. (41) by two orders of magnitude, and thus drive convection.
Once the disk becomes optically thin, however, and κ decreases, ∇ rad increases and excess heat can be radiated away. The flow approaches thermodynamic equilibrium with the stellar radiation field, and loses the ability to drive convection ultimately leading to a radial outflow (C in Fig. 7), which mechanically maintains the structure. We expect the vortex center -in analogy to the "eye" of the terrestrial hurricane -to be relatively quiet, as rotational velocities, as well as corresponding latent heat fluxes, are minimal. Part A, B, and C in Fig. 7 correspond to the first three legs of the Carnot cycle that characterizes a mature hurricane (Emanuel 1991). Note that the convective instability we invoked to mechanically maintain the vortex can be produced by a variety of mechanisms in protoplanetary disks. Accretion, for example driven by the magneto-rotational instability (Balbus & Hawley 1991), can also render optically thick disks convective (Garaud & Lin 2007;Jankovic et al. 2021), which would relax the additional requirement that Eq. (39) poses on the latent heat flux.

Mixing model
In our parameterized model of the latent heat flux in Eq. (1), the mixing efficiency is described by the turbulent exchange coefficient C E , which describes the efficiency with which water vapor is mixed into the flow thereby enabling it to transfer its latent heat. Past experimental measurements of C E (e.g., Verma et al. 1978) while difficult to transfer to the protoplanetary disk context, have been accompanied by theoretical considerations in the Earth context. In this section, we will build upon this framework to calculate an order-of-magnitude analytical prediction of C E in the protoplanetary disk context, thereby motivating our choice of C E = 0.1.
For a latent heat flux from Earth's ocean to the atmosphere, the latent heat transfer coefficient above the interface is known to depend on wind speed and sea roughness, i.e. wave height, period, steepness, as well as the relative orientation of the wind velocity vector to the direction of wave propagation (e.g., Cronin et al. 2019). In addition, the transfer efficiency is expected to decrease with height above the boundary layer as turbulent eddies cascade to decay in the presence of gas viscosity. This decrease is typically modeled using a semi-empirical logarithmic law, and modified by functions ψ developed in Monin-Obukhov Similarity (MOS) theory (Monin & Obukhov 1954), which handle complexities stabilizing (stratifying) or destabilizing (convective) conditions. The wind speed at height z above the boundary layer then can be written as where κ f ≈ 0.4 is the von Kármán constant, z 0 is the interface's roughness length, and v f the friction velocity which in turn depends on shear stress τ and fluid density ρ via v 2 f = τ /ρ. Given a roughness length z 0 and MOS correction term, the turbulent exchange coefficient for latent heat can be expressed in the form (see Garratt 1994) In a protoplanetary disk, the boundary between dust grains and gas fluid is not distinct as it is with the oceanair interface on Earth. Instead, the pressure gradient, η, prescribes the intrinsic scale height of the particle layer as regulated by KHI (e.g., Chiang 2008;Gerbig et al. 2020). As a consequence, even in equilibrium, gas and dust are not well-separated flows, but form a combined fluid, with a vertical gradient in dust-to-gas ratio and resulting drift velocities. Our goal is to use the established theory of KHI in disks to construct a function ψ(z, z 0 ) that satisfies the semi-logarithmic wind profile in Eq. (42) allowing us to calculate an exchange coefficient via Eq. (43).
A flow with velocity component v(z) = (v 2 φ + v 2 R ) 1/2 perpendicular to a gravitational acceleration g = gẑ is subject to KHI if its Richardson number falls below a critical number, for Cartesian flows given by Ri = 1/4 (Chandrasekhar 1961). In equilibrium, and For the purposes of our estimate, we evaluate the exchange coefficient C E at z = ηR, which is a typical length scale for the particle layer (e.g., Gerbig et al. 2020), and notice that for Richardson number of order unity (Johansen et al. 2006), f 1/4 too is of order unity. Thus, In this form, the (equilibrium) turbulent exchange coefficient depends on the local disk properties pressure gradient η, aspect ratio H/R, viscosity α, the flow Richardson number Ri, and the ratio of roughness scale z 0 to characteristic length scale ηR. Since we did not consider potential vortex modifications to the velocity profile but used the equilibrium velocity profile v 0 , Eq. (53) likewise is a prescription for the equilibrium exchange coefficient. However, we can use the roughness length z 0 to incorporate the effect of faster velocities, as faster wind speeds are expected to increase the roughness length. Fig. 8 shows the turbulent exchange coefficient at a height z = ηR above the mid-plane vs roughness length z 0 for multiple values of α. Our choice of C E = 0.1 is achieved for high z 0 (of order ηR/2) and moderate α-values. Note that while high α-values lead to more efficient mixing, and thus to an increase in latent heat flux and energy input, they also would cause the vortex to decay at a faster rate (compare to e.g., Godon & Livio 1999). The overall impact of an external α is therefore not obvious, and will be studied in more detail in future work.

Dust trapping and planetesimal formation
Anticyclonic vortices are associated with maxima in gas pressure and are thus commonly invoked as a dust trap (e.g. Heng & Kenyon 2010;Raettig et al. 2015;Owen & Kollmeier 2017) that can break the meter barrier and lead to planetesimal formation (see Chiang & Youdin 2010, for a review). In models of planetesimal formation that rely on such a dust trapping mechanism, the lifetime of traps can be a bottleneck (e.g., Lenz et al. 2019;Gerbig et al. 2019), since vortices in accretion disks have finite lifetimes in the absence of any support mechanisms (Lesur & Papaloizou 2009) Therefore, if the heat engine successfully operates in disks and elongates vortex life times, we expect the framework discussed here to provide a zeroth-order favoring of planetesimal formation. Planetesimal formation feedback, however, would remove dust grains from the disk, in the process decreasing the surface available for water vapor sublimation and water ice deposition. This may shut off the heat engine, much as a terrestrial Hurricane is shut off when it hits land. In addition, if dust feedback is considered, the inertia of the captured dust stream can also destroy a vortex (Fu et al. 2014). Similarly, the convection we invoked to maintain the heat engine may displace particles vertically, which can disrupt regions of gravitational collapse, similar to the simulations studied in Schäfer et al. (2020); Gole et al. (2020), where gas turbulence hinders sedimentation. Another point of consideration is that the phase disequilibrium driving the vortex requires that vapor and water ice co-exist in the immediate environment of the vortex. This constrains the ideal location to just outside the ice line. Promisingly, Pfeil & Klahr (2021) found that vortices can be found throughout the disk, including the regions around the iceline. On the other hand, in global simulations by Flock et al. (2020), a large scale vortex forms at 30 AU. Such a structure would be wholly unaffected by latent heat flux.
The overall impact to planetesimal formation in the context of vortices as dust trap location is therefore not obvious, and must be investigated in more detail -ideally with numerical simulation that include dust particles.

Conclusions
Terrestrial hurricanes have been studied intensively (in part as a consequence of their fearsome economic importance (e.g., Emanuel 2005)) and they constitute truly remarkable self-organized systems powered by a flux of latent heat and maintained by the thermodynamic disequilibrium between the ocean surface and the atmosphere. Moreover, it seems doubtful that even the existence of these cyclonic storms would have been proposed in the absence of direct observations and on the basis of purely theoretical investigations of idealized models of Earth's ocean-atmospheric structure. It is thus of interest to use terrestrial hurricanes as a guide to explore whether a similar situation can arise in protoplanetary disks. Our conclusion is if disk gas flows are undersaturated, particularly just outside the water ice line, then the essence of the hurricane mechanism may spring into operation.
Under-saturation incites sublimation of water from ice-coated dust-grains in the mid-plane, leading to an energy increase in the flow. Our heuristic two-dimensional simulations suggest that modest under-saturations are sufficient to maintain pre-existing vortices that would otherwise dissipate (Sect. 4.2). Our results also suggest that the shear inherent to gas and particle streams in protoplanetary disks can enhance the enthalpy transfer mechanism, and potentially cause formation of vortices (Sect. 4.3).
Three-dimensional simulations, ideally including the dust-particles independently, are necessary to confirm the feasibility of hurrican-like storms in prostellar disks. As with hurricanes on Earth, a carnot engine in a disk will be inherently three-dimensional. Our work here merely suggests that latent heat fluxes can strongly affect the flow dynamics near the water ice line. Further investigation thus seems warranted.