{'fraction_non_alphanumeric': 0.03538631076282601, 'fraction_numerical': 0.013223950042855395, 'mean_word_length': 4.49372197309417, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability of thermal structures with an internal heating source

We study the thermal equilibrium and stability of isobaric, spherical structures having a radiation source located at its center. The thermal conduction coefficient, external heating and cooling rates are represented as power laws of the temperature. The internal heating decreases with distance from the source r approximately as exp(-tau)/(r**2), being tau the optical depth. We find that the influence of the radiation source is important only in the central region, but its effect is enough to make the system thermally unstable above a certain threshold central temperature. This threshold temperature decreases as the internal heating efficiency increases, but, otherwise, it does not depend on the structure size. Our results suggest that a solar-like star migrating into a diffuse interstellar region may destabilize the surrounding medium.


I. INTRODUCTION
Most of the structures observed in Astrophysics can be fully explained only if one takes into account several physical processes simultaneously, such as dynamics, self-gravitation, magnetic fields, radiation transfer, thermal conduction, complex chemical networks, etc.
Since the seminal paper by Field [1], a lot of work has been devoted to the study of the thermal and mechanical equilibrium of structures under the action of different physical mechanisms [2]. One of the motivations behind this kind of studies has been to explain the formation of structures via thermal instabilities. Nowadays, it is clearly established that thermal instabilities play a major role in several astrophysical contexts, such as solar chromosphere and corona, interstellar medium, planetary nebulae, and intergalactic medium [2,3,4,5,6]. A full modeling of these systems may be a relatively difficult task requiring large computational resources (e.g. see Refs. [7,8,9]). On the contrary, the development of simple models including few physical processes allows to clarify the specific role of each process, setting up the framework for more realistic models.
When only thermal conduction and a net heat-loss function are considered, a variety of (isobaric) thermal structures are obtained and stability criteria can be derived analytically [10]. These criteria can be extended analytically up to the second order [11]. In spite of the simplicity of this approach (many physical processes are not considered) these results can be applied successfully to the atomic phase of the ISM [12], molecular clouds [13], and the solar corona [14]. Usually, it is assumed that the heating mechanisms are "local" functions, in the sense that they depend only on the local thermodynamical state (temperature, density).
However, two factors may complicate this picture. On the one hand, the heating mechanism could be an explicit function of the position in the structure. This situation applies, for example, if a heating source is located close to the structure, such as an interstellar region with a nearby star (or star cluster), or the solar corona heated by the chromosphere. For coronal loops, it has been shown that the existence or not of thermal equilibrium solutions and their stability is subordinated to the spatial dependence of the heating [14]. On the other hand, if the heating mechanism is associated with radiative processes, the opacity adds a non-local character to the heating function making it difficult to perform an analytical study.
In Ref. [15] the stability of structures heated by radiative processes was analyzed. Generally, the results showed that the effect of increasing the opacity is to increase the stability of the thermal structures. However, in order to achieve certain analytical results that work was limited to the unrealistic case of a slab-like configuration irradiated simultaneously on both sides.
In the present paper, our attention is focused on investigating the thermal equilibrium and the stability of spherical structures heated by a radiation source located at its center.
Besides the pure academic interest of generalizing the previous results [10,11,12,13,15] by including an internal heat source, this kind of study could be very useful because this is a common configuration, in particular in Astrophysics. We try to keep the approach as general as possible and, therefore, our results can in principle be directly applied to any situation in which an isobaric, spherical structure has an internal radiation source. However, we use values for constants and parameters that are appropriate for typical interstellar conditions because in the future we plan to generalize this work in order to be applicable to interstellar regions in the Galaxy. It is known, for example, that the diffuse ISM is linearly thermally stable but nonlinearly unstable (e.g. [16]). Moreover, it is expected that stars systematically pass through the ISM. Thus, it is relevant to study whether these stars can sufficiently perturb an interstellar region as to induce thermal instability. Obviously, many more physical mechanisms have to be considered to describe this situation adequately but, as mentioned before, in a first approximation this study can help to understand the role played by this type of realistic heating. In Sec. II we present the basic equations, and in Sec. III we briefly explain the stability criterion we are using. The effect of the internal source on the stability is discussed in Sec. IV and, finally, the main conclusions are summarized in Sec. V.

II. BASIC EQUATIONS
The equation of energy conservation for a static, thermally conducting fluid, with pressure p constant in time t can be written in the form [17]: where T is the temperature, ρ the density, c p the specific heat per unit mass, and κ(ρ, T ) is the thermal conduction coefficient. The functions Γ(ρ, T ) and Λ(ρ, T ) represent the heating and cooling rates per unit volume, respectively. In many astrophysical situations the quantities κ, Γ and Λ can be approximated as being proportional to power laws of the form ∼ ρ α T β [18]. Under the assumption of isobaricity, we can eliminate the density as an independent variable and write where the constants κ 0 , Γ 0 , Λ 0 and the indices k, m and n are given by the physical processes involved. For instance, a typical interstellar region in the range of temperature 10 2 < ∼ T < ∼ 10 4 K with thermal conduction by neutral particles, heating proportional to the density, and cooled by collisions between particles, may be characterized by k = 1/2, m = −1 and n = −3/2 [12,18].
In this work we consider an additional heating mechanism. Let us assume a spherically symmetric region with radius R b that has a radiation source located at its center. The photons coming from this source can heat the surrounding gas. In the ISM this heating is mainly because the UV radiation field is able to photoeject electrons from atoms and dust which subsequently thermalize through collisions [3]. The energy absorbed per unit volume is [19]: where F (r) is the radiation flux coming from the central source, n is the number of atoms per cm 3 , σ is the photon capture cross section, and the efficiency parameter γ is the fraction of incident energy which is actually absorbed by the gas. The radiation flux decreases with distance from the center (r) as 1/r 2 . In order to avoid dealing with divergences at r = 0 when integrating Eq. 1, we impose a spatial smoothing scale R s . Thus, the radiation flux can be represented in the form [20]: where τ is the total optical depth from the central source to r: The value of the arbitrary constant a determines how quickly the behavior changes from ∼ F (0) when r ≪ R s to ∼ 1/r 2 when r ≫ R s . We have kept a unchanged (a = 10), but several tests showed that our results do not depend on a as long as the condition R s ≪ R b remains fulfilled. The flux at the center can be calculated by assuming emission from a black body of effective temperature T 0 : B ν being the Planck's function. The integration is done in the far ultraviolet (900 − 1100Å) because this is the range in which photons heat the ISM efficiently [21,22].
With these considerations, Eq. 1 can be rewritten as This is the integro-differential equation whose stationary solutions (and their stability) we want to study in this work. It is convenient, however, to express it in terms of some nondimensional quantities. If T eq is the temperature to which Γ = Λ, i.e.: then we can define a non-dimensional temperature as: Given the (external) heating and cooling mechanisms, T eq is a constant representing the temperature at the thermal equilibrium when there exists no internal heating (Γ s = 0). We also write the position normalized to the region size: We have kept the smoothing length fixed at the arbitrary value z s = 10 −8 , but, as mentioned before, our results do not depend on the exact z s value as long as z s ≪ 1 is satisfied. The total optical depth from the center to the boundary for the particular case T = T eq = constant is: where the ideal gas equation of state p/k B = nT (being k B the Boltzmann's constant) has been used. Therefore, including Eq. 9 to 11 and normalizing the flux at the center to the equilibrium flux, i.e.F = F (0)/F eq where F eq is the black body flux calculated using T = T eq , we see that the stationary solutions of Eq. 7 are the solutions of: where and To solve Eq. 12 we need two boundary conditions for the temperature. This paper is aimed to quantify the effects of an internal heat source on the stability of the steady structures taking as reference the previous papers for which the stability criteria are known. Then, for convenience, we assume exactly the same boundary conditions used in Refs. [10,11,12,13,15,18]: The first condition follows from the symmetry of the configuration: the resulting temperature distribution from the center z = 0 to the z = +1 side should be identical to the z = −1 side. This is the only condition compatible with this property if we impose the additional physically motivated requirements that both the temperature and its derivative exist and be finite. The second condition, i.e. fixed boundary temperature θ b , is necessary in order to integrate the energy equation. We write down explicitly this condition (fixed temperature) at z = 1 only to be consistent with the previous papers, but actually a fixed value at any other position would be equally valid as long as the integration over z can be performed. In fact, we are assuming a fixed central temperature and the integration of Eq. 12 is performed outwards from this central position (see Sec. IV). This is equivalent to converting the two point boundary value problem (which requires considerably more numerical effort) into an initial value problem. In practice there is no difference between both approaches because the result of interst for this work is the entire family of solutions rather than a single solution and, therefore, we have to span the whole range of possible central (and boundary) temperatures.

III. STABILITY CRITERION
When η = 0 (no internal heating source) Eq. 12 reduces to the case studied in Ref. [10]. In spite of the simplicity of the case η = 0, its study has already led to a better understanding of the general problem of stability of structures, even allowing to extend the study analytically up to the second order [11]. We see that for the case η = 0 the trivial solution θ = 1, i.e. the thermal equilibrium solution T = T eq , is solution of Eq. 12. It has been shown that this solution is linearly unstable if λ fulfills the condition [10]: The dimensionless parameter λ (Eq. 13) measures the ratio between heating (or cooling) at equilibrium and thermal difussion through the scale length R b . Moreover, thermal diffusion is a stabilizing factor because it tends to diminish any temperature fluctuation. Thus, physically, the condition λ > λ cri means that the generation of heat is so large that it cannot be removed efficiently by thermal difussion, and therefore a thermal instability develops. The analysis of the rest of the stationary solutions (different from the trivial one) cannot be done analytically but the stability can be inferred from the position on the curve θ 0 vs. θ b . In curve for the case λ < λ cri . This case corresponds to the stable solution according to Eq. 17.
On the contrary, the trivial solution is on the negative slope branch for the case λ > λ cri (unstable solution). Moreover, for λ = λ cri (marginally stable solution) the trivial solution separates the positive and negative slope branches. According to Ref. [10] (see also Refs. [17,18,23]), the positive slope branch can be identified as the stable branch. Therefore, stationary solutions located on the positive (negative) slopes branches are expected to be linearly stable (unstable). This is the stability criterion we are using in this work.
For the more general case (η = 0) we see that θ = 1 is not a stationary solution of Eq. 12. Both the spatial dependence of the internal heating and its non-local character (i.e., the local heating depends on material state in other points through the optical depth) complicate the mathematical treatment in such a way that it prevents us from finding an analytical stationary solution. Thus, a stability analysis is not possible but nevertheless we expect the above stability criterion to remain valid. The argument for this claim comes from the relative behaviors of thermal conduction and net heating/cooling with respect to the temperature. When two stationary solutions are obtained for a given boundary temperature (for example when θ b > 1 for the case λ = λ cri in Figure 1), only one of them is stable (see Refs. [17,23]). Let us assume, for example, that both thermal conduction Q(T ) and net heating rate L(T ) are increasing functions of temperature. Two stationary solutions means that Q = L at two temperature values (T = T eq ). Necessarily, we have dQ/dT > dL/dT at one point and the opposite in the other. The point at which dQ/dT > dL/dT is stable because any small temperature variation will be eventually damped by thermal conduction.
On the other hand, if the heating rate increases with temperature faster than the thermal conduction then any small temperature rise will grow. This argument should remain true as long as the temperature dependences remain the same. For a fixed position r, the dependence of the internal heating on T is the same as in the previous papers (∼ T −1 , i.e. proportional to density). This is the reason why we expect the stability criterion to be valid for the case under study.
In any case, we performed some numerical tests in order to check the validity of this criterion for the more general and complicated case of a non-local heating mechanism. For a given set of free parameters (λ and η) we first find numerically the corresponding stationary solutions for central temperatures around the "turning point temperature", i.e. around the central temperature value that separates the positive and negative slope branches. Then we impose different temperature perturbations (we use gaussian, exponential, and sinusoidal profiles) around randomly chosen positions in the structure. What we did was to check that ∂T /∂t changes its sign exactly at the turning point, in such a way that on the positive slope branch the sign is opposite to the perturbation (stable branch) and viceversa. We performed many tests with different combinations of parameters and the results obtained were the same. Notwithstanding this is not a demonstration, these tests allow us to be much more confident about the validity of the stability criterion for the case under consideration.

IV. EFFECT OF THE INTERNAL HEATING
To numerically solve Eq. 12 with the boundary conditions (15)(16) we use the Runge-Kutta method of order 4 with adaptive step size control [24]. The distinctive feature of Eq. 12 is the fact that to evaluate the internal heating term we need to know the solution we are looking for. The density is linked to the temperature (because of the isobaric assumption) and therefore the integral in Eq. 5 ultimately depends on the temperature distribution from the source (as we have indicated explicitly in Eq. 12). Thus, we have to proceed by successively iterating Eq. 12 until θ(r) converges to the stationary distribution. Then, given a central temperature θ 0 , in the first step we assume an initial distribution θ(z) = θ 0 to obtain an initial estimate of the optical depth τ (z). Then we use this result to numerically integrate Eq. 12 and obtain a new temperature profile. Now this temperature distribution is used to recalculate τ (z), and so on. The initial conditions are always satisfied, thus the central temperature remains constant during these steps. The procedure is then repeated until a convergence criterion is satisfied. In particular, we require that temperature values computed at two successive steps differ less than a prefixed small tolerance. Furthermore, we require this condition to be met at several points inside the structure. The final result is the whole temperature distribution from the center to the boundary. The central temperature can then be changed to search for other solutions with their corresponding boundary temperatures.
The coefficients λ and η are determined by the physical processes involved (pressure, external heating, cooling, thermal conduction, chemical composition, etc.). Here we consider typical ISM conditions. The mean pressure of the warm gas in the ISM is p/k B ≃ 10 3 K cm −3 [25], which is heated mainly (in absence of nearby stars) by the galactic radiation field. Under these conditions, we have Γ 0 = 2.7 × 10 −23 , Λ 0 = 10 −21 , m = −1, and n = −1.5 [3,26] (in this work Γ and Λ are in units of erg cm −3 s −1 ). This yields an equilibrium temperature T eq ≃ 1370 K. Under these conditions, the thermal conduction is mainly due to neutral particles and therefore κ 0 = 2.5 × 10 3 and k = 1/2 [27]. The cross section is assumed to be σ = 6.3 × 10 −18 cm 2 [3]. We keep as constants these fiducial values whereas the region size (R b ) and the internal heating efficiency (γ) are treated as free parameters. From Eqs. 13 and 17 we get that the case in which the trivial solution is marginally stable corresponds to structures with radius R cri ≡ λ cri κ 0 T k−m+1 eq /Γ 0 1/2 ≃ 0.12 pc. The total optical depth (in thermal equilibrium) would be τ 0 ≃ 1.7. As mentioned before, λ measures the ratio between net external heating at equilibrium and thermal conduction through the structure.
The trivial solutions is stable when λ < λ cri because thermal conduction smooths the energy gradient more efficiently in relatively small structures (R b < R cri ). On the other hand, η gives in some way the relative importance of the internal heating and thermal conduction terms. Thus, given three different mechanisms of energy transport (thermal conduction, net external heating, and internal heating) there are only two free parameters (R b and γ, related to λ and η through Eqs. 13 and 14) determining the relative importance among them.
The temperature distributions inside the structure for the case R b = R cri and γ = 10 −3 [19] are shown in Fig. 2. Each line in this figure is one solution θ(z) of Eq. 12 for a given value of central temperature θ 0 and, therefore, of internal heating. We see that there are no solutions with boundary temperature (θ b ) below 1. Then, all the solutions with θ 0 < 1 have positive temperature gradients. As we will see later, these solutions are thermally unstable.
For θ 0 > 1, we have that θ 0 increses as θ b increases, but for θ 0 > ∼ 4 this behavior changes notoriously. For the highest central temperature shown, we can see a sign change in the temperature gradient around z ≃ 0.55. To understand this behavior, we have calculated the relative contribution of each term in Eq. 12 and we have seen that at central temperatures θ 0 < ∼ 2 the internal heating is always much smaller than the (external)heating/cooling term. Thus, the situation is almost the same as if γ = 0, i.e., (external)heating/cooling balanced only by thermal conduction. In contrast, at higher θ 0 values the internal heating dominates over the (external)heating/cooling, but this internal heating decreases very quickly with the distance from the center and then a change in regime occurs at some point (depending on the central temperature). For the highest central temperature shown in Fig. 2, this change in regime occurs precisely at z ≃ 0.55. We will see in next figure that this central temperature is very close to the maximum central temperature for which there exists a solution. Above this maximum, the temperature gradient becomes so large that the temperature goes down to unphysical negative values (i.e. there is no physical solution). For central temperatures between ∼ 2 and this maximum temperature the internal heating term always dominates at z < ∼ 0.5. Therefore, the influence of the radiation source seems to be important only around the central region (for the values of constants and parameters considered here), but its effect can be enough to make the system unstable, as we will show next.
The θ 0 vs. θ b diagram for R b = R cri and various values of γ is shown in Fig. 3 occurs at θ 0,max ≃ 3.9, but clearly θ 0,max decreases as γ increases.
In Fig. 4 we have fixed the internal heating efficiency (γ = 10 −3 ) and we have varied the radius R b . At low central temperatures the behavior is similar to the case γ = 0, i.e., the trivial solution is stable for R b < R cri and unstable for R b > R cri . At relatively high central temperatures, we obtain that the threshold temperature θ 0,max for thermal instability (indicated as open circles in Fig. 4) is approximately constant. This interesting result is a direct consequence of the "local" character of the internal heating: this heating mechanism decreases too quickly with z that its destabilizing effect does not depend on the structure size.
It can be seen that there is a small variation in θ 0,max for the extreme case R b = 10 −6 R cri , but this value is so small (∼ 5R sun for our fiducial values) that its validity is questionable.
In the case of a star that migrates into a diffuse interstellar region, the quantity θ 0,max would represent the photospheric (effective) temperature of a star above which the surrounding medium becomes thermally unstable. We have to point out that this simple interpreta- with γ = 10 −3 and R b = R cri ≃ 0.12 pc we obtain θ 0,max ≃ 3.9. An interesting result is that θ 0,max does not depend on, or depends very weakly on, R b (at least under the conditions and assumptions considered in this work). In other words, our results suggest that a solar-like star with surface temperature ∼ 5300 K may possibly destabilize the medium in its vecinity, almost independently of the region size. After this, the surrounding medium will undergo a transition to another state (or phase). Differences in the chemical composition and/or the amount of dust from region to region in the Galaxy (or in other galaxies) may change the value of γ [19], and consequently the maximum central temperature. Here we obtain θ 0,max = 3.6 for γ = 10 −2 and θ 0,max = 3.2 for γ = 10 −1 .
If stars eventually migrate into interstellar regions, then thermal instabilities caused by these stars may play a non-negligible role in stimulating phase transitions in the ISM. However, the present model is clearly too simple to be directly applicable to the ISM. The strongest restriction is perhaps isobaricity, which is a difficult condition to be fulfilled even in the diffuse ISM [25]. Notwithstanding these limitations, our results suggest that thermal instabilities driven by stars could be a key mechanism in the ISM.
We have studied the equilibrium and stability of isobaric, spherical structures with an internal radiation source located at its center. We find that there exists a threshold central temperature above which the structure becomes thermally unstable, and this temperature decreases as the internal heating efficiency increases. Interestingly, the threshold temperature does not depend on the region size, because the destabilizing mechanism has only a local effect. We have shown that thermal instability triggered by an embedded source may have interesting consequences in the ISM evolution which should be addressed in future, more complete, studies.