Upon linearizing the equations of non-adiabatic gas dynamics with perturbations of the form exp(ω t + i k ⋅ x ), where ω = ω _R + iω _I is the complex-valued frequency of a mode with wavenumber k = | k |, two speeds emerge: the phase velocity v _p = ω _I /k that gives the propagation speed of the mode and the ‘condensation velocity’ v _c = ω _R /k that characterizes the speed of the advective flow that grows or damps the amplitude of the density perturbation. Condensation modes are non-propagating (v _p = 0) linear modes with v _c ≠ 0, whereas adiabatic acoustic modes are non-condensating (v _c = 0) with v _p ≠ 0 . This statement holds so long as all non-adiabatic source terms aside from a heat flux due to thermal conduction can be written as a volumetric term L = L ( ρ , T ) .

As discussed in detail in Section 2.3, the classifier ‘isochoric’ refers to the type of derivative of L / T that governs the stability of the isochoric modes; it does not, however, imply evolution that differs from the entropy mode during the linear growth phase. In a homogeneous gas (or when the local approximation holds), all three modes grow or damp according to the analytic solution: ρ x , t = ρ 0 + A ρ 0 e ω R t ⁡ cos k ⋅ x v x , t = v 0 − A v c e ω R t ⁡ sin k ⋅ x p x , t = p 0 − A ρ 0 v c 2 e ω R t ⁡ cos k ⋅ x , (1) where ρ ₀, v ₀, and p ₀ are the density, velocity, and pressure, respectively, of the uniform background flow; A is the perturbation amplitude; and v _c is the solution to the cubic dispersion relation, which can be expressed as v c = − N ρ k t cool R λ + v c / c s , 0 2 1 + v c / c s , 0 2 . (2) Here, c s , 0 = γ p 0 / ρ 0 and t _cool is the characteristic cooling time defined as t cool ≡ E 0 Λ 0 , (3) where E 0 = c v T 0 is the gas internal energy and Λ₀ is the cooling rate (in units of erg s⁻¹ g⁻¹) evaluated in the background flow. The quantity R _λ is the ratio R λ ≡ N p γ N ρ , (4) where N _p and N _ρ are dimensionless measures of how the net cooling rate varies with temperature: ¹ N p = T 0 Λ 0 T ∂ L / T ∂ T p T = T 0 + λ F λ 2 , (5) N ρ = T 0 Λ 0 ∂ L ∂ T ρ T = T 0 + λ F λ 2 . (6) The subscript on R _λ emphasizes the wavelength dependence of this ratio; in Eq. 5 and Eq. 6, λ _F is the Field length defined as λ F = 2 π κ 0 T 0 ρ 0 Λ 0 , (7) where κ ₀ is the (isotropic) thermal conductivity evaluated in the background flow.

Condensation modes can be classified by the solution branch that they occupy when solving Eq. 2. Equivalently, they can be classified by their stability criteria. The stability of the entropy mode is always governed by the sign of N _p . In particular, when thermal conduction is negligible (requiring either λ _F → 0 or λ ≫ λ _F in Eq. 5 and Eq. 6), the stability of the entropy mode depends only on the sign of the isobaric temperature derivative of L / T , i.e., they are unstable if the inequality ∂ L / T ∂ T p < 0 , (8) known as Balbus’ instability criterion (after Balbus, 1986), is satisfied. Likewise, the stability of the fast/slow isochoric modes is determined by the sign of N _ρ ; the instability criterion for them is ∂ L ∂ T ρ < 0 , (9) when neglecting thermal conduction. Note the temperature derivative is not of L / T in the isochoric case; this is the discrepancy mentioned above that we will examine further in Section 2.4.

Typically, only one of the condensation modes, namely, the entropy mode or the fast isochoric mode, will exist in the plasma for any given value of R≡ R_λ(λ_F= 0), the other two modes being acoustic. The main exception is the long-wavelength regime for a gas with R < 0; there is a critical wavelength λ _crit such that for λ > λ _crit, the acoustic modes transition into the fast and slow isochoric condensation modes, so all three modes can exist simultaneously. This critical wavelength is, to a good approximation, λ crit = 2 π | N ρ | − B R λ cool , (10) where B = 27 ( R − 1 / 3 ) 2 / 4 − 1 and λ _cool is the so-called cooling length that is defined as λ cool = c s , 0 t cool . (11) For λ < λ _crit, only the entropy mode exists; when the isochoric modes appear at λ > λ _crit, they are stable when the entropy mode is unstable (i.e., when N _p < 0) and vice versa. The R > 0 regime is very different, as the fast isochoric mode takes the place of the entropy mode as the main condensation mode accompanying acoustic modes, and these acoustic modes are also unstable (or to use proper terminology, they are overstable), their stability being governed by the isochoric criterion provided R < 1. Note by Eq. 4 that instability when R > 0 requires N _p < 0 and N _ρ < 0—both the isobaric and isochoric instability criteria are satisfied. For R > 1, the only condensation mode that can exist is the fast isochoric one; the overstability of acoustic modes is no longer determined by the sign of N _ρ but rather by the sign of the isentropic derivative ( ∂ L / ∂ T ) s , where s is the specific entropy.

A number of properties can be inferred from Eq. 1 and Eq. 2, making it clear under what circumstances a transition away from isobaric growth rates implies nonisobaric evolution in the nonlinear regime, regardless of which stability criterion is satisfied. Notice that the quantity A v c e ω R t is the maximum magnitude of the instantaneous velocity at which plasma is supplied to the slightly cooler gas by the slightly warmer gas to grow the density perturbation with time. At the end of the linear regime when A e ω R t ∼ 1 , v _c is the maximum velocity reached. It is therefore the characteristic velocity of advective flows upon entering the saturation phase of nonlinear growth. Stated differently, not only is the quantity v _c the solution to the dispersion relation of the linearized equations of gas dynamics but also its value reveals what type of dynamics to expect in the nonlinear regime of TI.

