# UFL Dual Spaces, a proposal

David Ham

(January 2021)

###### Abstract

This white paper highlights current limitations in the algebraic closure
Unified Form Language (UFL). UFL currently represents forms over finite
element spaces, however finite element problems naturally result in objects in
the dual to a finite element space, and operators mapping between primal and
dual finite element spaces. This document sketches the relevant mathematical
areas and proposes changes to the UFL language to support dual spaces as first
class types in UFL.

##  1 Introduction

UFL is the Unified Form Language. It is a language used to describe
multilinear forms. Typically these forms are to be used in a finite-
dimensional variational problem, specifically a finite element problem. UFL
has not, hitherto, concerned itself with the spaces to which forms belong or
into which they map: this has been the concern of the assembler, such as
Firedrake or FEniCS and, to some extent, the linear algebra backend (PETSc in
the relevant cases, though other choices would be possible). These spaces
involve the duals to the spaces in which UFL objects such as Coefficients and
Arguments are defined. The absence of dual spaces in UFL means that the
language is not closed: there are operations visible to the UFL user which
result in objects that are not UFL objects, or which rely on bespoke functions
where a more closed system would have a first class operator. Examples of this
include:

  1. 1.

Interpolation, which is a form whose range space is the dual of the space
being interpolated into.

  2. 2.

Nonlinear constitutive laws for materials such as ice. In these cases, UFL-
based systems currently rely on an unnecessary Galerkin projection to
represent an operation which is not a differential operator.

  3. 3.

The adjoint, which currently has to rely on intricate reasoning to work out
which objects are primal or dual, where this should be obvious from their
type.

This document is an outline proposal for the inclusion of dual spaces and
related objects in UFL. It has particular implications for ongoing work in
Firedrake on interpolation and operators external to the UFL system.

##  2 Preliminaries

We will assume that there is a domain $\Omega\subset\mathbb{R}^{n}$ for some
$n\geq 0$ which is the region on which we will calculate.

######  Definition 2.1 (Function space)

A _function space_ is a vector space over a field K comprising functions
$\Omega\rightarrow K^{m}$.

That is to say, if $V$ is a function space and $f\in V$ then $f$ associates a
value in $K^{m}$ with each point in the domain $\Omega$. The field $K$ will be
either $\mathbb{R}$ or $\mathbb{C}$.

######  Definition 2.2 (linear $k$-form)

A _linear $k$-form_ is a linear mapping from $k\geq 0$ vector spaces to $K$.
That is, a linear $k$-form is a function:

$V_{k-1}\times V_{k-2}\times\ldots\times V_{0}\rightarrow K$ (1)

for suitable vector spaces $V_{0}\ldots V_{k-1}$.

###### Note 2.1

The argument numbering in definition 2.2 is deliberately from right to left.
This reflects the numbering of arguments in UFL. It also caters for a
difference in notation between forms and linear algebra: the assembled tensors
are the transpose of the form, so the function spaces number from left to
right when we consider the eventual linear algebra operations.

######  Definition 2.3 (sesquilinear $k$-form)

A _sesquilinear $k$-form_ is a mapping from $k\geq 0$ vector spaces to $K$
which is conjugate linear (also termed antilinear) in the final argument (for
$k>0$), and linear in the others.

######  Definition 2.4 (Inner product)

An inner product defined on a vector space $V$ is a sesquilinear 2-form
$V\times V\rightarrow K$. Which is Hermitian symmetric and positive definite.
That is, for any $u,v\in V$

  1. 1.

$\langle u,v\rangle=\overline{\langle v,u\rangle}$

  2. 2.

$\langle u,u\rangle\geq 0$

######  Definition 2.5 (Hilbert space)

A Hilbert space is a complete normed vector space equipped with an inner
product. All of the function spaces we are interested in are Hilbert spaces.

######  Definition 2.6 ($L^{2}$ inner product)

The $L^{2}$ inner product between two functions $u,v:\Omega\rightarrow K$ is
given by:

$\langle u,v\rangle_{L^{2}}=\int_{\Omega}u\cdot\overline{v}\ \mathrm{d}x$ (2)

######  Definition 2.7 ($L^{2}$)

The function space $L^{2}(\Omega)^{n}:\Omega\rightarrow K^{n}$ is the Hilbert
space of functions where the norm induced by the $L^{2}$ inner product is
finite:

$u\in L^{2}(\Omega)^{n}\Leftrightarrow\sqrt{\langle
u,u\rangle_{L^{2}}}\in\mathbb{R}$ (3)

######  Definition 2.8 (Dual space)

If $V$ is a vector space then $V^{*}$, the dual space to $V$, is the space of
bounded linear functionals $V\rightarrow K$.

For the finite-dimensional spaces with which we are concerned, it is easy to
establish the following:

###### Theorem 2.1

If $V$ is a finite-dimensional vector space, then $|V|=|V^{*}|$

######  Theorem 2.2 (Dual basis)

If $\\{\phi_{i}\\}$ is a basis for a finite vector space $V$ then there exists
$\\{\phi_{i}^{*}\\}$ a basis for $V^{*}$ such that:

$\phi_{i}^{*}(\phi_{j})=\delta_{ij}$ (4)

The basis $\\{\phi_{i}^{*}\\}$ is termed the _dual basis_.

Where it is necessary to make the distinction, we will refer to the space to
which a dual space is dual as the _primal space_ and its basis as the _primal
basis_.

###### Definition 2.9

Reflexive vector space A reflexive vector space is one for which $(V^{*})^{*}$
is isomorphic to $V$ under the canonical map. That is, we can identify
$(V^{*})^{*}$ and $V$.

###### Theorem 2.3

All finite-dimensional vector spaces are reflexive.

This is a completely intuitive result. If $u\in V$ and $v\in V^{*}$ then this
induces a $u^{\prime}$ in $(V^{*})^{*}$ given by:

$u^{\prime}(v)=v(u)$ (5)

We simply identify $u^{\prime}$ and $u$. Since the dimensions of the spaces
match and the functionals are linear, this identification is an isometric
isomorphism. In other words we can’t distinguish the spaces so the
identification is valid.

###### Theorem 2.4

All Hilbert spaces are reflexive.

##  3 Functional and tensor product notation

In definition 2.2 and elsewhere, we use a notation for multilinear forms which
takes a Cartesian product of spaces and returns a scalar. For example for a
2-form:

$V_{1}\times V_{0}\rightarrow K$ (6)

This notation emphasises that the form is a function which takes two vectors
and returns a scalar. However we could just as well Curry this form and write:

$V_{1}\rightarrow V_{0}\rightarrow K$ (7)

This notation is left-associative, so this means:

$V_{1}\rightarrow(V_{0}\rightarrow K)$ (8)

In other words, a 2-form is a function which takes a vector and returns a
1-form.

##  4 Dual spaces and $k$-forms

In this section we restrict our attention to real-valued function spaces so as
not to have to consider where all the conjugates should be placed. We will
come back to the complex valued case.

It follows immediately from definitions 2.2 and 2.8 that the space of linear
1-forms on vector space $V$ is the same thing as $V^{*}$. This is nothing else
than writing:

$V^{*}=V\rightarrow\mathbb{R}$ (9)

A direct consequence of this is the following equivalence:

$V_{1}\times V_{0}\rightarrow\mathbb{R}=V_{1}\rightarrow V_{0}^{*}$ (10)

This says that every form can equivalently be thought of as an operator
mapping into the dual of its 0-number space.

##  5 The interpolation operator

Usually we are concerned with vector spaces which are finite-dimensional
subspaces of larger, typically infinite dimensional, spaces. In the finite
element method, we often define the basis of these finite spaces by defining
the space and choosing the dual basis. The primal basis can then be derived by
solving a linear system (see, for example, the finite element course). Let’s
briefly examine the dual basis functions, also known as nodes. These are
functionals that take in a function and (linearly) produce a scalar. What do
such objects look like? The simplest case is probably the point evaluation
functional, or Dirac delta distribution. For a function
$v:\Omega\rightarrow\mathbb{R}$, the point evaluation functional at the point
$x$ is given by:

$\delta_{x}(v)=v(x)$ (11)

Other common functionals include the integral of a function over one mesh
cell, or the component of a vector value in a particular direction at a
particular point. All of these have in common that they are trivially
extensible to a much larger space of functions than the finite element
function space whose dual they define. For example, point evaluation
functionals can be evaluated for any space of functions whose members are
well-defined at the points in question. This is useful because int is a very
common requirement to approximate some function from a more general space in a
finite element function. This most commonly occurs when a forcing function or
reference solution is known either from external (physical) data, or from an
analytic expression. Suppose $v\in V$ and that $\\{\phi_{i}\\}$ is a basis for
$V$, with $\\{\phi_{i}^{*}\\}$ the corresponding dual basis. The for a
function $f:\Omega\rightarrow\mathbb{R}$ such that the dual functionals are
defined, the interpolation of $f$ into $v$ is given by:

$\hat{f}=\hat{f}_{i}\phi_{i}$ (12)

where:

$\hat{f}_{i}=\phi_{i}^{*}(f)$ (13)

In other words, the $i$-th basis coefficient is the evaluation of the $i$-th
dual basis functional on the function to be interpolated. We can call this
operator $I_{V}(f)$.

###### Note 5.1

If $f\in V$ then it is straightforward to show that $I_{V}(f)=f$. That is to
say, interpolation is a projection.

In fact, the interpolation operator is nothing other than the evaluation of
the dual basis. We could instead define an operator, a generalisation of the
delta distribution, which takes a functional $v\in V^{*}$ and applies it to a
function $u\in U$ where $U$ is a function space on which the functionals in
$V^{*}$ are defined. That is to say:

$\delta(u,v)=v(u)$ (14)

Notice that $\delta:U\times V^{*}\rightarrow\mathbb{R}$. Since $V^{*}$ is
itself a vector space, this makes $\delta$ a 2-form, albeit a slightly unusual
one. We would be more used to forms in $U\times V\rightarrow\mathbb{R}$.

##  6 UFL Forms

UFL is the Unified Form Language [1] also described in [2]. As one might
expect, it provides support for $k$-linear forms. UFL is a symbolic algebra
language. In common with other such languages, it provides for symbols to
stand for unknown values, including the arguments to $k$-linear forms. UFL is
a purely symbolic language which problem solving environments such as
Firedrake and FEniCS extend by subclassing in order to attach numerical data
and actually solve problems. In understanding how UFL objects do (or should)
behave, it is sometimes necessary to consider how they would be integrated
into the problem solving environment, and especially what should happen when a
form is assembled.

In an attempt to avoid replicating the documentation of UFL, we will start
here from UFL `FunctionSpace` and gloss over the cells, meshes, and finite
elements from which these are constructed.

