Eugenio Cinquemani
email: eugenio.cinquemani@inria.fr

Inference of the statistics of a modulated promoter process from population snapshot gene expression data

Keywords: Reporter gene systems, Controlled Markov chains, Regularized estimation, Splines

In previous work, we have developed mathematical tools for the analysis of single-cell gene expression data from population snapshots, and an inference algorithm for the estimation of stationary statistics of promoter activation. In this work, we address the inference problem in the nonstationary case of modulated processes. This is of special relevance to control scenarios, where an exogenous input modulates the time evolution of promoter activation. We provide an effective method for the computation of the output statistics of a reaction network with a nonstationary, causal input process of modulated form. Based on this we devise and demonstrate an algorithm for the reconstruction of the promoter (input) process statistics from snapshot data.

INTRODUCTION

Noise in gene expression is at the roots of important phenomena in cell life and evolution, such as bet-hedging and the emergence of specific gene regulatory patterns [START_REF] Raj | Nature, nurture, or chance: Stochastic gene expression and its consequences[END_REF][START_REF] Rao | Control, exploitation and tolerance of intracellular noise[END_REF]. Its characterization is thus key for understanding cellular dynamics, and it has been the object of intense investigation [START_REF] Thattai | Intrinsic noise in gene regulatory networks[END_REF][START_REF] Paulsson | Models of stochastic gene expression[END_REF][START_REF] Swain | Intrinsic and extrinsic contributions to stochasticity in gene expression[END_REF][START_REF] Kaern | Stochasticity in gene expression: From theories to phenotypes[END_REF] . Among the frontiers of systems and synthetic biology is control of gene expression [START_REF] Uhlendorf | Long-term model predictive control of gene expression at the population and single-cell levels[END_REF][START_REF] Chait | Shaping bacterial population behavior through computer-interfaced control of individual cells[END_REF][START_REF] Milias-Argeitis | Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth[END_REF]. A primary objective of control is the exploration of cellular dynamics. In particular, control enables a deeper investigation of gene expression noise. This calls for the development of statistical analysis and inference methods [START_REF] Komorowski | Bayesian inference of biochemical kinetic parameters using the linear noise approximation[END_REF][START_REF] Munsky | Listening to the noise: Random fluctuations reveal gene network parameters[END_REF][START_REF] Neuert | Systematic identification of signal-activated stochastic gene regulation[END_REF][START_REF] Ocone | Reconstructing gene regulatory dynamics from high-dimensional single-cell snapshot data[END_REF] where the presence of control is explicitly accounted for.

In previous work [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF], we have addressed inference of promoter activity statistics from gene expression population-snapshot data, that is, measurements of the statistical distribution of the gene expression product in samples from the cellular population collected at different times [START_REF] Hasenauer | Identification of models of heterogeneous cell populations from population snapshot data[END_REF]. We framed the problem in the context of (first-order) stochastic reaction networks with a causal input on reaction rates, and developed Generalized Moment Equations (GMEs), a marginalized form of the well-known Moment Equations (MEs, see e.g. [START_REF] Lestas | Noise in gene regulatory networks[END_REF]) that describes the statistics of the network state process as a function of the input process statistics. We then used the GMEs to devise a method to reconstruct, in particular, the autocovariance function of an arbitrary This work was funded in part by the French national research agency (ANR) via project MEMIP (ANR-16-CE33-0018), and by the Inria IPL CoSy. The author thanks Jakob Ruess for useful discussions.

stationary promoter process. Unfortunately, the method does not apply to controlled gene expression processes, due to the typical nonstationarity of the control action.

In this work, we propose a generalization of the methods in [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF] to modulated input processes [START_REF] Zhao | Inference for modulated stationary processes[END_REF]. With this terminology we refer to stationary processes reshaped by a time-varying signal. In their matrix formulation, modulated processes account for a variety of scenarios of special relevance to gene expression, where a control input acts as an additive or multipicative factor on random gene expression noise (see e.g. [START_REF] Munsky | Listening to the noise: Random fluctuations reveal gene network parameters[END_REF]; [START_REF] Uhlendorf | Long-term model predictive control of gene expression at the population and single-cell levels[END_REF][START_REF] Milias-Argeitis | Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth[END_REF]). Based on a very flexible characterization of modulating functions via splines, we develop a strategy for the efficient numerical solution of GMEs in presence of modulated inputs. Then, we extend the inference method in [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF] to the estimation of modulated process statistics, and discuss reconstruction performance and the role of the modulating input via numerical simulations.

GMEs are reviewed in Sec. 2. The inference problem and method developed for stationary promoter processes in discussed is Sec. 3. Modulated input processes as well as an effective solution of the corresponding GMEs are presented in Sec. 4. The novel, generalized inference method is then presented in the same section. Numerical simulations demonstrating the method are in Sec. 5. Finally, conclulsions are drawn in Sec. 6. Notation: N, Z, R and R + denote natural, integer, real and nonnegative real numbers, respectively. For a set T ⊂ R, 1 T (•) is the indicator function of T , and

1(•) is the unit step function 1 [0,+∞) (•). || • || denotes Euclidean norm. For two random vectors X and Y , E[X] denotes expectation of X, Cov(X, Y ) = E[(X -E[X])(Y -E[Y ]
) T ] (superscript " T " denoting transposition) and Var(X) = Cov(X, X). P[ • ] denotes probability.

GENERALIZED MOMENT EQUATIONS

A reaction network is a family of n ∈ N chemical species and m ∈ N reactions that may occur among them in a given reaction volume. At a time t, let X i (t), with i = 1, . . . , n, be the number of molecules of the ith species. Let S j , with j = 1, . . . , m, be the column vector whose ith row is the net change in X i when reaction j occurs.

S = [S 1 • • • S m ] ∈ Z n×m is the network stoichiometry matrix. Assuming a well-stirred reaction volume, X = [X 1 • • • X n ]
T can be modelled as a Markov process, with

P[X(t + dt) = x + S j |X(t) = x] = w j (x, t)dt + o(dt),
where the functional form of the reaction rates (propensities

) w(x, t) = [w 1 (x, t) • • • w m (x, t)]
T is typically fixed by the mass-action laws [START_REF] Gillespie | A rigorous derivation of the chemical master equation[END_REF], and o(dt) tends to 0 faster than dt. We consider the case where rates are affine functions of the state, w(X, t) = W X(t) + F (t), (1) with W ∈ R m×n + a constant matrix and F (t) some exogenous input taking values in R m + . This captures zerothand first-order reactions and suffices to describe many biochemical processes of interest (notably gene expression) in an exact or approximate manner [START_REF] Kaern | Stochasticity in gene expression: From theories to phenotypes[END_REF][START_REF] Thattai | Intrinsic noise in gene regulatory networks[END_REF]. Depending on the context, F may represent environmental perturbations, intracellular regulation, or control.

Let µ(t) = E[X(t)], Σ(t) = Var X(t) . Given S, functions w(•, •) and a deterministic function F (•), the time dynamics of µ and Σ (and of higher-order moments) are described by a well-known set of Ordinary Differential Equations (ODEs) called the Moment Equations (MEs), see e.g. [START_REF] Lestas | Noise in gene regulatory networks[END_REF]. For the more general case considered in this paper, where F (t) is a stochastic (vector) process with well-defined second-order moments, let us further define µ F (t) = E[F (t)], ρ F (t, τ ) = Cov F (t), F (τ ) , ξ F (t) = Cov X(0), F (t) . In [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF], we proved in particular the following result. Proposition 1. Assume that first-and second-order moments of F are uniformly bounded. Then, for t ≥ 0,

μ(t) = SW µ(t) + Sµ F (t), (2) Σ(t) = SW Σ(t) + Σ(t)W T S T + Q(t)+ (3) V ξ (t) + V T ξ (t) + V ρ (t) + V T ρ (t), where, denoting (t) = exp(SW t)1(t), Q(t) = Sdiag W µ(t) + µ F (t) S T , (4)
V ξ (t) = Sξ F (t) T (t) T , (5)
V ρ (t) = +∞ 0 dτ Sρ F (t, τ )S T (t -τ ) T . (6)
We call Eq. ( 2)-(3) the Generalized Moment Equations (GMEs). For F a deterministic profile (that is, F = µ F , ρ F = 0, ξ F = 0), they reduce to the standard MEs. Note that the result holds irrespective of F being Markovian or stationary, however, a form of stochastic causality between the input process F and the state process X is subsumed by the definition of the rate functions (1) (see [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF]). Crucially, GMEs are linear in the input statistics µ F , ρ F and ξ F . As discussed in the sequel, this fact enables the development of efficient methods for practical applications.

INFERENCE OF PROMOTER STATISTICS

At the single-cell level, gene expression is well described by the so-called random telegraph model, i.e. the reaction network defined by the four stochastic reactions [START_REF] Kaern | Stochasticity in gene expression: From theories to phenotypes[END_REF][START_REF] Paulsson | Models of stochastic gene expression[END_REF], where M and P denote mRNA and protein species, respectively, while θ = (k M , d M , k P , d P ) are positive rate parameters. The binary random process U describes the state of the promoter that controls gene expression. Let X 1 and X 2 be the number of mRNA and protein molecules, respectively. Ordering reactions ( 7) from left to right, then top to bottom,

∅ k M •U ---→ M, M d M --→ ∅, M k P --→ M + P, P d P --→ ∅ (7)
S = 1 -1 0 0 0 0 1 -1 .
Moreover, the rate equations are as in (1), with

W =    0 0 d M 0 k P 0 0 d P    , F (t) = KU (t), K =    k M 0 0 0    .
The random outcomes of the joint process (U, X) may be thought of as different gene expression time profiles in different cells of a genetically homogeneous population. Randomness of U accounts in particular for different activation profiles and can be interpreted as extrinsic noise [START_REF] Swain | Intrinsic and extrinsic contributions to stochasticity in gene expression[END_REF].

In order to study the expression of a gene of interest, it is possible to engineer cells such that P is a fluorescent protein (see e.g. de Jong et al. ( 2010)). Then, individual-cell fluorescent protein levels X 2 (t) can be quantified by means of videomicroscopy or flow-cytometry. In particular, this allows one to get population snapshots of gene expression distribution at different times. We focus on the empirical mean and variance of gene expression snapshots, which provide noisy measurements of the mean and variance of X 2 . Denoting these measurements with μk and σ2 k , for a sequence of measurement times

t k with k = 1, . . . , M , μk = Cµ(t k ) + e µ k , (8)
σ2 k = C T Σ(t k )C + e σ k , (9)
where C = [0 1] T selects from µ and Σ the mean and variance of the protein count X 2 (the proportionality factor relating X 2 with observed fluorescence is ignored for simplicity). Measurement errors e µ k and e σ k can be statistically characterized [START_REF] Zechner | Moment-based inference predicts bimodality in transient gene expression[END_REF].

Inference of the statistics of a stationary process

On the basis of the gene expression model above (with known parameters θ), in Cinquemani (2019), a method for the reconstruction of the second-order statistics of the promoter state process U from snapshot data was developed for the case where U is stationary. In particular, under a few assumptions (notably X(0) = 0, as met by appropriate experiments, which implies V ξ = 0 in (3)) the method leverages the linearity of the GMEs (2)-(3) in µ F and ρ F to infer the (scalar-valued) stationary autocovariance function ρU (•) = Cov U (t + •)U (t) from measurements (9), without any mechanistic model for U . The interest of this problem is that the autocovariance function carries information about hidden regulatory processes, and its reconstruction from data is thus an important first step in characterization of gene expression. Notice that µ