Our reference point is therefore the isobaric value of v _c that can be derived by taking the limit |v _c | ≪ c _s,0 in Eq. 2, namely, v _c = −N _p /(γt _cool k). Since v _c ≡ ω _R /k, this corresponds to the growth rate ω _R = −N _p /(γt _cool) derived by Field (1965) for the short-wavelength limit (when neglecting thermal conduction) that defines the isobaric regime. Unstable growth requires ω _R > 0, giving N _p < 0 as the instability criterion. This derivation of the maximum growth rate and Balbus’ instability criterion from taking the limit |v _c | ≪ c _s,0 thus reveals the dynamics associated with isobaric evolution; at the end of the linear regime, only tiny pressure gradients are necessary to give rise to flows with |v _c | ≪ c _s,0 that regulate the growth and saturation of entropy modes.

Thus, the nonisobaric regime can be defined as the wavelengths for which solutions to Eq. 2 do not satisfy |v _c | ≪ c _s,0. As mentioned above, the entropy mode is the only condensation mode that can exist for wavelengths λ < λ _crit and R < 0. It occupies the solution branch of Eq. 2 with the property that v _c rises smoothly from 0 in the short-wavelength (k → ∞) limit to | v c | = c s , 0 − R in the long-wavelength (k → 0) limit. Hence, as the wavelength increases, the entropy modes undergo a transition from isobaric to nonisobaric dynamics unless |R| ≪ 1. The two isochoric condensation mode branches only appear for λ > λ _crit and asymptote to | v c | = c s , 0 − R and | v c | = ( k t cool ) − 1 N ρ in the k → 0 limit. We named the mode with v _c ∝ k ⁻¹ the fast isochoric mode due to this divergent property. Because k can always be chosen small enough to make |v _c | > c _s,0, this mode is associated with an explosive saturation regime at sufficiently long wavelengths, the expected outcome being a self-fragmentation of the cloud after the saturation phase. With | v c | / c s , 0 = − R for both the entropy mode and slow isochoric mode in the k = 0 limit, the advective flows of long enough wavelength modes will become supersonic during saturation provided R ≲ − 1. In WP19, we speculated that this could also lead to self-fragmentation, referring to either instance of saturation leading to self-fragmentation as ‘splattering’. This has since been demonstrated by Gronke and Oh (2020), although they did not specify if the unstable condensation modes in their simulations were fast isochoric ones or entropy modes, and they incorrectly attributed the self-fragmentation seen to the now discredited (see WP19; Das et al., 2021; Farber and Gronke, 2022) ‘shattering’ hypothesis of McCourt et al. (2018) that is discussed below.

The nonisobaric behavior at wavelengths giving |v _c |≳ 0.1c _s,0, but before reaching the splattering regime at |v _c |∼ c _s,0, has been confirmed by other researchers and is the following. The response to the halting of the compression phase once the cooling gas reaches its equilibrium temperature is for the core of the condensation to ‘bounce’ and enter an expansion phase. Jennings and Li (2021) coined the term ‘contraction rebound’ to conceptualize this saturation dynamics. Further compressive/expansive motions recur as damped oscillations, causing the cloud to ‘pulsate’—the descriptor used by Gronke and Oh (2020). If this rich dynamics is studied using a single entropy mode for the initial conditions, the frequency of the pulsations is the linear theory growth rate, ω _R (see WP19).

In the literature on nonlinear TI, there are a number of conflicting claims about the differences between the short- and long-wavelength regimes and their relationships to isobaric and isochoric instability (e.g., Meerson, 1996; Burkert and Lin, 2000; Vázquez-Semadeni et al., 2003; McCourt et al., 2018; Mandelker et al., 2021). The most controversial claim came from McCourt et al. (2018), who presented simulations that appear to show an isochorically cooling background flow containing unstable perturbations spontaneously shatter during the nonlinear saturation phase of TI. Their study was aimed at establishing λ _cool evaluated in the cold phase gas after TI saturates as the more relevant length scale characterizing cloud sizes when compared with the already established value, namely, λ _cool evaluated in the background flow out of which the clouds had condensed (e.g., Perry and Dyson, 1985; Burkert and Lin, 2000). We showed in WP19 that λ _cool in the cold phase gas will generically be smaller than the Field length of the background flow, making McCourt et al. (2018)’s hypothetical ‘cloudlets’ subject to immediate evaporation (see Begelman and McKee, 1990). Additionally, we interpreted their simulations as being a demonstration not of ‘shattering’ but rather of ‘isobaric takeover’, a process described succinctly by Burkert and Lin (2000): “the fluctuations that can first reach non-linearity would dominate the growth of all perturbations with longer wavelengths and homogenize disturbances with smaller wavelengths. Thus, they determine the characteristic size and mass of the cold dense clumps that would emerge from the cooling of an initially nearly homogeneous cloud.” In other words, since McCourt et al. (2018) introduced a spectrum of perturbations into their initial conditions (a pre-existing cloud with a density contrast χ = 10), the short-wavelength condensation modes with the highest growth rates saturated while the ones with the longest wavelengths were still in the linear regime, making it appear as though the cloud as a whole underwent fragmentation. We return to this point in Section 4.1.