If `V` is a `ufl.FunctionSpace` then `c = ufl.Coefficient(V)` defines a
_known_ function in `V`. Conversely, `a = ufl.Argument(V, 0)`111Users more
frequently write ``a = ufl.TestFunction(V)—, but that is simply syntactic
sugar. A test function is just argument 0 to a form. defines a placeholder
symbol for an _unknown_ function in `V`. Consequently, `f_0 = c * dx` is a
symbolic expression for the integral of `c` over the domain and is a real
value. UFL calls this a 0-form. In contrast, `f_1 = a * dx` represents the
integration of the unknown function `a` over the domain. It’s therefore a
1-form, or a function in $V\rightarrow\mathbb{R}$.

In order to make the argument order clear, we introduce a second function
space, `W`. We define `b = ufl.Argument(W, 1)`222Equivalently, ``b =
ufl.TrialFunction(W)—. and create the simplest 2-form: `f_2 = c * b * dx` This
is a function in $W\times V\rightarrow\mathbb{R}$. Notice that, by convention,
argument 1 comes before argument 0.

###  6.1 Assembly

The problem solving environment takes symbolic UFL forms and evaluates them.
In Firedrake and FEniCS, the operation which does this is called assembly. The
problem solving environment encodes in its subclass of `ufl.FunctionSpace` a
basis for the function space. We will write ${\phi_{i}}$ for the basis for $V$
and ${\psi_{i}}$ for the basis for $W$. The problem solving environment will
also implement its own subclass of `ufl.Coefficient`, in Firedrake and FEniCS
this is called a `Function`. If `c` is a `Coefficient`, then mathematically:

$c=c_{i}\phi_{i}$ (15)

The consequence of this is that the core of the implementation of a `Function`
is a reference to a `FunctionSpace` (which encodes $\\{\phi_{i}\\}$) and a
vector of values $c_{i}$. The latter give the `Function` its numerical value,
while the former gives that value a mathematical meaning.

If we execute:

⬇

g = assemble(f_0)

Then the integral is evaluated and $t_{0}$ is a scalar value. Mathematically,
nothing has really happened: $g$ and $f_{0}$ represent the same scalar value.
Computationally, however, we have gone from a symbolic object which evaluates
to a scalar, to the scalar itself.

So, what happens in the following case?

⬇

h = assemble(f_1)

`f_1` contains the argument, or unknown `Function`, `a`. In other words, we
need to evaluate

$\begin{split}h(a)&=\int a\,\mathrm{d}x\\\ &=\int
a_{i}\phi_{i}\,\mathrm{d}x\end{split}$ (16)

for unknown values $a_{i}$. However the $a_{i}$ are scalar so they can be
lifted out of the integral:

$\int a\,\mathrm{d}x=\int\phi_{i}\,\mathrm{d}x\ a_{i}$ (17)

Observe that we can now evaluate all the integrals, as they no longer contain
unknown values. Next, we use the fact that
$\phi^{*}_{i}(\phi_{j})=\delta_{ij}=\mathrm{I}$, the identity matrix. Hence:

$\begin{split}h(a)&=\int\phi_{i}\,\mathrm{d}x\ \mathrm{I}_{ij}a_{j}\\\
&=\int\phi_{i}\,\mathrm{d}x\ \phi^{*}_{i}(\phi_{j})a_{j}\\\
&=\int\phi_{i}\,\mathrm{d}x\ \phi^{*}_{i}(a_{j}\phi_{j})\\\
&=\int\phi_{i}\,\mathrm{d}x\ \phi^{*}_{i}(a)\\\ \end{split}$ (18)

We used the linearity of the dual basis functionals for the penultimate line.
We have now expressed $h$ in terms of the basis for $V^{*}$:
$h=h_{i}\phi^{*}_{i}$ where:

$h_{i}=\int\phi_{i}\ \mathrm{d}x$ (19)

The consequence of this is that the value of $h$ is encoded in a vector of
numbers, of the same size as that which encodes $c$. Now one could say that
the only purpose of assembly is to create a linear system to be solved, so the
vector of values is the only thing which is required. This is the current
FEniCS approach. This approach is problematic because it neglects the fact
that assembled operators can and are used by users in other contexts.

What is missing, of course, is a reference to the dual basis
$\\{\phi^{*}_{i}\\}$. The challenge is that UFL does not have a class encoding
$\\{\phi^{*}_{i}\\}$, and hence also does not have a suitable class for $h$.
The Firedrake approach is to appeal to the natural isomorphism between $V$ and
$V^{*}$ and simply make `h` a `Function`. This has the advantage of preserving
the link between the numerical values and the function space, but at the
expense of incorrectly labelling $t_{1}$ as being in $V$. What we should
really do is associate the list of coefficients $h_{i}$ with an object
encoding the dual basis, in a manner analogous to the association between a
list of coefficients and the basis for a function space in a `Function`.

To illustrate the manner in which the assembled form encodes the same
information as the symbolic form, observe that:

$h(c)=f_{1}(c)=\int c\,\mathrm{d}x=f_{0}=g$ (20)

However the evaluation of $h(c)$ is actually accomplished like this:

$\begin{split}h(c)&=h_{i}\phi^{*}_{i}(c_{j}\phi_{j})\\\
&=h_{i}\phi^{*}_{i}(\phi_{j})c_{j}\\\ &=h_{i}\delta_{ij}c_{j}\\\
&=h_{i}c_{i}\end{split}$ (21)

So the assembly has turned the evaluation of an integral into the dot product
between two vectors.

###  6.2 UFL dual spaces and cofunctions

The previous section suggests an outline of at least some missing UFL
functionality. There should be an object representing the dual to a function
space (`DualFunctionSpace` or maybe just `DualSpace`?) Presumably one would
add a property to `ufl.FunctionSpace` to access the dual, and of course
`DualSpace` would also have this property but it would point back to the
primal space. Maybe `ufl.FunctionSpace` and `ufl.DualSpace` have a common
parent class `ufl.VectorSpace`.

We also need to be able to define a concrete member of $V^{*}$. Consistency
with the existing UFL usage would suggest that this class would be
`ufl.Cocoefficient` but this is a ridiculous and cumbersome name so perhaps
`ufl.Cofunction` would be better. Firedrake and FEniCS subclasses of
`ufl.Cofunction` would, in a manner analogous to `Function` have both a vector
of values and a reference to the `ufl.DualSpace`. The question of which
operations should be valid for a cofunction requires more thought, though
these would certainly include the basic vector operations of addition and
multiplication by a scalar.

###  6.3 Cofunctions are forms

A cofunction is mathematically just a 1-form. The distinction between this and
an existing UFL 1-form is that the latter is a symbolic expression that
represents a 1-form, while the former is a numerical 1-form represented
directly in the dual basis. Since they are mathematically objects of the same
type, they should be arithmetically interchangeable. In particular, if `f =
ufl.Cofunction(V.dual)` and `v = ufl.Argument(V, 0)` then the following makes
mathematical sense as a form and should be allowed:

⬇

g = f + v*dx

This possibly implies that cofunctions intrinsically have an argument drawn
from the relevant primal space. Assembly of `g` can be achieved by simply
defining assembly of a cofunction as the identity operation.

###  6.4 2-forms

What do we get from:

⬇

A = assemble(f_2)

Mathematically this is defined as follows:

$A_{ij}=\int\phi_{i}\psi_{j}\ \mathrm{d}x$ (22)

That is to say, $A$ is a matrix. Recall that $\\{\phi_{i}\\}$ is the basis for
$V$, which is the function space of argument 0 in this case, and
$\\{\psi_{j}\\}$ is the basis for $W$, which is argument 1. This means that
the assembled operator has its arguments in ascending order, while the form
notation has them in descending order. This is a frequent cause of confusion
about which we will have to remain aware. FEniCS and Firedrake both agree that
assembling a 2-form produces a matrix, though UFL does not have a
corresponding class, so currently matrices are outside its scope. However, a
matrix is an assembled 2-form in just the same way that a cofuncion is an
assembled 1-form, so the logical conclusion would be the creation of a
`ufl.Matrix` class which would behave as an assembled 2-form, and would have
algebraic properties analogous to a cofunction. A matrix would presumably have
two intrinsic arguments.

##  7 The delta operator and dual space arguments

Let’s now return to the interpolation operator. If we wish to incorporate it
into UFL as a first class citizen then we need to take our extensions of the
language just a little further. Let’s suppose we can define arguments in dual
spaces. These are mathematical objects of a different kind to the arguments we
have already met: they are unknown members of $V^{*}$ rather than $V$. We
suppose, therefore, that they are of a new type: `ufl.Coargument`. In order to
minimise syntactic clutter, one might imagine overloading the object
constructor such that `ufl.Argument(V.dual, 0)` is valid and returns a
`ufl.Coargument`.

This enables us to define a new form `ufl.Delta` such that if `e` is a UFL
expression and `v_ = ufl.Argument(V.dual, 0)` then `delta(e, v_)` is a
symbolic expression for the interpolation of `e` into $V$. If we assume that
`e` does not contain any arguments and execute:

⬇

f = assemble(Delta(e, v_))

Then `f` will be a `Function` in $V$ representing the result of this
interpolation. Conversely, if `u = Argument (W, 1)` then:

⬇

A = assemble(Delta(u, v_))

results in a matrix which is the interpolation operator from $W$ into $V$. One
can also construct the adjoint of this operator by reversing the argument
positions:

⬇

v = ufl.Argument(W, 0)

u_ = ufl.Argument(V.dual, 1)

AT = assemble(ufl.Delta(v, u_))

As a particular case of interest, if we write:

⬇

u = ufl.Argument(V, 1)

v_ = ufl.Argument(V.dual, 0)

I = assemble(ufl.Delta(u, v_))

Then `I` is the identity matrix for $V$

For completeness, we can also address the other cases. If we write:

⬇

f = firedrake.Cofunction(V.dual)

c = firedrake.Function(W)

t = ufl.Delta(c, f)

s = assemble(t)

then $t$ is a symbolic expression for $f(c)$ and $s$ is the scalar which
results from evaluating it. Indeed, one might want to implement the `__call__`
special method on `ufl.Cofunction` and `ufl.Coargument` so that the syntax
`f(c)` works directly. Note that this would simply return `ufl.Delta(c, f)`.
`ufl.Delta` is still needed as the object which encodes the symbolic operation
of dual evaluation.

Finally, if we write:

⬇

f = firedrake.Cofunction(V.dual)

u = ufl.Argument(W, 0)

g = assemble(ufl.Delta(u, f))

Then $g\in W^{*}$ is the interpolation of $f$ into $W^{*}$. (Note that
Cofunction interpolation is the transpose of function interpolation).

##  8 Adjoint

The dual or adjoint of a 2-form is the conjugate of that form with its
arguments relabelled in reversed order. This carries over to forms which map
from or to dual spaces. If we take the example of the delta form again, then:

$\delta:V\times W^{*}\rightarrow\mathbb{R}.$ (23)