F (t) = K μU , ρ F (t + δ, t) = K ρU (δ)K T , (10
) so that linearity of ( 2)-(3) in µ F and ρ F implies their linearity in μU and ρU . Assuming identification of the unknown constant μU from the mean data (8) and Eq. ( 2) in a preliminary step, the method seeks the autocovariance function ρU such that the corresponding solution of (3) best fits the measurements (9) in a Least-Squares (LS) sense. This is done via a finite expansion of the unknown autocovariance function. Let R(δ) = [r 1 (δ), . . . , r N (δ)] T be N appropriately chosen functions such that the convex cone generated by R approximates well the space of stationary autocovariance functions. Then, for some vector

of nonnegative coefficients c = [c 1 , . . . , c N ] T , the equation ρU (δ) = c 1 r 1 (δ) + . . . + c N r N (δ) = R(δ) • c
(11) holds with good approximation at all time lags δ. For inference purposes, R is fixed a priori, and c is the quantity to be determined from the data. For any t ≥ 0, let Σ(t|ρ U , μU ) denote the solution Σ(t) of the GMEs (2)-(3), with µ F and ρ F as in (10), as a function of μU and ρU . Define

V (•) = [v 1 (•), . . . , v N (•)] where, for every l, v l (•) = C T Σ(• |r l , 0)C, and v 0 (•) = C T Σ(• |0, μU )C. Using (11), in view of the linearity of the GMEs, it holds that C T Σ(t|ρ U , μU )C = v 0 + V (t) • c.
Therefore, the LS estimate ĉ of c from data (9) can be written as

inf c∈R N + M k=1 α 2 k σ2 Y (t k ) -v 0 (t k ) -V (t k )c 2 , (12)
where the weights α k = std(e σ k ) -1 are known [START_REF] Zechner | Moment-based inference predicts bimodality in transient gene expression[END_REF]. The estimate of ρU is then defined as ρU = R • ĉ. In [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF], additional constraints and a regularization term are introduced to ensure that the problem is well-posed [START_REF] De Nicolao | Nonparametric input estimation in physiological systems: Problems, methods, and case studies[END_REF] and that ρU is a well-defined autocovariance function, while preserving convexity of the optimization problem.

The viability of the method crucially depends on the ability to compute the transforms V of the approximating functions R. For a convenient choice of univariate, scalarvalued functions r l (notably, indicator functions), the integral term V r l (t) (see Eq. ( 6)) can be written out as an explicit function of t. In turn, for every l, this allows one to calculate the basis function transform v l (t k ) at all times t k by standard numerical integration of (3) (with r l in place of ρU ), leading to an efficient numerical implementation of the inference method (see [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF] for details).

EXTENSION TO MODULATED PROCESSES

We now want to pursue the efficient solution of the GMEs and a method for reconstruction of input process statistics from snapshot data in the case where the input process F is subjected to a control signal. We look at a specific class of nonstationary processes, which we refer to as modulated processes [START_REF] Zhao | Inference for modulated stationary processes[END_REF]:

F (t) = G(t)E(t), (13)
where E(t) is a second-order stationary (vector) process taking values in R q and G(t) is a deterministic (matrix) function of time taking values in R m×q . If μE and ρE (•) are respectively the stationary mean and stationary autocovariance function of E, then

µ F (t) = G(t)μ E , ρ F (t, τ ) = G(t)ρ E (t -τ )G(τ ) T . (14)
This case is of relevance, in particular, for gene expression models. Indeed, thanks to its general matrix-vector form, Eq. ( 13) includes as special cases additional and multiplicative forms of promoter regulation and control, acting at the level of basal expression rate or in terms of expression strength (see e.g. [START_REF] Chait | Shaping bacterial population behavior through computer-interfaced control of individual cells[END_REF]; [START_REF] Berthoumieux | Shared control of gene expression in bacteria by transcription factors and global physiology of the cell[END_REF]; [START_REF] Fiore | In vivo real-time control of gene expression: A comparative analysis of feedback control strategies in yeast[END_REF]).

Here, we limit ourselves to the case where the control signal G(t) is known. For simplicity we focus on the case where E is a scalar process (q = 1). The approach we take is analogous to the one reviewed in the previous section. For a given G(t), we aim at reconstructing the stationary autocovariance ρE from variance measurements (9) by a LS approach, leveraging the GMEs (2)-( 3) and ( 14). (As above, we assume that the constant μE is identified from mean data in a preliminary step). From an estimate of ρE , an estimate of the nonstationary ρ F follows via (14).

Given functions R for the finite expansion of ρE , the main challenge in the new scenario is the computation of the transformed functions V . Formally, these are still given by v l = C T Σ(• |r l , 0)C, where, in the light of ( 14), Σ(• |r, m) now expresses the solution of ( 2)-(3) relative to a hypothetical mean m and autocovariance function r of E. In practice, though, calculating Σ(• |r, m) by numerical integration of ( 2)-( 3) requires evauation at all times t of the integral term (6), which is now

V ρ (t) = +∞ 0 dτ SG(t)ρ E (t -τ )G(τ ) T S T (t -τ ) T . (15)
For a generic G(t), this term cannot be written as an explicit function of t. To solve the GMEs, the idea is then to express V ρ itself as the solution of a system of differential equations to be numerically solved jointly with (2)-(3). For convenient classes of functions G(t) (possibly used as approximations of more general functions) this is the subject of the section that follows.

Computation of the approximating function transforms

For a modulated process (13), eliminating τ from (15) by the change of variable δ = t -τ , one finds that, for a generic function r in place of ρE ,

V r (t) = SG(t) t 0 dδr(δ)G(t -δ) T S T (δ) T = SG(t)H 0,r (t), ( 16
)
where H 0,r (t) is defined as the integral in the first line.

We aim at giving a differential expression to H 0,r (t). For later use, for indices i such that the ith derivative

G (i) (•) of G(•) exists, define H i,r (t) = t 0 dδr(δ)G (i) (t -δ) T S T (δ) T , (17)
which is also consistent with the above definition of H 0,r . For ease of notation, we will now temporarily write H i in place of H i,r .

Proposition 2. Suppose that G(t) is a matrix of polynomials of order up to d -1. Then, for t ≥ 0, Ḣi (t) = r(t)G (i) (0) T S T (t) T + H i+1 (t), H i (0) = 0, (18) with i = 0, 1, . . . , d -2, and Ḣd-1 (t) = r(t)G (d-1) (0) T S T (t) T , H d-1 (0) = 0. (19)

Proof: For any i, the derivative of ( 17) is

[r(δ)G (i) (t - δ) T S T (δ) T ] δ=t + t 0 dδ d dt [r(δ)G (i) (t -δ) T S T (δ) T ]
, which is (18). For i = d -1, since G (d) (t) = 0 for all t, the above gives ( 19). The initial conditions H i (0) = 0 are an immediate consequence of (17).

Eq. ( 18)-( 19) constitute an ODE system for H 0 (t), as desired. In practice, though, the utility of this result is limited. In general G(t) may not be made of polynomials, and a polynomial approximation of a function over an interval generally requires polynomials of high degree. This would lead to a blow-up of the differential system size, which increases with d. It turns out that the above result can be extended to splines, a class of piecewise polynomial functions that are well suited to function approximation. In particular, suppose that the entries of G are piecewise polynomials of order up to d -1, with continuous derivatives up to G (d-2) [START_REF] Wahba | Spline models for observational data[END_REF]. Let the knots (points of junction of the different polynomials) be placed at T 1 , . . . , T p , with T 1 < T 2 < . . . < T p (without loss of generality, we assume all entries of G to have equally placed knots). Notice that, for some constant matrices G j ,

G (d-1) (τ ) =    G 0 , τ < T 1 , G j , τ ∈ [T j , T j+1 ), j = 1, . . . , p -1, G p , τ ≥ T p .
(20) Proposition 3. Suppose that the entries of G(t) are piecewise polynomials as defined above. Then, for t ≥ 0 and i = 0, 1, . . . , d -2, Eq. ( 18) holds. Moreover,

Ḣd-1 (t) = H + 0 (t) + p-1 j=1 H + j (t) -H - j-1 (t) , (21)
with H d-1 (0) = 0, where, for all relevant j, H

+ j (t) = r(t -T j )G T j S T (t -T j ) T , H - j (t) = r(t -T j+1 )G T j S T (t -T j+1 ) T . ( 22
)
Proof: The validity of Eq. ( 18) for i = 0, 1, . . . , d -2 is proven as in Prop. 2, also in view of the continuity of the derivatives of G up to G (d-2) . Next, using Eq. ( 20), simple manipulations of Eq. ( 17) allow one to expand the expression of H d-1 (t) as

t t-T1 dδr(δ)G T 0 S T (δ) T + p-1 j=1 t-Tj t-Tj+1 dδr(δ)G T j S T (δ) T + t-Tp 0 dδr(δ)G T p S T (δ) T
(since (δ) = 0 for δ < 0, this holds for any value of t).

Taking the derivative in t term by term yields

H + 0 (t) -H - 0 (t) + p-1 j=1 H + j (t) -H - j (t) + H + p (t).
Eq. ( 21) is a straightforward rearrangement of the above.

Notice that the generic jth term of the summation in ( 21) is nonzero only for t > T j . The generalization of ( 21) to an unbounded number of knots over an infinite time horizon is thus straightforward. In conclusion, for an arbitrary (measurable) function r, Eq. ( 18) (with i = 0, . . . , d -2) and ( 21) jointly give a system of ODEs for the computation of H 0,r (•). For any time t, V r (t) then follows from (16).

Inference method

We are now ready to state our method for the inference of ρE and ρ F from measurements (9), for a known function G. For the sake of this paper, we assume that E is a scalar (stationary) process, and

R = [r 1 • • • r N ] is a collection of indicator functions, though generalizations are immediate.
As in [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF], the method assumes X(0) = 0. Parameters S and W are assumed known and, as explained above, μE is identified in advance.

(i) (Computation of the v l ) Perform numerical integra- tion of μ(t) = SW µ(t) + SG(t)μ E , Σ(t) = SW Σ(t) + Σ(t)W T S T + Sdiag W µ(t) + G(t)μ E S T
, with µ(0) = 0 and Σ(0) = 0 to get the solution Σ(t k |0, μE ) at all measurement times, and set v 0 (t k ) = C T Σ(t k |0, μE )C, with k = 1, . . . , M . Next, for every function r l , with l = 1, . . . , N , perform numerical integration of the augmented ODE system