To quote again from Burkert and Lin (2000), when discussing the nonisobaric regime, they state: “Because of its long sound-crossing timescale, the perturbation cannot be compressed significantly while cooling; it cools almost isochorically (Parker, 1953).” We view this wording as confusing on account of the fact that as already mentioned, all condensation modes undergo compression according to Eq. 1 during the linear growth phase, regardless of the sound-crossing timescale. To clarify matters, we note that there are two common usages of the word isochoric in the literature on TI: i) a reference to the isochoric instability criterion or to the two condensation modes obeying this criterion and ii) a reference to long timescales for any small-amplitude thermal fluctuations in the background flow (which can be decomposed into a superposition of sinusoidal condensation modes of various wavelengths) to lead to changes in the density, i.e., for t _ρ ≡ |∂ ln ρ/∂t|⁻¹ ≫ t _cool. Because t _ρ can only be as large as the relevant dynamical time—the perturbation sound crossing time, t _cross, in local simulations—usage ii) entails having t _cool ≪ t _cross. This in turn requires distinguishing between the short- and long-wavelength condensation modes, but the concept of ‘long wavelength’ only makes sense when the modes become nonlinear so that they can excite and interact with sound waves to ‘be informed’ of their lengths. Once this interaction occurs, if t _cool ≪ t _cross, the sound waves cannot quickly communicate the pressure changes accompanying the thermal fluctuations. As we demonstrated in WP19, what happens instead of isochoric evolution (the density remaining approximately constant in this nonlinear interaction phase with sound waves) is the characteristic nonisobaric behavior described above: oscillations commence throughout the condensation. Thus, ‘isochoric evolution’ is not a meaningful concept unless the oscillations fully dampen, this mechanical equilibrium state being reached only if the condensation is free of thermal disturbances.

In summary, in the linear phase of TI, the descriptors ‘isobaric’ and ‘isochoric’ only have meaning as classifiers of the instability criteria. When discussing evolution in the nonlinear phase of TI, ‘isobaric’, ‘nonisobaric’, and ‘isochoric’ are associated with t _cross ≲ t _cool, t _cross ≳ t _cool, and t _cross ≫ t _cool, respectively. However, it is clear that the regime t _cross ≫ t _cool represents ‘extreme nonisobaric behavior’ that will not allow t _ρ ∼ t _cross. For this reason, we feel that in studies discussing the long-wavelength regime of TI, phrases such as ‘the evolution should be nearly isochoric’ should be interpreted as ‘the evolution should be highly nonisobaric’.

All of the possibilities for TI can be assessed graphically given the function L . In general, this function can only be computed numerically using, for example, a photoionization code, and so it is impossible to arrive at analytic expressions for the instability thresholds—the zero contours of N _p and N _ρ . Recently, however, in extending the theory of TI to radiation hydrodynamics (RHD), we arrived at a simple function L that does have analytic expressions for N _p and N _ρ , and the curves defining their zero contours (see Proga et al., 2022). In the Appendix, we take the optically thin limit of these equations to facilitate their use in hydro/MHD codes and to make our results here easily reproducible. In Figure 1, we plot the contour L = 0 on a pressure–density phase plane (referred to as just the ‘phase diagram’ hereafter). In textbook presentations of TI, typically only the regions where L > 0 and L < 0 are marked, masking any connection to nonlinear dynamics. We will see below how overplotting TI zones makes the phase diagram useful for understanding nonisobaric evolution. In Section 3, we employ a combined graphical and analytical analysis to uncover the saturation process.

We can use the cooling function defined in Eq. A2 in the Appendix to highlight the significance of the discrepancy mentioned at the beginning of Section 2, namely the existence of an isochoric TI zone according to Balbus’ criterion and the lack of unstable parameter space for isochoric modes according to Field’s criterion. Taking the temperature derivative at fixed density of Eq. A2 gives ∂ L ∂ T ρ = 1 2 c A ff ρ T − 4.5 a T 4 + 7 E r , s + c A C E r . (12) Here, E _r and E _r,s are both constant values of the radiation energy density (see the Appendix) and A _ff and A _C are positive constants. We see that this derivative is always positive and thus cannot satisfy Field’s instability criterion as given by Eq. 9. However, Balbus’ instability criterion for isochoric modes is [ ∂ ( L / T ) / ∂ T ] ρ < 0 , which expands to ∂ L ∂ T ρ < L T . (13) Substituting Eq. 12 into Eq. 13 yields an expression that can be satisfied, demonstrating how allowing for L > 0 increases the parameter space for TI. We define any circumstance where the background flow has L ≠ 0 as it being in a state of thermal non-equilibrium (TNE). There has been confusion in the literature regarding the onset of TI in a TNE state, which we address directly in Section 4.2.

At this point, it is useful to recognize that the contour where ( ∂ L / ∂ T ) ρ = L / T is equivalent to the zero contour of the dimensionless quantity N ρ ′ = T 0 Λ 0 T ∂ L / T ∂ T ρ . (14) The region on the phase diagram satisfying Eq. 13, or equivalently N ρ ′ < 0 , is what we refer to as the isochoric TI zone; the gray region in Figure 1 is therefore where the fast/slow isochoric modes are unstable and where the acoustic modes can be overstable. The TI zone marking where entropy modes are unstable is likewise defined by the zero contour of the quantity: N p ′ = T 0 Λ 0 T ∂ L / T ∂ T p . (15) We avoided terming this the ‘isobaric TI zone’ because the long-wavelength entropy modes that give rise to nonisobaric evolution in the nonlinear regime of TI are still nevertheless governed by the ‘isobaric’ instability criterion N p ′ < 0 . Note that, upon neglecting thermal conduction, Eq. 15 is identical to Eq. 5 (after evaluating it at a given T ₀), whereas Eq. 14 is equivalent to Eq. 6 only for L = 0 .

From Figure 1, notice that the gray region is a subset of the red region, implying both N p ′ < 0 and N ρ ′ < 0 there; in terms of the discussion from Section 2.1, along the equilibrium curve, this net cooling function lacks parameter space where R _λ < 0 due to N _p > 0 and N _ρ < 0. Again, the only way the isochoric TI zone will be entered is if the gas undergoes cooling with L > 0 , i.e., if the background flow is in a TNE state.