Equivalently:

$\delta:V\rightarrow W.$ (24)

If we write $\delta^{*}$ for the adjoint to $\delta$ then:

$\delta^{*}:W^{*}\times V\rightarrow\mathbb{R},$ (25)

which is to say:

$\delta^{*}:W^{*}\rightarrow V^{*}.$ (26)

Which is simply another way of saying that cofunction interpolation is the
dual of coefficient interpolation.

##  9 External operators

The delta form is simply one example of an form whose first argument is in a
dual space. In fact any operator which outputs a `Function` and which is
linear in any other arguments it has is a form of similar type. Note that such
an operator could have known inputs (e.g. coefficients) with respect to which
it is nonlinear, it is only the unknown arguments with respect to which it has
to be linear. In this context it is important to distinguish between the
inputs to such an operator, which we will call its _operands_ and arguments.
Arguments may appear in operands, but an operand also contain other UFL
expressions. For example an operator $N$ might take some UFL expression as an
input, evaluate a neural net and return the output as a `Function`. That is
$N(e,m)\in V$ for a UFL expression $e$, $m\in\mathbb{R}^{k}$ and $V$ some
function space. Of course saying that $N(e,m)\rightarrow V$ is the same as
saying $N(e,m)\times V^{*}\rightarrow\mathbb{R}$, so an external operator is
also a form.

##  10 Form composition

Recall that a $k$-form of type $V_{k-1}\times\ldots\times
V_{0}\rightarrow\mathbb{R}$ can be interpreted as an operator
$V_{k-1}\times\ldots\times V_{1}\rightarrow V_{0}^{*}$. In other words, it can
mathematically be used anywhere were an object in $V_{0}^{*}$ is expected.
Similarly, the reflexivity of the vector spaces we work with means that a
$k$-form of type $V_{k-1}\times\ldots\times V_{0}^{*}\rightarrow\mathbb{R}$
can be interpreted as an operator $V_{k-1}\times\ldots\times V_{1}\rightarrow
V_{0}$. That is, it can be used where a value in $V_{0}$ is expected.

Hitherto, UFL forms have only taken inputs from primal spaces and produced
outputs in dual spaces, so composing forms by using the output of one as an
input to another was not an available option. However, now that there are
forms which produce primal space outputs, and forms which accept dual space
inputs, we need to examine the consequences of this.

Consider this simple case:

⬇

v_ = ufl.Argument(V.dual, 0)

f = firedrake.Function(W)

v = ufl.Argument(V, 0)

F_1 = ufl.Delta(f, v_) * v * dx

l = assemble(F_1)

This is a slightly contrived case, though it could really occur if, for
example, $W$ is defined on a different mesh from $V$. The assembly algorithm
for `l` would be:

⬇

tmp = firedrake.Function(V)

l = assemble(ufl.replace(F_1, {ufl.Delta(f, w_): tmp}))

I.e., we first evaluate the delta, and then substitute the result into the
surrounding form. The type of the temporary variable can be inferred from the
arguments to the delta. Consider now the 2-form case:

⬇

v_ = ufl.Argument(V.dual, 0)

w = ufl.Argument(W, 1)

v = ufl.Argument(V, 0)

F_2 = ufl.Delta(w, v_) * v * dx

A = assemble(F_2)

This results in the following assembly algorithm:

⬇

tmp_1 = assemble(Delta(w, v_)) # tmp_1 is a |V| x |W| matrix.

tmp_2 = assemble(ufl.replace(F_2, {ufl.Delta(w, v_): ufl.Argument(V, 1)}))

A = tmp_2 @ tmp_1

Again the, the types of the temporary variables can be inferred from the
arguments of the delta. The delta has 2 arguments and they are in the
canonical order, so it is replaced by an argument numbered 1 but in the
function space dual to the argument 0 function space. The sequence of the
matrix multiplication is similarly inferable from the arguments and the
sequence of composition (the delta occurs inside the other form, and not vice
versa). The general rule is that the form is replaced by an argument in the
dual space to its zero argument, and with an argument number which is the next
available argument number in the form into which it is being substituted.

###  10.1 Arguments of composite forms

What are the arguments of `F_2`? One might be tempted to say `v_`, `w`, and
`v`, but this would characterise `F_2` as a 3-form rather than a 2-form. It
would also have two different arguments numbered 0 (`v_`, and `v`), so it
would not even be a well-formed form. The answer to this puzzle is that the 0
argument of a form nested into another form is consumed by the substitution.
In this case, this means that the arguments of the composite form are `v` and
`w`.

###  10.2 No substitution for argument 0

Consider instead the case where the interpolation happens for the test
function. By symmetry with `F_2`, one might expect to be able to write:

⬇

w_ = ufl.Argument(W.dual, 0)

v = ufl.Argument(V, 0)

u = ufl.Argument(V, 1)

F_3 = u * ufl.Delta(v, w_) * dx

However this is ill-formed since the delta contains two arguments numbered 0.
This difficulty is the result of the fact that form composition only makes
sense if we think of forms as operators into the dual space of argument 0. In
this context, argument 0 is not a form input at all, but is instead a
placeholder for the result of evaluating the form for known values of all of
its higher-numbered arguments. However, composition of linear operators works
in both directions so we can instead write:

⬇

w = ufl.Argument(W, 0)

v = ufl.Argument(W, 0)

u = ufl.Argument(V, 1)

F_4 = ufl.Delta(w, u * v * dx)

Which is well-formed UFL, and mathematically equivalent to what `F_3`
attempted to achieve. As one might expect, $F_{4}:V\rightarrow W^{*}$.

##  11 Action

The action of a form on a coefficient or argument simply replaces the highest
numbered argument with that coefficient. This carries over directly to the
dual case. If the highest numbered argument of a form lies in a dual space,
then the action is only well defined on a cofunction or coargument from the
appropriate space.

Form composition means that it should also be possible to take the action of a
form on another form, so long as the space of the highest numbered argument of
the first form is dual to the space of the 0 argument to the second form. That
is, an alternative definition of `F_4` is given by:

⬇

w = ufl.Argument(W, 0)

omega = ufl.Argument(W.dual, 1)

v = ufl.Argument(V, 0)

u = ufl.Argument(V, 1)

F_4 = ufl.action(ufl.Delta(w, omega), u * v * dx)

Of course this lines up with the assembly strategy for `F_4`, which would be:

⬇

