Sustained oscillations in Interstellar chemistry models

Non-linear behavior in interstellar chemical models has been recognized for 25 years now. Different mechanisms account for the possibility of multiple fixed-points at steady state, characterized by the ionization degree of the gas. Chemical oscillations are also a natural behaviour of non-linear chemical models. We study under which conditions spontaneous sustained chemical oscillations are possible, and what kind of biffurcations lead to, or quench, the occurrence of such oscillations. Methods. The well known Ordinary Differential Equations (ODE) integrator VODE is used to explore initial conditions and parameter space in a gas phase chemical model of a dark interstellar cloud. We recall that the time evolution of the various chemical abundances under fixed temperature conditions depends on the density over cosmic ionization rate nH/{\zeta} ratio. We also report the occurrence of naturally sustained oscillations for a limited but well defined range of control parameters. The period of oscillations is within the range of characteristic time scales of interstellar processes and could lead to spectacular resonances in time dependent models. Reservoir species (C, CO, NH3, ...) oscillation amplitudes are generally less than a factor two. However, these amplitudes reach a factor ten to thousand for low abundance species, e.g. HCN, ND3, that may play a key role for diagnostic purposes.The mechanism responsible for oscillations is tightly linked to the chemistry of nitrogen, and requires long chains of reactions such as found in multi-deuteration processes.