The bottom panel of Figure 1 illustrates our point about it being necessary to replace the concept of ‘isochoric evolution’ with ‘highly nonisobaric evolution’. The black dot marks an unstable position on the equilibrium curve. If the background flow is assigned the initial conditions at this location, then the linear growth phase of TI takes place on this dot because any noticeable deviation away from the initial conditions implies nonlinear amplitudes. The two sets of ‘tracks’ shown therefore represent the paths that can be taken by an unstable condensation mode once its amplitude becomes nonlinear. The horizontal set corresponds to isobaric evolution. A vertical set would correspond to isochoric evolution. What is shown instead is what actually happens in the long-wavelength regime: a condensation mode follows the equilibrium curve. This was first revealed by numerical simulations (see WP19), but it is clear in hindsight that this is what the nonlinear dynamics requires; when t _cool ≪ t _cross, the gas will be driven to equilibrium as it evolves. ²

The red solid line in Figure 1 is the boundary of the TI zone that we refer to as the Balbus contour. We terminate the tracks once they cross this contour to emphasize that this stage of evolution marks the beginning of the saturation process, as we show in the following section.

Our starting point is the evolution equation for the specific entropy (neglecting thermal conduction) D s D t = − L T , (16) where D/Dt = ∂/∂t + v ⋅∇ is the advective derivative. It is convenient to make this equation dimensionless. If we introduce the primed variables s′ ≡ s/c _v, t′ ≡ t/t _cool, L ′ ≡ L / Λ 0 , and T′ ≡ T/T ₀, this equation becomes D s ′ / D t ′ = − L ′ / T ′ for t _cool = c _v T ₀/Λ₀ as defined in Eq. 3. ³ To proceed without an excessive use of ‘primes’, all equations from here are understood to be written in terms of these primed variables.

As mentioned by Balbus (1986), applying a Lagrangian perturbation operator to Eq. 16 is the key to understanding TI anytime the flow is in a TNE state, having departed from the equilibrium curve. This hints at it being useful to apply a second advective derivative instead, to give D 2 s D t 2 = − D Dt L T . (17) In the Appendix, we show that for net cooling functions of the form L = L ( ρ , T ) , Eq. 17 is equivalent to D 2 s D t 2 = N p ′ D ⁡ ln ⁡ ρ D t − N ρ ′ D ⁡ ln ⁡ p D t , (18) where N p ′ and N ρ ′ are just the quantities from Eq. 14 and Eq. 15 expressed in dimensionless variables: N p ′ = T ∂ L / T ∂ T p ; N ρ ′ = T ∂ L / T ∂ T ρ . (19) Eq. 18 is therefore a dynamical relationship that connects the gas density and pressure gradients with the quantities determining the stability criteria and growth rates of the condensation modes. Whereas in linear theory, N _p and N _ρ are considered constant parameters describing the background flow, here N p ′ and N ρ ′ are time dependent, and this equation can be applied to the modes themselves to understand their nonlinear evolution. It becomes immediately clear that once the Balbus contour is crossed (i.e., once N p ′ changes the sign), D ln ρ/Dt and D ln p/Dt must also change, implying that the linear growth regime given by Eq. 1 has ended.

Less elegant versions of Eq. 18 are required for understanding the saturation dynamics of TI in detail. In going from Eq. 17 to Eq. 18, we have already assumed an ideal gas, making Eq. 18 not yet fully simplified. The specific entropy is s = ln(p/ρ ^γ ) up to an additive constant. Hence, either the left-hand side of Eq. 18 can be written in terms of D ² p/Dt ² and D ² ρ/Dt ², or either D ln p/Dt or D ln ρ/Dt can be eliminated in favor of Ds/Dt on the right-hand side. Retaining both forms of the second option is the key to understanding the full saturation mechanism. By analogy with Eq. 4, we define the ratio R ′ ≡ N p ′ γ N ρ ′ , (20) leading to, after a bit of algebra, D 2 s D t 2 + N p ′ γ D s D t = − N ρ ′ 1 − R ′ D ⁡ ln ⁡ p D t ; D 2 s D t 2 + N ρ ′ D s D t = − γ N ρ ′ 1 − R ′ D ⁡ ln ⁡ ρ D t . (21)

By eliminating Ds/Dt and D ² s/Dt ² using Eq. 16 and Eq. 17, respectively, and then solving for D ln p/Dt and D ln ρ/Dt, we arrive at D ⁡ ln ⁡ p D t = 1 N ρ ′ 1 − R ′ D D t L T + N p ′ γ L T ; D ⁡ ln ⁡ ρ D t = 1 γ N ρ ′ 1 − R ′ D D t L T + N ρ ′ L T . (22) The first of these final expressions reveals that the sign changes of D ln p/Dt are determined by the location on the phase diagram (the value of N p ′ , N ρ ′ , R′, and L / T ), as well as by D ( L / T ) / D t , the rate of change of the instantaneous entropy production rate as measured by tracking a fluid element. The latter quantity accounts for the kinematic motion of the fluid element that is determined by the continuity and force equations.

The second expression is interpreted similarly, but it is distinct in that N p ′ no longer appears in the expression in brackets. As we show in Section 3.3, for the dynamics taking place once the gas lands on the equilibrium curve, sign changes in D ln ρ/Dt are accompanied by sign changes in D ln p/Dt, and this is associated with ‘pulsations’—oscillations in the size of the cloud.