Σ(t) = SW Σ(t) + Σ(t)W T S T + SG(t)H 0,r l (t) + H 0,r l (t) T G(t) T S T , Ḣ0,r l (t) = r l (t)G (0) (0) T S T (t) T + H 1,r l (t), . . . Ḣd-2,r l (t) = r l (t)G (d-2) (0) T S T (t) T + H d-1,r l (t), Ḣd-1,r l (t) = H + 0,r l (t) + p-1 j=1 H + j,r l (t) -H - j-1,r l (t) ,
with H 0,r l (0) = . . . = H d-1,r l (0) = 0 and Σ(0) = 0, to get the solution Σ(t k |r l , 0) at all measurement times. Set 1) of U are as in [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF]. To explore the role of the (known) control input g in the estimation problem, we consider two different modulating signals, G and G . For G we take g(t) = (1 + γ • t) -1 , while for G we take g(t) = 1 -(1 + γ • t) -1 (γ = 0.02). These control inputs are qualitatively different: G is a monotonically decreasing signal with unit value at zero, while G is monotonically increasing and null at the origin (from now on, with some abuse of terminology, we refer to G and G instead of their first entry as the control inputs). Toward the application of the methods of Sec. 4, we approximate G and G by cubic splines (d = 3), with knots at T j = (j -1) • T , with T = 1, j = 1, . . . , p and p = 100. Plots of (the spline approximations of) G and G ) are visualized in the insets of Fig. 1.

V (t k ) = [v 1 (t k ), . . . , v N (t k )], where v l (t k ) = C T Σ(t k |r l , 0)C, with l = 1, . . . , N and k = 1, . . . , M . (ii) (Estimation of c) Compute ĉ as a solution to the convex optimization problem min c∈R N M k=1 α 2 k σ2 Y (t k ) -v 0 (t k ) -V (t k )c 2 + λ • c T Qc s.t. T (
For this stochastic gene expression system, we simulate data (8) and ( 9) analogous to [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF]. That is, we consider empirical mean and variance measurements collected from samples of 10 6 cells at times t k = 5 • k, with k = 1, . . . M and M = 20, with initial state X(0) = 0. For inference, we assume that the gene expression parameters θ are known, and so is the control input. Additionally, given the simplicity of estimating μE from measurements (8), we assume the latter to be known, and focus on testing the reconstruction of autocovariance functions from measurements (9). For both G and G , we simulate 100 datasets, and draw estimates ρE and ρF for every dataset with the method of Sec. 4.2. Approximating (indicator) functions R, regularization parameter λ and roughness penalty terms are defined as in [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF]. The statistics of the estimation results are reported in Fig. 1 and Fig. 2, for a fixed, educated choice of λ = 10 6 as well as for λ chosen from the data in every estimation run. The whole simulation was implemented and tested in Matlab on a 3GHz Intel Xeon (Ubuntu) laptop, using standard ODE solvers and CVX (CVX Research, Inc., 2012) for the solution of the optimization problems.

The first observation is that the prodecure is numerically viable. The computation of the function transforms V takes about 3 minutes, while the solution of every optimization problem (with N = 156 unknown parameters, as many as the size of R) takes about 3 seconds (for the automated choice of λ, optimization is repeated several times in the search of its best value). From the statistics of the estimates ρE in Fig. 1, the second observation is that estimates for both the cases G = G and G = G have limited bias (no bias at most lags δ), with increasing uncertainty toward the tail of ρE . Here, the automated choice of the regularization parameter λ reduces estimation error relative to a fixed λ. The increasing error in the tails has to be understood as the result of the memory of the reaction network vanishing for large lag times δ, and the smaller amount of measurements effectively entering the computation of ρE at large lags. All this is in perfect analogy with the results in absence of a modulat-ing signal of [START_REF] Cinquemani | Stochastic reaction networks with input processes: Analysis and application to gene expression inference[END_REF]. Interestingly, estimation performance is clearly worse for G than for G in the tails of ρE (larger estimation error variability) and also around 0 (larger bias). In our interpretation, this is strictly related with the null value of G at time zero, which reduces the sensitivity of the data to the information-rich (transient) system dynamics at early times. On the other hand, the estimation error of ρE does not equally reflect in the estimation error of ρF . Indeed, comparing Fig. 1 with Fig. 2, one finds that the error strengh in the estimation of ρ F is smaller (in a relative sense) than for ρE for large differences |z -t|. Thus, in general, performance should be assessed based on the context, depending whether the input (E(t)) or the controlled (G(t)E(t)) process is the object of primary interest. In summary, the procedure is viable and capable to estimate the statistics of a modulated noise process, with performance depending on the modulating (control) signal.

CONCLUSIONS

Starting from own previous work, in this paper we have developed an approach for the effective numerical solution of so-called GMEs for a class of modulated input processes where the modulating input is in the form of (or can be approximated by) splines. We have used this result to develop an algorithm for the estimation of the statistics of a controlled promoter activation process from gene expression population-snapshot data. Demonstrated via simulation, the method successfully reconstructed the promoter activity autocovariance function in the limits of a realistic dataset. The solution showed interesting dependencies on the choice of the control input that deserve future quantitative investigation. While developed with specific reference to gene expression reporter systems, many results of this paper are applicable to more general reaction networks and in other fields, where spline approximation of control signals and analysis of controlled continuous-time Markov chains is of relevance.

  c) ∈ C N where T (c) denotes the symmetric Toeplitz matrix with first column equal to c, C N is the convex cone of positive semi-definite matrices of order N , and Q ∈ C N . (iii) (Calculation of ρE and ρF ) Return the estimate of ρE as ρE (•) = R(•)ĉ and the estimate of ρ F as ρF (z, t) = G(z)ρ E (z -t)G(t) T . In (ii), the penalty term c T Qc is used to enforce regularity of the solution. For a suitable choice of Q, c T Qc accounts in particular for the average curvature of the candidate solution Rc. The optimization constraint ensures that Rc is a well-defined (positive semi-definite) autocovariance function. Factor λ ≥ 0 is a regularization weight that can be chosen automatically via repeated executions of step (ii) with different candidate values of λ. See more details in Cinquemani (2019). 5. NUMERICAL DEMONSTRATION We now demonstrate by simulation the inference procedure of Sec. 4. We consider the gene expression model of Sec. 3, modified for the presence of a modulating signal. That is, in place of F (t) = KU (t), we assume that F (t) = G(t)E(t), where G(t) = [k M g(t) 0 0 0] T and E(t) = U (t). The parameter values θ = (k M , d M , k P , d P ) = (0.5, 0.1, 0.2, 0.01), as well as the definition of the mean μE = μU = 0.5 and autocovariance function ρE = ρU (shown in Fig.

Fig. 1 .

 1 Fig. 1. Results from 100 runs of the estimation of ρE . Solid black line: True profile of ρE ; Dotted lines and shaded regions: At all values of δ, median of 100 estimates ρE (δ) and region between the 10% and 90% quantiles of the estimates, for fixed (λ = 10 6 , blue and magenta) and automated (grey and red) choices of the regularization factor, in the case G = G (magenta and red) and G = G (blue and grey). Plots are bottom-cropped for better scaling. The inset shows the modulating functions G (red) and G (black).

Fig. 2 .

 2 Fig. 2. Statistics of the estimation error ρF (z, t) = ρF (z, t) -ρ F (z, t) (case G = G , automated choice of λ). At every point (z, t), true autocovariance function ρ F (z, t) (colored shading) vs. maximal (absolute) value among the 10% and 90% quantile of ρF (z, t) over 100 simulation runs (grey mesh).


Gonzalo Santoro

Lidia Martínez

Koen Lauwaet

Mario Accolla

Guillermo Tajuelo-Castilla

Pablo Merino

Jesús M Sobrado

Ramón J Peláez

J Víctor

Isabel Herrero

Tanarro

Víctor J Herrero

Isabel Tanarro

Álvaro Mayoral

Marcelino Agúndez

Hassan Sabbah

Christine Joblin

José Cernicharo

José Ángel Martín-Gago

THE CHEMISTRY OF COSMIC DUST ANALOGUES FROM C, C 2 , AND C 2 H 2 IN C-RICH CIRCUMSTELLAR ENVELOPES

Keywords: stars, AGB and post-AGB -circumstellar matter -dust, extinction -ISM, lines and bands -methods, laboratory, molecular -methods, laboratory, solid state

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Corresponding author: Gonzalo Santoro, José Ángel Martín-Gago gonzalo.santoro@icmm.csic.es, gago@icmm.csic.es Carbonaceous dust is ubiquitously found in space, from the interstellar medium (ISM) to nova ejecta and circumstellar shells (Chiar et al. 2013). Astronomical observations in the Mid-Infrared (MIR) range reveal the presence of both aromatic and aliphatic dust components. In the diffuse ISM, the 3.4 µm absorption band (assigned to the CH vibrational stretching modes of CH 2 and CH 3 aliphatic moieties), along with the weaker absorption features at 6.8 and 7.3 µm (which correspond to the CH bending modes of aliphatic groups) have been attributed to hydrogenated amorphous carbon (Pendleton & Allamandola 2002;Dartois et al. 2004). On the other hand, the so-called aromatic infrared emission bands (AIBs) are widely observed in environments of our Galaxy that are submitted to ultraviolet photons (Peeters et al. 2002(Peeters et al. , 2004)). The AIBs, which fall in the spectral range from 3 to 20 µm (main bands at 3.3, 6.2, 7.7, 8.6 and 11.2 µm) have generally been assigned to polyaromatic carriers that are small enough to be stochastically heated by the absorption of a single ultraviolet (UV) photon, which constitutes the polycyclic aromatic hydrocarbon (PAH) hypothesis (Leger & Puget 1984;Allamandola et al. 1985Allamandola et al. , 1989;;Puget & Leger 1989). There are also observational evidences for an evolutionary scenario from aliphatics to aromatics calling for very small grains of mixed aliphatic-aromatic composition both in evolved star environments [START_REF] Goto | [END_REF]Kwok & Zhang 2011) and in molecular clouds (Pilleri et al. 2015). Overall, these observations call for a better understanding on the formation of PAHs and of carbon dust in general in astrophysical environments.

In our Galaxy, carbon-rich (C-rich) asymptotic giant branch (AGB) stars are the major sources of carbonaceous dust [START_REF] Gehrz | IAU Symposium[END_REF], thus chemical models have been developed to account for the formation of PAHs and hydrogenated carbon clusters in these environments (e.g. [START_REF] Pascoli | [END_REF]Cherchneff 2011). According to observations, acetylene (C 2 H 2 ) is one of the most abundant molecules in the inner regions of circumstellar envelopes (CSEs) where dust is formed, presenting a relative abundance to H 2 of 8 × 10 -5 (Fonfría et al. 2008), while carbon atoms are also expected to be highly abundant. Chemical equilibrium calculations predict the concentration of carbon atoms to be 1-2 orders of magnitude below that of C 2 H 2 in the region where carbon molecules are expected to condense into dust (Agúndez et al. 2020).

It is well known that the polymerization of C 2 H 2 gives rise to polyacetylenic chains (Cernicharo 2004), whose formation mechanism has been recently mapped with ALMA in the outer regions of the C-star envelope IRC +10216 (Agúndez et al. 2017). Furthermore, diacetylene (C 4 H 2 ) has been recently detected in this star through high spectral resolution MIR observations (Fonfría et al. 2018). Acetylene is also believed to play a crucial role in the formation of benzene and PAH compounds (Frenklach & Feigelson 1989). The formation of benzene in evolved stars has been rationalized from experiments of C 2 H 2 pyrolysis by the reaction of two propargyl (C 3 H 3 ) radicals (Miller & Melius 1992), a mechanism that requires high temperatures and high acetylene concentrations, which are not encountered in the dust nucleation zone of AGB stars. The temperature in the region between 2 and 4 stellar radii is estimated to range from 1500 K to 1000 K whereas a density of acetylene of 10 7 -10 5 molecules/cm 3 has been derived by Fonfría et al. (2008). This abundance for acetylene is in very good agreement with recent model predictions for for C-rich AGB stars (Agúndez et al. 2020).