Introduction
The evolution of interstellar chemical abundances through formation/destruction chemical processes is described by a set of nonlinear differential equations that fall under the frame of dynamical systems. Then, it is not surprising to find some of the standard features characteristic of these systems. In particular multiple fixed-points for a given set of control parameters have been demonstrated by Pineau des Forets et al. (1992) and first analyzed by Le Bourlot et al. (1993Bourlot et al. ( , 1995a. The two stable solutions correspond to very different chemical signatures and have been labelled as LIP, corresponding to a Low Ionization Phase and HIP for an High Ionization Phase. LIP features are described by abundant saturated molecules and molecular ions whereas HIP signatures rather point to abundant radicals, carbon chains and atomic ions. Although these results led to some debate (e.g. Shalabiea & Greenberg 1995), they were later confirmed by various groups (Lee et al. 1998;Charnley & Markwick 2003;Boger & Sternberg 2006;Dufour & Charnley 2019). Some puzzling observations were even suggested to be representative of the HIP phase by, e.g., Gerin et al. (1997) or Ceccarelli et al. (2011).
However, apart from a passing mention in Le Bourlot et al. (1995a), no steady chemical oscillations have been reported to date. This is rather surprising, as such a behavior is a normal feature of nonlinear sets of chemical equations (see Gray & Scott 1994) and results probably from the limited range of integration time and parameter values considered in the current timedependent chemical studies applied to interstellar clouds (Wakelam et al. 2015). Various time dependent behaviors including oscillations are explicitly quoted in Maffucci et al. (2018). This paper will be discussed below, as we believe their study to be not fully complete. Furthermore, our analysis explains naturally a large fraction of their findings.
In this paper, we describe in some details naturally occurring oscillations in a model of dark cloud chemistry. The underlying formalism is presented in Sect. 2, a fiducial oscillatory solution is introduced in Sect. 3. We discuss the effect of various control parameters in Sect. 4 and present our analysis and a discussion in Sect. 5.

Time dependent chemical equations
We restrict ourself to an idealized dark cloud. That is, photo processes come only from secondary photons following the creation of high energy electrons by cosmic-ray. These electrons excite the different electronic states of H 2 that eventually decay by emitting ultraviolet photons in the 750 − 1750 Å window (Gredel et al. 1989).
We also make two very strong restrictions by assuming a constant temperature and discarding gas-grain chemical interactions. We are fully aware that, as a consequence, the resulting models cannot be readily applied to "real" astrophysical objects, nor be compared to observations. But this is a necessary preliminary step to try and understand the chemical mechanisms at work. Extended models that overcome one or both of these approximations will be introduced later (note however a brief discussion in the last section).
On the other hand, our chemical model includes an extended and detailed chemical network including 254 species, coupled by 4872 reactions taken from our previous studies (Roueff et al. Article number, page 1 of 13 arXiv:2010.01348v2 [astro-ph.GA] 22 Oct 2020 A&A proofs: manuscript no. AA_39085 2015; Cernicharo et al. 2018). In particular it includes a detailed deuteration scheme with multiple deuterated species from NH 3 to ND 3 , H 2 CO to D 2 CO, H 2 CS to D 2 CS (Marcelino et al. 2005), CH 4 to CD 4 . We introduce the value of the ortho-to-para ratio of molecular hydrogen, (O/P), which drives the deuteration efficiency through the H + 3 + HD H 2 D + + H 2 reaction (Pagani et al. 1992) as an external fixed parameter, following Roueff et al. (2015). This parameter is also critical for the initial step of the nitrogen hydrogenation reaction N + + H 2 → NH + + H which is slightly endothermic when H 2 is in its para form as first emphasized by Le Bourlot (1991). We follow the prescription of Dislaire et al. (2012b) based on the experimental studies of Marquette et al. (1988) to account for this reaction. We also include the experimental reaction rate coefficients to account for N + + HD and N + + D 2 reactions.
Parameters used for various models are displayed in Table 1. The range shown is typical of the values explored, although we did not compute all possible models. 1

Nondimensional equations
One major advantage of using a constant temperature is that it leads to constant chemical reaction rates once the other external parameters are chosen. In a low density gas only two body collisions are possible. So the time evolution of any species X i abundance n i = n (X i ) writes: Two bodies rates k i l,m are in cm 3 s −1 and one body rates k i j are in s −1 . These condensed equations account both for formation/destruction processes with appropriate signs. Given our hypothesis, all k i j are proportional to the molecular hydrogen cosmic-ray ionization rate ζ. This suggests to use 1/ζ as a characteristic timescale, and to introduce a nondimensional time variable τ = ζ t. Writing now n i = n H x i , with x i the relative abundance of species X i , we can introduce nondimensional rates for all chemical reactions by writing: k i l,m = ζ n H r i l,m ; k i j = ζ q i j . Subsequently, Eq. (1) reduces to: This analysis shows that we can lower the number of control parameters by defining an "ionization efficiency parameter" leads to values of I ep from 1 to 10 for typical interstellar conditions. This property has been previously recognized for steadystate results (Lepp & Dalgarno 1996;Boger & Sternberg 2006) or for the maxima corresponding to so-called early-time results (Lee et al. 1998). However, to our knowledge, this occurrence was never taken full advantage of for interstellar applications.
Let us also point out that this property is maintained for gasgrain interactions as long as a single grain temperature is considered and a standard rate equations formalism is used to describe surface chemistry (Hasegawa et al. 1992). 1 The code is made available at http://ism.obspm.fr.

Dynamical system properties
Equation (2) is in the standard form of a dynamical system: where X is the vector of (unknown) variables, α a set of (fixed) control parameters, and F a vector of functions. Below, we summarize the main properties of dynamical systems as described in the literature; for example Gray & Scott (1994). Fixed points of this system are found by setting F (X, α) = 0. The (possibly multiple) solutions may or may not be stable. The stability can be assessed by computing the eigenvalues of the Jacobian matrix J at the fixed point. IfX is a fixed point, then: Close toX, the system (2) can be linearized. Its evolution is then given by where λ i are the eigenvalues of J X . We see that, if all real parts (λ i ) < 0, then the fixed point is stable, and evolution may lead to it, provided that all control parameters remain constant during the evolution. If a single λ i has a positive real part, then the fixed point is unstable, and the system leaves its vicinity. As soon as X (τ) −X is no longer infinitely small, the linearized approximation does not apply and higher order terms control the evolution. This is independent of whether λ i eigenvalues have an imaginary part or not. If they do, then they occur in pairs of complex conjugates. Writing λ i = µ i + i ω i , all X (τ) include a contribution of the form exp (µ i τ) cos (ω i τ + φ i ), which is an oscillatory contribution.
-If the fixed point is stable, then the system behaves as a damped oscillator, with a damping characteristic time of where µ j is the smallest eigenvalue (in absolute value). If the associated λ i has an imaginary part, then there is an oscillatory damping. -If the fixed point is unstable, the system leaves its vicinity with a timescale associated to the most positive eigenvalue.
Here again, an oscillatory component may or may not be present.
The subsequent evolution depends on the topology of the phase space. Evolution may lead to another (stable) fixed point (steadystate), or to sustained oscillations if a stable limit cycle exists. In the following, we show how a stable fixed point undergoes a Hopf bifurcation, giving birth to a stable limit cycle at a specific value of a control parameter. This occurs when the real part of an eigenvalue of the Jacobian matrix crosses a value of zero. We note that this also explains the phenomenon of "critical slowdown" observed in the vicinity of a bifurcation point. As the evolution depends on exp (µ i τ), the damping or growing time diverges as µ i → 0. An example is given below in Sect. 4.2.

The role of initial conditions
If several attractors exist (e.g., one or more stable fixed points, and a stable limit cycle) then the evolution depends critically Hydrogen nucleus number density ζ 10 −17 s −1 1 − 300 Hydrogen cosmic-ray ionization rate T K 5 − 20 Gas temperature δ C -10 −5 − 10 −4 Carbon gas-phase abundance δ N -10 −5 − 10 −4 Nitrogen gas-phase abundance δ O -10 −5 − 10 −4 Oxygen gas-phase abundance δ S -8 10 −8 Sulfur gas-phase abundance on the chosen initial conditions. Different initial conditions may lead to one or the other attractor as previously shown in Le Bourlot et al. (1993). The initial part of the evolution has thus no particular meaning, depending solely on an arbitrary choice of initial conditions. Relaxation to the attractor depends mainly on its stability character and its distance -in parameter spacefrom a possible bifurcation point. We would like to stress here that the so-called early-time behavior is entirely dependent on an arbitrary choice of initial conditions; it should be considered with great care as it provides almost no information on the intrinsic behavior of the system. When application to interstellar medium conditions is examined, one should remember that interstellar clouds do not form out of nothing at a given time, but rather result from fluctuations within a turbulent medium that is permanently stirred by various dynamical process (e.g., supernovae or the differential rotation of the galaxy). Part of the gas disappears in new stars and fresh material is constantly added to the medium by the evolution of stars (stellar winds, loss of envelopes, supernova) or cooling flows from intergalactic material.
Thus, when a dense cloud separates from its environment, its composition is typical of the diffuse medium and quite different from the arbitrary simple values often found in the literature (e.g., H 2 , C + , O, N, ...). Model results obtained with such initial conditions may therefore be significantly misleading.

The oscillatory fiducial model
We report and discuss now our finding of oscillatory solutions obtained through a mixture of serendipity and systematic exploration of the model parameters displayed in Table 1. Table 2 shows the values of the physical parameters leading to an oscillatory behavior of the reduced abundances, and defines our fiducial model.

Typical oscillations
11 K 3 10 −5 6 10 −5 3 10 −5 8 10 −8 1.5 10 −8 Examples of oscillations are shown in Fig. 1. All quantities are nondimensional (x relative to n H and time in units of τ). Panel (a) shows how the reservoir of carbon oscillates from C to CO and back. Oscillations are simple and in perfect phase opposition. We note that the C/CO ratio is of the order of 1, which corresponds to unusual but observationally found situations with a significant neutral carbon abundance such as TMC1 (Schilke et al. 1995). The oscillatory behavior has indeed been found for a C/O elemental ratio of 1 corresponding to carbon-rich environments. We also point out that the electronic fraction displayed in panel (b) of Fig. 1 oscillates around 4 10 −7 , with an amplitude of ∼ 2 10 −7 , corresponding to a relatively high value for dense interstellar cloud conditions. The chemical properties of the corresponding solutions are similar to those of the HIP phase (Le Bourlot et al. 1993). Panel (c) shows how minor species react to this balance. NH 3 oscillation amplitudes are mild, but the full isotopic substituted ND 3 shows spectacular amplitude variations, which is the case for all fully deuterated substituted species. The maximum variation is over three orders of magnitude for CD + 5 . Panel (d) shows, for SH, that more complex oscillations are achievable that probably emerge from intricate coupling terms in the chemical evolution.
We can also display the dependence between different evolving variables in phase space. After dissipation of the transitory regime, the phase trajectory converges towards a limit cycle as displayed for NH 3 , CH 4 , and CO in Fig. 2.

Initial conditions effects
Most if not all time-dependent studies use initial conditions where carbon is totally ionized as C + , sulfur as S + , whereas oxygen and nitrogen are atomic. Hydrogen is most often taken as fully molecular but atomic initial conditions are also used sometimes. These assumptions are obviously far from reality as a dense cloud results from the previous evolution of diffuse clouds which are found to be chemically diverse from absorption observations in front of randomly available continuum radio sources; see for example Liszt et al. (2019). The partitioning of elemental carbon amongst the different potential reservoirs, C + , C, and CO is much more likely to represent plausible initial conditions. If some amount of CO is present as an initial condition, this impacts the sharing of oxygen as well in order to maintain the full elemental abundances. Figure 3 demonstrates this effect on the time evolution of the fractional abundance of a typical complex molecule, D 2 CO 2 . All initial abundances are identical except for C + , O, and CO. Initial C is kept to a vanishingly low value here, but qualitatively identical results are obtained if we use C in place of C + . For these, the total amount of carbon and oxygen is distributed with a fraction α on C + and O and a fraction 1 − α on CO (recall that C/O = 1 here). We see that D 2 CO oscillations are exactly the same after about two periods, exhibiting only a phase shift. The transitory regime occurring before the oscillations take place shows amplified response to initial conditions that are far from the asymptotic behavior. But that initial variation can either be towards a stronger peak (the so-called earlytime regime) or a very low abundance, with all intermediate values possible. We see that the resulting abundances are entirely a consequence of the choice of initial conditions. It is therefore very important to select physically plausible conditions, which excludes the usual choice of a single species found in theoretical papers. A reasonable solution is for example to run a first case until a permanent regime is reached, and then to change the control parameters (e.g., the density) and start from the solution of the previous run.

Timescales
The period of the oscillations can be deduced from the position of successive maxima. We find a period of τ 0 = 1.1562 10 −3 with I ep = 4, which translates to 0.73274 Myr for ζ = 5 10 −17 s −1 . To convert from a nondimensional period τ 0 to "real" time, one needs only to divide by the value of ζ in s −1 . If I ep is given, then: and so For a given value of I ep , the oscillation period decreases as 1/n H . Figure 4 shows the relation between n H and t osc for the range of I ep found in Sect. 4.1. Using initial conditions leading to a high first maximum of the time evolution (e.g., α = 1 in Fig. 3), we can also study how these maxima converge towards a steady-state behavior by plotting the differences between successive maxima ∆x and the asymptotic value. Figure 5 shows a typical example 3 . A clear exponential decrease emerges from which we deduce a characteristic damping time for the transitory regime, τ r 6.5 10 −4 , slightly above a half period.
The oscillations are also clearly visible when performing a Fourier transform of the abundances. The equations were integrated over 1400 periods, which allows us to reach a resolution of δν = 0.7736 (unit-less). The corresponding spectra are shown in Fig. 6, where the Fourier transform power P of x (C) and x (SH) evolutions are displayed for comparison.
Although very close, the peak intensities are slightly different at low frequency. First, the main peak at ν 0 = 864.93 confirms the value of the period of τ 0 = 1.1562 10 −3 reported above. All the harmonics at n × ν 0 are very strong. Additional weaker peaks arise at two other frequencies, ν 1 = 298.63 and ν 2 = 566.27. However, subsequent peaks at ν 2 + n × ν 0 get increasingly stronger when those at ν 1 + n × ν 0 disappear. This is particularly perceptible for SH. We also point out that ν 1 is close to ν 0 /3 and ν 2 is close to 2ν 0 /3 but the difference between these fractions is markedly greater than the frequency resolution, and is therefore significant. These additional peaks are therefore not sub-harmonics of the main frequency.

Dependence of the oscillatory behavior on the control parameters
Such a chemical evolution could affect the dynamics of the environment if its occurrence complies with likely observable conditions, as discussed in Sect. 5. Here, we study within which range of control parameters oscillations are obtained.

Effect of the ionization efficiency parameter I ep
The standard method to find the extent of an oscillatory zone in parameter space is the following: 1. Fix a value of the parameter I ep and find the exact value of T at the bifurcation. This should be achieved by computing the eigenvalues of the Jacobian matrix and solving for λ m (T ) = 0, where λ m is the highest eigenvalue real part. 2. Then use a continuation scheme to follow the bifurcation point in parameter space. To this end, we require that a variation δI ep of I ep leads to a variation δT of T to still be at a bifurcation point.
Then, to first order: If ∂λ m ∂T 0, we can solve the 1D ordinary differential equation (EDO): Article number, page 5 of 13 A&A proofs: manuscript no. AA_39085 Unfortunately, this scheme breaks down in the present case. The Jacobian matrix happens to be very ill-conditioned, and the computation of the eigenvalues becomes unstable close to the bifurcation. This has been confirmed using two different numerical routine libraries: the latest LAPACK version of DGEEVX, and the oldest but quadruple precision version of EISPACK routine "rg". Both routines agree very well on all eigenvalues, except in the vicinity of a bifurcation point where the smallest eigenvalues in norm are unreliable 4 .
We therefore resorted to a less accurate but robust method to explore parameter space: brute force. Figure 7 shows the domain of oscillations in the I p − T parameter space obtained from systematic sampling, together with the oscillation period represented by the color coding. The domain of oscillation spans a narrow range of temperatures between about 7 and 15 K with ionization efficiency parameters between 1.3 and 10.4. For a typical cosmic-ray ionization parameter ζ = 5 10 −17 s −1 , the range of density is 5 10 3 −5 10 4 cm −3 . We stress that all these values are compatible with physical conditions found in interstellar clouds. The oscillatory period varies by less than a factor of two, between 450 and 900 kyr for this value of cosmic-ray ionization rate. Such values are also similar to the expected evolution time of interstellar clouds. Small ionization efficiency parameters lead to larger values of the oscillation period. We can therefore consider that these findings are much more than a numerical curiosity and could significantly impact our knowledge of the interstellar medium.

Effect of the temperature
Keeping other parameters from Table 2 fixed, we now explore the range of temperatures over which oscillations are found. Figure 8 is a bifurcation diagram displaying the steady-state solutions when T is varied from 8.5 K to 13 K, keeping I ep = 4. A single fixed point is obtained for T c1 < 8.752 K and T c2 > 12.347 K, which is exhibited in green. Stable oscillations are present in between and the corresponding maxima and minima values of the fractional abundances are also reported in green. The fixed point itself does not disappear, but becomes unstable, and this is displayed in red in Fig. 8. This behavior is standard, 10 -1 8.5 9 9.5 10 10.5 11 11.5 12 12.5 13 τ r T (K) Fig. 9. Relaxation time of initial conditions. We highlight the large uncertainties far from the bifurcation points. Conversely, close to T c1 and T c2 the large relaxation time allows for a very accurate determination. and representative of a Hopf bifurcation, as detailed for example in Gray & Scott (1994). This effect of temperature is illustrated using HCN, but the behavior is generic and identical for all species. Below T c1 = 8.752 K, there is a stable fixed point. Starting from various (but not necessarily any) initial conditions, the system converges towards it. However, convergence is slower and slower corresponding to an increasing relaxation time of the initial conditions, as seen in Fig. 9, which displays the relaxation time deduced from the decline of successive maxima as a function of temperature.
This is the critical slow-down phenomenon, characteristic of bifurcations. Close to T c1 , the fixed point becomes unstable, and oscillations appear through a Hopf bifurcation. Their amplitudes are initially vanishingly small, but they grow fast as T is increased. The time evolution for two values of the temperature that bracket the bifurcation point is shown in Fig. 10. In the close vicinity of T c1 and for T < T c1 , a nice damped oscillatory structure is shown for the ND 3 fractional abundance, which converges slowly towards the stable fixed point (in green on Fig. 10). The red curve displays the sustained oscillations taking place immediately after T c1 . One can see that many oscillations are necessary before reaching the permanent regime.
The amplitude of the oscillations remains significant until the vicinity of T c2 = 12.347 K where a second Hopf bifurcation occurs. Beyond this second bifurcation point, the fixed point is stable again. A fixed point is still present between the two bifurcation points but the corresponding solution is unstable. These points are displayed in red in Fig. 8 and were computed using a Newton-Raphson scheme by directly solving the d/dt = 0 coupled chemical equations. The critical slow-down in the vicinity of a bifurcation point is important, because it leads to timescales that are totally unrelated to any other characteristic times of the system, either chemical or dynamical. Hence, the system may seem to behave strangely given the chosen parameters. This is simply the direct consequence of a close-by bifurcation point in parameter space. It is also remarkable that the range of temperatures where the oscillatory behavior takes place is within the values derived from interstellar cloud observations.
The value T c2 is driven by our hypothesis on (O/P), the H 2 ortho-to-para ratio. Here, we assume that it is the thermal equilibrium value, where T is the gas temperature. Increasing T increases (O/P), which has a dramatic effect on the N + + H 2 → NH + + H reaction rate coefficient (see discussion on the chemistry in Sect. 5.1 and Appendix A and on the (O/P) value in Sect. 5.3).

Effect of the elemental abundances
Elemental abundances in the gas phase constitute another set of control parameters. We focus our study on the most abundant elements after hydrogen and helium, namely C, N, and O, and introduce a relation between δ C and δ O . Instead of varying the latter independently, we explore the influence of the δ C /δ O ratio at constant total δ C + δ O (total metallicity) or the influence of metallicity for a constant δ C /δ O = 1. Figure 11 shows a bifurcation diagram for δ N . All other parameters are kept at the values of our fiducial model (cf Table 2). As the density n H is kept constant, we chose the ionization degree of the gas as a representative variable, as it does not require a rescaling by elemental abundances. Two Hopf bifurcations are present, with oscillations maintained up to a very high (and somewhat unrealistic) abundance value for nitrogen (over 7 10 −4 ).
Keeping δ N at its fiducial value (6 10 −5 ), Fig. 12 shows that oscillations encompass typical values of δ C /δ O and not only the selected value of 1 we used up to here. Figure 13 shows in a similar way that the total metallicity can be varied up to values found in many environments, although somewhat lower than in the solar vicinity.  Table 2.
We conclude that the occurrence of oscillations does not require very specific values of the elemental abundances, although our studies are more representative of low-metallicity or depleted environments.
The chance occurrence of a spectacular critical slow-down found during the exploration of this parameter space is presented in Appendix B.

Chemistry
Elucidating the chemical mechanisms responsible for oscillations is much more challenging than for steady-state results. There is indeed no set of equations to solve that would provide the answer. One has to follow the time variations of leading species to point to the relevant interactions. Here, we do not know which are the leading species. We can however derive some chemical properties that comply with the oscillatory behavior: 1. The full evolution remains within HIP chemistry. This corresponds to a significant electronic fraction (∼ a few 10 −7 ) which further implies the presence of C at a level close to that of CO. 2. Full deuteration appears compelling (within our set of chemical reactions). This implies that long chains of reactions exist leading from some very simple species (say, e.g., N + ) to fully deuterated species (say, e.g., ND 3 , ND + 4 ). It may be possible that other long chains of reactions also lead to oscillations, for example the formation of cyanopolyines, but we did not explore this point further. 3. Nitrogen chemistry plays a central role. The slightly endothermic N + + H 2 → NH + + H reaction, which is subject to large uncertainties (Dislaire et al. 2012b;Marquette et al. 1988), has a huge impact on the collective behavior.
A first clue as to the underlying mechanism is seen from the successive maxima of chemically linked species. Figure 14 shows the sequence of different deuterated compounds of NH 3 as a function of time. The abundances are rescaled to the (dimensionless) interval [0 : 1] through (n(t)−n min ) (n max −n min) to enable comparison. We first observe the peak of NH 3 itself. Then, offset by nearly half a period, we find the peak of NH 2 D, which is closely followed by that of ND 2 H, and finally by that of ND 3 . Meanwhile, NH 3 has time to decrease to its minimum, which quenches the deuteration process. The deuterated species are then destroyed more efficiently than they are formed, and they fall to very low levels. The process then starts again. Figure 15 shows the ratio of the maximum to minimum abundances of all species as a function of their peak abundance time with respect to that of NH 3 . The peak of NH 3 is well separated from all other species (by at least one-tenth of a period), which is a strong indication that it plays a particular role. This latter peak is followed by a small number of "lead" molecules, half of which are hydrogenated nitrogen with no or one deuterium. Then, at about 0.6 τ 0 , most species peak, with an amplitude that may exceed a factor of 10 4 . This is where all multiply deuterated molecules are found, with the single exception of D + 2 . This peak extends over a limited time interval, which expresses the trigger of successive interconnected maxima. We interpret these features as a signature of the importance of nitrogen deuteration in the mechanism responsible for oscillations.

NH 3 /CN connection
Formation of NH 3 and its deuterated substitutes proceeds progressively in time as shown in Fig. 14. The initial step is provided by the reaction of N + with H 2 , HD, and D 2 as shown in Fig. 16.
In the HIP environment considered here, atomic carbon C Hickson et al. 2015) and ionized carbon C + are the main destruction agents (see Fig. 17) which finally Others multi-D Lead N Lead others Fig. 15. Ratio of maximum to minimum abundances as a function of phase shift with respect to NH 3 . "multi-D" are molecules with two or more D. The two vertical lines delineate the time range where multideuterated species peak.
lead to the formation of a CN bond within a large variety of species (HCN, HNC, HCNH + ,..), as displayed in Fig. 18. Focusing on these H, C, and N compounds as displayed on the left of Fig. 18, we see that the CN radical emerges. Reactions of CN with C and O then lead back to atomic N whereas the reaction with N produces N 2 as found experimentally by Daranlot et al. (2012). N + is subsequently restored via reactions of N 2 with He + and cosmic-ray ionization of atomic N, and the full cycle starts again.
During the whole process, minor reactions may create leaks towards other branches of the full chemical set (e.g., through NO for a coupling to oxygen chemistry), and so it is very sensitive, which explains why we find oscillations only within a restricted part of parameter space. We further study nitrogen chemistry in Appendix A.

Thermal balance
Our current model does not include a sophisticated treatment of thermal balance, as found for example in the Meudon PDR code (Le Petit et al. 2006). However, simple analytic formulae of the main heating and cooling processes are included, following Flower & Pineau des Forêts (2015), that allow the derivation of the temperature within a given equation of state,  Fig. 16. NH 3 and deuterated isotopolog formation. namely n = Cte here, with n the density of the gas, that is, n = n (H) + n (H 2 ) + n (He) + ....
Solving for thermal balance removes one control parameter from the system, as temperature is now computed by solving its evolution equation and is not set arbitrarily. However, the cooling and heating terms do not scale in the same way with ζ and n H , and so the ionization efficiency I ep cannot be used as a single parameter. The number of control parameters is therefore unchanged.
We did not conduct a full exploration of that parameter space, but we checked that thermal oscillations can be found and are fully compatible with the mechanism described here. For I ep = 4, Fig. 19 shows a Hopf bifurcation close to n H = 9.5 10 3 cm −3 . We kept a constant value for I ep , to help comparisons with Fig. 7. One can see that the range of temperatures spanned within the oscillations is compatible with the values found at constant temperature.
This example shows that oscillations can be expected within a detailed model using a more sophisticated thermal balance. However, it is not yet possible to determine the precise domain of parameters involved.  Fig. 18. CN links destruction and nitrogen recycling. However, this latter process is very poorly understood. (O/P) is usually supposed to be three at the formation process on grains, but this is true only if a significant fraction of the formation enthalpy is used to that end. Some mechanisms may also lead to a significant transfer of energy to the grain, or favor a different statistical equilibrium. For example, a strong interaction with the surface could lead to independent spins for the two protons, resulting in a (O/P) value of one at formation. However, our choice is certainly a lower limit, and higher values are expected, as found by Flower et al. (2006). The thermal value of (O/P) for the fiducial model at T = 11 K is 1.67 10 −6 . We explored the influence of using a higher value of (O/P) and found that oscillations persist up to (O/P) = 1.65 10 −4 in this case. We therefore find that oscillations are entirely compatible with typical values of (O/P).
We checked that the impact of the fraction of ortho-H 2 on the occurrence of oscillations is entirely due to N + +H 2 → NH + +H. The other (O/P)-sensitive reaction, namely H 2 D + + H 2 → H + 3 + H, has no impact on whether or not oscillations occur.

Grain chemistry
The present study is restricted to pure gas-phase chemistry. However, recombination of atomic ions on the grain surfaces reduces the ionization fraction as emphasized by Shalabiea & Greenberg (1995) in their discussion of the effects of grains on bistability properties. We recall that Le Bourlot et al. (1995b) showed that the bistability domain is shifted, but not suppressed, when only neutralization on grains is introduced.
Following the treatment reported in Le Bourlot et al. (1995b) and summarized in Cuppen et al. (2017), we obtain where a is the grain radius, v X + is the velocity of X + , f C is the Coulomb factor 2 e 2 3 k T a , and x gr is the dust-to-gas density ratio. For a typical grain radius a of 0.1 µ, a dust-to-gas mass ratio G = 10 −2 , a volumic mass density ρ of 3 g cm −3 , and a temperature T = 10 K, (1 + f C ) = 12, x gr = 3 × 1.4 m H G 4π ρ a 3 = 1.85 10 −12 and the equation reduces to: where M X + is the atomic mass of X + in amu. We find that such a value is compatible with the oscillatory behavior. In addition, the neutralization rates are reduced significantly in low-metallicity regions where the dust-to-gas mass ratio is small, and in dense cores where coagulation leads to a large mean grain radius 5 . 5.5. Discussion of the oscillatory behavior reported in Maffucci et al. (2018) In a paper devoted to cyanopolyynes and other oxygenated complex organic molecules in TMC1, Maffucci et al. (2018) mention oscillatory solutions for their time-dependent model. Unfortunately, the authors do not extend their calculations further than two millions years. Thus, one cannot be sure that these oscillations are maintained over a large period of time. There are nevertheless striking similarities between our findings: -The oscillatory behavior takes place for specific physical conditions, corresponding to a large C/O ratio, low densities, and a large cosmic ionization rate, which probably would correspond to a chemical HIP phase, with large abundance of carbon chains. -The largest amplitude of the oscillations is obtained for the large cyanopolyynes whereas we find a similar feature for the highly substituted deuterated species. These species correspond to the extremities of successive complex chemical reactions.
These authors did not try to further explore such a behavior, and the origin of these oscillations could be more thoroughly investigated. However, we note that this model includes grain surface chemistry, which shows that this is not restricted to a pure gasphase chemical network. We also anticipate that the oscillatory behavior concerns all chemical species, as this is a basic mathematical property of the chemical equations, as mentioned above. The observed fluctuations of long chains in TMC1 could come from oscillations triggered with different initial conditions from place to place. Therefore, variations in abundances would result from a common origin and common environment. We also note that these latter authors use initial conditions with all abundances set to zero except for a single atom or atomic ion for each element (except for hydrogen). This is not likely to happen in a real cloud that condensates from diffuse conditions. Hence, their "early-time" results are entirely a consequence of this unfortunate choice and tell us nothing about the age of the cloud (see Sect. 3.2 above). In particular, we can expect that starting from conditions where most carbon is locked in CO would quench the formation of long cyanopolyyne chains entirely if oscillations do not restore some C to the gas phase.

Summary and conclusions
The present study reveals for the first time a well-known feature of nonlinear effects of chemical kinetics in the context of interstellar chemistry, namely the occurrence of oscillatory solutions of the chemical abundances over a restricted but representative 5 The overall dependence of the neutralization rate coefficient on grains is inversely proportional to a, the grain radius. physical parameter space. The oscillations are still present in an isochoric model where the thermal balance is solved concomitantly with the chemistry. We obtained different bifurcation diagrams by varying several control parameters, such as n H /ζ, temperature, and elemental depletions. The spectacular outcome, in our point of view, lies in the range of physical parameters where the oscillatory behavior is obtained, which corresponds to the values derived from interstellar observations in terms of density, temperature, cosmic ionization rate, and elemental abundances. In addition, the timescales obtained for the oscillation periods are between some hundred thousand years and some million years, which are also within the same order of magnitude as other important timescale values relevant to the interstellar medium, such as the free-fall time or the crossing time (Klessen & Glover 2016).
The oscillatory behavior occurs within the so-called HIP phase, corresponding to electronic fractions which are typically of the order of several 10 −7 , significantly larger than those corresponding to the standard LIP chemical cloud conditions. These conditions are favored by a relatively high C/O elemental ratio which in turn leads to a significant number of carbon chains, as found in several interstellar conditions.
The triggering mechanism of the oscillations is not precisely identified but the detailed analysis of the differential behavior of the various species points to the role of long chains of reactions. In our study, these are provided by the sequential deuteration processes which end up with fully deuterated species, such as CD 4 and ND 3 , and the corresponding molecular ions CD + 5 and ND + 4 . The amplitude of the oscillations increases dramatically with the number of substituted deuterium nuclei. We also find that nitrogen chemistry plays a critical role through the link between carbon and oxygen. Finally, we want to emphasize the crucial importance of realistic initial conditions for deriving sensible conclusions on early-time behaviors of the chemical abundances.
We can expect that the (magneto-)hydrodynamical behavior of a fluid subject to chemical oscillations could be quite distinct from that of a more inert fluid. Various timescales might couple and lead to unexpected consequences. In particular, we propose that some characteristic length scales could emerge, which would be a natural explanation for variations in the star formation rate and in the mass distribution of young stars in a giant molecular complex, depending on whether or not chemical oscillations, and the induced oscillation of the cooling function, are triggered during the collapse.