tmp_1 = assemble(Delta(w_, ufl.Argument(V.dual, 1))

tmp_2 = assemble(u * v * dx)

A = tmp_1 @ tmp_2

## References

  * [1] Martin S Alnæs et al.  “Unified form language: A domain-specific language for weak formulations of partial differential equations”  In _ACM Transactions on Mathematical Software (TOMS)_ 40.2 ACM New York, NY, USA, 2014, pp. 1–37
  * [2] Miklós Homolya, Lawrence Mitchell, Fabio Luporini and David A Ham  “TSFC: a structure-preserving form compiler”  In _SIAM Journal on Scientific Computing_ 40.3 SIAM, 2018, pp. C401–C428


11institutetext: Max Planck Institute for Extraterrestrial Physics,
Giessenbachstrasse 1, D-85748 Garching, Germany
11email<EMAIL_ADDRESS>22institutetext: KU Leuven, Department of
Chemistry, Celestijnenlaan 200F, B-3001 Leuven, Belgium 33institutetext:
Laboratoire Aimé-Cotton, Bâtiment 505, CNRS, Université Paris-Saclay, F-91405
Orsay, France 44institutetext: Radboud University, Institute for Molecules and
Materials, FELIX Laboratory, Toernooiveld 7, 6525ED Nijmegen, the Netherlands
55institutetext: Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa,
Italy 66institutetext: Dipartimento di Chimica “Giacomo Ciamician”, Università
di Bologna, Via Selmi 2, 40126 Bologna, Italy

# Deuteration of c-C3H2 towards the pre-stellar core L1544

K. Giers 11 S. Spezzano 11 F. Alves 11 P. Caselli 11 E. Redaelli 11 O. Sipilä
11 M. Ben Khalifa 22 L. Wiesenfeld 33 S. Brünken 44 L. Bizzocchi 5566

###### Abstract

Context.  In the centre of pre-stellar cores, the deuterium fractionation is
enhanced due to the cold temperatures and high densities. Therefore, the
chemistry of deuterated molecules can be used to probe the evolution and the
kinematics in the earliest stages of star formation.

Aims. We analyse emission maps of cyclopropenylidene, c-C3H2, to study the
distribution of the deuteration throughout the prototypical pre-stellar core
L1544.

Methods. We use single-dish observations of c-C3H2, c-H13CC2H, c-C3HD, and
c-C3D2 towards the pre-stellar core L1544, performed at the IRAM 30 m
telescope. We derive the column density and deuterium fraction maps, and
compare these observations with non-LTE radiative transfer simulations.

Results. The highest deuterium fractions are found close to the dust peak at
the centre of L1544, where the increased abundance of H2D+ ions drives the
deuteration process. The peak values are N(c-C3HD)/N(c-C3H${}_{2})=0.17\pm
0.01$, N(c-C3D2)/N(c-C3H${}_{2})=0.025\pm 0.003$ and
N(c-C3D2)/N(c-C3HD$)=0.16\pm 0.03$, which is consistent with previous single
point observations. The distributions of c-C3HD and c-C3D2 indicate that the
deuterated forms of c-C3H2 in fact trace the dust peak and not the c-C3H2
peak.

Conclusions. The N(c-C3D2)/N(c-C3HD) map confirms that the process of
deuteration is more efficient towards the centre of the core and demonstrates
that carbon-chain molecules are still present at high densities. This is
likely caused by an increased abundance of He+ ions destroying CO, which
increases the amount of carbon atoms in the gas phase.

###### Key Words.:

astrochemistry – ISM: clouds – ISM: molecules – ISM: abundances – stars:
formation – ISM: individual objects: L1544

##  1 Introduction

Pre-stellar cores represent the initial conditions of low-mass star formation.
The self-gravitating cores (Ward-Thompson et al., 1999) are dense ($n>10^{5}$
cm-3) and cold ($T<10\,K$) towards their centre (Crapsi et al., 2007; Caselli
et al., 2012), showing signs of contraction motion and chemical evolution
(Crapsi et al., 2005). Therefore, pre-stellar cores are ideal laboratories to
study the physical and chemical processes taking place in the early stages of
low-mass star formation.

In the cold and dense conditions of pre-stellar cores, the deuteration of
molecules is a favoured chemical process (Ceccarelli et al., 2014), increasing
the abundance of deuterated species. The deuterium fractionation is driven by
the gas phase reaction (only true with respect to para states of reactants and
products; Pagani et al. 1992):

$\mathrm{H}_{3}^{+}+\mathrm{HD}\rightleftharpoons\mathrm{H}_{2}\mathrm{D}^{+}+\mathrm{H}_{2}+232\,\mathrm{K}.$
(1)

This reaction is exothermic. As the temperature decreases, the abundance of
H2D+ ions is enhanced, which drives the deuteration of other molecules. As the
reaction goes more efficiently from left to right at low temperatures, H2D+
will no longer be destroyed when meeting an H2 molecule, but mostly destroyed
by CO. However, in cold and dense regions like pre-stellar cores, CO is highly
depleted from the gas phase and frozen out onto dust grains (Willacy et al.,
1998; Caselli et al., 1999; Bacmann et al., 2003; Crapsi et al., 2005), which
further increases the deuterium fractionation (Dalgarno & Lepp, 1984). The
para-H2D+ is further deuterated by successive reactions with HD, leading to
noticeable densities of both D2H+ and D${}_{3}^{+}$ (Hily-Blant et al., 2018;
Caselli et al., 2019). All three deuterium substituted H${}_{3}^{+}$ ions are
driving the deuteration of other molecules, by transfer of a D+ atom.
Furthermore, the level of deuteration also depends on the ortho-to-para ratio
of H2. The higher the ratio, the more H2D+ is destroyed by ortho-H2, as the
energy of the ground state of ortho-H2 compared to para-H2, $\sim 170$ K, is
close to the exothermicity of reaction (1). However, the ortho-to-para ratio
of H2 in pre-stellar cores is usually $\sim$10-3 (Kong et al., 2015), further
promoting the deuterium fractionation. Deuterated molecules hence have proven
to be very good tracers of the central parts of dense cores that are on the
verge of star formation (Crapsi et al., 2007; Redaelli et al., 2019), as a
large deuterium enrichment is characteristic for these high density nuclei
(Crapsi et al., 2005). The deuterium fraction is measured by the ratio $R_{\rm
D/H}$ of the molecular tracer, dividing the column density of the D-bearing
molecule by the column density of the hydrogenated isotopologue. In addition,
the level of deuteration is sensitive to the evolutionary stage of starless
cores (Crapsi et al., 2005; Chantzos et al., 2018). Therefore, it can work as
a chemical clock and can be used as a diagnostic tool.

Cyclopropenylidene (c-C3H2) is one of the most abundant and widespread
molecules in the Galaxy. Its first detection in space was made by Thaddeus et
al. (1981), which was later confirmed by the first laboratory detection by
Thaddeus et al. (1985). Since then, c-C3H2 has been observed in diffuse gas,
cold dark clouds, giant molecular clouds, photodissociation regions,
circumstellar envelopes, planetary nebulae, circumstellar disks, and low-mass
protostars (Thaddeus et al., 1985; Matthews & Irvine, 1985; Seaquist & Bell,
1986; Vrtilek et al., 1987; Cox et al., 1987; Madden et al., 1989; Lucas &
Liszt, 2000; Teyssier et al., 2005; Qi et al., 2013; Majumdar et al., 2017).
As the formation of c-C3H2 is believed to happen solely in the gas phase (Park
et al., 2006), it preferentially traces dense and chemically young gas, which
is gas where the C atoms have not yet been mainly locked into CO (Spezzano et
al., 2016). However, it also suffers from depletion, and therefore it does not
trace the central regions of pre-stellar cores but more external layers
(Spezzano et al., 2017). Following the detection of its singly and doubly
deuterated isotopologues (Gerin et al., 1987; Spezzano et al., 2013), c-C3H2
is a unique probe to study the deuteration processes that take place in the
gas phase in pre-stellar cores, as it is only efficiently deuterated there
(Spezzano et al., 2013).

The low-mass pre-stellar core L1544 is located in the Taurus Molecular Cloud
complex, one of the nearest star forming regions, and has a distance of 170 pc
(Galli et al., 2019). The dense core is a well-known object; numerous and
extensive studies have led to a profound knowledge of its physical and
chemical structure (Crapsi et al., 2007; Keto et al., 2004; Keto & Caselli,
2010; Keto et al., 2015; Caselli et al., 2002a). It is centrally concentrated
(Ward-Thompson et al., 1999), with a high central density of around 106 cm-3
within the central 2000 AU (Crapsi et al., 2007; Caselli et al., 2019), and
very cold, with a temperature ranging from 12 K in the outskirts down to
$\sim$ 6 K towards the centre (Crapsi et al., 2007). The core has an elongated
shape and shows signs of contraction motion, suggesting that it is on the
verge of gravitational collapse (Williams et al., 1999; Ohashi et al., 1999;
Lee et al., 2001; Caselli et al., 2002a, 2012). Towards its centre, L1544
exhibits a high degree of CO freeze-out (93%; Caselli et al. 2002a; Crapsi et
al. 2005), and a high level of deuteration (Crapsi et al., 2005; Redaelli et
al., 2019). It is chemically rich, showing spatial inhomogeneities in the
distribution of molecular emission (Spezzano et al., 2017; Redaelli et al.,
2019; Chacón-Tanarro et al., 2019).

In this work we present emission maps of c-C3H2, the singly and doubly
deuterated isotopologues, c-C3HD and c-C3D2, as well as the isotopologue
c-H13CC2H towards L1544. In Sect. 2, we describe the observations, followed by
the results in Sect. 3. The analysis in Sect. 4 is divided into the derivation
of the c-C3H2 deuteration maps and a comparison of the dust emission peak
column density with non-LTE simulations. We discuss the results in Sect. 5 and
present our conclusions in Sect. 6.

##  2 Observations

The emission maps presented in this paper were obtained using the IRAM 30 m
single dish telescope located in Pico Veleta, Spain, in October 2013. The
observed $2.5^{\prime}\times 2.5^{\prime}$ on-the-fly (OTF) maps are centred
on the source dust emission peak
($\alpha_{2000}=05^{\mathrm{h}}04^{\mathrm{m}}17^{\mathrm{s}}.21$,
$\delta_{2000}=+25^{\circ}10^{\prime}42^{\prime\prime}.8$, Ward-Thompson et
al. 1999). We applied position switching, with the reference position set at
(-180′′,180′′) offset with respect to the map centre. The EMIR E090 receiver
was used along with the Fourier transform spectrometer backend (FTS), with a
spectral resolution of 50 kHz. The antenna moved along an orthogonal pattern
of linear paths separated by $8^{\prime\prime}$ intervals, corresponding
roughly to one third of the beam full width at half maximum (FWHM). The
weather conditions during the mapping were good ($\tau\sim 0.03$), with a
typical system temperature of $T_{\mathrm{sys}}\sim$ 90$-$100 K.

The observed transitions are summarised in Table 1. Figure 1 shows the
corresponding integrated intensity maps for each molecule. The data processing
was done using the GILDAS software (Pety, 2005) and Python. All emission maps
presented in this paper have been gridded to a pixel size of
$8^{\prime\prime}$ with the CLASS software in the GILDAS packages; this
corresponds to 1/3 -1/4 of the actual beam size, depending on the frequency.
The antenna temperature $T_{A}^{*}$ was converted to the main beam temperature
$T_{\mathrm{mb}}$ using the relation
$T_{\mathrm{mb}}=F_{\mathrm{eff}}/B_{\mathrm{eff}}\cdot T_{A}^{*}$. The
corresponding values for the 30 m forward ($F_{\mathrm{eff}}$) and main beam
efficiency ($B_{\mathrm{eff}}$) ratios are given in Table 1.

Table 1: Observed lines and observation details.

Molecule | Transition | Frequency aa$a$Extracted from Cologne database for molecular spectroscopy (Müller et al., 2001) | rms | HPBW | $F_{\mathrm{eff}}/B_{\mathrm{eff}}$ | Velocity range bb$b$Velocity ranges where the integrated emission has been computed. | Reference
---|---|---|---|---|---|---|---
| $(J_{K^{\prime}_{a}K^{\prime}_{c}}^{\prime}-J_{K_{a}K_{c}})$ | (GHz) | (mK) | (”) |  | (km s-1) |
c-C3H2 | $3_{22}-3_{13}$ | 84.728 | 22 | 30.6 | 0.95/0.838 | 6.90$-$7.42 | 1
c-H13CC2H | $2_{12}-1_{01}$ | 84.186 | 22 | 30.8 | 0.95/0.839 | 6.83$-$7.35 | 2
c-C3HD | $3_{03}-2_{12}$ | 104.187 | 27 | 24.9 | 0.95/0.802 | 6.78$-$7.48 | 2
c-C3D2 | $3_{13}-2_{02}$ | 97.762 | 21 | 26.5 | 0.95/0.814 | 6.89$-$7.49 | 3

111

(1) Thaddeus et al. (1985); (2) Bogey et al. (1987); (3) Spezzano et al.
(2012).

##  3 Results

The integrated intensities were derived calculating the zeroth moment of the
emission maps, integrating over the velocity ranges given in Table 1. The
ranges contain all velocity channels whose map showed a signal-to-noise ratio
larger than or equal to 3. To enable the comparison between the molecules, all
maps are convolved to an angular resolution of $31^{\prime\prime}$, which
corresponds to the half power beam width (HPBW) of the largest observed beam
(c-H13CC2H2).

##### c-C3H2.

The integrated intensity map of para-c-C3H2 ($3_{22}-3_{13}$) is shown in the
top left panel of Fig. 1. The emission peak is located to the south-east of
the dust peak. This is driven by a non-uniform illumination from the
interstellar radiation field, which dominates towards the southern part of the
core and leads to more C atoms in the gas phase and subsequently an increased
formation of c-C3H2, as discussed in Spezzano et al. (2016).

##### c-H13CC2H2.

The distribution of the c-H13CC2H2 integrated emission (Fig. 1, top right
panel) follows the distribution of c-C3H2, with an emission peak located at
the c-C3H2 peak and similar values of integrated intensity.

##### c-C3HD.

The integrated intensity map of c-C3HD, bottom left panel of Fig. 1, peaks in
an elongated region passing through the dust peak, where temperature and
density are suitable for the deuteration process to take place. The emission
peak follows the direction of the major axis and extends more to the north-
west than to the south-east, suggesting a flattened central structure,
predicted by simulations of magnetized core contraction (see e.g. Caselli et
al. 2019).

##### c-C3D2.

The integrated intensity map of c-C3D2 is shown in the bottom right panel of
Fig. 1. The emission peak is located very close to the dust peak, confirming
that closer to the centre of the core the level of deuteration is higher. In
addition, the emission is less extended than for c-C3H2 and c-C3HD, probably
tracing a higher density and more central shell of the core.

Figure 1: Integrated intensity maps of the observed transitions. The black
dashed lines represent the 90%, 50% and 30% of the H2 column density peak
value derived from Herschel maps (Spezzano et al., 2016), 2.8$\times$1022
cm-2. The solid lines represent contours of the integrated intensity starting
at $3\sigma$ with steps of $3\sigma$. From top to bottom and left to right,
the average errors on the integrated intensity are 9, 9, 12, 8, in units of mK
km s-1. The dust peak is indicated by the black triangle. The beam size of the
30 m telescope, $\rm HPBW=31^{\prime\prime}$, is shown by the red circle in
the bottom left of each panel, and the scalebar is shown in the bottom right
corners.

##  4 Analysis

###  4.1 Column densities

In this work, all column densities are calculated using the assumption of
optically thin emission, as presented in Mangum & Shirley (2015). Furthermore,
we apply the approximation of a constant excitation temperature throughout the
core (CTex), following Caselli et al. (2002b) and Redaelli et al. (2019):

$N=\frac{8\pi\nu^{3}}{c^{3}}\frac{Q_{\mathrm{rot}}(T_{\mathrm{ex}})}{g_{u}A_{\mathrm{ul}}}\left[J_{\nu}(T_{\mathrm{ex}})-J_{\nu}(T_{\mathrm{bg}})\right]^{-1}\frac{\mathrm{e}^{\frac{E_{u}}{kT_{\mathrm{ex}}}}}{\mathrm{e}^{\frac{h\nu}{kT_{\mathrm{ex}}}}-1}\int
T_{\mathrm{mb}}\mathrm{d}v\;,$ (2)

where $Q_{\mathrm{rot}}(T_{\mathrm{ex}})$ is the partition function of the
molecule at an excitation temperature $T_{\mathrm{ex}}$, $g_{u}$ and $E_{u}$
are the degeneracy and energy of the upper level of the transition,
respectively, $A_{\mathrm{ul}}$ the Einstein coefficient for spontaneous
emission, $T_{\mathrm{bg}}=2.73$ K the temperature of the cosmic mircowave
background, $J(T)$ the Rayleigh-Jeans equivalent temperature and
$T_{\mathrm{mb}}$ the main beam temperature. The corresponding parameters for
each transition used in the derivation of the column density are listed in
Table 2.

As excitation temperature, we use 6 K for the main species and for the
13C-bearing isotopologue (see Gerin et al. 1987). For the deuterated
isotopologues, we choose an excitation temperature of 5 K, following Spezzano
et al. (2013). The effect of the excitation temperature on the derived column
density ratios was found to be small, with a change of few percent upon a
variation of $\pm$ 1 K (as stated in Spezzano et al. 2013).
Our assumption of a constant excitation temperature applied across the maps
follows the analysis presented in the appendix of Redaelli et al. (2019).
There, the authors prove that their results are not sensitive to the
variations of the excitation temperature across the core and thus the analysis
can be done assuming a constant value. As our observations have similar sized
maps and angular resolution, we apply this assumption.

The derived column density maps are shown in Fig. 7 in the Appendix. The
morphologies of the maps follow the integrated intensity maps shown in Fig. 1.
The main species shows the highest column densities, peaking at
2.0$\times$1013 cm-2. The 13C isotopologue and the singly deuterated form have
a lower peak column density than the main isotopologue, with peak values
around 1.3$\times$1012 cm-2 and 4.0$\times$1012 cm-2. Furthermore, the doubly
deuterated isotopologue only peaks at $\sim$ 6.3$\times$1011 cm-2, also one
order of magnitude less in column density than its precursor molecule, c-C3HD.
On average, the uncertainties on the column density maps are 1$\times 10^{12}$
cm-2, 9$\times 10^{10}$ cm-2, 2$\times 10^{11}$ cm-2 and 1$\times 10^{11}$
cm-2 for c-C3H2, c-H13CC2H2, c-C3HD and c-C3D2, respectively. The error maps
are derived through standard error propagation from the rms maps, which are
calculated from the velocity channel maps that do not show any signal.

Furthermore, we extracted spectra at the dust peak and the c-C3H2 peak of
L1544, respectively, which are shown in Fig. 8 in Appendix B. Using Gaussian
fitting in Python, we derive the column densities at the two peaks. The
resulting best-fit parameters along with the corresponding column densities
and optical depths are summarised in Table 2. In general, the derived column
densities of the deuterated isotopologues are higher towards the dust peak
than the c-C3H2 peak. The derived optical depths show that all observed lines
are optically thin. However, by multiplying the column density of the 13C
isotopologue with the isotopic ratio for the local interstellar medium,
12C/13C = 68 (Milam et al., 2005), divided by two (to account for the
degeneracy of c-H13CC2H regarding the position of the 13C atom), we obtain an
additional way to determine the column density of the main isotopologue.
Despite the large uncertainty of the resulting value, $N(^{13}$C,CTex$)=(44\pm
10)\times 10^{12}$ cm-2, it might suggest that the c-C3H2 ($3_{33}-3_{13})$
transition is moderately optically thick.

Table 2: Parameters used in the derivation of the column densities. Molecule | Tex (K) | Q(T${}_{\mathrm{ex}})$ | Eu (K) | A $(10^{-5}$s${}^{-1})$ | gu
---|---|---|---|---|---
c-C3H2 ($3_{22}-3_{13}$) | 6 | 37.68 | 16.14 | 1.04 | 7
c-H13CC2H ($2_{12}-1_{01}$) | 6 | 38.69 | 6.33 | 2.17 | 10
c-C3HD ($3_{03}-2_{12}$) | 5 | 50.10 | 10.86 | 3.96 | 21
c-C3D2 ($3_{13}-2_{02}$) | 5 | 86.37 | 9.89 | 3.88 | 42
Table 3: Properties of the spectral lines extracted at the dust peak and the
c-C3H2 peak.

Molecule | $T_{\mathrm{mb}}$ |  $\tau$ aa$a$Derived using the excitation temperatures listed in Table 2. | rms | $W$ | $\rm{v}_{\mathrm{LSR}}$ |  $\Delta\rm{v}$ bb$b$Observed linewidth. |  $N$ aa$a$Derived using the excitation temperatures listed in Table 2.
---|---|---|---|---|---|---|---
| (K) |  | (mK) | (K km s-1) | (km s-1) | (km s-1) | (1012 cm-2)
Dust peak
c-C3H2 | 0.28(4) | 0.10(1) | 35 | 0.10(2) | 7.20(2) | 0.33(5) | 17(3)
c-H13CC2H | 0.19(3) | 0.07(1) | 32 | 0.12(3) | 7.15(4) | 0.6(1) | 1.3(3)
c-C3HD | 0.59(3) | 0.36(2) | 39 | 0.27(2) | 7.20(1) | 0.43(3) | 5.3(4)
c-C3D2 | 0.14(2) | 0.07(1) | 27 | 0.07(2) | 7.19(4) | 0.47(9) | 0.88(3)
c-C3H2 peak
c-C3H2 | 0.49(4) | 0.18(1) | 28 | 0.14(2) | 7.21(1) | 0.27(3) | 23(3)
c-H13CC2H | 0.32(3) | 0.10(1) | 32 | 0.13(2) | 7.21(2) | 0.40(4) | 1.4(2)
c-C3HD | 0.74(4) | 0.47(2) | 39 | 0.24(2) | 7.29(1) | 0.31(2) | 4.7(4)
c-C3D2 | 0.16(3) | 0.08(2) | 27 | 0.04(1) | 7.30(2) | 0.25(6) | 0.5(1)
222

###  4.2 Deuterium fraction

The deuterium fraction maps are derived by dividing the column density maps of
the deuterated isotopologues by that of the main species, pixel by pixel. To
avoid issues with the optical depth of the main isotopologue, we use the
column density map for c-C3H2 derived from the 13C isotopologue multiplied by
the isotopic ratio for the local interstellar medium divided by two,
${}^{12}\mathrm{C}/^{13}\mathrm{C}$ = 68/2 = 34. Colzi et al. (2020) showed
that in the local interstellar medium the ratio is molecule-dependent and
varies with time, volume density, and temperature, and can deviate
significantly from 68. For example, for a fixed density and temperature, it
may change by a factor of 2 between 105 and 106 years of cloud evolution,
depending on the molecule. It is important to note that the emission of the
different isotopologues of c-C3H2 across the core comes mostly from a layer of
the core with volume density similar to the critical density of the
transitions (a few $\times$105 cm-3). Due to this, we only expect minor
variation of the 12C/13C ratio of c-C3H2 across L1544. While variations of the
12C/13C ratio with respect to the ISM value will affect the absolute values of
the column densities and column density ratios we derive, they will not affect
the overall morphology of the maps and the conclusions that we draw here.

The resulting deuterium fraction maps are shown in Fig. 2, where we show
pixels with signal above $3\sigma$. The dust continuum peak is indicated by a
black triangle, while the white cross shows the c-C3H2 peak. The map of
$R_{\rm D/H}$(c-C3H${}_{2})=N($c-C3HD)/$N$(c-C3H2) shows a peak located
$\sim$20′′ to the east of the dust peak, reaching a maximum of
$0.17\,\pm\,0.01$. In the map of $R_{\rm
D2/H2}$(c-C3H${}_{2})=N$(c-C3D2)/N(c-C3H2), the peak is located $\sim$20′′
east of the dust peak, with a maximum of $0.025\,\pm\,0.003$. The level of
deuteration in the map of $R_{\rm D2/HD}$(c-C3H${}_{2})=N$(c-C3D2)/$N$(c-C3HD)
ranges between 0.10 and 0.16, while the peak is located at the dust peak, with
a maximum of $0.16\,\pm\,0.03$. The distribution shows a concentration of
pixels with values of 0.15 - 0.16 spreading between the dust peak and the
c-C3H2 peak.

Figure 2: Deuteration maps of c-C3H2. Only pixels above $3\sigma$ are plotted,
pixels below are grey. The black dashed lines represent 90%, 50% and 30% of
the H2 column density peak value derived from Herschel maps (Spezzano et al.,
2016), 2.8$\times$1022 cm-2. The black triangle indicates the dust continuum
peak (Ward-Thompson et al., 1999), the white cross shows the c-C3H2 peak. The
HPBW is indicated by a red circle in each panel. The dashed line in the left
panel shows the cut used to produce Figure 6.

###  4.3 Non-LTE modelling

We test the assumption of local thermodynamic equilibrium (LTE) used in the
derivation of the column densities by modelling the spectra at the dust peak,
using non-LTE simulations. For this, we run the radiative transfer codes RADEX
(van der Tak et al., 2007) and MOLLIE (Keto, 1990; Keto et al., 2004).

The radiative transfer code RADEX is one-dimensional, using the escape
probability formulation to simplify the calculations. For this work, we apply
the geometry of a static, spherically symmetric and uniform medium for the
escape probability used to derive the column densities.

In contrast to RADEX, MOLLIE solves the radiative transfer equation on a
three-dimensional grid, using a spherically symmetric model of the source.
This physical model is taken from Keto et al. (2015), who describe an unstable
quasi-equilibrium Bonnor-Ebert sphere with a peak central H2 volume density of
$n_{0}\approx 10^{7}$ cm-3 and a central gas temperature of 6 K. The model
provides the infall velocity, density and gas temperature structure for MOLLIE
(see Fig. 3). The simulated molecular abundance profiles are derived using a
state-of-the-art gas-grain chemical model (Sipilä et al., 2019), assuming an
external visual extinction of $A_{\mathrm{V}}=2\,$mag to account for the fact
that L1544 is embedded in a molecular cloud. We employ the same initial
abundances and chemical networks as in Sipilä et al. (2019). The initial
abundances are reproduced here in Table 4. We assume ‘standard’ values for
various model parameters; we set a constant cosmic-ray ionization rate
$\zeta=1.3\times 10^{-17}\,\rm s^{-1}$, a grain radius of 0.1 $\mu$m
(monodisperse grains), and a grain material density of $2.5\,\rm g\,cm^{-3}$.
Redaelli et al. (2021) found a cosmic-ray ionization rate of $\zeta=\sim
3\times 10^{-17}\,\rm s^{-1}$, which is a factor of $\sim 2$ different to the
value we apply. In Appendix C, we discuss the effect of an increased cosmic-
ray ionization rate on the chemistry. The inital ortho-to-para ratio of H2 is
set to $10^{-3}$ and evolves with the rest of the chemistry. We adopt a two-
phase description of gas-grain chemistry, that is, the ice on the grain
surfaces is treated as a single reactive layer. Molecules may be desorbed
thermally, or via cosmic ray induced desorption, photodesorption, or chemical
desorption (with 1% efficiency).

With the two different codes we intend to test how a rather ‘simple’ and often
used code like RADEX, which only applies a general spherical core geometry
with constant molecular abundance, volume density and temperature, performs
compared to a more sophisticated but time-consuming simulation like MOLLIE,
which considers the specific physical structure of a source and molecular
abundance profiles across it.

Table 4: Initial abundances (with respect to $n_{\rm H}\approx 2\,n({\rm
H_{2}})$) used in the chemical modelling.

Species | Abundance
---|---
$\rm H_{2}$ | $5.00\times 10^{-1}\,{}^{(a)}$
$\rm He$ | $9.00\times 10^{-2}$
$\rm C^{+}$ | $1.20\times 10^{-4}$
$\rm N$ | $7.60\times 10^{-5}$
$\rm O$ | $2.56\times 10^{-4}$
$\rm S^{+}$ | $8.00\times 10^{-8}$
$\rm Si^{+}$ | $8.00\times 10^{-9}$
$\rm Na^{+}$ | $2.00\times 10^{-9}$
$\rm Mg^{+}$ | $7.00\times 10^{-9}$
$\rm Fe^{+}$ | $3.00\times 10^{-9}$
$\rm P^{+}$ | $2.00\times 10^{-10}$
$\rm Cl^{+}$ | $1.00\times 10^{-9}$
$\rm HD$ | $1.60\times 10^{-5}$
333(a) The initial $\rm H_{2}$ ortho/para ratio is $1\times 10^{-3}$.

Figure 3: Profiles of the gas temperature (red), H2 number density (black, in
logarithmic scale) and infall velocity (blue, in units of 0.1 km s-1) for the
L1544 model developed by Keto et al. (2015). The velocity in the model is
negative, but here it is shown as positive for better readability.

####  4.3.1 RADEX: Results

The results are obtained by manually fitting the synthetic spectra derived by
RADEX to the observed spectra at the dust peak given in Figure 8 in the
appendix. This was done using the CASSIS software (Vastel et al., 2015),
assuming $n(\mathrm{H}_{2})=10^{6}$ cm-3. We set $T_{\mathrm{kin}}=5$ K for
c-C3D2 ($3_{13}-2_{02}$) and $T_{\mathrm{kin}}=8$ K for c-C3H2
($3_{22}-3_{13}$) and c-C3HD ($3_{03}-2_{12}$), corresponding to the lowest
temperatures with available collision rates. The collision rates of c-C3H2,
c-C3HD and c-C3D2 with H2 were provided by Ben Khalifa et al. (2019) for the
main species and Ben Khalifa et al. (in prep.) for the isotopologues. Because
of its different symmetry, the c-H13CC2H isotopologue collisional rates cannot
be deduced from the other rates, and they have not been computed. Hence, we do
not model the 13C isotopologue in the non-LTE formalisms.

The total column densities corresponding to the best fit solutions of RADEX
are listed in column two of Table 4. The corresponding uncertainties of RADEX
are estimated by varying the column density to result in a change of the peak
temperature by 25%. To derive the total column densities of c-C3H2 and c-C3D2,
the RADEX output column densities are multiplied by the ortho-to-para ratio,
which is 3:1 and 1:2 for the two molecules, respectively. This is done because
the collision rates used in the calculations treat ortho and para as two
different species, so they are modelled separately.

Using RADEX, we can, within the errorbars, reproduce the column density of
c-C3H2 observed at the dust peak. The value for c-C3HD is only slightly
underestimated; c-C3D2, however, is overestimated by a factor of $\sim 2$.

####  4.3.2 MOLLIE: Results

As output, MOLLIE produces a synthetic spectral line profile. The molecular
column density is obtained by integrating the abundance profile multiplied by
the gas density, and convolved to a desired beam size. The collisional
coefficients needed for the line radiative transfer are provided by Ben
Khalifa et al. (2019) and Ben Khalifa et al. (in prep).

The total column densities corresponding to the best fit solutions of MOLLIE
are listed in columns three and four of Table 4. The values for c-C3H2 and
c-C3D2 are derived by multiplying the MOLLIE output column densities by the
ortho-to-para ratio (3:1 and 1:2 for c-C3H2 and c-C3D2, respectively) to
account for the separate treatment of ortho and para species during the
modelling. The errors of the column densities derived with MOLLIE lie roughly
within a factor of 2, and are constrained by the resolution of the time steps
of the chemical model.

Figure 4 shows the total molecular abundance profiles of c-C3H2, c-C3HD and
c-C3D2 with respect to H2 at the two best fitting time steps of the chemical
model, $t_{1}$ = 1.0$\times$105 yrs (solid lines) and $t_{2}$ = 1.4$\times$105
yrs (dashed lines). At the first time step, c-C3H2 peaks at $\sim$
5$\times$103 AU, while c-C3HD and c-C3D2 peak at $\sim$ 3$\times$103 AU,
slightly shifted towards the centre of the core. At the second time step, all
peaks are shifted to larger radii, as the molecules get depleted in the centre
of the core with the chemistry evolving. Furthermore, c-C3H2 shows a rather
flat profile in the outer radii, being abundant also in less dense parts,
while the D-bearing species quickly drop off after they peak, staying more
concentrated within the central parts.

We only vary the time step in the chemical model in our analysis but not any
other input parameters of the model. Therefore, we tested the effect of the
various parameters by running four single-point simulations. The results of
these test simulations are presented in Appendix C.

Figure 5 compares the spectra modelled with MOLLIE (coloured) at the two time
steps $t_{1}$ = 1.0$\times$105 yrs (solid lines) and $t_{2}$ = 1.4$\times$105
yrs (dashed lines) with the observed spectra (black). The later time step
shows a very good fit with the observed spectrum of c-C3H2. The deuterated
isotopologues, on the other hand, are both underestimated by a factor of 3.
They are better reproduced by $t_{1}$, where the main species is overestimated
by a factor of 7. The same applies for the derived column densities; at
$t_{2}$, the observed values of c-C3HD and c-C3D2 are reproduced, while,
within the error bars, c-C3H2 is reproduced at $t_{1}$. One would expect the
main species to match with the earlier time step, as the deuteration process
is expected to happen late in the core evolution, with the deuterated
isotopologues peaking at a later time. However, it is important to note that
the chemical code used to produce the abundance profiles is static (i.e., the
physical structure does not evolve with time). Thus, the timesteps of the
chemical model are not linked to the dynamical evolution of the cloud. This
behaviour can be explained by the fact that the chemical time scale is
approximately inversely proportional to the volume density (more reactants
imply faster rates), so the central regions have faster chemical time scales
compared to the outer less dense regions. When using a static core, this could
lead to inconsistencies, as, in reality, a contracting core spends longer
times at low densities, while the denser regions form only at a later stage of
its evolution. Therefore, the small difference between the time steps and the
large difference between the resulting spectra demonstrates the sensitivity of
the chemical model to the chosen time step and points out its limitations.

Figure 4: Total molecular fractional abundance profiles with respect to H2
used to obtain the MOLLIE best fit solutions. Plotted are c-C3H2 (red, c-C3HD
(pink) and c-C3D2 (blue), at two different time steps of the chemical model,
$t_{1}$ = 1.0$\times$105 yrs (solid lines) and $t_{2}$ = 1.4$\times$105 yrs
(dashed lines). Figure 5: Comparison of the best fit solutions of spectra
(coloured) modelled with MOLLIE at two different time steps of the chemical
model, $t_{1}$ = 1.0$\times$105 yrs (solid) and $t_{2}$ = 1.4$\times$105 yrs
(dashed), with the observed spectra (black). Plotted are the synthetic spectra
of c-C3H2 (left), c-C3HD (middle), and c-C3D2 (right).  Table 5: Comparison of
total column densities, in units of $\times 10^{12}$ cm-2

Molecule |  $N$(RADEX) |  $N_{t_{1}}$(MOLLIE) |  $N_{t_{2}}$(MOLLIE)
---|---|---|---
c-C3H2 | 14(4) | 160 | 30
c-C3HD | 3.8(9) | 3.8 | 0.8
c-C3D2 | 1.8(3) | 0.79 | 0.23
444The column densities derived with MOLLIE (NM) are given at the chemical
times $t_{1}$ = 1.0$\times$105 yrs and $t_{2}$ = 1.4$\times$105 yrs. The best
fit for the deuterated species is $t_{1}$, for the main species $t_{2}$. The
error on the derived values lies within a factor of 2.

##  5 Discussion

The deuterium fraction maps presented in this work suggest that the
deuteration of c-C3H2 is more efficient towards the centre of L1544. This
implies that there is still c-C3H2 available to be deuterated also at high
densities, despite the freeze-out. However, carbon chains are not expected to
trace the dense central regions of pre-stellar cores, as they are expected to
be frozen out. The detection of c-C3D2 peaking close to the dust peak
demonstrates that a significant fraction of this species is still present
towards the core centre. This might be due to the large amount of CO freeze-
out, as our chemical models show that with increasing CO freeze-out, the
formation of He+ is enhanced. With increasing level of He+ in the gas phase,
also the amount of carbon atoms needed to produce carbon chains increases, as
He+ efficiently destroys the CO molecules left in the gas phase via the
reaction $\mathrm{He}^{+}+\mathrm{CO}=\mathrm{C}^{+}+\mathrm{O}+\mathrm{He}$.
This was already predicted by Ruffle et al. (1999), whose models show an
increase of C+ at $3\times 10^{6}$ yrs, when CO starts its dramatic freeze-
out, corresponding to a plateau in the HC3N abundance (Figure 1 and 2 in
Ruffle et al. 1999)

The deuteration process of c-C3H2 consists of two reactions, the proton-
deuteron transfer with H2D+ and the subsequent dissociative recombination with
electrons (see Spezzano et al. 2013). Therefore, the D/H ratio of c-C3H2 is
directly related to that of H${}_{3}^{+}$. The distribution of H2D+ across
L1544 was shown to be extended and an excellent tracer of the dust continuum
(Vastel et al., 2006; Koumpia et al., 2020), with the emission peak at the
dust peak. The correlation between the degree of deuteration in molecules and
the abundance of H2D+, found by Vastel et al. (2006), is now further confirmed
by our map of c-C3D2.

The deuterium fraction maps of c-C3H2 reveal how the efficiency of the
deuteration varies across the core. However, it has to be taken into account
that the deuterated isotopologues are much less abundant at larger radii than
the main species, and more concentrated at higher densities (see abundance
profiles in Fig. 4). Therefore, the column density of the normal species also
contains contributions from the outer layers of the core, where c-C3H2 is
still abundant, while c-C3HD and c-C3D2 are not present. Due to this, the
observed emission of the normal isotopologue does not necessarily come from
the same layers where the deuterated isotopologues are located and the
corresponding deuterium fraction maps have to be interpreted with caution. By
analysing $R_{\rm D2/HD}$(c-C3H2), this problem is avoided, as the two
molecules should trace more similar regions of the core (as also shown by the
abundace profiles in Figure 4). The $R_{\rm D2/HD}$(c-C3H2) map shows a peak
slightly shifted to the south-east of the dust peak, in between the dust and
c-C3H2 peaks. This confirms that the deuteration process, and especially
higher order deuteration, is more efficient towards the centre of the core.
The peak values of the deuteration maps derived in this work are
$R_{\mathrm{HD/H2}}($c-C3H${}_{2})=0.17$, $R_{\rm D2/H2}$(c-C3H${}_{2})=0.025$
and $R_{\rm D2/HD}$(c-C3H${}_{2})=0.16$, in good agreement with previous
measurements towards the dust peak of L1544 (see e.g. Spezzano et al. 2013,
Chantzos et al. 2018).

Using the non-LTE models RADEX and MOLLIE, we are able to reproduce the
observed spectra towards the dust peak of L1544 in both cases. The derived
column densities agree with the LTE calculations within a factor of 2. This
shows that for our observed transitions, we can safely use RADEX for non-LTE
calculations instead of a costly and time-consuming simulation like MOLLIE. We
further used the results of RADEX and MOLLIE to estimate the optical depth of
the possibly optically thick transition of c-C3H2, ($3_{22}-3_{13}$). We found
values ranging from 0.07 up to 0.3, indicating optically thin emission.
Nevertheless, the non-LTE simulations confirm that the assumption of CTex is a
useful approximation to derive the column density of optically thin lines.

Recently, Redaelli et al. (2019) and Chacón-Tanarro et al. (2019) studied the
deuteration throughout L1544, using emission maps of the molecules N2H+, HCO+
and CH3OH, H2CO, respectively. The molecules discussed in Redaelli et al.
(2019), N2H+ and HCO+, show a similar behaviour to c-C3H2: the D-bearing
species peak closer to the dust peak than their corresponding normal species,
and the deuteration maps show a compact morphology around the dust peak. N2H+
is more deuterated than c-C3H2, as N2, and therefore N2H+, is a late-type
molecule and suffers less of depletion in the dense central parts of the core.
Thus, N2H+ reaches higher levels of deuteration (26%) than c-C3H2 (17%),
tracing the densest regions. HCO+, on the other hand, is formed from CO, which
is highly depleted in the centre, so it traces outer layers of the core,
similar to c-C3H2, only showing low levels of deuteration ($\approx 3.5$%).

Using methanol (CH3OH) and formaldehyde (H2CO), which form on dust grain
surfaces, Chacón-Tanarro et al. (2019) focus on the deuteration process in the
solid phase. A comparison of the deuteration maps of methanol to the $R_{\rm
HD/H2}$(c-C3H2) map and the $R_{\rm D2/H2}$(c-C3H2) map shows an interesting
similarity. The deuteration peaks of both molecules are shifted from the dust
peak towards the opposite direction of their respective molecular peaks, with
respect to the major axis of L1544. Following this, the morphology of the maps
might not be the consequence of an increased deuteration, but rather a steep
decrease of the column density of the main species in the outer layers of the
core. The emission of the deuterated isotopologues comes instead only from the
inner layers, and is not affected by the inhomogeneous interstellar radiation
field around the core (that drives the segregation among methanol and c-C3H2).
This is confirmed by the fact that the $R_{\rm D2/HD}$(c-C3H2) ratio peaks
just next to the dust peak. The distribution of the deuterium fraction of H2CO
shows an opposite behaviour with respect to c-C3H2. The $R_{\rm
HD/H2}$(c-C3H2) map is extended through the core, while for $R_{\rm
HD/H2}$($\rm H_{2}CO$) the map is more centred around its peak located towards
the north-west of the dust peak. For the double-deuteration maps the same
applies, however, this time the distribution of $\rm H_{2}CO$ is more extended
and c-C3H2 is more centered. The more extended $R_{\rm D2/H2}$($\rm H_{2}CO$)
map might be linked to the fact that the $\rm H_{2}CO$ isotopologues have both
grain-surface and gas-phase formation mechanisms, while c-C3H2 isotopologues
solely form in the gas-phase. The $R_{\rm D2/HD}$($\rm H_{2}CO$) map is not
discussed in Chacón-Tanarro et al. (2019) due to high uncertainties and hence
can not be used for a comparison. A major difference between c-C3H2 and H2CO
is their deuteration level. For H2CO, both $R_{\rm HD/H2}$($\rm H_{2}CO$) and
$R_{\rm D2/H2}$($\rm H_{2}CO$) lie below 0.05. Contrary to this, for c-C3H2,
the singly deuterated species is much more abundant than the doubly deuterated
one, reaching $R_{\rm HD/H2}$(c-C3H${}_{2})=0.17$ and $R_{\rm
D2/H2}$(c-C3H${}_{2})=0.025$, respectively. This indicates that for $\rm
H_{2}CO$, higher order deuteration happens much more efficiently.

By cutting across the dust peak and the $R_{\rm HD/H2}(\mathrm{HCO}^{+})$
maximum, Redaelli et al. (2019) compared the trends of the deuteration maps of
their molecular tracers. Figure 6 applies the same cut to the deuteration maps
of c-C3H2 derived in this work, plotting the deuteration fractions as a
function of the impact parameter (=projected distance to the dust peak) and
comparing them to HCO+. The cut is visualized as a dashed line in the left
panel of Fig. 2. Along the cut, the deuteration levels of $R_{\rm
HD/H2}$(c-C3H${}_{2})$ and $R_{\rm D2/HD}$(c-C3H2) decrease faster than
$R_{\rm D2/H2}$(c-C3H${}_{2})$ and $R_{\rm HD/H2}(\mathrm{HCO}^{+})$. The
slope of $R_{\rm HD/H2}$(c-C3H${}_{2})$ shows a similar behaviour as HCO+, as
both peak between 10$"$-20$"$ from the dust peak and then drop off, though the
one of c-C3H2 is steeper. The fast decrease of $R_{\rm D2/HD}$(c-C3H2) with
increasing distance from the dust peak confirms once more that the deuteration
process is most efficient in the central region. This enhanced c-C3H2
deuteration towards the centre of a pre-stellar core could explain the large
$R_{\rm HD/H2}$(c-C3H${}_{2})$ values found by Chantzos et al. (2018) towards
the protostar HH211, one of the youngest class 0 sources.

Figure 6: Comparison of the trends of the HCO+ deuteration map presented in
Redaelli et al. (2019) and the deuteration maps of c-C3H2 derived in this
work, along the cut shown in Figure 10 in Redaelli et al. (2019). The data
points of $R_{\rm HD/H2}(\mathrm{HCO}^{+})$ and $R_{\rm D2/H2}$(c-C3H${}_{2})$
are shifted upwards by 0.10 to allow an easier comparison.

##  6 Conclusions

In this work, we studied the deuteration process of cyclopropenylidene towards
the pre-stellar core L1544. We presented the first emission map of doubly
deuterated c-C3H2, along with the emission maps of its singly deuterated and
non-deuterated isotopologues. Peaking at the centre of the core, the column
density map of c-C3D2 proves that carbon-chain molecules are still present at
the high central densities. This is possible due to an increased abundance of
He+ that destroys CO molecules and therefore increases the amount of free
carbon atoms available to form c-C3H2 (Ruffle et al., 1999). Furthermore, the
distributions of c-C3HD and c-C3D2 indicate that the deuterated forms of
c-C3H2 are preferentially located towards the centre of L1544, where the
increased abundance of H2D+ drives the deuteration process, and do not trace
the c-C3H2 peak located in the south-east, which is caused by external
illumination of the interstellar radiation field. Using the column density
maps calculated assuming CTex, we derived the first deuteration maps of
c-C3H2. Showing a maximum at the dust peak, the maps confirm that the
deuteration process of c-C3H2 is more efficient towards the centre of the
core. The maximum values in the maps are $R_{\rm HD/H2}$(c-C3H${}_{2})=0.17\pm
0.01$, $R_{\rm D2/H2}$(c-C3H${}_{2})=0.025\pm 0.003$ and $R_{\rm
D2/HD}$(c-C3H${}_{2})=0.16\pm 0.03$, which are in good agreement with previous
single point measurements (Spezzano et al., 2013; Chantzos et al., 2018).

Using the non-LTE radiative transfer models RADEX and MOLLIE, we were able to
reproduce the spectra of c-C3H2, c-C3HD and c-C3D2 at the dust peak of L1544.
The derived non-LTE column densities agree with the CTex calculations within a
factor of 2. The comparison of our methods shows that for our observed
transitions, we can safely use RADEX for non-LTE calculations instead of a
time-consuming simulation like MOLLIE.

Finally, we compared the deuterium fraction maps derived in this work to the
deuteration of other molecular tracers. The discussion showed how studying the
deuteration of different species can bring important information on the
chemistry of different parts of the source and helps to build a comprehensive
picture of L1544. However, so far only a limited number of deuteration maps is
available in the literature. Future projects will be able to address this
deficit and lead to a deeper understanding of the deuteration process of
individual molecules. In particular, future studies can further investigate
the detection of carbon-chain molecules at unexpectedly high densities to put
tighter constraints on gas-grain chemical models. Suitable candidates for
observations of other abundant carbon-chains are for example CH3CCH and HC3N,
which show a not yet understood segregation towards the pre-stellar core L1544
(see Spezzano et al. 2017).

###### Acknowledgements.

S.S. and K.G. wish to thank the Max Planck Society for the Max Planck Research
Group funding. All others authors affiliated to the MPE wish to thank the Max
Planck Society for financial support.

## References

  * Bacmann et al. (2003) Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2003, ApJ, 585, L55
  * Ben Khalifa et al. (2019) Ben Khalifa, M., Sahnoun, E., Wiesenfeld, L., et al. 2019, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 21, 1443
  * Bogey et al. (1987) Bogey, M., Demuynck, C., Destombes, J. L., & Dubus, H. 1987, Journal of Molecular Spectroscopy, 122, 313
  * Caselli et al. (2012) Caselli, P., Keto, E., Bergin, E. A., et al. 2012, ApJ, 759, L37
  * Caselli et al. (2019) Caselli, P., Pineda, J. E., Zhao, B., et al. 2019, ApJ, 874, 89
  * Caselli et al. (1999) Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C. 1999, ApJ, 523, L165
  * Caselli et al. (2002a) Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002a, ApJ, 565, 331
  * Caselli et al. (2002b) Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002b, ApJ, 565, 344
  * Ceccarelli et al. (2014) Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 859
  * Chacón-Tanarro et al. (2019) Chacón-Tanarro, A., Caselli, P., Bizzocchi, L., et al. 2019, A&A, 622, A141
  * Chantzos et al. (2018) Chantzos, J., Spezzano, S., Caselli, P., et al. 2018, ApJ, 863, 126
  * Colzi et al. (2020) Colzi, L., Sipilä, O., Roueff, E., Caselli, P., & Fontani, F. 2020, A&A, 640, A51
  * Cox et al. (1987) Cox, P., Guesten, R., & Henkel, C. 1987, A&A, 181, L19
  * Crapsi et al. (2005) Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2005, ApJ, 619, 379
  * Crapsi et al. (2007) Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, A&A, 470, 221
  * Dalgarno & Lepp (1984) Dalgarno, A. & Lepp, S. 1984, ApJ, 287, L47
  * Galli et al. (2019) Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A, 630, A137
  * Gerin et al. (1987) Gerin, M., Wootten, H. A., Combes, F., et al. 1987, A&A, 173, L1
  * Hily-Blant et al. (2018) Hily-Blant, P., Faure, A., Rist, C., Pineau des Forêts, G., & Flower, D. R. 2018, MNRAS, 477, 4454
  * Keto & Caselli (2010) Keto, E. & Caselli, P. 2010, MNRAS, 402, 1625
  * Keto et al. (2015) Keto, E., Caselli, P., & Rawlings, J. 2015, MNRAS, 446, 3731
  * Keto et al. (2004) Keto, E., Rybicki, G. B., Bergin, E. A., & Plume, R. 2004, ApJ, 613, 355
  * Keto (1990) Keto, E. R. 1990, ApJ, 355, 190
  * Kong et al. (2015) Kong, S., Caselli, P., Tan, J. C., Wakelam, V., & Sipilä, O. 2015, ApJ, 804, 98
  * Koumpia et al. (2020) Koumpia, E., Evans, L., Di Francesco, J., van der Tak, F. F. S., & Oudmaijer, R. D. 2020, A&A, 643, A61
  * Lee et al. (2001) Lee, C. W., Myers, P. C., & Tafalla, M. 2001, ApJS, 136, 703
  * Lucas & Liszt (2000) Lucas, R. & Liszt, H. S. 2000, A&A, 358, 1069
  * Madden et al. (1989) Madden, S. C., Irvine, W. M., Matthews, H. E., Friberg, P., & Swade, D. A. 1989, AJ, 97, 1403
  * Majumdar et al. (2017) Majumdar, L., Gratier, P., Andron, I., Wakelam, V., & Caux, E. 2017, MNRAS, 467, 3525
  * Mangum & Shirley (2015) Mangum, J. G. & Shirley, Y. L. 2015, PASP, 127, 266
  * Matthews & Irvine (1985) Matthews, H. E. & Irvine, W. M. 1985, ApJ, 298, L61
  * Milam et al. (2005) Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126
  * Müller et al. (2001) Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49
  * Ohashi et al. (1999) Ohashi, N., Lee, S. W., Wilner, D. J., & Hayashi, M. 1999, ApJ, 518, L41
  * Pagani et al. (1992) Pagani, L., Salez, M., & Wannier, P. G. 1992, A&A, 258, 479
  * Park et al. (2006) Park, I. H., Wakelam, V., & Herbst, E. 2006, A&A, 449, 631
  * Pety (2005) Pety, J. 2005, in SF2A-2005: Semaine de l’Astrophysique Francaise, ed. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani, 721
  * Qi et al. (2013) Qi, C., Öberg, K. I., Wilner, D. J., & Rosenfeld, K. A. 2013, ApJ, 765, L14
  * Redaelli et al. (2019) Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15
  * Redaelli et al. (2021) Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109
  * Ruffle et al. (1999) Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691
  * Seaquist & Bell (1986) Seaquist, E. R. & Bell, M. B. 1986, ApJ, 303, L67
  * Sipilä et al. (2019) Sipilä, O., Caselli, P., Redaelli, E., Juvela, M., & Bizzocchi, L. 2019, MNRAS, 487, 1269
  * Spezzano et al. (2016) Spezzano, S., Bizzocchi, L., Caselli, P., Harju, J., & Brünken, S. 2016, A&A, 592, L11
  * Spezzano et al. (2013) Spezzano, S., Brünken, S., Schilke, P., et al. 2013, ApJ, 769, L19
  * Spezzano et al. (2017) Spezzano, S., Caselli, P., Bizzocchi, L., Giuliano, B. M., & Lattanzi, V. 2017, A&A, 606, A82
  * Spezzano et al. (2012) Spezzano, S., Tamassia, F., Thorwirth, S., et al. 2012, ApJS, 200, 1
  * Teyssier et al. (2005) Teyssier, D., Hily-Blant, P., Gerin, M., et al. 2005, in ESA Special Publication, Vol. 577, ESA Special Publication, ed. A. Wilson, 423–424
  * Thaddeus et al. (1981) Thaddeus, P., Guelin, M., & Linke, R. A. 1981, ApJ, 246, L41
  * Thaddeus et al. (1985) Thaddeus, P., Vrtilek, J. M., & Gottlieb, C. A. 1985, ApJ, 299, L63
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
  * Vastel et al. (2015) Vastel, C., Bottinelli, S., Caux, E., Glorian, J. M., & Boiziot, M. 2015, in SF2A-2015: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 313–316
  * Vastel et al. (2006) Vastel, C., Caselli, P., Ceccarelli, C., et al. 2006, ApJ, 645, 1198
  * Vrtilek et al. (1987) Vrtilek, J. M., Gottlieb, C. A., & Thaddeus, P. 1987, ApJ, 314, 716
  * Ward-Thompson et al. (1999) Ward-Thompson, D., Motte, F., & Andre, P. 1999, MNRAS, 305, 143
  * Willacy et al. (1998) Willacy, K., Langer, W. D., & Velusamy, T. 1998, ApJ, 507, L171
  * Williams et al. (1999) Williams, J. P., Myers, P. C., Wilner, D. J., & Di Francesco, J. 1999, ApJ, 513, L61

##  Appendix A Column density maps

In Figure 7 we show the derived column density maps of the observed molecules,
using the parameters listed in Table 2.

Figure 7: Column density maps of the molecules, used to calculate the
deuteration fraction maps. The black dashed lines represent the 90%, 50% and
30% of the H2 column density peak value derived from Herschel maps (Spezzano
et al. 2016), 2.8$\times$1022 cm-2. The solid lines represent contours of the
column density, starting at $3\sigma$ with steps of $3\sigma$. From top to
bottom and left to right, the average errors on the column density are
1$\times 10^{12}$, 9$\times 10^{10}$, 2$\times 10^{11}$ and 1$\times 10^{11}$,
in units of cm-2. The dust peak is indicated by the black triangle. The beam
size of the 30 m telescope, $\rm HPBW=31^{\prime\prime}$, is shown by the red
circle, and the scalebar is shown in the bottom right corners.

##  Appendix B Averaged spectra

Figure 8 shows the spectra extracted at the dust peak (Ward-Thompson et al.
1999) and the c-C3H2 peak (Spezzano et al. 2016) of L1544. Overlaid are fitted
Gaussian profiles (red).

Figure 8: Comparison of extracted spectra (black) of the transitions listed in
Table 1, at the dust peak (left) and the c-C3H2 peak (right) of L1544.
Overlaid are the fitted Gaussian profiles (red). The vertical lines indicate
the systemic velocity of the source, $7.2$ km s-1.

##  Appendix C Variation of input parameters of the chemical model

To test the effect of the various parameters, we run four single-point
simulations, assuming physical conditions roughly representative of the outer
core region where the c-C3H2 abundance peaks: $n(H_{2})=10^{4}$ cm-3, $T=10$
K, $A_{\mathrm{V}}=5$ mag. In addition to the fiducial parameter set, we test
three variations with respect to the fiducial model: 1) decreased C/O ratio
(from $\sim$.47 to $\sim$.35); 2) higher cosmic-ray ionization rate
($\zeta=10^{-16}$ s-1); 3) higher initial H2 ortho-to-para ratio (3). Figure 9
shows the resulting D/H ratios as a function of time in these four cases. The
solid lines represent the fiducial case, while the dashed, dotted, and dash-
dotted lines represent cases 1 to 3, respectively, as defined above.

