The vertical shear instability VSI is widely believed to be effective in driving turbulence in protoplanetary discs PPDs. Prior studies on VSI exclusively exploit the reflecting boundary conditions BCs at the disc surfaces. VSI depends critically on the boundary behaviours of waves at the disc surfaces. We extend earlier studies by performing a comprehensive numerical analysis of VSI with partially reflecting BCs for both the axisymmetric and nonaxisymmetric unstable VSI modes. We find that the growth rates of the unstable modes diminish when the outgoing component of the flow is greater than the incoming one for highorder body modes. When the outgoing wave component dominates the growth of VSI is notably suppressed. We find that the nonaxisymmetric modes are unstable and they grow at a rate that decreases with the azimuthal wavenumber. The different BCs at the lower and upper disc surfaces naturally lead to nonsymmetric modes relative to the disc midplane. The potential implications of our studies for further understanding planetary formation and evolution in PPDs are also briefly discussed.
['Astrophysics - Earth and Planetary Astrophysics']
Vertical shear instability with partially reflecting boundary conditions
{'Vertical Shear Instability with Partially Reflecting Boundary Conditions': 'Yuzi Wu ( 吴 玉 孜 ), 1 , 2 Cong Yu ( 余 聪 ), 1 , 2 ★ Can Cui ( 崔 灿 ) 3  \n1 School of Physics and Astronomy, Sun Yat-Sen University, Zhuhai 519082, China  \n- 2 CSST Science Center for the Guangdong-Hong Kong-Macau Greater Bay Area, Zhuhai 519082, China\n- 3 Department of Astronomy and Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'The vertical shear instability (VSI) is widely believed to be effective in driving turbulence in protoplanetary disks. Prior studies on VSI exclusively exploit the reflecting boundary conditions (BCs) at the disk surfaces. VSI depends critically on the boundary behaviors of waves at the disk surfaces. We extend earlier studies by performing a comprehensive numerical analysis of VSI with partially reflecting BCs for both the axisymmetric and non-axisymmetric unstable VSI modes. We find that the growth rates of the unstable modes diminish when the outgoing component of the flow is greater than the incoming one for high-order body modes. When the outgoing wave component dominates, the growth of VSI is notably suppressed. We find that the non-axisymmetric modes are unstable and they grow at a rate that decreases with the azimuthal wavenumber. The different BCs at the lower and upper disk surfaces naturally lead to non-symmetric modes relative to the disk midplane. The potential implications of our studies for further understanding planetary formation and evolution in protoplanetary disks (PPDs) are also briefly discussed.  \nKey words: hydrodynamics - instabilities - protoplanetary discs.', '1 INTRODUCTION': 'Protoplanetary disks serve as the nurseries for planet formation. Turbulence is a key factor driving angular momentum transport within the disks, and further exploration of the physical factors that drive turbulence would enhance our understanding on disk mass and angular momentum transfer, which is of great significance in the processes of planet formation and evolution (Youdin & Lithwick 2007; Lesur et al. 2023). With the development of observational technology, the finer structures of protoplanetary disks are gradually emerging. The high-resolution observations by the Atacama Large Millimeter/submillimeter Array (ALMA) reveal the distribution of dust and gas within protoplanetary disks, as well as clear substructures such as rings and gaps (Teague et al. 2016; Yamaguchi et al. 2024). Concurrently, ALMA has detected a certain level of turbulence within protoplanetary disks, suggesting that turbulence may play a significant role in various stages of planet formation (Manger & Klahr 2018; Latter & Papaloizou 2018; Flaherty et al. 2019, 2024). The investigation into the origin of disk turbulence will provide an indispensable theoretical framework for a more profound comprehension of the formation and evolution processes of exoplanets.  \nIn different regions of the protoplanetary disk, the dominant mechanisms for generating turbulence are diverse (Pfeil & Klahr 2019). In the inner disk region, with radii smaller than about 0.3 astronomical units and the outermost disk region, beyond 100 astronomical units, turbulence is primarily driven by magnetorotational instability (MRI). In the intermediate region, however, the weak ionization results in a poor coupling between the magnetic field and gas, suppressing the MRI (Bai 2011; Cui & Bai 2021). Within this "dead  \nzone", purely hydrodynamic instabilities become the primary driver of turbulence (Yellin-Bergovoy et al. 2021), such as the VSI (Urpin & Brandenburg 1998; Urpin 2003), zombie vortex instability (ZVI) (Marcus et al. 2015), and convective overstability (COS) (Klahr & Hubbard 2014). The operation of these instabilities depends on specific environmental conditions. VSI, due to its extensive influence in protoplanetary disks, has received widespread attention.  \nThe VSI arises in the presence of angular velocity shear 𝜕 Ω / 𝜕𝑧 within protoplanetary disks, serving as an effective source of turbulence. The VSI is similar in form to the Goldreich-Schubert-Fricke instability (GSF), which was initially discovered in slowly rotating stars (Goldreich & Schubert 1967). It does not require strong radial stratification and is induced solely by vertical shear (Urpin 2003).  \nFor protoplanetary disks, the dependence of Ω on 𝑧 significantly alters the stability characteristics of the disk. However, fluid elements also experience stabilizing effects due to vertical buoyancy, hindering the growth of VSI. Effective rapid cooling can overcome this stabilizing effect caused by buoyancy. Under rapid cooling with 𝑡 c < Ω -1 K ℎ | 𝑞 |/( 𝛾 -1 ) , VSI will be effective (Lin & Youdin 2015), where 𝑡 c is the cooling timescale, ℎ is the disk aspect ratio, 𝑞 is the power-law index of temperature profile and 𝛾 is the adiabatic index. For cooling times longer than the critical cooling time of VSI, COS modes will emerge, experiencing similar growth rates as GSF (Klahr 2024). Realistic protoplanetary disks have complex density and temperature structures, with thermal relaxation depending on the intricate spatial distribution of local disk stratification. Threedimensional simulations have studied the thermal coupling of gas and dust particles in optically thick and optically thin regions within protoplanetary disks. The decoupling of oscillations of dust and gas particles in the upper atmosphere leads to ineffective thermal relaxation, thereby suppressing VSI turbulence (Pfeil & Klahr 2021).  \nThe VSI is classified into surface modes and body modes (Nelson et al. 2013; Lin & Youdin 2015). Surface modes concentrate on the disk surface. Body modes extend vertically within the disk, exhibiting prominent vertical oscillations with radial wavenumbers much larger than vertical wavenumbers. Based on the symmetry of body modes, they can be further classified into breathing modes and corrugation modes (Nelson et al. 2013).  \nNelson et al. (2013) studied the nonlinear evolution of VSI through numerical simulations. They showed that, although surface modes grow rapidly, unstable modes are eventually dominated by body modes, especially corrugation modes, which are characterized by a narrow band of large-scale vertical motions (Nelson et al. 2013). This can lead to turbulence with an effective parameter 𝛼 ≈ 10 -4 ∼ 10 -3 , causing non-negligible angular momentum transport (Nelson et al. 2013; Stoll & Kley 2014).  \nThe VSI mainly occurs in the "dead zone" with weak magnetic coupling within protoplanetary disks. The magnetic field provides a stabilizing effect through magnetic tension or magnetorotational turbulence, suppressing the VSI. With increasing magnetization, surface modes will disappear before body modes (Cui & Lin 2021; Latter & Kunz 2022). The Ohmic resistivity in non-ideal MHD can overcome magnetic stability, restoring the instability (Cui & Lin 2021). Disk winds could coexist and interact with the VSI, enhancing the accretion rate of disk winds and triggering strong VSI turbulence. Weak ambipolar diffusion strength or enhanced coupling between gas and magnetic fields will also suppress shear instability (Cui & Bai 2020).  \nPrevious studies on the VSI adopted a vertically global treatment with the reflecting BCs (e.g. Nelson et al. 2013; Lin & Youdin 2015; Cui & Lin 2021; Svanberg et al. 2022). The effects of the disk\'s vertical extent on the VSI were studied by Barker & Latter (2015). Lin & Youdin (2015) applied free and rigid BCs to the VSI. Free BCs establish a zero Lagrangian pressure perturbation at the boundary, resembling a reflecting BC with a phase difference, as we detailed in Section A. The rigid BC means no vertical velocity perturbation at the boundary, consistent with the one used in Barker & Latter (2015). It is also a reflecting BC, albeit with half-wave loss, as we explained in Section 4.3. This work investigates the VSI with partially reflecting BCs, introducing a new factor worth further exploration to enhance our understanding of the VSI\'s behavior in more realistic astrophysical scenarios.  \nThe complex environment of protoplanetary disks at their surfaces can indeed give rise to the outgoing flows. Zhang et al. (2024) considered that full disks and transition disks with small cavities have a superheated atmosphere and cool midplane, with this temperature structure producing an outgoing flow layer at optical depth 𝜏 ∗ < 1 on top of an ingoing flow layer at 𝜏 ∗ ∼ 1. Notably, magnetohydrodynamic (MHD) disk winds are one such phenomenon that is effectively launched and driven from the disk surface, providing a physical foundation for understanding and justifying the existence of outgoing BCs. In addition, the presence of anti-symmetric magnetic field modes in protoplanetary disks results in one-sided characteristics for disk winds, meaning that these winds occur solely on one side of the disk relative to its midplane, which can potentially lead to asymmetry in the outgoing flow (Wang et al. 2024).  \nConsequently, it is essential to fill the research gap regarding the behavior of VSI under conditions of partial reflection at the boundaries. We focus on the effect of partially reflecting BCs, considering the axisymmetric case and further including a nonzero azimuthal wavenumber to account for non-axisymmetric VSI effects.  \nThe structure of the paper is as follows: In § 2, we present the dynamical equations. In § 3, we introduce the vertical shear of the equilibrium disk model. In § 4, we show the linearized equations for  \nboth axisymmetric and non-axisymmetric perturbations and describe the partially reflecting BCs. In § 5, we analyze the growth rates and symmetry of VSI. Finally, we summarize our results in § 6.', '2 BASIC EQUATIONS': 'The continuity, momentum and energy conservation equations of hydrodynamics for a protoplanetary disk are as follows:  \n𝜕𝜌 𝜕𝑡 + ∇ · ( 𝜌 𝒖 ) = 0 , (1)  \n𝜕 𝒖 𝜕𝑡 + ( 𝒖 · ∇) 𝒖 = -1 𝜌 ∇ 𝑃 - ∇ Φ , (2)  \n𝜕𝑃 𝜕𝑡 + ( 𝒖 · ∇) 𝑃 = -𝛾𝑃 ∇ · 𝒖 -Λ , (3)  \nwhere 𝜌, 𝒖 , 𝑃 are the density, three-dimensional velocity and pressure, respectively. To study the protoplanetary disk, we adopt the cylindrical coordinates( 𝑟, 𝜙, 𝑧 ). The gravitational potential by the central star with mass 𝑀 is Φ = -𝐺𝑀 /( 𝑟 2 + 𝑧 2 ) 1 / 2 , where 𝐺 is the gravitational constant. In equation (3), 𝛾 is the adiabatic index, and the sink term Λ on the right side represents the radiative cooling.', '3 EQUILIBRIUM DISK MODEL': 'The radial power-law distribution for both density and temperature in the cylindrical coordinate are:  \n𝜌 ( 𝑟 ) = 𝜌 0 GLYPH<18> 𝑟 𝑟 0 GLYPH<19> 𝑝 , 𝑇 ( 𝑟 ) = 𝑇 0 GLYPH<18> 𝑟 𝑟 0 GLYPH<19> 𝑞 , (4)  \nwhere 𝑝 is the power-law index of the mass density and 𝑞 is the power-low index of the temperature, 𝜌 0 and 𝑇 0 are the density and temperature at the fiducial radius 𝑟 0 , respectively. We adopt the ideal equation of state, 𝑃 = R 𝜇 𝜌𝑇 . Here R is the gas constant and 𝜇 is the mean molecular mass of gas. For a characteristic timescale 𝑡 c, the radiative cooling takes the form:  \nΛ = 𝜌 R 𝜇 𝑇 -𝑇 eq 𝑡 c = 𝑃 -𝑐 2 s 𝜌 eq 𝑡 c , (5)  \nwhere 𝑇 eq and 𝜌 eq are the equilibrium temperature and density, respectively. In the absence of cooling, the disk will be stable against the Solberg-Hoiland criteria. Sufficiently rapid cooling is an indispensable factor to trigger VSI since it requires rapid thermal relaxation to inhibit buoyancy.  \nThe equilibrium state of the disk at radial and vertical directions meet:  \n𝜕 Φ 𝜕𝑟 + 1 𝜌 𝜕𝑃 𝜕𝑟 = 𝑟 Ω 2 , 𝜕 Φ 𝜕𝑧 + 1 𝜌 𝜕𝑃 𝜕𝑧 = 0 , (6)  \nwhere Ω = Ω ( 𝑟, 𝑧 ) is the equilibrium rotation field. The density and pressure are related by the isothermal sound speed 𝑐 s, which is 𝑐 s ( 𝑟 ) = √︁ 𝑃 / 𝜌 . The density profile is:  \n𝜌 ( 𝑟, 𝑧 ) = 𝜌 ( 𝑟 ) exp " 1 𝑐 s 2 𝐺𝑀 √︁ 𝑟 2 + 𝑧 2 -𝐺𝑀 𝑟 !# , (7)  \nThe equilibrium angular velocity is  \nΩ 2 ( 𝑟, 𝑧 ) = Ω 2 K ( 𝑟 ) " 1 + 𝑞 + ( 𝑝 + 𝑞 ) ℎ 2 -𝑞 𝑟 √︁ 𝑟 2 + 𝑧 2 # , (8)  \nwhere Ω K = √︁ 𝐺𝑀 / 𝑟 3 is the Kepler frequency. We assume thin disk models with a small disc aspect ratio ℎ , i.e. ℎ = 𝐻 / 𝑟 ≪ 1, where 𝐻 = 𝑐 𝑠 / Ω 𝐾 is the scale height of the disc. In this case, the density profile in equation (7) is simplified as:  \n𝜌 ( 𝑟, 𝑧 ) = 𝜌 ( 𝑟 ) exp GLYPH<18> -𝑧 2 2 𝐻 2 GLYPH<19> . (9)  \nRadial variations in temperature inevitably lead to a vertical gradient in the angular velocity, or vertical shear:  \n𝜕 Ω 2 𝜕𝑧 = 𝑞𝑧 𝑟 2 Ω 2 K GLYPH<0> 1 + 𝑧 2 / 𝑟 2 GLYPH<1> 3 / 2 . (10)  \nThe radial gradient of the angular velocity squared is:  \n𝜕 Ω 2 𝜕𝑟 = -3 Ω 2 𝑟 -𝑞𝑧 2 𝑟 3 Ω 2 K GLYPH<0> 1 + 𝑧 2 / 𝑟 2 GLYPH<1> 3 / 2 . (11)', '4 LINEAR PERTURBATION EQUATIONS': 'In this work, the linear perturbations for both the axisymmetric and non-axisymmetric cases with partially reflecting BCs will be discussed. The Eulerian perturbations ( 𝛿𝑢 𝑟 , 𝛿𝑢 𝜙 , 𝛿𝑢 𝑧 , 𝛿𝜌, 𝛿𝑃 ) have space and time dependence of the form exp ( 𝑖𝑘 𝑟 𝑟 + 𝑖𝑚𝜙 -𝑖𝜐𝑡 ) , where 𝑘 𝑟 is the radial wavenumber, 𝑚 is the azimuthal wavenumber and 𝜐 is the complex frequency. The linearized perturbation equations for the velocity, density and pressure perturbations are:  \n𝑖 ( 𝜐 -𝑚 Ω ) 𝛿𝜌 𝜌 = 𝑖𝑘 𝑟 𝛿𝑢 𝑟 + 𝛿𝑢 𝑟 𝑟 + 𝑖𝑚 𝑟 𝛿𝑢 𝜙 + 𝜕𝛿𝑢 𝑧 𝜕𝑧 + 𝜕 ln 𝜌 𝜕𝑧 𝛿𝑢 𝑧 , (12)  \n𝑖 ( 𝜐 -𝑚 Ω ) 𝛿𝑢 𝑟 = 𝑖𝑘 𝑟 𝛿𝑃 𝜌 -2 Ω 𝛿𝑢 𝜙 , (13)  \n𝑖 ( 𝜐 -𝑚 Ω ) 𝛿𝑢 𝜙 = 𝑖𝑚 𝑟 𝛿𝑃 𝜌 + 𝑟 𝜕 Ω 𝜕𝑧 𝛿𝑢 𝑧 + 𝜅 2 2 Ω 𝛿𝑢 𝑟 , (14)  \n𝑖 ( 𝜐 -𝑚 Ω ) 𝛿𝑢 𝑧 = 𝑑 𝑑𝑧 GLYPH<18> 𝛿𝑃 𝜌 GLYPH<19> + 𝜕 ln 𝜌 𝜕𝑧 GLYPH<18> 𝛿𝑃 𝜌 -𝑐 s 2 𝛿𝜌 𝜌 GLYPH<19> , (15)  \n𝑖 ( 𝜐 -𝑚 Ω ) 𝛿𝑃 𝜌 = 𝑐 s 2 𝛾 GLYPH<18> 𝑖𝑘 𝑟 𝛿𝑢 𝑟 + 𝛿𝑢 𝑟 𝑟 + 𝑖𝑚 𝑟 𝛿𝑢 𝜙 + 𝜕𝛿𝑢 𝑧 𝜕𝑧 GLYPH<19> + 𝑐 s 2 𝜕 ln 𝜌 𝜕𝑧 𝛿𝑢 𝑧 + 1 𝑡 c GLYPH<18> 𝛿𝑃 𝜌 -𝑐 s 2 𝛿𝜌 𝜌 GLYPH<19> , (16)  \nwhere 𝜅 is the epicyclic frequency and 𝜅 2 ≡ 4 Ω 2 + 𝑟 𝜕 Ω 2 𝜕𝑟 . Note that our equations retain the curvature terms neglected in Lin & Youdin (2015). The impact caused by the curvature terms is weak. We have ignored the radial gradients of the equilibrium density and pressure to simplify the equation. This type of radially local approximation was initially proposed in McNally & Pessah (2014) and has also been adopted in Lin & Youdin (2015). However, we retain the radial temperature gradient to ensure the VSI operates.', '4.1 Wave Equations for the Axisymmetric VSI': 'For the axisymmetric VSI, we take 𝑚 = 0. In terms of the variables 𝑦 1 ≡ 𝛿𝑢 𝑧 , and 𝑦 2 ≡ 𝛿𝑃 / 𝜌 , Equations (12)-(16) can be cast into two  \nfirst-order ordinary differential equations (ODEs) 1 :  \n𝑑𝑦 1 𝑑 ˆ 𝑧 = 𝐴 11 𝑦 1 + 𝐴 12 𝑦 2 , 𝑑𝑦 2 𝑑 ˆ 𝑧 = 𝐴 21 𝑦 1 + 𝐴 22 𝑦 2 . (17)  \nThe relevant coefficients in the above linear ODEs are  \n𝐴 11 = ( 𝜒 + 𝑄 ) ˆ 𝑧 , 𝐴 12 = 𝑖 ˆ 𝜐 ( 𝜒 + 𝑈 ) , (18)  \n𝐴 21 = 𝑖 ˆ 𝜐 -ˆ 𝑧 2 ( 𝜒 -1 ) 𝑖 ˆ 𝜐 , 𝐴 22 = - ( 𝜒 -1 ) ˆ 𝑧 , (19)  \nwhere 𝜒 = ( 1 -𝑖 ˆ 𝜐𝛽 )/( 1 -𝑖 ˆ 𝜐𝛽𝛾 ) , 𝑄 = 𝑖ℎ𝑞 ˆ 𝑘 𝑟 + 𝑞ℎ 2 , and 𝑈 = ˆ 𝑘 2 𝑟 -𝑖 ˆ 𝑘 𝑟 ℎ . Note that the equations above have been written in a dimensionless form adopting the following rules: 𝑧 = ˆ 𝑧𝐻 , 𝑘 𝑟 = ˆ 𝑘 𝑟 / 𝐻 , 𝜐 = ˆ 𝜐 Ω K , 𝛽 = 𝑡 c Ω K . The eigenvalue ˆ 𝜐 is a complex number. Its real part ˆ 𝜔 is the frequency and the imaginary part ˆ 𝜎 means the growth rate. With our convention, the sound speed 𝑐 𝑠 = 𝐻 Ω 𝐾 is normalized to unity.  \nWe are interested in the frequency range in which 𝜐 is small compared to the epicyclic frequency 𝜅 . This low-frequency approximation filters out the acoustic waves and retains the inertial-gravity waves (Lubow & Pringle 1993). With this approximation and rapid cooling 𝛽 = 10 -3 , the dispersion relation (B1) in Appendix B can be written as  \nˆ 𝜐 2 ≃ ˆ 𝑘 2 𝑧 + 𝑖 GLYPH<16> 𝑖ℎ𝑞 ˆ 𝑘 𝑟 + 𝑞ℎ 2 + 1 GLYPH<17> ˆ 𝑧 ˆ 𝑘 𝑧 1 + ˆ 𝑘 2 𝑟 -𝑖 ˆ 𝑘 𝑟 ℎ , (20)  \nSince ℎ ≪ 1 and ˆ 𝑘 𝑟 ≫ 1 , ˆ 𝑘 𝑧 , the equation shows the dispersion relation for inertial waves (Goodman & Lackner 2009)  \n𝜐 2 ≈ Ω 2 𝐾 ˆ 𝑘 2 𝑧 + 𝑖 GLYPH<16> 𝑖ℎ𝑞 ˆ 𝑘 𝑟 + 1 GLYPH<17> ˆ 𝑧 ˆ 𝑘 𝑧 ˆ 𝑘 2 𝑟 ∼ Ω 2 𝐾 ˆ 𝑘 2 𝑧 ˆ 𝑘 2 𝑟 . (21)  \nIt is worthwhile to note that the VSI actually corresponds to the inertial wave. The vertical component of the group velocity and phase velocity for inertial waves are always in the same direction.', '4.2 Wave Equations for the Non-axisymmetric VSI': "In the non-axisymmetric case, the azimuthal wavenumber 𝑚 ≠ 0. Westill adopt the low-frequency approximation to simplify our analysis, i.e., the Doppler shifted frequency is much smaller than the epicyclic frequency, | 𝜐 -𝑚 Ω | 2 ≪ 𝜅 2 . The governing ODEs for the non-axisymmetric case read:  \n𝑑𝑦 1 𝑑 ˆ 𝑧 = 𝐴 ' 11 𝑦 1 + 𝐴 ' 12 𝑦 2 , 𝑑𝑦 2 𝑑 ˆ 𝑧 = 𝐴 ' 21 𝑦 1 + 𝐴 ' 22 𝑦 2 , (22)  \nwhere  \n𝐴 ' 11 = GLYPH<0> 𝜒 ' + 𝑄 ' GLYPH<1> ˆ 𝑧 , (23)  \n𝐴 ' 12 = 𝑖 ( ˆ 𝜐 -𝑚 ) GLYPH<0> 𝜒 ' + 𝑈 ' GLYPH<1> + 𝑊 2 -ˆ 𝑘 𝑟 𝑚𝑞ℎ 3 2 ˆ 𝑧 2 , (24)  \n1 In a more compact matrix form,  \n𝑑 𝒀 𝑑 ˆ 𝑧 = 𝑨𝒀 , where 𝑨 = GLYPH<18> 𝐴 11 𝐴 12 𝐴 21 𝐴 22 GLYPH<19> , 𝒀 ≡ ( 𝑦 1 , 𝑦 2 ) 𝑇 = GLYPH<18> 𝑦 1 𝑦 2 GLYPH<19> .  \nIt should be stressed that all the wave propagation properties are hidden in the matrix 𝑨 . The two eigenvalues ( 𝑖 ˆ 𝑘 𝑧 1, 𝑖 ˆ 𝑘 𝑧 2) and two eigenvectors ( 𝒓 1, 𝒓 2) of the matrix 𝑨 are important for our analysis of BCs. (Detailed derivation of the eigenvalues and eigenvectors are given in Appendix A. Further explanation can be found in Appendix B.)  \n𝐴 ' 21 = 𝑖 ( ˆ 𝜐 -𝑚 ) -ˆ 𝑧 2 ( 𝜒 ' -1 ) 𝑖 ( ˆ 𝜐 -𝑚 ) , 𝐴 ' 22 = -GLYPH<0> 𝜒 ' -1 GLYPH<1> ˆ 𝑧 . (25)  \nRelevant variables in the above matrix coefficients are  \n𝜒 ' = [ 1 -𝑖 ( ˆ 𝜐 -𝑚 ) 𝛽 ] /[ 1 -𝑖 ( ˆ 𝜐 -𝑚 ) 𝛽𝛾 ] , 𝑈 ' = 𝑈 + 𝑚 2 ℎ 2 ,  \n𝑄 ' = 𝑄 + ( ˆ 𝜐 -𝑚 ) 𝑚𝑞ℎ 2 / 2 , 𝑊 = 4 𝑖𝑚ℎ 2 -3 𝑚ℎ ˆ 𝑘 𝑟 ,  \nrespectively. Compared to equation (17), the additional terms in the above equations are due to the non-axisymmetric effect.", '4.3 Boundary Conditions': "We adopt the WKB approximation to specify the BCs. Note that it is a local approximation of the linear second-order ODEs. The WKBapproximation is only exploited locally at the boundary, which provides certain constraints on the global integration. Our global solutions are obtained by numerical integration rather than the WKB approximation. According to Appendix A, the WKB approximation can be understood in terms of the eigenstructures of the coefficient matrix of the wave equations. The variables 𝒀 ≡ ( 𝑦 1 , 𝑦 2 ) 𝑇 can be viewed as the sum of two eigenvectors 𝒓 1 and 𝒓 2 :  \n𝒀 = 𝑎 1 𝒓 1 + 𝑎 2 𝒓 2 . (26)  \nAs mentioned in Appendix B 2 , we know that the wave group velocity associated with 𝒓 1 is incoming and 𝒓 2 outgoing. The coefficients 𝑎 1 and 𝑎 2 represent the strength of waves in these two directions, respectively. The boundary parameter 𝑅 ≡ 𝑎 2 / 𝑎 1 is determined by 𝑦 1 and 𝑦 2 in the following way (Please refer to Appendix A for the derivation):  \n𝑅 = h -𝐴 1 ˆ 𝑧 + 𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 i 𝑦 1 -𝐵 1 𝑦 2 h 𝐴 1 ˆ 𝑧 + 𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 i 𝑦 1 + 𝐵 1 𝑦 2 . (27)  \nThe parameters, 𝐴 1 , 𝐴 2 , 𝐴 3 and 𝐵 1 in equation (27) are explicitly given in Appendix A. Note that equation (27) can be written alternatively as a linear combination of 𝑦 1 and 𝑦 2 :  \nGLYPH<18> 𝐴 1 ˆ 𝑧 + 𝑅 -1 𝑅 + 1 𝑖 √︃ 𝐴 3 -𝐴 2 ˆ 𝑧 2 GLYPH<19> 𝑦 1 + 𝐵 1 𝑦 2 = 0 . (28)  \nThe above equation is the actual BC adopted in our numerical calculations. Perfect reflection occurs when 𝑅 = 1,  \n𝐴 1 ˆ 𝑧𝑦 1 + 𝐵 1 𝑦 2 = 0 . (29)  \nFor a reflection with a half-wave loss, the BC can be expressed as 𝑦 1 = 0 with 𝑅 = -1, i.e. 𝛿𝑢 𝑧 = 0 at ˆ 𝑧 = ± ˆ 𝑧 max. As 𝑅 approaches infinity, the mode solely exhibits an outgoing inertial wave. Purely outgoing boundary conditions are widely used in the studies of Rossby wave instability (Li et al. 2000; Yu & Lai 2013; Huang & Yu 2022).  \nTheBCforthenon-axisymmetric case could be similarly obtained, which is:  \nGLYPH<18> 𝐴 ' 1 ˆ 𝑧 + 𝑅 -1 𝑅 + 1 𝑖 √︃ 𝐴 3 ' + 𝐴 ' 4 ˆ 𝑧 4 -𝐴 2 ' ˆ 𝑧 2 GLYPH<19> 𝑦 1 + GLYPH<16> 𝐵 1 ' -𝐵 2 ' ˆ 𝑧 2 GLYPH<17> 𝑦 2 = 0 , (30)  \nwhere  \n𝐵 ' 1 = 2 𝑖 ( ˆ 𝜐 -𝑚 ) GLYPH<0> 𝜒 ' + 𝑈 ' GLYPH<1> + 𝑊 , 𝐵 ' 2 = ˆ 𝑘 𝑟 𝑚𝑞ℎ 3 , (31)  \n𝐴 ' 1 = 2 𝜒 ' + 𝑄 ' -1 , (32)  \n𝐴 ' 2 = GLYPH<0> 𝑄 ' + 1 GLYPH<1> 2 -2 𝑖 ( ˆ 𝜐 -𝑚 ) 𝐵 ' 2 + 4 GLYPH<0> 𝜒 ' -1 GLYPH<1> GLYPH<0> 𝑄 ' -𝑉 GLYPH<1> , (33)  \n𝐴 ' 3 = 4 ( ˆ 𝜐 -𝑚 ) 2 GLYPH<0> 𝜒 ' + 𝑉 GLYPH<1> , 𝐴 ' 4 = 2 𝑖 ( 𝜒 ' -1 ) ˆ 𝜐 -𝑚 𝐵 ' 2 , (34)  \nrespectively. The parameter 𝑉 is given by 𝑉 = 𝑈 ' -𝑖𝑊 /[ 2 ( ˆ 𝜐 -𝑚 )] . Following Lubow & Pringle (1993), it should be noted that when applying the low-frequency approximation to non-axisymmetric analyses, a critical requirement exists that the modes must be tightly wrapped, leading to a condition 𝑚 ≪ ˆ 𝑘 𝑟 / ℎ .", '5 NUMERICAL RESULT': 'To investigate the axisymmetric or non-axisymmetric linear perturbations of VSI, we need to solve Equation (17) or Equation (22). The eigenvalues of these equations are closely related to the BCs, which we adopt as Equation (28) or Equation (30). These BCs represent a linear combination of outgoing and incoming wave components of different strengths. The numerical solution of this kind of complex two-point boundary eigenvalue problem can be efficiently obtained using the relaxation method (Press et al. 1992). In this method, we discretize the ODEs into finite difference equations (FDEs) across a mesh that spans the domain of interest. Through iterative trial and error, we progressively determine eigenvalues and eigenfunctions by a method analogous to the Newton -Raphson method. Throughout these iterations, both the ODEs and BCs are satisfied simultaneously. In this study, we employ 5001 uniform mesh points to achieve a high-accuracy numerical solution with an average error of less than 10 -8 . The basic parameters of the disk are set as 𝑞 = -1, 𝑝 = -1 . 5, ℎ = 0 . 05, 𝛾 = 1 . 4, respectively.', '5.1 Axisymmetric VSI unstable modes': 'Usually, VSI modes exhibit odd or even symmetry. Our calculations shown in Fig. 2-3 with symmetric BCs at 𝑧 = -𝑧 max and 𝑧 = 𝑧 max obtain similar results to Nelson et al. (2013). However, there exist highly asymmetric cases, as demonstrated in the bottom right panel of Figure 20 in Nelson et al. (2013). There is still a lack of a clear explanation for the asymmetry. For simplicity, we choose the domain 𝑧 ∈ [ 0 , 𝑧 max ] to avoid such asymmetric surface modes. At the midplane 𝑧 = 0, we set the BCs according to the parity of the eigenfunction. While at 𝑧 = 𝑧 max, we set partially reflecting BCs.  \nOur numerical results with the cooling parameter 𝛽 = 10 -3 are shown in Fig. 1, which depicts the variation of growth rate with the frequency for partially reflecting BCs. Each symbol denotes a specific order mode. For a general set of parameters, lower-order modes possess a smaller eigenfrequency | 𝜔 | . In the upper panel, we present the result for smaller radial wavenumber ˆ 𝑘 𝑟 = 10. Only the body modes are observed in this case. The results for larger radial wavenumber ˆ 𝑘 𝑟 = 30 are displayed in the lower panel. In this case, both the body and surface modes show up. The growth rate varies with the frequency in a non-monotonic way. With the gradual increase of frequency, the growth rate first increases with the frequency (i.e., lower-order modes). Once the growth rate reaches a maximum, it would decrease with the frequency (i.e., higher-order modes). The growth rate of the low order body mode is almost not affected by  \nFigure 1. Variations of the growth rate with frequency for the radial wavenumber ˆ 𝑘 𝑟 = 10 and ˆ 𝑘 𝑟 = 30, respectively. The cooling parameter is 𝛽 = 10 -3 . The disk respect ratio ℎ = 0 . 05 and maximum disk height 𝑧 max = 5 𝐻 . The boundary condition parameter is set to 𝑅 = 2 (magenta pluses), 𝑅 = 0 . 5 (blue pluses) and 𝑅 = 1 (black pluses), respectively.  \n<!-- image -->  \nthe BCs. However, for the higher-order body modes, the influence of the BCs becomes more obvious. However, for surface modes, the boundary has a relatively minor effect on the growth rate. As the parameter 𝑅 decreases, the growth rates become larger. On the other hand, as the parameter 𝑅 increases, the growth rates decrease. Our result indicate that the surface of the disk actually acts as a trap of the inertial waves. When the trapping effects become weaker ( 𝑅 becomes larger), the VSI would become weaker as well.  \nFig. 2 shows the vertical velocity perturbation, 𝛿𝑢 𝑧 , of the loworder corrugation mode (symmetric with respect to midplane) with growth rate 𝜎 = 0 . 73 ℎ Ω K . In the upper panel, perfectly reflecting BC is set at both the upper and lower boundaries. The resulting eigenfunction is symmetric relative to the midplane. In the lower panel, the perfectly reflecting BC ( 𝑅 = 1) is adopted at the lower boundary 𝑧 = -5 𝐻 and a partially reflecting BC ( 𝑅 = 2) at the upper boundary 𝑧 = 5 𝐻 . The strength of the outgoing wave is twice that of the incoming wave, breaking the symmetry with respect to the midplane. At the upper boundary, where the flow is predominantly outgoing, the vertical velocity perturbations are stronger than those observed at the lower boundary. Owing to its nature as a low-order mode, the growth rate is barely affected by BCs. For higher order modes, such as the breathing mode (anti-symmetric with respect to the midplane) with a frequency | 𝜔 | = 2 . 15 ℎ Ω 𝐾 shown in Fig. 3, the growth rate is clearly affected by the BCs. For the partially reflecting BC with 𝑅 = 2, the growth rate has been suppressed by 0 . 02 ℎ Ω 𝐾 .  \nFigure 2. Vertical velocity perturbation 𝛿𝑢 𝑧 of the low order corrugation mode. Top: with same boundary parameter 𝑅 = 1 at 𝑧 = ± 𝑧 max. Bottom: with boundary parameter 𝑅 = 1 on the lower boundary and boundary parameter 𝑅 = 2 on the upper boundary, respectively.  \n<!-- image -->  \nFigure 3. Same as Figure 2 but for higher order breathing mode.  \n<!-- image -->  \nMeanwhile, the anti-symmetric configuration originally relative to the midplane can no longer be maintained.  \nWecompare three cases of 𝑅 = 0 . 01, 𝑅 = 1, and 𝑅 = 100 in Fig. 4. When the strength of the outgoing wave is much greater than that of the incoming wave, the growth rate is significantly suppressed. Especially, the growth rate of the higher order modes with frequency | 𝜔 | ≳ 5 ℎ Ω K are preferentially damped. An increase in the strength of the outgoing component also results in a shift of the wave frequency  \nFigure 4. Growth rate and intrinsic frequency of the unstable modes for ˆ 𝑘 𝑟 = 30andextremeboundaryparameters 𝑅 = 0 . 01(yellow pluses), 𝑅 = 100 (green pluses) and the case of 𝑅 = 1 (black pluses) as a reference.  \n<!-- image -->  \nFigure 5. Growth rate of unstable modes vs. vertical extent. Blue and yellow pluses symbolize boundary parameters 𝑅 = 1 and 𝑅 = 2 for 𝑧 max = 5 𝐻 , while black and magenta pluses represent the same parameters for 𝑧 max = 8 𝐻 .  \n<!-- image -->  \nwhere the maximum growth rate is achieved towards the lower frequency. We find that the VSI is quite similar to the Papaloizou-Pringle instability (PPI), which requires a reflecting boundary for PPI to operate (Papaloizou & Pringle 1984). As a comparison, there exists more robust instability, such as Rossby wave instability (Lovelace et al. 1999; Huang & Yu 2022). The reflecting boundary is not necessary for RWI. The fully outgoing boundary conditions (i.e., 𝑅 = ∞ ) are often adopted in the RWI calculations.  \nIn Fig. 5, we note that the growth rate varies non-monotonically with the frequency. Below a certain threshold frequency, the growth rate increases with the frequency. Above the threshold, the growth rate decreases with the frequency. The disk vertical extent has an important effect on the threshold frequency and corresponding growth rate. This phenomenon applies not only to reflecting BCs but also to partially reflecting BCs, as shown in Fig. 5. The threshold frequency for 𝑧 max = 8 𝐻 is larger than the threshold for 𝑧 max = 5 𝐻 . A larger vertical extent has a higher growth rate at the threshold frequency as well. Comparing the changes in growth rates due to contrasting BCs across the two vertical extents, the partially reflecting BCs have less effect on the growth rate in the larger one. Even so, as the ver-  \nFigure 6. Growth rate and wave frequency of the unstable modes with ˆ 𝑘 𝑟 = 30, disk respect ratio ℎ = 0 . 05, and thermal relaxation timescale 𝛽 = 0 . 001 (black), 𝛽 = 0 . 01 (magenta), 𝛽 = 0 . 02 (blue), 𝛽 = 0 . 04 (yellow) for different boundary parameters 𝑅 = 2 (colored pluses) and 𝑅 = 1 (colored dots).  \n<!-- image -->  \nFigure 7. Growth rate and frequency of the unstable modes with radial wavenumber ˆ 𝑘 𝑟 = 10, disk aspect ratio ℎ = 0 . 05, cooling time 𝛽 = 10 -3 . We set perfectly reflecting BCs 𝑅 = -1 and azimuthal wavenumbers 𝑚 = 0 (black pluses), 𝑚 = 1 (magenta pluses), 𝑚 = 2 (blue pluses), 𝑚 = 3 (cyan pluses).  \n<!-- image -->  \nincreases, a very high 𝑅 still significantly suppresses the unstable modes, similar to the findings in Fig. 4.  \nSince the cooling is essential for VSI to work, it would be interesting to study how the cooling affects VSI with our new BCs. We compare in Fig. 6 the mode growth rate and frequency for different values of cooling parameter 𝛽 with 𝑅 = 2 (partially reflecting) and 𝑅 = 1 (perfectly reflecting). It shows that, for either partially reflecting or perfectly reflecting BCs, the growth rates increase as 𝛽 varies from 0.04 to 0.001. At sluggish thermal relaxation, the growth of unstable modes is notably reduced. In this case, only the lower order unstable modes can survive. Furthermore, even with slower cooling, VSI is still affected by BCs. The augmentation of outgoing wave leads to a reduction in the growth rate of higher order modes, consistent with the scenario observed when cooling is set at 0.001 (as illustrated in Fig. 1).  \n̂  \nFigure 8. unstable modes with different BCs and same azimuthal wavenumber (top), same BC but different azimuthal wavenumbers (bottom).  \n<!-- image -->', '5.2 Non-axisymmetric VSI unstable modes': 'In Fig. 7, we show our results of non-axisymmetric VSI. The results for ˆ 𝑘 𝑟 = 10 with both 𝑚 ≠ 0 (non-axisymmetric) and 𝑚 = 0 (axisymmetric) are displayed. The parameter 𝑅 = -1 means perfect reflection with a half-wave loss. With the increase of 𝑚 , there is an overall decrease in the growth rate. Note that the azimuthal wavenumber 𝑚 ≪ ˆ 𝑘 𝑟 / ℎ , thus the low-frequency approximation could be used to deal with the non-axisymmetric perturbations.  \nFig. 8 shows non-axisymmetric VSI with larger radial wavenumber for both the body and surface modes. The top panel shows the effect of BCs on non-axisymmetric unstable modes with 𝑚 = 2. Consistent with the trend in the axisymmetric scenario, boundaries characterized by a stronger outgoing component suppress the growth rate of nonaxisymmetric instability as well. The bottom panel elucidates the variation in the growth rate caused by non-axisymmetric effects. In the lower panel, we fix the parameter 𝑅 = 2 at the upper boundary 𝑧 = 5 𝐻 . For instance, at the frequency of | 𝜔 | = 3 . 22 ℎ Ω 𝐾 , the growth rate decreases from 0 . 381 ℎ Ω 𝐾 to 0 . 307 ℎ Ω 𝐾 as the wavenumber 𝑚 increases from 0 to 3.  \nComparing the two panels of Fig. 8, the evident changes of unstable modes caused by BCs occur primarily for higher-order body modes, whereas non-axisymmetric effects demonstrate a global impact across the entire frequency range for both surface and body modes. The non-axisymmetric effect is significant and can not be neglected. Though we do not yet have a satisfactory physical explanation for the 𝑚 -dependence of VSI, our calculations indicate the term -3 𝑚ℎ ˆ 𝑘 𝑟 in the variable 𝑊 = 4 𝑖𝑚ℎ 2 -3 𝑚ℎ ˆ 𝑘 𝑟 plays an essential role for this dependence.', '6 CONCLUSION': 'In this paper, we analyze the VSI with partially reflecting BCs for vertically global and radially local protoplanetary disks. The partially reflecting BCs are viewed as a linear combination of the incoming and outgoing wave components. By adjusting the relative magnitudes of the two wave components, we can investigate the response of the VSI to different BCs.  \nWe find that the growth rate of VSI is closely related to the BCs, especially for higher-order breathing and corrugation modes whose frequency exceeds the threshold. The growth rate decreases when the boundary parameter, 𝑅 , increases. In other words, the growth rate of the instability diminishes when the outgoing wave component dominates. As the outgoing component increases, the location of the maximum growth rate shifts toward a lower frequency.  \nHow waves propagate away from the surface of realistic protoplanetary disks is determined by the atmosphere of the disk (Ogilvie & Lubow 1999). The fraction of the reflected and transmitted wave depends critically on the transition region between the disk and the atmosphere. In this paper, we use a free boundary parameter, 𝑅 , to describe the wave behavior at this disk-atmosphere transition region. Whether the disk atmosphere enhances or suppresses the VSI depends on the characteristics of this transition region. Our current understanding of the physical conditions of protoplanetary disks and their atmosphere is still limited. It would be interesting to perform further studies on wave propagation in a more realistic disk-atmosphere transition region.  \nThe vertical velocity perturbation exhibits symmetric or antisymmetric properties in the case of reflective BC at both the upper and lower boundaries (Nelson et al. 2013). It is possible that the upper and lower disk-atmosphere transition regions have different physical properties. In this case, the waves would behave differently at the upper and lower boundaries, and the inherent symmetry (or anti-symmetry) would be broken.  \nWe examine the growth rate of VSI with different vertical extents. As the vertical extent increases, the maximum growth rate of VSI increases as well. The finding is consistent with the discoveries made by Barker & Latter (2015). Under partially reflecting BCs, this property remains unchanged. In Lin & Youdin (2015), the free BCs establish a zero Lagrangian pressure perturbation, consistent with a reflecting boundary with a phase difference. Moreover, the rigid BC ( 𝛿𝑢 𝑧 = 0) they used aligns with our calculations with 𝑅 = -1. We have also performed calculations with boundary parameters 𝑅 = 1 (reflecting BC) and 𝑅 = -1. In both of these cases, the growth rates of VSI are found to be consistent.  \n. We find that the non-axisymmetric modes are also unstable and they grow at a rate that decreases with the azimuthal wavenumber. In a non-axisymmetric scenario, the BCs favoring the outgoing wave component still have a suppression effect on the VSI.  \nMHDdisk winds can potentially initiate outgoing flow in realistic PPDs. The presence of anti-symmetric magnetic field modes leads to one-sided characteristics in disk winds, which can further cause asymmetry in the outgoing flow of the disk. These factors will affect the development of VSI and may significantly influence planet formation and evolution. Realistic protoplanetary disks are magnetized. It would be interesting to take into account the effects of magnetic fields, incorporating our new treatment of BCs to explore the VSI in magnetized protoplanetary disks.', 'ACKNOWLEDGEMENTS': "We are grateful for the anonymous referee's valuable comments and suggestions that improve the manuscript. This work has been supported by the National SKA Program of China (Grant No. 2022SKA 0120101) and the National Key R&D Program of China (No. 2020YFC2201200), the science research grants from the China Manned Space Project (No. CMS-CSST-2021-B09, CMSCSST2021-B12 and CMS-CSST-2021-A10), and opening fund of State Key Laboratory of Lunar and Planetary Sciences (Macau University of Science and Technology) (Macau FDCT Grant No. SKLLPS(MUST)-2021-2023). C.Y. has been supported by the National Natural Science Foundation of China (Grant Nos. 11521303, 11733010, 11873103, and 12373071).", 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable request to the corresponding author.', 'REFERENCES': 'Bai X.-N., 2011, ApJ, 739, 51 Barker A. J., Latter H. N., 2015, MNRAS, 450, 21 Cui C., Bai X.-N., 2020, ApJ, 891, 30 Cui C., Bai X.-N., 2021, MNRAS, 507, 1106 Cui C., Lin M.-K., 2021, MNRAS, 505, 2983 Cui C., Tripathi A., Yu C., Lin M.-K., Youdin A., 2024, arXiv e-prints, p. arXiv:2407.02103 Flaherty K., Hughes A. M., Mamajek E. E., Murphy S. J., 2019, ApJ, 872, 92 Flaherty K., et al., 2024, in American Astronomical Society Meeting Abstracts. p. 106.04 Goldreich P., Schubert G., 1967, ApJ, 150, 571 Goodman J., Lackner C., 2009, ApJ, 696, 2054 Huang S., Yu C., 2022, MNRAS, 514, 1733 Klahr H., 2024, arXiv e-prints, p. arXiv:2404.15933 Klahr H., Hubbard A., 2014, ApJ, 788, 21 Latter H. N., Kunz M. W., 2022, MNRAS, 511, 1182 Latter H. N., Papaloizou J., 2018, MNRAS, 474, 3110 Lesur G., et al., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 465 Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023 Lin M.-K., Youdin A. N., 2015, ApJ, 811, 17 Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805 Lubow S. H., Pringle J. E., 1993, ApJ, 409, 360 Manger N., Klahr H., 2018, MNRAS, 480, 2125 Marcus P. S., Pei S., Jiang C.-H., Barranco J. A., Hassanzadeh P., Lecoanet D., 2015, ApJ, 808, 87 McNally C., Pessah M., 2014, ApJ, 811 Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610 Ogilvie G. I., Lubow S. H., 1999, ApJ, 515, 767 Papaloizou J. C. B., Pringle J. E., 1984, MNRAS, 208, 721 Pfeil T., Klahr H., 2019, ApJ, 871, 150 Pfeil T., Klahr H., 2021, ApJ, 915, 130 Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in FORTRAN. The art of scientific computing. Press Syndicate of the University of Cambridge Schiff L. I., 1968, Quantum Mechanics, 3rd edn. McGraw-Hill Education Stoll M. H. R., Kley W., 2014, A&A, 572, A77 Svanberg E., Cui C., Latter H. N., 2022, MNRAS, 514, 4581 Teague R., et al., 2016, A&A, 592, A49 Urpin V., 2003, A&A, 404, 397 Urpin V., Brandenburg A., 1998, MNRAS, 294, 399 Wang L., Xu S., Wang Z., Fang M., Goodman J., 2024, ApJ, 972, 142  \nYamaguchi M., et al., 2024, PASJ, 76, 437  \nYellin-Bergovoy R., Umurhan O. M., Heifetz E., 2021, Geophysical and Astrophysical Fluid Dynamics, 115, 674 Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588 Yu C., Lai D., 2013, MNRAS, 429, 2748 Zhang S., Zhu Z., Jiang Y.-F., 2024, ApJ, 968, 29', 'APPENDIX A: DERIVATION OF THE BOUNDARY CONDITIONS': 'In this appendix, we show the derivation of Equation (28). According to the WKB approximation 𝑦 1 , 𝑦 2 ∝ 𝑒 𝑖 ˆ 𝑘 𝑧 ˆ 𝑧 , we can readily know that 𝑑𝑦 1 / 𝑑 ˆ 𝑧 = 𝑖 ˆ 𝑘 𝑧 𝑦 1 , 𝑑𝑦 2 / 𝑑 ˆ 𝑧 = 𝑖 ˆ 𝑘 𝑧 𝑦 2 . Note that ˆ 𝑘 𝑧 here is a local constant associated with the boundaries 3 . Then Eq. (17) becomes  \n( 𝐴 11 -𝑖 ˆ 𝑘 𝑧 ) 𝑦 1 + 𝐴 12 𝑦 2 = 0 , 𝐴 21 𝑦 1 + ( 𝐴 22 -𝑖 ˆ 𝑘 𝑧 ) 𝑦 2 = 0 . (A1)  \nTo ensure the existence of nontrivial solutions of 𝑦 1 and 𝑦 2 , the matrix 𝑨 must satisfy  \nGLYPH<12> GLYPH<12> 𝑨 -𝑖 ˆ 𝑘 𝑧 𝑰 GLYPH<12> GLYPH<12> = GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 𝐴 11 -𝑖 ˆ 𝑘 𝑧 𝐴 12 𝐴 21 𝐴 22 -𝑖 ˆ 𝑘 𝑧 GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> = 0 . (A2)  \nThe equation (A2) can be rewritten as a quadratic equation for 𝑖 ˆ 𝑘 𝑧 :  \nGLYPH<16> 𝑖 ˆ 𝑘 𝑧 GLYPH<17> 2 -𝑖 ˆ 𝑘 𝑧 ( 1 + 𝑄 ) ˆ 𝑧 -Δ = 0 , (A3)  \nwhere Δ = -ˆ 𝜐 2 ( 𝜒 + 𝑈 ) + ( 𝜒 -1 ) ˆ 𝑧 2 [ 𝑄 -𝑈 ] . This quadratic equation for 𝑖 ˆ 𝑘 𝑧 can be readily solved. The two roots are  \n𝑖 ˆ 𝑘 𝑧 1 = 1 2 GLYPH<20> ( 1 + 𝑄 ) ˆ 𝑧 + 𝑖 √︃ 𝐴 3 -𝐴 2 ˆ 𝑧 2 GLYPH<21> , (A4)  \n𝑖 ˆ 𝑘 𝑧 2 = 1 2 GLYPH<20> ( 1 + 𝑄 ) ˆ 𝑧 -𝑖 √︃ 𝐴 3 -𝐴 2 ˆ 𝑧 2 GLYPH<21> . (A5)  \nNote that there are two eigenvectors associated with the two eigenvalues of the matrix 𝑨 4 . According to the standard algorithm in linear algebra, the two eigenvectors associated with 𝑖 ˆ 𝑘 𝑧 1 and 𝑖 ˆ 𝑘 𝑧 2 can be written as ：  \n𝒓 1 = GLYPH<18> 𝐵 1 -𝐴 1 ˆ 𝑧 + 𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 GLYPH<19> , 𝒓 2 = GLYPH<18> 𝐵 1 -𝐴 1 ˆ 𝑧 -𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 GLYPH<19> , (A6)  \nwhere  \n𝐴 1 = 2 𝜒 + 𝑄 -1 , 𝐴 2 = ( 𝑄 + 1 ) 2 + 4 ( 𝜒 -1 ) ( 𝑄 -𝑈 ) , (A7)  \n3 We adopt the Wentzel-Kramers-Brillouin (WKB) approximation, which is typically used for a semi-classical calculation in quantum mechanics in which the wave function is recast as an exponential function (Schiff 1968). That is the theme of our treatment for the boundary condition. This method aims to find approximate solutions to linear ODEs with spatially varying coefficients. For the global vertical shear instability analysis, the relevant disk physical variables are supposed to be viewed as wave components with vertical wave number that varies with the vertical positions. But when it comes to the boundaries, we actually focus on the very regions that are spatially local. As a result, the variations of wave number are small and could be approximated as a local constant. In addition to quantum mechanics, such WKB approximations are also widely used in the astrophysical wave dynamics (Li et al. 2000; Cui et al. 2024).  \n𝐴 3 = 4ˆ 𝜐 2 ( 𝜒 + 𝑈 ) , 𝐵 1 = 2 𝑖 ˆ 𝜐 ( 𝜒 + 𝑈 ) , (A8)  \nrespectively. Note that these two eigenvectors represent waves with different group velocities.  \nThe equation (26) can be written explicitly as:  \n( 𝑎 1 + 𝑎 2 ) = 𝑦 1 / 𝐵 1 , (A9)  \n( 𝑎 1 -𝑎 2 ) = 𝑦 2 + 𝐴 1 ˆ 𝑧 ( 𝑎 1 + 𝑎 2 ) 𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 = 𝑦 2 + 𝐴 1 ˆ 𝑧𝑦 1 / 𝐵 1 𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 . (A10)  \nWith some simple manipulations, we can arrive at:  \n𝑅 ≡ 𝑎 2 𝑎 1 = h -𝐴 1 ˆ 𝑧 + 𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 i 𝑦 1 -𝐵 1 𝑦 2 h 𝐴 1 ˆ 𝑧 + 𝑖 √︁ 𝐴 3 -𝐴 2 ˆ 𝑧 2 i 𝑦 1 + 𝐵 1 𝑦 2 , (A11)  \nwhich could be recast as a linear combination of 𝑦 1 and 𝑦 2 , namely:  \nGLYPH<18> 𝐴 1 ˆ 𝑧 + 𝑅 -1 𝑅 + 1 𝑖 √︃ 𝐴 3 -𝐴 2 ˆ 𝑧 2 GLYPH<19> 𝑦 1 + 𝐵 1 𝑦 2 = 0 . (A12)  \nThis is exactly Equation (28) in the main text. A larger 𝑅 indicates a greater outgoing component. Its non-axisymmetric counterpart, the Equation (30) could be obtained in a similar manner.  \nNote that our BCs are quite general and the free BC (zero Lagrangian pressure perturbation, i.e., Δ 𝑃 = 0) in Lin & Youdin (2015) can be recovered with an appropriate value of 𝑅 free . Numerical calculation shows that the modulus of | 𝑅 free | = 1. The free BC actually means a reflecting boundary with a phase difference between the incoming and outgoing waves.', 'APPENDIX B: DISPERSION RELATION AND IDENTIFICATION OF WAVE PROPAGATION DIRECTION': 'With the WKB approximation, the local dispersion relation is:  \n𝑐 4 ˆ 𝜐 4 + 𝑐 3 ˆ 𝜐 3 + 𝑐 2 ˆ 𝜐 2 + 𝑐 1 ˆ 𝜐 + 𝑐 0 = 0 , (B1)  \nwhere coefficients 𝑐 0 , 𝑐 1 , 𝑐 2 , 𝑐 3 and 𝑐 4 are  \n𝑐 0 = ˆ 𝑘 2 𝑧 + 𝑖 ( 𝑄 + 1 ) ˆ 𝑧 ˆ 𝑘 𝑧 , 𝑐 1 = 𝑖 𝛽 ( 𝛾 -1 ) ( 𝑄 -𝑈 ) ˆ 𝑧 2 , (B2)  \n𝑐 2 = 𝛽 2 𝛾 2 𝑐 0 + 𝑖 𝛽𝛾𝑐 1 -𝑈 -1 , 𝑐 3 = -𝑖 𝛽 ( 𝛾 -1 ) , (B3)  \n𝑐 4 = -𝛽 2 𝛾 [ 1 + 𝛾𝑈 ] , (B4)  \nrespectively. Note that this dispersion relation is equivalent to the Equation (A3). According to the definition of group velocity 𝑣 𝑔 = 𝜕𝜐 𝜕 ˆ 𝑘 𝑧 , we know that the group velocity is ：  \n𝑣 𝑔 = -( ˆ 𝜐 2 𝛽 2 𝛾 2 + 1 ) GLYPH<2> 2 ˆ 𝑘 𝑧 + 𝑖 ( 𝑄 + 1 ) ˆ 𝑧 GLYPH<3> / 𝐷 . (B5)  \nwhere 𝐷 = 4 𝑐 4 ˆ 𝜐 3 + 3 𝑐 3 ˆ 𝜐 2 + 2 𝑐 2 ˆ 𝜐 + 𝑐 1 . The inverse function rule is applied during the calculation because ˆ 𝑘 𝑧 depends on 𝜐 and not the other way around.  \nNote that there are two roots of the quadratic equation (A3) for ˆ 𝑘 𝑧 . Each root of ˆ 𝑘 𝑧 has a corresponding eigenvector. The specific form of each eigenvector is shown in Equation (A6). It is clear that, for the wave associated with 𝒓 2 , the group velocity 𝑣 𝑔 > 0 (outgoing) at the outer boundary where 𝑧 = 5 𝐻 . For the wave associated with 𝒓 1 , the group velocity 𝑣 𝑔 < 0 (incoming) at outer boundary 𝑧 = 5 𝐻 . The  \nwave assignment is general for all eigenfunctions and independent of the chosen parameters. It also holds in general for the case of 𝑚 ≠ 0.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}