We now apply Eq. 22 to a single condensation mode in the nonlinear phase of TI to infer the sequence of dynamical events that must take place during saturation. We focus our analysis on the blue horizontal tracks in Figure 1, corresponding to nearly isobaric evolution within a TI zone (where N ρ ′ > 0 , N p ′ < 0 , and R′ < 0). It is clear from Figure 1 that N p ′ = 0 will occur well before the cooling gas reaches the equilibrium curve (at which point n ≈ 9 × 10¹¹ cm⁻³). Evaluating the D ln p/Dt equation at N p ′ = 0 gives D ⁡ ln ⁡ p D t = 1 N ρ ′ 1 − R ′ D D t L T . (23) This equation states that provided R′ < 1, there can be only two possibilities once the density has increased enough to place the cool gas at the Balbus contour: i) the sign of D ln p/Dt is the same as that of D ( L / T ) / D t or ii) both D ln p/Dt and D ( L / T ) / D t are zero. We will henceforth refer to the event of D ln p/Dt passing through zero as a ‘pressure reversal’. The proviso that R′ < 1 is confirmed in Figure 2.

Possibility ii) would imply that a pressure reversal accompanies a Balbus crossing and also coincides with the time that L / T reaches its maximum value. Our claim is that possibility i) would imply that all three events are not simultaneous but that a pressure reversal is nevertheless associated with the event of the gas reaching N p ′ = 0 . If the evolution were to exactly follow an isobar, then using the chain rule, D D t L T = ∂ L / T ∂ T p D T D t . (24) The Balbus contour is by definition where the first term in parenthesis vanishes, so in this instance, there is only possibility ii) above. Because perfect isobaric evolution is in violation of Eq. 1, we conclude that these events cannot be simultaneous. However, the causal link has been established and we expect the interval between a Balbus crossing and a pressure reversal to be a measure of (non)isobaricity.

To sort out the order in which the events occur, we can follow similar reasoning and express Balbus’ instability criterion from Eq. 8 as ∂ T ∂ t − 1 ∂ ∂ t L T p < 0 . (25) This expression reveals that unstable cooling gas has ∂ ( L / T ) / ∂ t > 0 at any fixed position within the overdense region of a condensation mode. It therefore tells us that the sign of D ( L / T ) / D t starts off positive during the linear growth phase. By Eq. 1, D ln p/Dt < 0 in the linear regime—physically, unstable cooling gas loses pressure support because the temperature drops faster than the density can rise. As the Balbus contour is approached, it is possible mathematically for D ln p/Dt to reach 0 before D ( L / T ) / D t does; by Eq. 22, D ln p/Dt = 0 corresponds to D D t L T = − N p ′ γ L T . (26) Above the Balbus contour N p ′ < 0 , giving D ( L / T ) / D t > 0 as required. However, this order of events implies that a pressure reversal was not preceded by a sign change in D ( L / T ) / D t or a change in stability conditions. On this basis, we can discount this outcome except in instances where the cooling gas passes through the isochoric TI zone (see Section 3.4).

The remaining outcome is the following sequence of events:

1. The cooling rate L / T reaches its maximum value prior to crossing the Balbus contour, i.e., D ( L / T ) / D t = 0 in a region where N p ′ < 0 . After this time, D ( L / T ) / D t < 0 .

Figure 2 shows that R′ < 1 during the entire path traced by the condensation mode as it isobarically approaches the equilibrium curve. This is the expected result based on linear theory (see Section 2.1). The sign changes in D ln p/Dt and D ln ρ/Dt are therefore controlled by the bracketed terms in Eq. 22. We will reserve the phrase ‘pressure reversal’ to apply exclusively to the circumstance that D ln p/Dt changes the sign due to crossing the Balbus contour. Other options for the sign changes will be referred to as ‘zero crossings’.

We have already mentioned that N p ′ does not appear in the bracketed term of the D ln ρ/Dt equation, and this means that contraction continues during and after a pressure reversal. Because D ( L / T ) / D t < 0 outside the TI zone after the Balbus contour is crossed, the first time that D ln ρ/Dt can reach 0 is, according to Eq. 22, when D D t L T = − N ρ ′ L T . (27) This event corresponds to the density reaching its maximum value, thus marking the transition from contraction to expansion. In other words, this marks the beginning of the contraction rebound process.

In an isobaric case, it is likely for Eq. 27 to be satisfied only once—this single contraction rebound event leading directly to the cloud reaching a steady state on the equilibrium curve. For nonisobaric evolution, the pulsations that occur after contraction rebound, as observed in numerical simulations (see WP19; Gronke and Oh, 2020), correspond to cycles between subsequent zero crossings of D ln p/Dt and D ln ρ/Dt, with each zero crossing satisfying Eq. 26 and Eq. 27. Notice that this will involve oscillations about the equilibrium curve L = 0 . A dedicated study of the D ln ρ/Dt equation in Eq. 22 is necessary to understand what parameters govern the pulsation damping rate.

Notice from Figure 1 that for the net cooling function given in the Appendix, it is possible to enter a regime of isochoric TI during isobaric evolution at just slightly higher pressures than that of the black dot. As mentioned in Section 3.2, this allows for a pressure reversal to occur inside the isochoric TI zone at a location where Eq. 26 can be satisfied. The only new constraint is found by setting N ρ ′ = 0 in Eq. 22 to give D ⁡ ln ⁡ ρ D t = − 1 N p ′ D D t L T . (28) This relationship must hold at the boundary of the isochoric TI zone, revealing that if density is to remain monotonically increasing upon both entering and exiting this region (throughout which N p ′ < 0 ), D ( L / T ) / D t must remain positive. Substituting Eq. 28 into the D ln p/Dt equation in Eq. 22 additionally shows that a pressure reversal will occur at this boundary if D ⁡ ln ⁡ ρ / D t = γ − 1 L / T there. The alternative is for L / T to reach its maximum value at the boundary of this zone, which would correspond to having D ln ρ/Dt = 0, thereby satisfying Eq. 27 also. This would imply that contraction rebound precedes the pressure reversal because D ⁡ ln ⁡ p / D t = − L / T by Eq. 22, i.e., the sign of D ln p/Dt is still that of the linear growth phase. The viability of either of these sequences of events has to be verified numerically.