Another chemical route for the formation of benzene, considering neutral-neutral reactions, and at much lower temperatures than those needed in the combustion theory, has been derived from acetylene discharges. It proceeds mainly by the addition of C 2 H 2 to the C 4 H 3 radical and to a lesser extent by the cyclization of C 6 H 4 (De Bleecker et al. 2006a). Further growth from benzene to PAHs is assumed to take place through the so-called Hydrogen Abstraction Acetylene Addition (HACA) mechanism (Frenklach & Warnatz 1987;Frenklach & Feigelson 1989;Shukla & Koshi 2012;Yang et al. 2016). Recently, a chemical route for the growth of large benzoid-PAHs (comprised only of 6membered ring structures) has been proposed through the so-called Hydrogen Abstraction Vinylacetylene Addition (HAVA) mechanism (Zhao et al. 2018). Thus, the complementary HACA/HAVA mechanisms might be needed for the molecular growth of PAHs in circumstellar environments.

Several laboratory experiments have dealt with the formation of stardust analogues from acetylene, either in dusty plasmas (Kovačević et al. 2005;[START_REF] Stefanović | [END_REF], in pyrolysis experiments (Jäger et al. 2009;Biennier et al. 2009) or in molecular jets exposed to plasma discharges (Contreras & Salama 2013). However, all of them lack of the addition of atomic carbon and therefore are not suitable to evaluate its impact on the growth of hydrocarbons with implications on, e.g., the formation of hydrocarbons with an odd number of carbon atoms (Cernicharo 2004).

Recently, we have shown that sputtering gas aggregation sources (SGAS) are particularly suited for studying dust formation in the CSE of AGBs (Martínez et al. 2020). In particular, we have demonstrated that the interaction of atomic carbon with molecular hydrogen promotes the formation of non-aromatic molecules, being acetylene one of the main gas phase products. Here, we expand our previous work to investigate the interaction of atomic carbon with acetylene due to its known importance in the formation and growth of hydrocarbons. In order to accelerate the chemistry between car-bon and acetylene so as to have access to most of the reaction pathways, as well as to investigate different conditions which would pertain to different regions of the CSE, we have intentionally increased the concentration of acetylene. In addition, our experimental conditions are also of interest for studying the envelopes of C-rich protoplanetary nebulae (PPNe) in which the UV radiation from the central star provokes the photodissociation of acetylene and methane providing a high abundance of carbon initiating a rich chemistry involving C, C 2 H 2 , CH 4 and other hydrocarbons (Cernicharo 2004).

We have found that, unlike what is commonly observed for dust analogues prepared in pyrolysis experiments or plasmas, a non-negligible amount of the produced material consists of carbon clusters (both pure and hydrogenated) with an odd number of carbon atoms, being a direct consequence of the interaction of atomic carbon with acetylene. Additionally, aromatics with aliphatic substitutions have been observed. We suggest that the formation of alkyl-substituted aromatics proceeds through the formation of CH 3 radicals, which is also a consequence of the addition of atomic carbon.

EXPERIMENTS

For the production and analysis of the dust analogues we employed the Stardust machine. The technical details of Stardust can be found elsewhere (Martínez et al. 2018(Martínez et al. , 2020)). The particular experimental setup used in this work is schematically depicted in Fig. 1. A pictorial view of the CSE region simulated by our experimental setup can be found in Martínez et al. (2020). The base pressure of the system is in the 10 -10 mbar range and the pressure during the production of the analogues at different positions of the machine is also indicated in Fig. 1.

Dust analogues were produced with a scaled-up Multiple Ion Cluster Source (MICS) (Martínez et al. 2012), a special type of sputtering gas aggregation source (SGAS), working in UHV conditions (Haberland et al. 1991). Atomic carbon was delivered by magnetron sputtering of a 2-inch graphite target (99.95 % purity) using Ar as sputtering gas with a flow rate of 150 sccm. The sputtered atoms coalesce inside the aggregation zone (aggregation length, i.e., distance from the magnetron target to the exit nozzle: 374 mm). C 2 H 2 (purity ≥ 99.6 % , diluted in acetone) was injected in the aggregation zone of the MICS through a lateral entrance at constant flow rates of 0, 4 × 10 -4 , 0.15 and 1 sccm. The DC magnetron is regulated in current and 0.2 A (corresponding to 100 W) were applied. Optical emission spectroscopy (OES) was performed inside the aggregation zone (at OES position in Fig. 1) and the light emitted by the sputtering plasma at a distance of around 1 cm from the magnetron target was collected through a fused silica window and a fused silica optical fibre. A 193 mm focal length, motorized Czerny-Turner spectrograph (Andor, model Shamrock SR-193i-A) equipped with a CCD camera (iDus DU420A-BVF) was employed. Two diffraction gratings with 1200 and 1800 grooves/mm, installed in a movable turret, provide spectral ranges of 300-1200 nm and 200-950 nm, respectively, and nominal spectral resolutions of 0.22 nm and 0.15 nm, respectively (for an input slit width of 20 µm). The relative spectral efficiencies of the spectroscopic equipment were quantified for both diffraction gratings with a calibrated tungsten lamp.

Mass spectrometry was used for in-situ analysis of the gaseous species produced during the production of the dust analogues. This was performed at positions labelled as QMS1 and QMS2 in Fig. 1. At QMS1 a quadrupole mass spectrometer (PrismaPlus, Pfeiffer) with a mass range of 0-100 amu and a Faraday cup detector was used. At QMS2, a Pfeiffer HiQuad QMG 700 with QMA 400 mass spectrometer (mass range of 0 to 512 amu) equipped with a CP 400 ion counter preamplifier was used, which enables a higher sensitivity. The pressure in the UHV chamber at position QMS1 during the analogue production precludes the use of a more sensitive Secondary Electron Multiplier (SEM) detector. Specific chemical compounds have been tentatively assigned from mass spectra considering the fragmentation pattern by electron impact ionization available at the NIST database (NIST Mass Spectrometry Data Center 2020).

Morphological characterization of the analogues was performed ex-situ by Atomic Force Microscopy (AFM) and Transmission Electron Microscopy (TEM). In both cases the deposits were collected through Sample entry 1 of Fig 1 . For the AFM measurements the dust analogues were collected on SiO x substrates whereas carbon grids were employed for TEM. The AFM measurements were carried out with Cervantes AFM System equipped with the Dulcinea electronics from Nanotec Electronica S.L. All images were recorded and analyzed using the WSxM software (Horcas et al. 2007). TEM was performed with a FEI-TITAN X-FEG transmission electron microscope used in scanning mode (STEM) and operated at 300 kV and at 120 kV. The images were acquired using a high-angle annular dark field (HAADF) detector. The microscope was equipped with a monochromator, Gatan Energy Filter Tridiem 866 ERS, a spherical aberration corrector (CEOS) for the electron probe (which allows for an effective spatial resolution of 0.08 nm) and an energy dispersive X-ray detector for EDS analysis.

The composition of the dust analogues was investigated in-situ by infrared spectroscopy in transmission geometry. The dust analogues were deposited on KBr substrates and transferred to a separated UHV chamber by means of an UHV suitcase (P < 5 × 10 -9 mbar). This protocol ensures that the inspected deposits were not air-contaminated. A VERTEX 70 V (Bruker) instrument was employed and the complete optical path was kept in vacuum. The spectral resolution was set to 2 cm -1 and 256 scans were coadded. A liquid nitrogen cooled Mercury Cadmium Telluride (MCT) detector was used. For the analysis of the IR spectra, spectral deconvolution was performed when peaks were not resolved due to overlapping of IR bands.

Finally, the molecular composition of the analogues was investigated ex-situ by laser desorption/ionization mass-spectrometry (LDI-MS) employing the AROMA (Aromatic Research of Organics with Molecular Analyzer) setup (Sabbah et al. 2017) using the so-called laser desorption-laser ionization (L2MS) scheme in which desorption and ionization are separated in time and space and performed with two different lasers. The desorption step was performed with the fundamental mode of a Nd:YAG laser (λ = 1064 nm) whereas the ionization step was performed employing the fourth harmonic of a Nd:YAG laser (λ = 266 nm). This scheme is particularly sensitive to PAHs due to the wavelength of the ionization laser. The ions generated are stored and thermalized in a linear quadrupole ion trap and subsequently monitored by time-of-flight mass spectrometry.

RESULTS

3.1. Gas-phase molecules produced during the formation of the dust analogues Fig. 2a shows the optical emission spectra during the dust analogue formation inside the aggregation zone of the MICS in the spectral regions corresponding to the CH (A-X) band, the C 2 Swan band (d 3 Π g -a 3 Π u ) and the H α line of the Balmer series for the different C 2 H 2 flow rates employed. When no C 2 H 2 was injected inside the aggregation zone, the molecular C 2 Swan band is clearly observed. Since the sputtering of graphite in SGAS produces mainly atomic carbon, we rationalize that C 2 is formed via the three body reaction C + C + Ar → C 2 + Ar (Martínez et al. 2020), which involves Ar atoms from the sputtering gas used.

As C 2 H 2 is injected, the C 2 Swan band vanishes and a clear signature of the presence of atomic hydrogen is observed through the appearance of the H α line. This suggests that the excited C 2 is consumed to some extent in the reaction C 2 + C 2 H 2 → C 4 H + H, which is known to be rapid even at low temperatures (Canosa et al. 2007;Daugey et al. 2008). In addition, as C 2 H 2 is injected, it can directly react with atomic carbon through, e.g., C + C 2 H 2 → C 3 H + H and/or C + C 2 H 2 → C 3 + H 2 (Clary et al. 2002), which diminishes the amount of atomic carbon available to form C 2 and releases atomic hydrogen as well. Furthermore, emission from the CH radical is detected for the two highest flow rates investigated. We have already shown that this radical can be formed by the interaction of atomic carbon with H 2 at very low densities, even with the residual H 2 in the UHV chamber (Martínez et al. 2020). The concurrent detection of the H α line suggests an excess of atomic hydrogen, which, if it is not participating in further chemical reactions, will recombine into H 2 by three-body reactions, thus increasing the overall H 2 density. The bimolecular reaction C + H 2 → CH + H is highly endothermic unless atomic carbon is electronically excited (Sato et al. 1998). The reaction of C with vibrationally excited H 2 could help in overcoming the barrier energy but under the physical conditions of the Stardust machine we do not expect high abundance of H 2 in the v=1, 2 vibrational levels. Therefore, CH is most probably formed by a three body reaction involving Ar (sputtering gas) as third body (C + H + Ar → CH + Ar), though we cannot completely rule out a contribution from the dissociation of C 2 H 2 in the magnetron plasma. From gas flow calculations, we have estimated that around 2% of the injected acetylene reaches the magnetron and is therefore susceptible to dissociate by electron impact.

In order to investigate the molecular species that are formed during the dust analogue production, we have performed in-situ mass spectrometry at position QMS1 in Fig. 1. Once at this position, no further growth of the analogues takes place. When the production of the analogues is carried out in the absence of C 2 H 2 , only Ar, which is used as sputtering gas, and the residual gases in the chamber (H 2 , H 2 O and CO) were detected (Fig. 3a). The peak at m/z = 80 corresponds to the formation of Ar dimers (Ar 2 ) during the gas expansion through the nozzle of the MICS.

On the contrary, when C 2 H 2 is injected we detected diacetylene (C 4 H 2 ) and for the highest flow rates employed also triacetylene (C 6 H 2 ) (Fig. 3a). The polymerization of acetylene to larger linear polyacetylenic chains is known to occur through the reaction with the ethynil radical

(C 2 H) via C 2 H + C 2n H 2 → C 2n+2 H 2 +
H. More interestingly, for the highest C 2 H 2 flow rate we have clearly detected C 6 H 4 and C 6 H 6 (Fig. 3c). C 6 H 4 has been suggested to be involved in the formation of benzene through a cyclization mechanism (De Bleecker et al. 2006a) and the concomitant detection of both species in our experiments seems to support this mechanism.

As shown in Fig. 3d, the formation of all the molecules detected by mass spectrometry is accompanied with the consumption of C 2 H 2 once atomic carbon is released from the magnetron target. In addition, we have observed the production of H 2 (m/z = 2). Atomic hydrogen is produced in high amounts by the polymerization reactions of C 2 H 2 , thus it is available to participating in further chemical reactions, including its recombination to molecular hydrogen by three body reactions.

The signals corresponding to H 2 O (m/z = 18) and CO 2 (m/z = 44) were stable, indicating that they were not involved in the chemical reactions. However, we observed slight changes in the signal associated to the main fragment of acetone (m/z = 43), which is present as an impurity in the C 2 H 2 (concentration lower than 0.3 %), along with a slight increase in the amount of CO (m/z = 28). Nevertheless, we have not detected any other gas-phase O-bearing molecule, implying that the acetone impurity in acetylene is only weakly (if any) contributing to the chemistry involved in the formation of the analogues. Moreover, the reaction between C and H 2 O has a very large barrier of ≃20000 K (Mayer et al. 1967) and that of C with CO 2 has been measured at 300 K with a rate below 10 -15 (Husain & Young 1975); thus, the formation of O-bearing species is very unlikely in our experiments.

Higher sensitivity gas-phase mass spectrometry was performed at position QMS2 in Fig. 1 for the highest C 2 H 2 flow rate (Fig. 3e). The differences observed in the relative intensities among C 6 H 2 , C 6 H 4 , C 6 H 6 and Ar 2 in Figures 3 (c) and (e) are due to the different sensitivity of the QMS as well as the different pressures at locations QMS1 and QMS2. Concerning polyacetylenic chains, we detected up to C 8 H 2 . Furthermore, apart from C 6 H 6 , we observed the formation of larger aromatic molecules including naphthalene (C 10 H 8 ), acenaphthylene (C 12 H 8 ) and biphenyl/ethenylnaphthalene (C 12 H 10 ). These molecules are important for the growth of larger PAHs through different mechanisms including the well-known HACA mechanism (Shukla & Koshi 2012). We have also detected phenylacetylene (C 8 H 6 ) which is a stable intermediate in the growth of naphthalene from benzene through the HACA mechanism. Interestingly, unlike what is observed in acetylene pyrolysis (Shukla & Koshi 2012) or acetylene discharges ). The HACA mechanism does not explain the formation of aromatic compounds with aliphatic substituents. However, the addition of atomic carbon and its interaction with molecular hydrogen (both residual in the UHV system and formed by recombination of the released atomic hydrogen) promotes the formation of alkyl radicals, which open up chemical routes for the formation of alkylsubstituted aromatics (see Sec. 4).