A variation of the applied C/O ratio from $\sim$0.47 down to $\sim$0.35 only
slightly changes the results. An increase of the initial ortho-to-para ratio
of H2 up to o:p $=$ 3 simply delays the evolution of the D/H ratios,
eventually arriving at the same end result as the in the case of a low initial
ratio. This behaviour was also observed by Kong et al. (2015). Only the
cosmic-ray ionization rate seems to have a larger effect, where an increase up
to 10-16 s-1 reduces both D/H ratios. Previous work done by Redaelli et al.
(2021) found the ionization rate to be $\sim 3\times 10^{-17}$ s-1, which is
less than a factor of 3 different from our fiducial value. The test
simulations show that the time is not the only variable that influences the
D/H ratios of c-C3H2. While we only explored a small part of the paramter
space, and the tests were limited to a single point and cannot be used to
conclude anything on the column densities, we do not expect variations of more
than a factor of a few in the c-C3H2 D/H ratios, as long as the model
parameters remain reasonably close to those in the fiducial model. A complete
parameter-space explorations is out of the scope of the present work.

Figure 9: Comparison of the evolution of the D/H ratios of c-C3H2 with varied
input parameters in the chemical model: fiducial model (solid lines),
decreased C/O ratio (dashed lines), higher CR ionization rate (dotted lines),
higher initial H2 ortho-to-para ratio (dashed-dotted lines).