Being a general identity of non-adiabatic gas dynamics, Eq. 18 may prove useful in a variety of settings. ⁴ In particular, Eqs 18–28 are unchanged in ideal MHD, hence so should be the overall saturation mechanism. However, the thermal evolution can differ in detail because in MHD, the paths traced on the phase diagram depend on the field strength (Bottorff et al., 2000, WP19). The linear theory of TI in MHD has recently been revisited by Claes and Keppens (2019), and it would be interesting to determine the plasma beta at which unstable condensation modes no longer carry out the sequence of events listed in Section 3.2.

Another obvious application is to the connection between TI and convective instability in stratified atmospheres (see Balbus and Soker, 1989; Binney et al., 2009; Balbus and Potter, 2016). To show that Eq. 18 is a generalization of results derived in that context, it is sufficient to examine one of the identities in Eq. 21 in a steady state. Taking D/Dt → v ⋅∇ gives, for the D ln p/Dt equation, v ⋅ ∇ v ⋅ ∇ s + N p ′ γ v ⋅ ∇ s = − N ρ ′ 1 − R ′ v ⋅ ∇ ⁡ ln ⁡ p . (29) Noting that N ρ ′ ( 1 − R ′ ) = N ρ ′ − N p ′ / γ by Eq. 20, we now also assume that the background state of the gas is in thermal equilibrium with L = 0 . This reduces the quantity N ρ ′ − N p ′ / γ to T [ ( ∂ L / ∂ T ) ρ − γ − 1 ( ∂ L / ∂ T ) p ] , which is simply related to the quantity ( ∂ L / ∂ T ) s by the thermodynamic identity (e.g., see Balbus, 1995): γ − 1 γ ∂ L ∂ T s = ∂ L ∂ T ρ − 1 γ ∂ L ∂ T p , (30) bringing us to v ⋅ ∇ v ⋅ ∇ s + 1 γ ∂ L ∂ T p ∇ s + γ − 1 γ ∂ L ∂ T s ∇ ⁡ ln ⁡ p = 0 . (31) If we now apply an Eulerian perturbation operator to this equation, the result is a product rule giving, schematically, δ v ⋅ […] + v ⋅ δ[…] = 0. As a last step, we assume the background atmosphere is static and spherically symmetric, leaving just the bracket attached to δ v , which simplifies to ∂ L ∂ T p ∂ s ∂ r = − γ − 1 ∂ L ∂ T s ∂ ⁡ ln ⁡ p ∂ r . (32) This equation is identical to Eq. 96 in Balbus and Potter (2016), who obtained it through a more involved Lagrangian perturbation analysis. As they explain, this relationship implies that TI will be accompanied by convective instability in an unstable static atmosphere.

In complicated flows, TI zone boundary crossings are likely common if not inevitable events during thermal evolution. Pressure reversals could therefore also play a role in mediating fragmentation caused by flow interactions. To provide an illustrative example, first note that a pre-existing cloud put in ‘by hand’ as initial conditions in non-adiabatic cloud–wind simulations and the fragments that appear as the cloud’s surface layers are disrupted, all occupy the same parameter space on a phase diagram as the cloud that is formed from TI. Referring to the tracks in Figure 1, as this gas changes its temperature upon interacting with the wind, the cold phase gas may repeatedly enter and exit the TI zone. A less disruptive example where crossing the Balbus contour may initiate the underlying fragmentation process is the ‘tearing’ mechanism identified by Jennings and Li (2021). By following individual fluid elements and applying Eq. 22, it should be possible to decipher the causal relationship between forces and the overall thermal evolution.

Yet, another topic concerns the dynamic role of pulsations in a regime of isochoric TI. Pulsations correspond to a change in the sign of D ln ρ/Dt, as discussed in Section 3.3. We have not yet studied nonisobaric evolution in a regime of isochoric TI using numerical simulations, and it is unclear if other groups have either. It is not standard practice to calculate the value of R _λ in Eq. 4 or R′ in Eq. 20 that characterizes any particular TI regime. This is especially relevant for a recent work on cloud coalescence instability; Gronke and Oh (2022) claim that this instability is aided by overstable acoustic modes, but such modes only exist when the isochoric or isentropic instability criterion is satisfied. Because they do not indicate which instability regime they consider, their claim of unstable acoustic modes being the agent that increases the coalescence rate is questionable. In our work on this instability (Waters and Proga, 2019a), we mentioned that coalescence can potentially be very fast; its rate can reach the rate set by the dynamical timescale provided the flow is continually subject to thermal disturbances, for this (re)excites pulsations and their accompanying advective flows that mediate the entire merger process.

In Section 2, we did not draw a distinction on ‘nonisobaric evolution’ as it relates to TI and nonisobaric dynamics more generally. The latter is simply the tendency for the gas with t _cool ≲ t _cross to undergo oscillations in response to a thermal disturbance, and the cloud coalescence simulations we presented in Waters and Proga (2019a) provide an example of this. The former refers to the nonisobaric dynamics accompanying the saturation of TI, and as the wavelength of a condensation mode increases, nonisobaric evolution involves increasingly strong pulsations following contraction rebound. Extreme nonisobaric behavior is characterized by ‘splattering’, a term we introduced to indicate a contraction rebound so strong that the cloud undergoes self-fragmentation.