Finally, we detected a number of hydrogenated carbon clusters (HC-clusters) of aliphatic nature with odd number of carbon atoms from C 5 to C 11 , again differing from what is found in acetylene discharges/pyrolysis in which predominantly molecules with even number of carbon atoms are formed (see Fig. 3e). The apparent absence of C 7 HC-clusters might be due to a blurring of the signal by the peaks associated to toluene and its main fragment (m/z = 91/92). For C 3 clusters, it is the Ar signal which precludes its detection. The production of these C-clusters is a clear signature of the interaction of atomic carbon with C 2 H 2 . The reaction of C and C 2 H 2 is known to be fast, yielding C 3 and/or C 3 H (Liao & Herbst 1995;Clary et al. 2002). Analogous reactions of C with larger polyacetylenic chains are as well an efficient way to produce hydrocarbons with an odd number of carbon atoms in chemical environments where C and C 2 H 2 are abundant (Cernicharo 2004). In addition, the reaction of C with other large hydrocarbons is also known to occur very fast and without barrier (Haider & Husain 1992;Liao & Herbst 1995;Husain & X. Ioannou 1997).

Overall, the results concerning in-situ mass spectrometry evidence the formation of three different families of gas-phase chemical compounds, namely, polyacetylenic chains, aromatic species, either pure and with aliphatic substitutions, as well as C-and HC-clusters of aliphatic nature. The identified masses are listed in Table 1 red along with their tentative assignments to chemical compounds (NIST Mass Spectrometry Data Center 2020). A detailed description of the possible formation mechanisms is provided in Sec. 4.

Morphology of the dust analogues

Apart from the gas phase species, most of the material produced during the formation of the dust analogues consists of particles with diameters in the nanometer range (Fig. 4). Due to the low kinetic energy of the particles produced using gas aggregation sources, the particles soft land on the substrate and thus retain the gas phase morphology and structure. Fig. 4a andb show typical AFM images of the dust analogues when no C 2 H 2 was injected in the aggregation zone (pure C particles) and for a C 2 H 2 flow rate of 1 sccm, respectively. In both cases nanoparticles (NPs) are observed with mean NP diameters of 9.3 nm and 7.0 nm for pure C NP and for a C 2 H 2 flow rate of 1 sccm, respectively. However, a broader size distribution is obtained when injecting C 2 H 2 (see Fig. 4 c andd). Moreover, the production rate increased enormously when C 2 H 2 was injected (by a factor around of 350 in terms of collected NPs/µm 2 s), which evidences the interaction of atomic C with C 2 H 2 , resulting in a manifest acceleration of the chemistry.

The TEM analysis of the NPs revealed the formation of amorphous NPs irrespective of the injection of C 2 H 2 (Fig. 4e-f). The main difference consists in the shape of the particles. Pure C NPs present a well-defined round shape whereas those formed after C 2 H 2 injection tend to agglomerate. This agglomeration is observed in the AFM images as well, but it is better seen in the TEM images. Diffusion of the NPs once deposited on the substrates is not expected, thus the agglomerates might be already formed in the gas-phase. For the highest C 2 H 2 flow rate (1 sccm), we measured the IR spectrum of the dust analogues (Fig. 5). The IR spectrum reveals that the sample consists in a complex mixture of hydrocarbons with sp, sp 2 and sp 3 carbon hybridizations, thus the gas-phase products detected by mass spectrometry are incorporated into the solid dust analogues. The observed infrared bands are listed in Table 2 along with their tentative assignments.

Composition of the dust analogues

The presence of sp carbon is evidenced by the absorption bands at 3322 and 3300 cm -1 corresponding to the ≡C-H stretching modes, and by the band at 2105 cm -1 , which corresponds to the -C≡C-stretching mode of monosubstituted acetylene. In addition, a second -C≡Cstretching mode is observed at 2200 cm -1 which is ascribed to the stretching mode of conjugated triple bonds -C≡C-C≡C- [START_REF] Socrates | Infrared characteristic group frequencies[END_REF]). This band is also found in dust from acetylene plasmas (Stoykov et al. 2001) but is often too weak to be observed (Kovačević et al. 2005;[START_REF] Stefanović | [END_REF]. In addition, it is not commonly observed in dust analogues from acetylene pyrolysis (Biennier et al. 2009). The fact that we observe this band might be related to the addition of atomic carbon and the formation of C 2 , which in the end increases the formation of linear polyacetylenic chains via the C 4 H radical route (see Sec. 4).

On the other hand, numerous bands indicate the presence of sp 3 carbon and aliphatic moieties. The bands at 1369 cm -1 , 1379 cm -1 and 1445 cm -1 , assigned to the CH 3 symmetric bending and CH 3 asymmetric bending / CH 2 scissoring modes, respectively (see Table 2), are a clear signature of alkyl moieties. These bands are accompanied by the bands in the region between 3000 cm -1 and 2800 cm -1 , which correspond to the CH 2 and CH 3 stretching modes. Moreover, the band at 968 cm -1 is ascribed to the CH 3 rocking mode. Concerning sp 2 carbon, the bands at 1600 cm -1 , 1578 cm -1 , 1492 cm -1 and 1437 cm -1 reveal the presence of aromatics in the dust analogues. All are assigned to aromatic C=C stretching modes. In the case of alkyl substituted aromatics, the band at 1578 cm -1 is usually a shoulder of that at 1600 cm -1 , as is our case. Moreover, the band at 1492 cm -1 is intense for substituted aromatics [START_REF] Socrates | Infrared characteristic group frequencies[END_REF]. These results are consistent with the detection of aromatics with aliphatic substituents by mass spectrometry and reveal that the aromatic compounds incorporated into the dust analogues are not only pure but present aliphatic substituents as well. We note that a similar result has been recently obtained by Gavilan et al. (2020) during the formation of dust analogues using small PAHs as precursors in a molecular jet exposed to an electrical discharge. Other bands related to the aromatic compounds incorporated into the dust analogues are the CH out-ofplane bending modes. The different bands observed at 890 cm -1 , 843 cm -1 , 776 cm -1 and 756 cm -1 are related to different number of peripheral adjacent H atoms (solo, duo, trio and quartet, respectively) (Hony et al. 2001;Carpentier et al. 2012) and are characteristic of PAHs though insensitive to PAH substitution. However, the band at 700 cm -1 can be assigned to CH out-ofplane bending mode of substituted benzenes [START_REF] Socrates | Infrared characteristic group frequencies[END_REF]. Finally, the band at 640 cm -1 is attributed to the CCC in plane bending of aromatic compounds (Gadallah et al. 2012) and the combination band at 1938 cm -1 is ascribed to the combination of out-of-plane and in-plane CH bending modes (combination and overtone bands in the spectral region 2000-1700 cm -1 are characteristic of aromatic compounds).

Apart from aromatic compounds, there are also bands related to olefinic sp 2 carbon. The sp 2 CH stretching bands in the region 3100-3000 cm -1 are due to both aromatic and olefinic sp 2 carbon, but the shoulder at 1630 cm -1 is unambiguously assigned to the C=C stretching modes of olefinic compounds (Fig. 5b) [START_REF] Socrates | Infrared characteristic group frequencies[END_REF]Gavilan et al. 2017) . The presence of olefinic moieties is further confirmed by the characteristic wagging modes of olefins at 988 cm -1 and 910 cm -1 [START_REF] Socrates | Infrared characteristic group frequencies[END_REF]. The band at 1074 cm -1 might be due to the CH 2 rocking mode characteristic of vinyl groups, in accordance with the observation of styrene (and maybe also vinylnaphthalene) by mass spectrometry, though the aromatic inplane CH bending modes overlap in this region.

The broad band that appears in the 1800-1000 cm -1 region is characteristic of amorphous carbon (aC) and its intensity is related to the relative amount of sp 2 to sp 3 carbon hydridization (Ferrari et al. 2003;Rodil 2005). Its presence in the IR spectrum is consistent with the amorphous morphology observed by TEM and suggests that part of the dust analogues is made of aC material.

Finally, weak C=O carbonyl stretching bands are observed which we consider as related to the acetone impurity in the C 2 H 2 bottle since the base pressure of the Stardust machine is 10 -10 mbar and the sample was transferred to the IR UHV chamber by means of an UHV suitcase at a pressure lower than 5 × 10 -9 mbar.

In addition, we have tested the stability of the analogues after air exposure and we have observed that they present a high reactivity (Fig. 6), whereas we did not observe any significant change in the IR spectrum of the analogues after 24 hours in UHV. Carbonaceous materials produced from unsaturated precursors are known to react swiftly with oxygen, leading to the formation of C-O-C moieties. Mainly, we have observed the incorporation of oxygen through the increase of the carbonyl C=O stretching modes and through the appearance and growth of a band at 1250 cm -1 , assigned to the C-O-C stretching mode. In addition we have observed a reduction in the sp carbon content, which is evident by the decrease of the sp CH stretching modes (3300 cm -1 ) and of the -C≡C-stretching mode of conjugated alkynes (2200 cm -1 ). An increase in sp 3 aliphatic moieties was also observed whereas the sp 2 CH stretching modes remain unaltered (see e.g. the normalized optical depth (OD) for the bands at 2925 cm -1 and at 3056 cm -1 in Fig. 6b). Further changes in the 1000-600 cm -1 region can also be observed, though difficult to interpret. Anyhow, these results evidence the extreme importance of both working in ultra-clean environments (such as UHV conditions) and of in-situ characterizing the dust analogues produced from acetylene.

Albeit the tested chemical reactivity of the dust analogues towards atmospheric air, we have also performed ex-situ laser desorption/ionization mass-spectrometry (LDI-MS) of the dust analogue produced with a C 2 H 2 flow rate of 1 sccm (Fig. 7), using the AROMA setup (Sabbah et al. 2017).The results are not quantitative in terms of the relative amount of the observed compounds because of the known chemical evolution of the analogues after air exposure, but they provide additional information on the chemical composition. Linear polyacetylenic chains (in cationic form) were detected up to C 8 H 2 , which is somehow surprising due to the high volatility of these compounds and indicates that these species are trapped inside the dust particles and released by laser desorption. Other volatile species detected by in-situ mass spectrometry, such as benzene or toluene, were not detected by LDI-MS despite the high sensitivity of the technique to aromatics (Sabbah et al. 2017). C 6 H 4 could be identified, which was also de-tected by in-situ mass spectrometry and is an important molecule for the formation of benzene. Moreover, C 2 H + 2 and C 3 H + 3 are believed to be produced during the laser desorption stage from the common C 3 H 2 precursor (cf. supplementary material in Martínez et al. 2020).

On the other hand, C-and HC-clusters with carbon atoms from 5 to 20 dominate the mass spectrum, being C 5 H and C 11 particularly abundant. Relatively large aromatics are also present in the dust analogue from C 12 H 8 (acenaphthylene) to C 20 H 12 . Finally, the peaks assigned to C 3 HO and C 7 HO evidence oxidation of the sample likely due to air exposure. The results from LDI-MS are summarized in Fig. 7b by using a Double Bond Equivalent (DBE) analysis. Briefly, in the case of hydrocarbons, DBE can be defined as the sum of the number of double bonds involving carbon plus the number of rings for a particular species. A description of the DBE analysis can be found in, e.g, Sabbah et al. (2017), Gavilan et al. (2020) and in the Supplementary Information of Martínez et al. (2020). Three molecular families, namely aromatics, HC-and C-clusters, were considered and their relative intensities in terms of summed peak intensities are shown. Linear polyacetylenic chains, C 2n H 2 , are observed in the outer layers around C-rich stars and their formation involves the ethynyl radical (C 2 H), which is formed by the photodissociation of C 2 H 2 in the outer layers of the CSE and triggers the chemistry of linear polyacetylenic chains (Cernicharo 2004;Agúndez et al. 2017). Once the C 2 H radical is formed, the formation of polyacetylenes proceeds mainly through reactions of the type

C 2 H + C 2n H 2 → C 2n+2 H 2 + H (R1)
which are assumed to occur fast at low temperatures (see scheme in Fig. 8) (Chastaing et al. 1998;Agúndez et al. 2017). A second important chemical route for the formation of triacetylene and larger polyacetylenic chains involves the radical C 4 H through reactions of the type (Berteloite et al. 2010;Agúndez et al. 2017)

C 4 H + C 2n H 2 → C 2n+4 H 2 +H (R2).
Besides the photodissociation of diacetylene, the C 4 H radical can be efficiently formed by the interaction with dicarbon through the bimolecular reaction

C 2 + C 2 H 2 → C 4 H +H (R3).
In CSEs around C-rich AGB stars, the interaction of C 2 H 2 with energetic radiation can lead to dissociation

C 2 H 2 C 4 H 2 C 6 H 2 C 8 H 2 C 2n H 2 C 2 H -H C 2 H -H C 2 H -H C 2 H -H C 4 H -H C 4 H -H C 4 H -H Figure 8
. Reaction pathway for the formation of linear polyacetylenic chains.

but also to ionization. Taking into account that ionneutral reactions are usually faster than neutral-neutral reactions (e.g., the rate constant of Anicich 1993;Chastaing et al. 1998), chemical reactions involving ions may also play a role in the synthesis of polyacetylenic chains. The polymerization of C 2 H 2 in the external layers is dominated by neutral-neutral reactions involving C 2 H, rather than by ion-molecule reactions involving C 2 H + 2 , because photodissociation of C 2 H 2 by the interstellar UV field is 5-10 times faster than photoionization (Heays et al. 2017;Agúndez et al. 2018). However, in the inner clumps of CSEs, where dust acts as an efficient shielding against UV photons but cosmic rays can penetrate, some C 2 H 2 may be chemically processed due to cosmic-ray induced ionization and further ion-molecule reactions. The Cstar envelope IRC +10216 is known to present a clumpy structure (Cernicharo et al. 2015) and, in fact, cationneutral reactions have been suggested to play a role for the formation of, e.g., CH 3 CN, SiH 3 CN and CH 3 SiH 3 (Agúndez et al. 2008;Cernicharo et al. 2017). However, the measured abundances of ions (cations and anions) is very low (X≤10 -10 ), precluding an important role in the chemistry of the most abundant carbon-bearing species (note that the anions of all C n H radicals have been detected in the external and cold layers of the envelope but not in the inner warm regions).

C 2 H + 2 + C 2 H 2 → C 4 H + 2 + H 2 is 2.5-5 times higher than that of C 2 H + C 2 H 2 → C 4 H 2 + H;
In our experiments, as mentioned in the previous section, we estimate that around 2% of the injected acetylene reaches the magnetron and can therefore be dissociated into C 2 H by electron impact (see Fig. 9 for a scheme on the formation mechanism of radicals from atomic carbon and acetylene in SGAS). The interaction of C 2 with H 2 (either residual in the chamber as well as produced by recombination of the atomic hydrogen in excess) contributes as well to the formation of the C 2 H radical via the bimolecular reaction C 2 + H 2 → C 2 H + H (Cernicharo 2004;Martínez et al. 2020). In addition, and due to the higher cross section for electron impact ionization over electron impact dissociation (Mao et al. 2008), chemical reactions involving cations such as C 2 H + 2 and C 2 H + might take place. However, the elec-

C 2 H 2 C 2 H C C 4 H C 2 C 3 CH C 3 H e - H 2 -H C, Ar C, Ar H 2 , Ar C 2 H 2 C 2 H 2 -H C 2 H 2 -H -H -H -H 2 Figure 9
. Formation mechanism of radicals from atomic carbon and acetylene in SGAS. Molecular hydrogen is both residual in the UHV system and formed by recombination of the released atomic hydrogen.

tron temperature of the plasma in magnetron sputtering discharges and in SGAS sources is lower than 10 eV (Ivanov et al. 1992;Kousal et al. 2017), which favours electron impact dissociation over direct ionization or dissociative ionization events. This is different to the case of dusty plasmas in which apart from the neutral dissociative production of C 2 H, there is a higher probability for the formation of C 2 H + 2 (from C 2 H 2 ) and C 2 H + (from C 2 H) by electron-knocking off and dissociative ionization processes, respectively. Thus, the chemistry of the formation of nanoparticles in C 2 H 2 dusty plasmas presents a non-negligible ion-neutral contribution (Jiménez-Redondo et al. 2019). In addition, anionic polymeryzation routes can also be relevant in acetylene plasmas (Deschenaux et al. 1999;De Bleecker et al. 2006b;Jiménez-Redondo et al. 2019) due to the trapping of negative ions in the plasma potential, despite the electron attachment probability being several orders of magnitude lower than those for the electron impact dissociation/ionization. However, the chemistry involving anions is very unlikely in SGAS, in which it has been shown that the chemistry proceeds predominantly by neutral-neutral interactions (Haberland et al. 1991;Martínez et al. 2020) .

On the other hand, the results from OES suggest that the excited C 2 is consumed to some extent in the reaction (R3) to produce the C 4 H radical. Therefore, during the production of the dust analogues the growth of polyacetylenic chains involves mainly the C 2 H and C 4 H radicals (Fig. 8). Moreover, the incorporation of conjugated polyalkynes to the dust analogues in higher amounts than those observed in acetylene discharges, as probed by IR spectroscopy, supports a significant contribution of the C 4 H chemical route, which is initiated by the addition of atomic carbon and subsequent formation of C 2 in substantial amounts.

All the above mentioned reactions leading to the growth of polyacetylenic chains release either atomic or molecular hydrogen, which is consistent with the observed increase in the peak at m/z = 2 (Fig. 3c) and with the increase of the H α line in the optical emission spectra (Fig. 2a), being therefore available in high amounts to participating in further chemical reactions. We assume that the atomic hydrogen that is not consumed in chemical reactions is recombined into H 2 by three body reactions involving both Ar and the solid dust analogues.

Formation of PAHs