The ‘shattering’ hypothesis by contrast, described by McCourt et al. (2018) as there being a tendency for highly nonisobaric clouds to undergo spontaneous self-fragmentation—as opposed to breakup caused by a dynamical response to contraction—has not held up to scrutiny (Gronke and Oh, 2020; Das et al., 2021; Jennings and Li, 2021; Farber and Gronke, 2022). As we have explained in WP19 and again here in more detail (Section 2.3), McCourt et al.'s (2018) simulations, with the initial conditions of having a spectrum of perturbations, can be understood as a clear example of ‘isobaric takeover’, where short-wavelength perturbations form many isobaric condensations within a much larger pre-existing cloud. In other words, those simulations do not follow self-fragmentation of one entity but rather the evolution of a nonisobaric cloud serving as the background flow within which the shortest wavelength condensation modes (that have the fastest growth rates) can form clouds and interact. Burkert and Lin (2000) described such an outcome, claiming that the end result is indistinguishable from a fragmentation process. While we disagree with this, especially on the grounds that isobaric takeover will lead to coalescence (the opposite of fragmentation), even accepting the claim does not invalidate our takeaway point that ‘splattering’ refers to a definite mechanism leading to self-fragmentation (namely, contraction rebound) while ‘shattering’ does not.

Despite our pointing out in WP19 that McCourt et al.'s (2018) simulations are a demonstration of isobaric takeover, the term ‘shattering’ is still regularly invoked when describing the appearance of small-scale cloud fragments in wind–cloud or shock–cloud interaction simulations that include radiative cooling (e.g., Sparre et al., 2020; Banda-Barragán et al., 2021; Bustard and Gronke, 2022; Jennings et al., 2023). The use of the term here is likely because the size of these fragments appears to be the cooling length evaluated in the cold phase gas, denoted by min(λ _cool), which is the isobaric length scale identified by McCourt et al. (2018). A cloud fragment with a characteristic size d _c has an associated crossing time t _cross = d _c /c _s and can be said to be isobaric if t _cross ≲ t _cool. Equivalently, because λ _cool = c _s t _cool is the length scale over which sound waves can effectively communicate thermal disturbances, isobaric cloud sizes satisfy d _c ≲ λ _cool, meaning that opposite sides of the cloud remain in sonic contact. The finding that interactions with the wind cause fragments to reach this length scale is not an unexpected outcome, hence this occurrence should not be confused with the regions of the cloud not interacting with the wind undergoing fragmentation, which is what ‘shattering’ would entail. Moreover, it is crucial to note that in Jennings et al.'s (2023) simulations, these isobaric fragments all evaporated in the runs with thermal conduction, consistent with our result that min(λ _cool) is generically smaller than the Field length (λ _F ) evaluated in the warm phase gas (see WP19). Thus, while min(λ _cool) is a characteristic size scale for cloud fragments, it is not a physically relevant one unless thermal conduction is highly suppressed.

While ‘splattering’ refers exclusively to a highly nonisobaric regime of TI, the pulsation behavior it relates to could permit the characteristic size of cloud fragments to significantly exceed the scale min(λ _cool) in non-adiabatic wind–cloud interaction simulations. These fragments are expected to become larger upon increasing the wind temperature (while keeping the wind pressure the same), based on the following reasoning. The descriptors ‘isobaric’ and ‘nonisobaric’ have thus far been used for both the size and behavior of the cold phase gas, but they also apply to the confining warm phase gas. A parcel of the warm phase gas is isobaric on scales d _w < max(λ _cool), where max(λ _cool) denotes evaluating the cooling length at the warm phase temperature. This gas can easily remain isobaric, depending on how much hotter it is compared to the cold phase, and can therefore be effective in dampening the pulsations of nonisobaric clouds when max(λ _cool) ≫ d _c . In other words, two equally sized clouds, with d _c ≫ min(λ _cool) (and with the same pressure) embedded in different temperature environments, will not exhibit the same nonisobaric behavior when subjected to the same thermal disturbance. ⁵ The pulsation damping rate is expected to be larger in a higher temperature environment, making that cloud effectively less nonisobaric and therefore less prone to fragmentation when interacting with a shearing flow. Even if this effect is significant, the fragments would still evaporate when they include isotropic thermal conduction because λ _F ∝ T ^11/4/p (under Spitzer conductivity and fixed gas pressure p), while max(λ _cool) ∝ T ^5/2/p. However, higher wind temperatures might allow larger fragments to survive when thermal conduction is anisotropic.

In the solar physics literature, the relationship between TNE and TI is explained with reference to the TNE-TI cycle that describes observed phenomenology (e.g., see Antolin and Froment, 2022), specifically the dynamics taking place in coronal loops (magnetic flux tubes that extend from the chromosphere into the corona but are rooted to the solar surface at both ends). TNE itself has been described as a process that occurs when the heating rate profile in a loop drops off sharply with height and is unchanging with time, while the cooling rate remains a local quantity, making it possible for there to be no ‘nearby’ thermal equilibrium state (in ρ-T phase space) where cooling can reach a balance with heating (Klimchuk, 2019). Both evaporative and bulk flows can occur in response to pressure changes during a TNE cycle.

Klimchuk (2019) claimed that it is meaningless to consider TI occurring under TNE conditions on the (false) premise that diagnosing TI requires the background flow to be in a steady state. We stress that Balbus's (1986) stability criterion supersedes Field’s criterion for the very purpose of applying to TNE conditions. For there to be a consensus on this issue, the concept of TNE as given in the preceding paragraph would have to be generalized to adopt the definition that we have used in Section 2.4, namely, TNE is simply the circumstance that flow departs from the equilibrium curve, meaning L ≠ 0 . We should furthermore specify that the bulk flow is subjected to L ≠ 0 on dynamical timescales to distinguish TNE from thermal misbalance, which simply refers to there being fluctuations about L = 0 on very short timescales due to acoustic compression and expansion cycles that result in the damping of MHD waves (Kolotkov et al., 2019; Kolotkov et al., 2021).