Benzene formation is considered as the bottleneck for the growth of PAHs, which have been proposed as the carriers of the AIBs and the seeds of carbonaceous cosmic dust. Benzene has been detected in protoplanetary nebulae (Cernicharo et al. 2001b;Malek et al. 2011) but not yet in the CSE of C-rich AGB stars, which constitute the main factories of cosmic dust. In addition, no specific PAH structure has been identified from observations. In particular, the apolar (or weakly polar) character of these molecules along with large rotational partition functions, hinder its detection by rotational emission spectroscopy (e.g. [START_REF] Joblin | EAS Publications Series[END_REF]. Recently, benzonitrile has been detected in the ISM (McGuire et al. 2018), whose formation might be linked cyanopolyynes in the same manner as benzene is believed to be related to polyynes.

Despite the lack of successful detection of benzene or PAHs in CSEs around C-rich AGB stars, a bottomup chemical pathway, derived from combustion experiments, has been proposed for the formation of benzene and PAHs in the warm inner regions of C-rich CSEs (Cherchneff et al. 1992;Cherchneff 2012). In this scheme, the key reaction is the direct recombination of two propargyl radicals to form both cyclic and linear C 6 H 6 as well as the phenyl radical (C 6 H 5 )

C 3 H 3 + C 3 H 3 → C 6 H 6 (R4a) C 3 H 3 + C 3 H 3 → C 6 H 5 + H (R4b)
being reaction (R4a) more efficient than (R4b). A different chemical scheme involving hydrocarbon ions has been proposed to explain the formation of benzene in the protoplanetary nebula CRL 618 (Woods et al. 2002), where C 6 H 6 has been detected. In this scheme, the synthesis is initiated with the proton transfer to C 2 H 2 to form C 2 H + 3 , followed by successive reactions with C 2 H 2 to form cyclic C 6 H + 5 , and ended with the reaction of C 6 H + 5 with H 2 to form C 6 H + 7 , which yield benzene upon dissociative recombination with electrons.

Apart from a top-down etching process (Merino et al. 2014) not considered here, another possible bottomup mechanism for the formation of benzene involving neutral-neutral reactions has been derived from plasma discharges (De Bleecker et al. 2006a). In this case, C 2 H, which is also crucial to the formation of polyacetylenic chains, initiates the chemistry mainly by the termolecular association reaction

C 2 H + C 2 H 2 + M → C 4 H 3 + M (R5).
where M stands for a third-body. In our case, due to the temperatures involved in the sputtering process in SGAS (in the order of 500 K) we rule out the formation of benzene via propargyl radicals (note that the C 3 H + 3 observed by LDI-MS is formed during the laser desorption step by protonation of C 3 H 2 ). In fact, the concurrent detection of C 6 H 4 and C 6 H 6 by in-situ mass spectrometry points towards a similar neutral-neutral route as that proposed for acetylene discharges, involving C 2 H and the formation of the butadienyl radical (C 4 H 3 ). Moreover, the formation of C 4 H 3 is less efficient than the formation of C 4 H 2 , which explains the observation of C 6 H 6 only at the higher C 2 H 2 flow rates employed in our experiments. Once C 4 H 3 is formed, the main reactions with C 2 H 2 produce C 6 H 5 and C 6 H 4 (Fig. 10), which account for about 75 % and 25% of the reaction mechanism towards C 6 H 6 , respectively (De Bleecker et al. 2006a). In our case, the detection of C 6 H 5 is hindered due to the location of the mass spectrometers, which only allows stable species to be detected.

Once C 6 H 6 is formed, the growth towards PAHs has been considered in astrophysical models to proceed by the so-called HACA mechanism in which acetylene is the main species propagating the growth of aromatic rings (Fig. 10). In addition, the HAVA mechanism (Zhao et al. 2018) as well as fusion of PAHs (Shukla & Koshi 2010) have been proposed for the growth of large benzoid PAHs. However, none of them explains the formation of aromatics with aliphatic substituents (except those of acetylenic nature), such as those observed in our experiments, both in the gas-phase (by in-situ mass spectrometry) and in the solid dust analogues (by in-situ IR spectroscopy).

Formation of alkyl-substituted aromatics

Alkyl substituted aromatics have been proposed to be present in space based on the analysis of the 3.3/3.4 µm band intensity ratio (Joblin et al. 1996). Li & Draine (2012) aliphatic fraction in PAHs derived from the intensities of the 3.4 µm and 6.85 µm emission features. The formation of alkyl substituted aromatics is expected to involve the reactivity of the methyl radical (CH 3 ) with aromatic units and/or aromatic radicals (Shukla et al. 2010).

In CSE around C-rich AGB stars, CH 4 has been detected (Keady & Ridgway 1993) and it is likely that methyl radicals are present, although they have not been detected and their abundance is unknown. Chemical equilibrium predicts a relatively low abundance of 10 -9 relative to H 2 in the stellar atmosphere (Agúndez et al. 2020), and unfortunately chemical kinetics models of the inner wind (e.g., Cherchneff 2012) have not provided abundance estimates for CH 3 . A bottom-up synthesis of CH 3 must involve the methylidyne radical (CH), which might arise from the interaction of atomic carbon with molecular hydrogen via C + H 2 → CH + H. However, this reaction is endothermic by about 100 kJ/mol and necessitates temperatures in the range of several thousand Kelvin to be efficient. Thus, shock-wave induced chemistry would be needed to efficiently form the methylidyne radical from a bottom-up approach. Once CH radicals are formed, the interaction with H 2 through the following reaction

CH + H 2 + M → CH 3 + M (R6)
where M stands for a third-body, can lead to CH 3 . If this is the case, cosmic dust might be needed as third body to efficiently produce CH 3 . On the other hand, CH 3 radicals can be formed in the outer layers of the CSE by photodissociation of CH 4 the same way as C 2 radicals are formed from C 2 H 2 .

Methylated species have been detected in the C-star AGB IRC+10216 (Cernicharo et al. 2000;Agúndez et al. 2015;Cernicharo et al. 2017). However, their low abundances suggest that they play a minor role in the chemistry of the envelope. A different situation has been found in the CRL618 PPN, where methylpolyynes (CH 3 C 2 H and CH 3 C 4 H) have been detected with abundances much larger than those observed in IRC+10216 (Cernicharo et al. 2001a) Aromatic compounds with aliphatic moieties are neither produced in acetylene pyrolysis (Shukla & Koshi 2012) nor detected (or only very marginally) as neutral species in the gas phase of acetylene discharges (Deschenaux et al. 1999;Consoli et al. 2008). However, this is drastically different when other hydrocarbons are used, particularly in ethylene (C 2 H 4 ) pyrolysis and ethylene plasmas (Jäger et al. 2006;[START_REF] Gehrz | IAU Symposium[END_REF]. This is related to the formation of alkyl radicals which enter into the chemistry, thus promoting the formation of alkyl moieties. Shukla et al. (2010) demonstrated that pyrolysis of aromatics in the presence of CH 3 radicals can lead to aromatic growth through a Methyl-Addition Cyclization process andGavilan et al. ( 2020) have recently proposed a route for PAH growth at low temperatures (150-200 K) by the addition of alkyl radicals.

In our case, the CH radical, which we observe by OES, is formed by the reaction C + H + Ar/dust → CH + Ar/dust (R7) involving Ar and the solid dust analogues as thirdbody.

Once the CH radical is formed, subsequent three body associations with atomic H would lead to CH 3 . This chemical route continues to the formation of CH 4 , but the acetone impurity in the C 2 H 2 bottle prevents its detection by in-situ mass spectrometry. Nevertheless, we have recently demonstrated the formation of CH 4 and C 2 H 4 by the interaction of atomic C and H 2 even at very low H 2 concentrations (Martínez et al. 2020).

The formation of toluene (C 7 H 8 ) proceeds then by the reaction of phenyl (C 6 H 5 ) and methyl (CH 3 ) radicals through

C 6 H 5 + CH 3 → C 7 H 8 + H (R8)
being the phenyl radical formed as shown in Fig. 10. The growth towards longer alkyl-chain substitutions to form e.g. ethylbenzene (C 8 H 10 ) proceeds through hydrogen abstraction and methyl addition . This process is not limited to the phenyl radical but applies to any other aromatic radical formed during the growth of PAHs through the HACA mechanism (e.g., naphtyl: C 10 H 7 , acenaphtyl: C 12 H 7 , . . . ) and provides a route for the growth of alkyl substituted PAHs (Fig. 11), closely related to the Methyl-Addition Cyclization mechanism (Shukla et al. 2010).

On the other hand, the likely detection of styrene (C 8 H 8 ) (maybe vinylnaphthalene (C 12 H 10 ) as well) by in-situ mass spectrometry along with the identification of olefins by IR spectroscopy (including a band likely related to vinyl moieties) points towards the production of the vinyl radical (C 2 H 3 ). Bimolecular reactions involving C, C 2 H 2 and H 2 to form C 2 H 3 are highly endothermic for all the possible chemical routes. The same holds for the addition of C 2 H 3 to PAHs and aromatic radicals. Therefore, we assume that the formation of aromatics with vinyl substitutions proceeds by three body reactions involving both Ar and the solid dust analogues.

Formation of C-clusters with odd C atoms

Molecules and carbon clusters with odd number of carbon atoms have been detected in the circumstellar environments of C-rich stars (Hinkle et al. 1988;Bernath et al. 1989) as well as in interstellar clouds (Maier et al. 2001). However, the growth of carbon clusters from C and C 2 at the temperatures found in the CSE is still not well characterized since most chemical reactions involving small carbon clusters are still unknown. Up to date only the C 2 +C 2 reaction has been measured at high temperature (Kruse & Roth 1997), which has been shown to present a high reaction rate. Therefore, although it potentially contributes to the chemistry, much is yet to be understood on the role of the C n +C m reactions in the formation of large carbon clusters. An alternative and efficient chemical route is provided by the interaction of atomic carbon with C 2 H 2 through the reactions (Cernicharo 2004)

C + C 2 H 2 → C 3 H + H (R9a) C + C 2 H 2 → C 3 + H 2 (R9b).
Additionally, C 2 H plays again an important role, as in the formation of polyacetylenic chains and benzene, as well as C 4 H by the interaction with atomic carbon via

C + C 2 H → C 3 + H (R10a) C + C 4 H → C 5 + H (R10b).
Once C 3 and C 5 are formed by (R9) and (R10), the towards larger molecules with odd number of carbon atoms proceeds through the following reactions with acetylene Therefore, the growth of C-clusters with an odd number of carbon atoms proceeds through two different chemical routes: an acetylenic route (R11) and a radical route (R12). In fact, these two chemical pathways do not apply only to C 2n+1 -clusters but also to C 2nclusters, providing a chemical way to the formation of hydrogenated C-clusters.

C 2n+1 + C 2 H 2 → C 2n+3 H + H (R11a)
The above-mentioned mechanisms necessitate atomic carbon to initiate the chemistry, which is not present in acetylene combustion experiments or in acetylene discharges. In fact, in particle-forming acetylene discharges the dominant species are by far hydrocarbon molecules with even number of carbon atoms (Deschenaux et al. 1999;Consoli et al. 2008;Kovačević et al. 2003;Benedikt 2010). On the contrary, the supply of atomic carbon in our case opens up both the acetylenic and the C 2 H/C 4 H radical routes described, resembling more closely the mechanism operating in the circumstellar environments of C-rich stars. Figure 12 summarizes the chemical network for the growth of C-clusters. Three body reactions involving Ar, which cannot be excluded in our experiments, are also shown.

CONCLUSIONS

In this article we present a thorough study of the gas-phase chemistry of acetylene in the presence of C and C 2 which can pertain to the outer layers of C-richAGB stars and protpplanetary nebulae. We have formed and studied in-situ cosmic dust analogues employing atomic carbon and acetylene as precursors, using a a well-controlled ultra-clean experimental set-up. Apart from the formation of linear polyacetylenic chains, we have observed the formation of polyaromatic hydrocarbons (PAHs). More importantly, aromatics with aliphatic substitutions as well as pure and hydrogenated carbon clusters were produced as a direct consequence of the addition of atomic carbon. The growth of carbon clusters proceeds from the interaction of atomic carbon with acetylene whereas the formation of alkylsubstituted aromatics is rationalized by the formation of methyl radicals and subsequent incorporation into aromatics by a hydrogen abstraction-methyl addition mechanism. In addition, the observed gas-phase species incorporate into the dust analogues, which consist in a complex mixture of sp, sp 2 and sp 3 hydrocarbons with an amorphous morphology. Our results are of particular interest for unveiling chemical routes leading to formation of acetylene-based molecular species in the external layers of AGB stars and in protoplanetary nebulae and foster the search for alkyl-substituted aromatics in these environments.

Figure 1 .

 1 Figure1. Scheme of the experimental setup used for the analog production and analysis. The pressure at each UHV chamber during the formation of the analogues is indicated.

Figure 2 .

 2 Figure 2. a) Optical emission spectra in the spectral regions of the CH (A-X) band, C2 Swan band (D 3 Π g -a 3 Π u ) and the Hα line of the Balmer series for the different C2H2 flow rates employed. The flow rates are indicated in the figure. The spectra are vertically shifted for clarity. b) Evolution of the CH, C2 and Hα line intensities with the C2H2 flow rate. The wavelength employed for each line is indicated in the figure.