It follows that under TNE conditions, the criterion for entropy modes to be thermally unstable is given by Eq. 8, which expands to ∂ L ∂ T p < L T . (33) This reduces to the isobaric instability criterion of Field (1965) when the gas reaches thermal equilibrium with L = 0 . In the astrophysics literature, the canonical example of a plasma that violates either criterion is one with (e.g., McCourt et al., 2012; Mościbrodzka and Proga, 2013; Balbus and Potter, 2016) L = A ρ T d − B , (34) where A, B, and d are constants satisfying A > 0, B ≥ 0, and d < 1. This is, however, the same example used by Klimchuk (2019) in an attempt to illustrate a condensation forming under TNE conditions in the absence of TI. ⁶

Because the gas can be in TNE due to a variety of causes, a further distinction has to be drawn between local TI and dynamical TI. The theory of local TI, as summarized in Sections 2.1–2.2, is not confined to a homogeneous gas; it applies whenever the so-called local approximation holds, i.e., on length scales Δx over which Δx ≪ λ _q , where λ _q ≡|∇ ln q|⁻¹ is the scale length for the gradient of the background flow quantity q to be significant. Dynamical TI is the situation when the growth rates predicted by local TI cease to be valid because Δx ∼ λ _q for at least one of the relevant flow variables among q = (ρ, v , p, T).

Given the above considerations, the background flow for local TI can be an initially uniform region that is in an evolving equilibrium state or in a TNE state. As discussed in Section 2.4, the various isochoric instability regimes are associated with the gas starting off in TNE, as the entire isochoric TI zone is far from the equilibrium curve. In Proga et al. (2022), we provided a qualitative example of a column of gas starting off thermally stable but becoming unstable to local TI upon following a time-dependent equilibrium curve. In dynamical TI, by contrast, a flow can attempt to evolve along the equilibrium curve to maintain L = 0 , but the timescale associated with adiabatic cooling can become shorter than t _cool, thereby causing the flow to tend toward a steady state at some position off the equilibrium curve where L ≠ 0 .

If the coronal rain observed in coronal loops is attributable to the saturation of condensation modes, then due to the presence of flows and a chromosphere–corona transition region in these loops, the condensations are a clear instance of dynamical TI. Rather than the black dot in Figure 1, which represents a constant density initial condition appropriate for local TI, let us consider instead a possible initial coronal loop TNE state having a density spanning the range n ≈ 10¹¹−10¹² cm⁻³ and a nearly flat pressure profile that places it below the equilibrium curve in a region of net heating (with L < 0 ). The plasma will approach the equilibrium curve with time, but the flow traversing a converging (diverging) portion of the loop will undergo slight adiabatic heating (cooling), making the phase diagram ‘tracks’ of this flow spread vertically at different rates even if the heating rate is uniform. The plasma with n ≳ 1.5 × 10¹¹ cm⁻³ is very close to the stable cold branch of the equilibrium curve and will therefore not be subject to runaway heating. However, the stable equilibrium state of the plasma with n < 1.5 × 10¹¹ cm⁻³ is the hot branch where n ≲ 10¹⁰ cm⁻³. If pressure equilibrium can be closely maintained, the bulk of the coronal region plasma may approach this branch without ever entering the TI zone.

This is, in essence, a TNE runaway heating scenario that leads to a multi-temperature plasma without ever invoking TI. We therefore agree with the overall premise of Klimchuk (2019), but our example above debunks his assertion that “the physics that governs the thermal runaway in a TNE loop is equivalent to the physics that governs the thermal runaway in a [thermally] unstable equilibrium loop”. To examine condensation formation occurring under TNE conditions because of TI, we could alternatively imagine initial conditions that place the gas above the equilibrium curve in Figure 1 in a region of net cooling. Here, we would have to confront both bulk runaway cooling and the saturation of unstable entropy modes, but we have hopefully made it clear that these are distinct processes.

Either coronal loop scenario is complicated by the fact that equilibrium curves are in general time dependent, i.e., the contour L = 0 can change with time according to how the radiation field or other background sources of heating evolve. If this evolution occurs on timescales shorter than t _cool, the plasma cannot maintain L = 0 even if it can settle on the equilibrium curve and will thus re-enter TNE. The heating and cooling rates used for illustrative purposes in this article allow studying this situation. Namely, the equilibrium curve in Figure 1 models a column of gas being irradiated by both a thermal and non-thermal source of X-rays, and the two parameters T _C,h and f _h (controlling the temperature and flux of the non-thermal photons, respectively) can be made time dependent to vary the effective Compton temperature, which controls the shape of the equilibrium curve (see the Appendix).

We have recently worked on several other applications of dynamical TI that reveal the interplay between TI and TNE. In Waters et al. (2021), we have shown that the multiphase radial outflow solutions discovered in 1D by Dannen et al. (2020) can reach a steady state, permitting a formal stability analysis of an outflow occupying parameter space with L ≠ 0 . There we also investigated the interconnection between TI and runaway heating in 2D simulations of accretion disk winds driven by external irradiation from X-rays, showing that it is possible to identify the presence of local TI even in a time-dependent, turbulent flow by plotting the ‘tracks’ of individual streamlines on the phase diagram to see if they pass through a TI zone. Finally, in Waters et al. (2022), we have addressed a different point that was raised by Klimchuk (2019). Referring to the stratified flow within a magnetic flux tube, Klimchuk posed the question: “if a perturbation grows, does it have time to reach a substantial amplitude before it is carried to the chromosphere by the flow?” We have shown how to calculate whether or not outflowing entropy modes have time to saturate when they sample a time-dependent growth rate as they ‘fly through’ a TI zone.