Figure 3 .

 3 Figure 3. a) In-situ mass spectra at position QMS1 in Fig. 1 for the different C2H2 flow rates employed (indicated in the figure). The peaks labelled by a star correspond to the acetone impurity of the C2H2 gas. The curves have been shifted for clarity. b) and c) Expanded regions for the 1 sccm flow rate. d) In-situ gas phase detection of m/z = 2, 18, 26, 28, 43, 44, 50, 74, 76 and 78, which are indicated in the figure along with the assignment to the main contributing chemical species for the 1 sccm flow rate. The dashed vertical lines indicate the switching on and off of the magnetron. e) Insitu mass spectra for the 1 sccm C2H2 flow rate at position QMS2 in Fig.1. Peaks labeled in red are assigned to polyacetylenic chains, in blue to aromatic molecules and in green to aliphatics. In a) C2nH2 label mass peaks associated to the fragmentation of polyacetylenic chains by electron impact in the quadrupole.

Figure 4 .

 4 Figure 4. Morphology of the dust analogues without (top) and with (bottom) C2H2 injection (flow rate of 1 sccm). (a,b) AFM images. (c,d) Size distribution extracted from AFM images. The blue lines correspond to the fitting to log-normal distributions. (e-f) TEM images.

Figure 5 .

 5 Figure 5. a) IR spectrum of the dust analogue. b) Spectral deconvolution in the range 1500-1750 cm -1 . For the deconvolution, a baseline has been subtracted to remove the contribution from the very broad amorphous carbon band.

Figure 6 .

 6 Figure 6. a) IR spectra of the dust analogues in UHV and after air exposure for 15 min and 24 hours. The spectra have been vertically shifted for clarity. b) Temporal evolution of the normalized optical depth for some selected IR bands after exposing to air the dust analogue.

Figure 7

 7 Figure 7. a) Ex-situ laser desorption/ionization massspectrometry (LDI-MS) of the dust analogue using the AROMA setup (Sabbah et al. 2017). b) Double Bond Equivalent (DBE) analysis. The stacked bar graph in the inset summarizes the composition of each molecular family.

  Formation of polyacetylenic chains (C 2n H 2 )

Figure 11 .

 11 Figure 11. Reaction pathway for the formation of alkyl substituted aromatics through hydrogen abstraction and subsequent methyl addition.

4 H 2 ) 74 (C 6 H 2 ) 78 (C 6 H 6 ) 76 (C 6 H 4 ) 44 (CO 2 ) 18 (H 2 O) 28 (CO) 43 (CH 3 COCH 3 )

	Ion current (A)	10 -12 10 -11 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5	H 2 a) C 2n H 2 * H 2 O Ar 2+ C 2 H 2	Ar	*	C 4 H 2	*	0 sccm 4×10 -4 sccm C 6 H 2 Ar 2 0.15 sccm 1 sccm
			0 10 20 30 40 50 60 70 80 90 100
	Ion current (pA) Ion current (pA)	1.0 2.0 3.0 4.0 5.0 1.2 1.4 1.6 1.8	45 b) c)	50 C 4 H 2 m/z C 6 H 6 C 6 H 4 C 6 H 2	m/z 55 Ar 2	d)	1.2 1.3 3 6 9 150 200 250 300 50 (C Ion current (pA) 26 (C 2 H 2 ) 2 (H 2 ) 2 4 6
	Intensity x 10 3 (cps)	20 40 60 80 70	C 5 H n e)	C 6 H 4 C 6 H 6 Ar 2 C 6 H 2 75 m/z	C 7 H 8	C 8 H 2 80	OFF C 8 H C 8 8 H 10 C 9 H n C 8 H 6	ON 10 H C 8	C 11 H n	OFF C 12 H 8 C 12 H 10	0 2 4 6 8	Intensity x 10 3 (cps)
			60 70 80	90 100 110 120 130 140 150 160
												m/z

  (Deschenaux et al. 1999) we have detected aromatic molecules with aliphatic substitutions (toluene (C 7 H 8 ), styrene (C 8 H 8 ) and xylene/ethylbenzene (C 8 H 10 )

Table 1 .

 1 Chemical species detected in the gas phase by mass spectrometry during the production of the dust analogues along with its tentative assignment to specific compounds.

	m/z	Chemical formula	Compound
		Polyacetilenic chains
	26	C2H2	acetylene
	50	C4H2	diacetylene
	74	C6H2	triacetylene
	98	C8H2	octatetrayne
		Aromatics
	78	C6H6	benzene
	91/92	C7H8	toluene
	102	C8H6	phenylacetylene
	104	C8H8	styrene
	106	C8H10	xylene/ethylbenzene
	128	C10H8	naphthalene
	152	C12H8	acenaphthylene
	153/154	C12H10	biphenyl/vinylnaphthalene
		Aliphatics
	62-68	C5Hn	C5-clusters
	76	C6H4	3-hexene-1,5-diyne
	115-123	C9Hn	C9-clusters
	140-144	C11Hn	C11-clusters

Table 2 .

 2 Vibrational band assignments for the IR bands with optical depth higher than 0.0025.

	ν	λ	Assignment	Intensity
	cm -1	µm
	3322	3.01	sp alkyne CH str.	sh
	3300	3.03	sp alkyne CH str.	s
	3100	3.23	sp 2 arom./olef. CH str.	sh
	3085	3.24	sp 2 arom./olef. CH str.	sh
	3056	3.27	sp 2 arom./olef. CH str.	m
	3029	3.30	sp 2 arom./olef. CH str.	m
	2973	3.36	sp 3 asym. CH3 str.	m-s
	2925	3.42	sp 3 asym. CH2 str.	s
	2870	3.48	sp 3 sym. CH3 str.	sh
	2845	3.51	sp 3 sym. CH2 str.	sh
	2200	4.55	sp alkyne CC str. (conjugated alkynes)	w
	2105	4.75	sp alkyne CC str	w
	1938	5.16	Combination of out-of-plane and in-plane CH bend. modes	w
	1800-1000 5.56-10.00 Amorphous carbon (π-band)	m-s
	1711	5.84	Carbonyl CO str.	w
	1686	5.93	Carbonyl CO str.	w
	1630	6.13	sp 2 olef. CC str.	sh
	1600	6.25	sp 2 arom. CC str.	m-s
	1578	6.34	sp 2 arom. CC str.	sh
	1492	6.70	sp 2 arom. CC str.	w
	1445	6.92	sp 3 CH3 asym. bend. / sp 3 CH2 sci.	m-s
	1437	6.96	sp 2 arom. CC str	sh
	1379	7.25	sp 3 CH3 sym. bend.	sh
	1369	7.30	sp 3 CH3 sym. bend.	w-m
	1074	9.31	sp 2 arom. CH in-plane bend. / sp 2 olef. CH2 rock.	w
	1029	9.72	sp 2 arom. CH in-plane bend.	w
	988	10.12	sp 2 olef. CH wag.	w
	968	10.33	sp 3 CH3 rock.	w
	910	10.99	sp 2 olef. CH wag.	sh
	890	11.24	sp 2 arom. CH out-of-plane bend. (solo)	m
	843	11.86	sp 2 arom. CH out-of-plane bend. (duo)	m
	776	12.89	sp 2 arom. CH out-of-plane bend. (trio)	sh
	756	13.23	sp 2 arom. CH out-of-plane bend. (quartet)	m-s
	700	14.29	sp 2 CH out-of-plane bend.	m
	640	15.63	sp 2 arom. CCC in-plane bend.	m
	Note-arom: aromatics; olef: olefins. The vibrational modes are abbreviated as str: stretching;
	bend: bending; sci: scissoring; rock: rocking; wag: wagging. Relative intensities are labelled as s:
	strong; m: medium; w: weak; sh: shoulder.

  estimated an upper limit of around 15 % for the Figure 10. Reaction pathway for the formation of PAHs.

				C T H 4
	C2H2		C 2 H 2			C2H2	C 10 H 7	C2H2
	C 2 H	C 4 H 3	C 2 H 2	H C T H 5	C 2 H 2	C 8 H 6	H	C 12 H 8
				H	C 6 H 6			HACA
							C 10 H 8
								PAHs
							HACA
				C 6 H 6	C 12 H 10

  2n+1 + C 2 H 2 → C 2n+3 + H 2 (R11b)as well as by the following interactions with C 2 H and C 4 HC 2n+1 + C 2 H → C 2n+3 + H (R12a) C 2n+1 + C 4 H → C 2n+5 + H (R12b).

		C 2	C 2 H 2 /C 2 H	C 4	C 2 H 2 /C 2 H	C 6	C 2 H 2 /C 2 H	C 2n
	C , A r	C , A r	C , A r	C , A r	C , A r	C , A r	C , A r	C , A r
	C	C 2 H 2 /C 2 H	C 3	C 2 H 2 /C 2 H	C 5	C 2 H 2 /C 2 H	C 7	C 2 H 2 /C 2 H	C 2n+1
			C 4 H			C 4 H			C 4 H
		C 2 H 2			C 2 H 2			C 2 H 2	C 2 H 2
			C 3 H			C 5 H			C 7 H	C 2n+1 H
			H, Ar		H, Ar		H, Ar	H, Ar
			C 3 H m		C 5 H m			C 7 H m	C 2n+1 H m
	Figure 12. Reaction pathway for the formation of C2n+1Hm
	clusters.

C

ACKNOWLEDGMENTS

We thank the European Research Council for funding support under Synergy Grant ERC-2013-SyG, G.A. 610256 (NANOCOSMOS). Also, partial support from the Spanish Research Agency (AEI) through grants MAT2017-85089-c2-1R, FIS2016-77578-R, FIS2016-77726-C3-1-P, AYA2016-75066-C2-1-P, and RyC-2014-16277 is acknowledged. Support from the FotoArt-CM Project (P2018/NMT 4367) through the Program of R&D activities between research groups in Technologies 2013, co-financed by European Structural Funds, is also acknowledged. G.T-C. acknowledges funding from the Comunidad Autónoma de Madrid (PEJD-2018-PRE/IND-9029).