11institutetext: Gifu University, Japan 22institutetext: Royal Holloway,
University of London, UK 33institutetext: Brunel University London, UK

#  Kmclib: Automated Inference and Verification of Session Types from OCaml
Programs

Keigo Imai 11 Julien Lange 22 Rumyana Neykova 33

###### Abstract

Theories and tools based on multiparty session types offer correctness
guarantees for concurrent programs that communicate using message-passing.
These guarantees usually come at the cost of an intrinsically top-down
approach, which requires the communication behaviour of the entire program to
be specified as a global type.

This paper introduces kmclib: an OCaml library that supports the development
of _correct_ message-passing programs without having to write any types. The
library utilises the meta-programming facilities of OCaml to automatically
infer the session types of concurrent programs and verify their compatibility
($k$-MC [13]). Well-typed programs, written with kmclib, do not lead to
communication errors and cannot get stuck.

###### Keywords:

Multiparty Session Types Concurrent Programming OCaml

##  1 Introduction

Multiparty session types (MPST) [6] are a popular type-driven technique to
ensure the correctness of concurrent programs that communicate using message-
passing. The key benefit of MPST is to guarantee statically that the
components of a program have compatible behaviours, and thus no components can
get permanently stuck. Many implementations of MPST in different programming
languages have been proposed in the last decade [16, 7, 10, 19, 15, 2, 8, 14,
5, 22], however, all suffer from a notable shortcoming: they require
programmers to adopt a top-down approach that does not fit well in modern
development practices. When changes are frequent and continual (e.g.,
continuous delivery), re-designing the program and its specification at every
change is not feasible.

Most MPST theories and tools advocate an intrinsically top-down approach. They
require programmers to specify the communication (often in the form of a
global type) of their programs before they can be type-checked. In practice,
type-checking programs against session types is very difficult. To circumvent
the problem, most implementations of MPST rely on _external_ toolings that
generate code from a global type, see e.g., all works based on the Scribble
toolchain [21].

In this paper, we present an OCaml library, called kmclib, that supports the
development of programs which enjoy all the benefits of MPST while avoiding
their main drawbacks. The kmclib library guarantees that threads in well-typed
programs will not get stuck. The library also enables _bottom-up development_
: programmers write message-passing programs in a natural way, without having
to write session types. Our library is built on top of Multicore OCaml [20]
that offers highly scalable and efficient concurrent programming, but does not
provide any static guarantees wrt. concurrency.

Figure 1 gives an overview of kmclib. Its implementation combines the power of
the type-aware macro system of OCaml (Typed PPX) with two recent advances in
the session types area: an encoding of MPST in OCaml (channel vector types
[8]) and a session type compatibility checker ($k$-MC checker [13]). To our
knowledge, this is the first implementation of type inference for MPST and the
first integration of compatibility checking in a programming language.

OCaml codeTyped PPXOCaml CompilerSession typed executable$k$-MC checker
[13]Session typesOCaml channel vector types [8]InferTranslateok: boundko:
counterexamples Figure 1: Workflow of the kmclib library.

The kmclib library offers several advantages compared to earlier MPST
implementations. (1) It is _flexible_ : programmers can implement
communication patterns (e.g., fire-and-forget patterns [13]) that are not
expressible in the synchrony-oriented syntax of global types. (2) It is
_lightweight_ as it piggybacks on OCaml’s type system to check and infer
session types, hence lifting the burden of writing session types off the
programmers. (3) It is _user-friendly_ thanks to its integration in Visual
Studio Code, e.g., compatibility violations are mapped to precise locations in
the code. (4) It is _well-integrated_ into the natural edit-compile-run cycle.
Although compatibility is checked by an external tool, this step is embedded
as a compilation step and thus hidden from the user.

##  2 Safe Concurrent Programming in Multicore OCaml

We give an overview of the features and usage of kmclib using the program in
Figure 2 (top) which calculates Fibonacci numbers. The program consists of
three concurrent _threads_ (user, master, and worker) that interact using
point-to-point message-passing. Initially, the user thread sends a request to
the master to start the calculation, then waits for the master to return a
work-in-progress message, or the final result. After receiving the result, the
user sends back a stop message. Upon receiving a new request, the master
splits the initial computation in two and sends two tasks to a worker. For
each task that the worker receives, it replies with a result. The master and
worker threads are recursive and terminate only upon receiving a stop message.

Figure 2 (bottom) gives a session type for each thread, i.e., the behaviour of
each thread wrt. communication. For clarity we represent session types as a
communicating finite state machine (CFSM [1]), where ! (resp. ?) denotes
sending (resp. receiving). For example, $\texttt{um}!\textit{compute}$ means
that the user is sending to the master a message compute, while
$\texttt{um}?\textit{compute}$ says that the master receives compute from the
user. Our library infers these CFSM representations from the OCaml code, in
Figure 2 (top), and verifies _statically_ that the three threads are
_compatible_ , hence no thread can get stuck due to communication errors. If
compatibility cannot be guaranteed, the compiler reports the kind of
violations (i.e., _progress_ or _eventual reception_ error) and their
locations in the code. Figure 3 shows how such semantic errors are reported
visually in Visual Studio Code.

Albeit simple, the common communication pattern used in Figure 2 cannot be
expressed as a global type, and thus cannot be implemented in previous MPST
implementations. Concretely, global types cannot express the intrinsic
asynchronous interactions between the master and worker threads (i.e., the
master may send a second task message, while the worker sends a result).

⬇

1let KMC (uch,mch,wch) = [%kmc.gen (u,m,w)]

2

3let user () =

4 let uch = send uch#m#compute 42 in

5 let rec loop uch : unit =

6 match receive uch#m with

7 | ‘wip(res, uch) ->

8 printf "in progress: %d\n" res;

9 loop uch

10 | ‘result(res, uch) ->

11 printf "result: %d\n" res;

12 send uch#m#stop ()

13 in loop uch

14

15let worker () =

16 let rec loop wch : unit =

17 match receive wch#m with

18 | ‘task(num, wch) ->

19 loop (send wch#m#result (fib num))

20 | ‘stop((), wch) -> wch

21 in loop wch

⬇

1let master () =

2 let rec loop (mch --!: [%kmc.check u]!\--) : unit =

3 match receive mch#u with

4 | ‘compute(x, mch) ->

5 let mch = send mch#w#task (x - 2) in

6 let mch = send mch#w#task (x - 1) in

7 let ‘result(r1, mch) = receive mch#w in

8 let mch = send mch#u#wip r1 in

9 let ‘result(r2, mch) = receive mch#w in

10 loop (send mch#u#result (r1 + r2))

11 | ‘stop((), mch) ->

12 send mch#w#stop ()

13 in loop mch

14

15let () =

16 let ut = Thread.create user () in

17 let mt = Thread.create master () in

18 let wt = Thread.create worker () in

19 List.iter Thread.join [ut;mt;wt]

u: $\texttt{um}!\textit{compute}$$\texttt{mu}?\textit{result}$$\texttt{um}!\textit{stop}$$\texttt{mu}?\textit{wip}$ |  m: $\texttt{um}?\textit{compute}$$\texttt{mw}!\textit{task}$$\texttt{mw}!\textit{task}$$\texttt{wm}?\textit{result}$$\texttt{mu}!\textit{wip}$$\texttt{wm}?\textit{result}$$\texttt{mu}!\textit{result}$$\texttt{um}?\textit{stop}$$\texttt{mw}!\textit{stop}$ |  w: $\texttt{mw}?\textit{task}$$\texttt{mw}?\textit{stop}$$\texttt{wm}!\textit{result}$
---|---|---

Figure 2: Example of kmclib program (top) and _inferred_ session types
(bottom).

Programming with kmclib. To enable safe message-passing programs, kmclib
provides two communication primitives, send and receive, and two primitives
for channel creation (KMC and %kmc.gen). We only give a user-oriented
description of these primitives here (see Appendix 0.A an overview of their
implementations).

The crux of kmclib is the session channel creation: [%kmc.gen (u,m,w)] at Line
1. This primitive takes a tuple of _role names_ as argument (i.e., (u,m,w))
and returns a tuple of communication channels, which are bound to
(uch,mch,wch). These channels will be used by the threads implementing roles
user (Lines 3-13), worker (Lines 15-21), and master (Lines 1-13). By default,
channels are implemented using concurrent queues from Multicore OCaml
(Domainslib.Chan.t) but other underlying transports can easily be provided.

Threads send and receive messages over these channels using the communication
primitives provided by kmclib. The send primitive requires three _arguments_ :
a channel, a destination role, and a message. For instance, the user sends a
request to the master with send uch#m#compute 20 where uch is the user’s
communication channel, m indicates the destination, and compute 20 is the
message (consisting of a label and a payload). Observe that a sending
operation returns a new channel which is to be used in the continuation of the
interactions, e.g., uch bound at Line 4. Receiving messages work in a similar
way to sending messages, e.g., see Line 6 where the user waits for a message
from the master with receive uch#m. We use OCaml’s pattern matching to match
messages against their labels and bind the payload and continuation channel.
See, e.g., Lines 7-10 where the user expects either wip or result message. The
receive primitive returns the payload res and a new communication channel uch.

New thread instances are spawned in the usual way; see Lines 15-18. The code
at Line 19 waits for them to terminate.

Figure 3: Examples of type errors.

Compatibility and error reporting. While the code in Figure 2 may appear
unremarkable, it hides a substantial machinery that guarantees that, if a
program type-checks, then its constituent threads are safe, i.e., no thread
gets permanently stuck and all messages that are sent are eventually received.
This property is ensured by kmclib using OCaml’s type inference and PPX
plugins to infer a session type from each thread then check whether these
session types are $k$-_multiparty compatible_ ($k$-MC) [13].

If a system of session types is $k$-MC, then it is safe [13, Theorem 1], i.e.,
it has the _progress_ property (no role gets permanently stuck in a receiving
state) and the _eventual reception_ property (all sent messages are eventually
received). Checking $k$-MC notably involves checking that all their executions
(where each channel contains at most $k$ messages) satisfy progress and
eventual reception.

The $k$-MC-checker [13] performs a bounded verification to discover the
_least_ $k$ for which a system is $k$-MC, up-to a specified upper bound $N$.
In the kmclib API, this bound can be optionally specified with [%kmclib.gen
$\mathit{roles}$ ~bound:$N$]. The $k$-MC-checker emits an error if the bound
is insufficient to guarantee safety.

The [%kmc.gen (u,m,w)] primitive also feeds the results of $k$-MC checking
back to the code. If the inferred session types are $k$-MC, then channels for
roles u, m and w can be generated. If $k$-MC cannot be guaranteed, then this
results in a type error. We have modified the $k$-MC-checker to return
counterexample traces when the verification fails. This helps give actionable
feedback to the programmer, as counterexample traces are translated to OCaml
types and inserted at the hole corresponding to [%kmc.gen]. This has the
effect of reporting the precise location of the errors.

To report errors in a function parameter, we provide an _optional_ macro for
types: [%kmc.check rolename] (see faded code in Line 2). Figure 3 shows
examples of such error reports. The left-hand-side shows the reported error
when Line 5 is commented out, i.e., the master sends one task, but expects two
result messages; hence progress is violated since the master gets stuck at
Line 9. The right-hand-side shows the reported error when Line 9 is commented
out. In this case, variable mch in Line 10 (master) is highlighted because the
master fails to consume a message from channel mch.

##  3 Inference of Session Types in kmclib

The kmclib API. The kmclib primitives allow the vanilla OCaml typechecker to
infer the session structure of a program, while simultaneously providing a
user-friendly communication API for the programmer. To enable inference of
session types from concurrent programs, we leverage OCaml’s structural typing
and row polymorphism. In particular, we reuse the encoding from [8] where
input and output session types are encoded as polymorphic variants and objects
in OCaml. In contrast to [8] which relies on programmers writing global types
prior to type-checking, kmclib infers and verifies local session types
automatically, without requiring any additional type or annotation.

Typed PPX Rewriter. To extract and verify session types from a piece of OCaml
code, the kmclib library makes use of OCaml PreProcessor eXtensions (PPX)
plugins which provide a powerful meta-programming facility. PPX plugins are
invoked during the compilation process to manipulate or translate the abstract
syntax tree (AST) of the program. This is often used to insert additional
definitions, e.g., pretty-printers, at compile-time.

A key novelty of kmclib is the combination of PPX with a form of type-aware
translation, whereas most PPX plugins typically perform purely syntactic
(type-unaware) translations. Figure 4 shows the workflow of the PPX rewriter,
overlayed on code snippets from Figure 2. The inference works as follows.

  1. 1.

The plugin reads the AST of the program code to replace the [%kmc.gen]
primitive with a hole, which can have any type.

  2. 2.

The plugin invokes the typechecker to get the typed AST of the program. In
this way, the type of the hole is inferred to be a tuple of channel object
types whose structure is derived from their usages (i.e., mch#u#compute).

To enable this propagation, we introduce the idiom “let (KMC $\ldots$) =
$\ldots$” which enforces the type of the hole to be monomorphic. Otherwise,
the type would be too general and this would spoil the type propagation (See §
0.B).

  3. 3.

The inferred type is translated to a system of (local) session types, which
are passed to the $k$-MC-checker.

  4. 4.

If the system is $k$-MC, then it is safe and the plugin instruments the code
to allocate a fresh channel tuple (i.e., concurrent queues) at the hole.

  5. 5.

If the system is unsafe, the $k$-MC-checker returns a violation trace which is
translated back to an OCaml type and inserted at the hole, to report a more
precise error location.

The translation is limited inside the [%kmc.gen] expression, retaining a clear
correspondence between the original and translated code. It can be understood
as a form of ad hoc polymorphism reminiscent of type classes in Haskell. Like
the Haskell typechecker verifies whether a type belongs to a class or not, the
kmclib verifies whether the set of session types belongs to the class of
$k$-MC systems.

Figure 4: Inferring session types from OCaml code.

##  4 Conclusion

We have developed a practical library for safe message-passing programming.
The library enables developers to program and verify arbitrary communication
patterns without the need for type annotations or user-operated external
tools. Our _automated verification_ approach can be applied to other general-
purpose programming languages. Indeed it mainly relies on two ingredients:
structural typing and metaprogramming facilities. Both are available, with a
varying degree of support, in, e.g., Scala, Haskell, TypeScript, and F#.

Our work is reminiscent of automated software model checking which has a long
history (see [9] for a survey). There are few works on inference and
verification of behavioural types, i.e., [18, 11, 12, 3]. However, Perera et
al. [18] only present a prototype research language, while Lange et al. [11,
12, 3] propose verification procedures for Go programs that rely on external
tools which are not integrated with the language nor its type system. To our
knowledge, ours is the first implementation of type inference for MPST and the
first integration of session types compatibility checking within a programming
language.

## References

  * [1] D. Brand and P. Zafiropulo.  On communicating finite-state machines.  J. ACM, 30(2):323–342, 1983.
  * [2] D. Castro, R. Hu, S. Jongmans, N. Ng, and N. Yoshida.  Distributed programming using role-parametric session types in Go: statically-typed endpoint apis for dynamically-instantiated communication structures.  PACMPL, 3(POPL):29:1–29:30, 2019.
  * [3] N. Dilley and Julien Lange.  Automated verification of Go programs via bounded model checking.  In International Conference on Automated Software Engineering (ASE). IEEE/ACM, July 2021.  To appear.
  * [4] J. Garrigue.  Relaxing the value restriction.  In Y. Kameyama and P. J. Stuckey, editors, Functional and Logic Programming, 7th International Symposium, FLOPS 2004, Nara, Japan, April 7-9, 2004, Proceedings, volume 2998 of Lecture Notes in Computer Science, pages 196–213. Springer, 2004.
  * [5] P. Harvey, S. Fowler, O. Dardha, and S. J. Gay.  Multiparty Session Types for Safe Runtime Adaptation in an Actor Language.  In A. Møller and M. Sridharan, editors, 35th European Conference on Object-Oriented Programming (ECOOP 2021), volume 194 of Leibniz International Proceedings in Informatics (LIPIcs), page 30, Dagstuhl, Germany, 2021. Schloss Dagstuhl – Leibniz-Zentrum für Informatik.
  * [6] K. Honda, N. Yoshida, and M. Carbone.  Multiparty asynchronous session types.  In POPL 2008, pages 273–284, 2008.
  * [7] R. Hu and N. Yoshida.  Hybrid session verification through endpoint API generation.  In FASE 2016, pages 401–418, 2016.
  * [8] K. Imai, R. Neykova, N. Yoshida, and S. Yuen.  Multiparty session programming with global protocol combinators.  In ECOOP, volume 166 of LIPIcs, pages 9:1–9:30. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2020.
  * [9] R. Jhala and R. Majumdar.  Software model checking.  ACM Comput. Surv., 41(4), Oct. 2009.
  * [10] D. Kouzapas, O. Dardha, R. Perera, and S. J. Gay.  Typechecking protocols with Mungo and StMungo.  In PPDP 2016, pages 146–159, 2016.
  * [11] J. Lange, N. Ng, B. Toninho, and N. Yoshida.  Fencing off Go: liveness and safety for channel-based programming.  In POPL 2017, pages 748–761, 2017.
  * [12] J. Lange, N. Ng, B. Toninho, and N. Yoshida.  A static verification framework for message passing in Go using behavioural types.  In ICSE 2018. ACM, 2018.
  * [13] J. Lange and N. Yoshida.  Verifying asynchronous interactions via communicating session automata.  In CAV (1), volume 11561 of Lecture Notes in Computer Science, pages 97–117. Springer, 2019.
  * [14] A. Miu, F. Ferreira, N. Yoshida, and F. Zhou.  Generating Interactive WebSocket Applications in TypeScript.  Electronic Proceedings in Theoretical Computer Science, 314:12–22, Apr. 2020.
  * [15] R. Neykova, R. Hu, N. Yoshida, and F. Abdeljallal.  A Session Type Provider: Compile-time API Generation for Distributed Protocols with Interaction Refinements in F$\sharp$.  In CC 2018. ACM, 2018.
  * [16] N. Ng, J. G. de Figueiredo Coutinho, and N. Yoshida.  Protocols by default - safe MPI code generation based on session types.  In CC 2015, pages 212–232, 2015.
  * [17] L. Padovani.  A simple library implementation of binary sessions.  J. Funct. Program., 27:e4, 2017.
  * [18] R. Perera, J. Lange, and S. J. Gay.  Multiparty compatibility for concurrent objects.  In PLACES 2016, pages 73–82, 2016.
  * [19] A. Scalas, O. Dardha, R. Hu, and N. Yoshida.  A linear decomposition of multiparty sessions for safe distributed programming.  In ECOOP 2017, pages 24:1–24:31, 2017.
  * [20] K. C. Sivaramakrishnan, S. Dolan, L. White, S. Jaffer, T. Kelly, A. Sahoo, S. Parimala, A. Dhiman, and A. Madhavapeddy.  Retrofitting parallelism onto ocaml.  Proc. ACM Program. Lang., 4(ICFP):113:1–113:30, 2020.
  * [21] N. Yoshida, R. Hu, R. Neykova, and N. Ng.  The Scribble protocol language.  In TGC 2013, volume 8358, pages 22–41. Springer, 2013.
  * [22] F. Zhou, F. Ferreira, R. Hu, R. Neykova, and N. Yoshida.  Statically verified refinements for multiparty protocols.  Proc. ACM Program. Lang., 4(OOPSLA):148:1–148:30, 2020.

##  Appendix 0.A Technical Details on the kmclib API

We explain the main communication primitives of kmclib and their translation
to session types. In particular, we reuse the encoding [8] where input and
output session types are encoded as polymorphic variants and objects, while
loops are naturally handled using equi-recursive types in OCaml.

Objects and variants in OCaml are structurally typed, which enables the
creation of ad-hoc types. This allows the channel structure to be derived from
the usage of channels in send and receive primitives.

Output types. Sending a message, e.g., in Line 4 of Figure 2, is parsed as
send (uch#m#compute) 42 where the two chained method calls yields a port for
sending compute label to role m, which in turn is passed to send (together
with a payload). This corresponds to an internal choice where the user
specifies a destination for its message and chooses a label from those offered
by the receiver.

The inferred type of channel object uch is a nested object type of the form
<m: <compute: (int, ch)> out> where m is a method that returns an object that
itself provides method compute (which returns a port for sending an int
payload and yielding a continuation channel ch). Note that the implementation
of these methods is not provided explicitly by the API nor the programmer,
instead they are constructed on-demand when invoking uch#m#compute; i.e.,
objects are generated automatically according to the method types that is
invoked on them (# denotes method invocation). Such object types correspond to
session types of the form m!compute<int>;ch (the translation is trivial).

Input types. To receive messages, as in Lines 6-12 of Figure 2, we use uch#m
to return a channel object which effectively corresponds to a port originating
from role m. This channel is then passed to the receive primitive, which
returns a polymorphic variant on which one needs to pattern match for expected
messages.

The inferred type of wch, which specifies the expected messages and their
respective continuation, is

<u: [‘compute of int * $\textit{ch}^{\prime}$ | ‘stop of unit * $\textit{ch}^{\prime\prime}$] inp>

This type corresponds to an external choice, i.e., session type of the form

{u?compute<int>;$\textit{ch}^{\prime}$} or
{u?stop<unit>;$\textit{ch}^{\prime\prime}$}.

Linearity. MPST require channels to be used linearly, i.e., each channel must
be used exactly once. If a channel is not used, this leads to a multiparty
compatibility issue (a message will not be sent/received), and hence our
implementation detects such issues statically via $k$-MC.

The idiomatic shadowing with the same variable names (e.g., re-binding of uch
in Line 4) in OCaml mitigates the risk of using a channel more than once. If
the program deviates from this best practice and a channel is used non-
linearly, an exception is raised at runtime.

Alternatively, kmclib provides an event-based alternative API (similar to that
of [22]), which eliminates the explicit need for linear channels, at the cost
of losing a direct-style API.111See
https://github.com/keigoi/kmclib/blob/tooldemo/test/paper/test˙handler.ml for
an example. We remark that there are other known ways to check linearity
statically [8], which can easily be adapted to our library.

⬇

1let KMC (uch,mch,wch) =

2 let um, mu, mw, wm = Chan.create_unbounded (), Chan.create_unbounded (), …
in

3 let uch = <m = <compute = Internal.make_out um (fun v -> ‘compute v) … > >
in

4 let mch = <u = Internal.make_inp mu … (fun v -> ‘compute v) > in

5 let wch = <m = Internal.make_inp mw … (fun v -> ‘stop v) > in

6 make_tuple (uch, mch, wch)

Figure 5: Code from Figure 2, instrumented at [kmc.gen]

##  Appendix 0.B Instrumented Code for Figure 2

Figure 5 shows the instrumented code for Figure 2. Line 2 allocates raw
channels using Chan.create_unbounded from Multicore OCaml, and Lines 3-5
create objects inhabiting the inferred type. We use shorthand <...> for in-
place objects object...end and abbreviate the continuations with an ellipsis.

The Internal functions make a channel from raw channels and a continuation. In
particular, make_out takes an extra function (fun v -> ‘label v), allocating a
variant tag representing the message label. Also, it uses type casts from Obj
module in OCaml, which is a common technique to implement session types in
OCaml (cf. [17]).

Line 6 (make_tuple) wraps the resulting tuple with the KMC constructor. As
mentioned in § 3, this makes the inferred hole type monomorphic. Normally, for
the top-level declarations, OCaml generalises the type to be polymorphic and
the hole type is inferred as $\forall\alpha.\alpha$ (can be instantiated with
any type at any site) if its occurrence is at the covariant position, spoiling
the propagation (cf. relaxed value restriction [4]). We avoid this by wrapping
the pattern with a type KMC : ’a -> ’a tuple declared explicitly as non-
covariant.

##  Appendix 0.C Error Reporting with Type Ascription

It is vital to show the location of the error to the programmer when an error
is found. To achieve this, the PPX plugin of kmclib instruments an extra
ascription of an incompatible type at the erroneous usage of a channel. For
example, see the error at Line 9 in the left of Figure 3, where the PPX plugin
assigns the variable mch a type [‘progress_violation] (a single variant
constructor type whose name is progress_violation), as the $k$-MC-checker
detects the input blocking forever at that point. Since it is used as <u:
[‘result of int * $\cdots$] inp > denoting an input of result with an int from
the user, the OCaml typechecker reports a type error.

For usability purposes, kmclib detects another kind of error, which we refer
to as format error. These errors happen when the inferred type of a channel is
not even in the form of output or input session type (channel misuse). For
example, if the one drops the role name (#m) writing send uch#compute 42, the
variable uch has the inferred type <compute: (int,$\cdots$) out>, which is not
a session type anymore. Figure 6 shows such an error. The highlighted part is
assigned a type [‘shoud_be_inp_or_out_object] type saying that the expression
needs another method call (or the expression should be used as an input). We
are planning to improve error messages to be more descriptive, e.g., as
[‘role_or_label_not_given].

Note that these format checks are all done within the [%kmc.gen] (or
[%kmc.check]) primitive. These errors could also be regarded as a “no
instance” error in the type class, as such ill-formatted types are not in the
class of $k$-MC systems (they are not even in the class of session type
syntax).

Figure 6: Example of a format error.


Vol.0 (20xx) No.0, 000–000

11institutetext: Yunnan Observatories, Chinese Academy of Sciences, Kunming
650216, P.R. China<EMAIL_ADDRESS>
22institutetext: Faculty of Information Engineering and Automation, Kunming
University of Science and Technology, No.727 South Jingming Rd., Chenggong
District, Kunming 650500, China
33institutetext: Center of Astrophysics, Guangzhou University, No.230,
Waihuanxi Rd, Panyu District, Guangzhou 510006, China
44institutetext: National Astronomical Observatory, Chinese Academy of
Sciences 20A Datun Rd., Chaoyang District, Beijing 100012, China
55institutetext: University of Chinese Academy of Sciences, 19A Yuquan Rd.,
Shijingshan District, Beijing 100049, China
Received 20xx month day; accepted 20xx month day

# A Study on Non-coplanar Baseline Effects for Mingantu Spectral
Radioheliograph

Qiu-ping Yang 112255 Feng Wang 331155 Hui Deng 33 Ying Mei 33 Wei Wang 44

###### Abstract

As a dedicated solar radioheliograph, the MingantU SpEctral RadioHeliograph
(MUSER) has a maximum baseline of more than 3000 meters and a frequency range
of 400 MHz – 15 GHz. According to the classical radio interferometry theory,
the non-coplanar baseline effect (i.e., w-term effect) would be considered and
calibrated for such a radio instrument. However, little previous literature
made the qualitative or quantitative analyses on w-term effects of solar
radioheliograph in-depth. This study proposes a complete quantitative analysis
of w-term effects for the MUSER. After a brief introduction of the MUSER, we
systematically investigate the baseline variations over a year and analyze the
corresponding variations of w-term. We further studied the effects of the
w-term in the imaging for the specified extended source, i.e., the Sun. We
discussed the possible effects of the w-term, such as image distortion and so
on. The simulated results show that the w-term is an essential and unavoidable
issue for solar radio imaging with high spatial resolution.

###### keywords:

Astronomical instrumentation: interferometric array – the Sun: interferometer
imaging – techniques: w-term correction

##  1 Introduction

Aperture-synthesis imaging technology (Ryle 1962) has been widely used in
radio astronomy because of its advantages, such as high resolution and quick
imaging. Radioheliograph is a type of radio interferometer specially designed
to observe the Sun. Several radioheliograph telescopes, such as Nobeyama
Radioheliograph (NoRH) (Nakajima et al. 1994), the Nancay Radioheliograph
(NRH) (Kerdraon and Delouis 1997; Kerdraon and Klein 2011), the Gauribidanur
Radioheliograph (Ramesh et al. 1998), the Siberian Solar Radio Telescope
(Grechnev et al. 2003)
(SSRT111http://en.iszf.irk.ru/The_Siberian_Solar_Radio_Telescope), and the
Frequency-Agile Solar Radiotelescope(FASR, (Bastian 2003; Gary 2003) had been
constructed or were under-constructed for monitoring solar activities.

The MUSER is a solar dedicated radio interferometric array with high temporal,
spatial, and spectral resolutions. It can simultaneously perform the spectral
and imaging observations of the full Sun in a wide frequency range(Yan et al.
2004, 2015, 2021). The MUSER has two sub-arrays: The MUSER-I contains 40
antennas with a 4.5m diameter for low frequency ranging from 400 MHz to 2 GHz.
The MUSER-II contains 60 antennas with a 2m diameter operating frequency
ranging from 2 GHz to 15 GHz. The specifications of the MUSER are listed in
Table 1 (Yan et al. 2011; Wang et al. 2015). All the 100 antennas are located
on three log-spiral arms with a maximum baseline length of about 3 km in
north-south and east-west directions (longitude=115.2505∘,
latitude=42.211833333∘, altitude=1365.0 meters).

Table 1: MUSER Characteristic and Performance Parameters | MUSER-I | MUSER-II
---|---|---
Antenna array | 40 (4.5 meters) | 60 (2 meters)
Frequency range | 0.4 GHz$-$2.0 GHz | 2.0 GHz$-$15 GHz
Frequency resolution | 64 channels | 520 channels
Channel bandwidth | 25 MHz
Time resolution | 25 ms | 206.25 ms
Spatial resolution | ${51.6}^{\prime\prime}-{10.3}^{\prime\prime}$ | ${10.3}^{\prime\prime}-{1.4}^{\prime\prime}$
Maximum Baseline length | more than 3000 m
Polarization | Dual circular left and right
Dynamic range | 25 db (snapshot)
Field of view | $0.6^{\circ}-7^{\circ}$

To obtain the high-quality image, a series of applications (Wang et al. 2015;
Wei et al. 2016; Mei et al. 2017; Wang and Yan 2019), such as data
preparation, calibration, gridding, and final image clean, had been developed
for the MUSER. Meanwhile, many possible errors were also calibrated and
corrected carefully. For example, the transport delay of the optical fiber and
the phase tracking error of calibrated source (i.e., geostationary satellite)
have been calibrated to guarantee the final quality image.

MUSER has a maximum baseline of more than 3000 meters and a frequency range of
400 MHz – 15 GHz. According to the classical radio interferometry theory, the
non-coplanar baseline effect (i.e., w-term effect) would be considered and
calibrated. However, there are few qualitative and quantitative analyses on
the non-coplanar baseline effect (i.e., w-term effect) for solar radio
observation. One possible reason is that the sun’s apparent diameter from the
earth is only about 32 arc minutes which is far from the concept of wide-field
imaging.

The w-term is a classical problem in radio interferometric imaging. There are
several methods to deal with the w-term effects of imaging. Perley (1989)
proposed the three-dimensional Fourier transform. However, because of its high
computational cost, this method is becoming more and more complicated when
dealing with the massive amount of observing data. The polyhedron imaging
method, or image plane faceting, was introduced by Perley (1989). An
alternative faceting method is $uv$ plane faceting (Sault 1999). The
w-projection algorithm is proposed based on the idea that $w$ could project
onto the $w=0$ plane (Cornwell et al. 2005; Cornwell 2008). The A-projection
algorithm pBhatnagar2008 extended upon this, also including instrumental
effects. W-stacking (Humphreys and Cornwell 2011; Offringa et al. 2014) is
another algorithm to deal with the w-term.

The rest of the paper is organized as follows. In section 2, a detailed
discussion of the w-term that arose from the non-coplanar baseline effects on
MUSER and the focus problem of w-term effects quantitative analyses. The
w-term effects on MUSER imaging performance are discussed in section 3.
Section 4 concludes the study and presents specific aspects of our future
work.

##  2 non-coplanar baseline effects Of MUSER

###  2.1 w-term issue of MUSER

For a radio interferometer, the received visibilities can be calculated
by(Thompson et al. 2008)

$\centering V(u,v,w)=\int\int\frac{I(l,m)}{\sqrt{1-l^{2}-m^{2}}}e^{-2\pi
i[ul+vm+w(\sqrt{1-l^{2}-m^{2}}-1)]}dldm\@add@centering$ (1)

where $(u,v,w)$ are the baseline coordinates, and the direction cosines
$(l,m,n)$ denote a position on the celestial sphere defined by
$l^{2}+m^{2}+n^{2}=1$. If the fields of view is sufficiently small, or the
observation time is short enough (snapshot observations), it can be simplified
by the approximation that $\sqrt{1-l^{2}-m^{2}}\approx 1$. Therefore, $2\pi
w(\sqrt{1-l^{2}-m^{2}}-1)$ is much less than unity. In such a situation, all
visibility data can be approximately considered to lie on a plane(Bhatnagar
1999). Equation 1 could be simplified as

$\centering V(u,v,w)=\int\int I(l,m)e^{-2\pi i(ul+vm)}dldm\@add@centering$ (2)

As observation fields of view becomes large, $2\pi w(\sqrt{1-l^{2}-m^{2}}-1)$
can no longer be ignored, thereby leading to w-term issues. Perley (1999)
proposed that w-term is the phase shift due to the dependence of direction
cosines and the value of $w$, which may be presented as

$\centering\phi(w)=2\pi w(\sqrt{1-{l}^{2}-{m}^{2}}-1)\@add@centering$ (3)

The maximum value of the term $2\pi w(\sqrt{1-l^{2}-m^{2}}-1)$ could be
expressed as $\frac{B\lambda}{D^{2}}$, where B is the longest baseline,
$\lambda$ is the observing wavelength, and $D$ is the diameter of the elements
of the array. The w-term can be ignored when

$|w(\sqrt{1-l^{2}-m^{2}}-1)|\ll 1$ (4)

Unlike ordinary radio interferometric arrays, MUSER’s observation target is
the Sun. Due to the revolution, the apparent position of the Sun changes
continuously during the year, thus changing the value of $(u,v,w)$. A
schematic diagram of the MUSER imaging result is given in Figure 1. The white
circle range is the solar disk. Since the diameter of the Sun is about 32 arc
minutes, the field of view of MUSER observation imaging is considered to be 1
square degree under normal conditions, which can effectively show the imaging
results at different heights of the Sun.

Figure 1: The observational clean image of MUSER-I, 1.7125GHz, 2015-11-01
04:08:49.354161240(UTC).

According to the study of Taylor et al. (1999), we can allow the w-term effect
to be as high as 0.1 radians without introducing severe errors in the image.
When the field of view is 32 arc minutes, therefore, a simple calculation
shows that the value of $w$ needs to be less than 1154.06 wavelength numbers
for Equation 4 to hold.

To quantify the w-term effects of MUSER, we calculated the w-term variation
for different frequencies of MUSER-I. Considering that the low frequencies are
more susceptible to the w-term, we calculated the $w$-value variation at noon
(China Standard Time (CST), UTC+8) for each day in a year(2020) using the 1
GHz observation of MUSER-I as an example. Figure 2 shows that the maximum
value of $w$ in 2020 was about 9412.12 wavelength numbers, which occurred on
December 11, and similarly, the minimum value was 3156.74 wavelength numbers
on June 18.

Figure 2: The maximum w value of MUSER at 12:00 noon (CST) every day in the
year 2020

We calculated the variation of the $w$ value during a day as well, giving the
distribution of the maximum $w$ value from 8:00 to 17:00 CST time with a
frequency of 1GHz(see Figure 3). In practice, observations are generally
considered to be made when the elevation angle of the Sun is greater than 15
degrees. Therefore, we can obtain that the observable time on December 11 was
from 9:45 to 14:45(CST), with the $w$ peak value of 10,141.84 at 13:45 and the
minimum value of 8,783.24 at 11:15; the $w$ peak value on June 18 was 8,408.73
at 17:00, and the minimum value was 3,156.74 at noon.

The results show that the $w$ value of MUSER was always greater than 1154.06
wavelength numbers during the year. Obviously, for the observation of the
MUSER-I with 400 MHz, the value of w-term would be much larger than unity.
W-term needs to be more explicitly dealt with, or the operation of the
instrument will be significantly restricted.

Figure 3:  Maximum $w$ value from 8:00(am) to 17:00(pm)(CST), w unit:number of
wavelength.

###  2.2 Quantitative analysis of Phase shift

The non-coplanar baseline effects will directly affect the phase of the
visibility data at different distances from the center of the solar disk. We
analyzed the phase shifts of MUSER using Equation 3 and two extremes of the
w-value already obtained, i.e., the maximum value of 9412.12 wavelength
numbers and the minimum value of 3156.74 wavelength numbers. The observational
frequency was 1GHz as well. In Figure4(a), the horizontal axis represents the
field of view. Moreover, the midpoint of the horizontal axis is the phase
center, and the value of the vertical axis represents the value of phase
shift. It can be seen that at the center of FOV, the phase shift is 0.
However, as the distance from the phase center increases, the phase shift will
keep increasing. And the larger the value of $w$, the larger the phase shift
increases with the distance. We could find that the phase shift of the maximum
$w$ value(9412.12) is almost three times that of the minimum $w$
value(3156.74) at the FOV of 32’. Furthermore, it can be derived that at the
location of the solar radius (16’), the phase shift increases by 0.2785
radians with respect to the location of 12’, while at the location of 32’, the
phase shift increases by 1.9215 radians compared to the location of 16’.

To make the phase variation more intuitive, we give the phase screen image
($e^{2\pi iw(\sqrt{1-l^{2}-m^{2}}-1)}$)for December 11, 2020 at 12:00:00
(CST), with $w$ value is 9412.12, as shown in Figure 4(b). When the value
phase shift is 0, the image value is 1.

Figure 4:  Diagrams of phase shifts with different baselines and phase screen
of maximum $w$. (a) shows the value of phase shifts. (b) shows the phase
screen at 12:00, December 11, 2020.

##  3 Imaging Quality Assessment

To analyze the impact of different w-term values on the final imaging and to
provide a reasonable basis for subsequent MUSER data processing, we adopted a
data simulation approach using the Radio Astronomy Simulation, Calibration and
Imaging Library (RASCIL222https://gitlab.com/ska-telescope/external/rascil) to
simulate the MUSER-I imaging observations. RASCIL is a fully open-source radio
interferometric array data processing package that has been widely used in SKA
data simulation and processing studies. Its processing results have been
compared with other software and confirmed to be reliable. The RASCIL is
developed in Python. The Gridder uses the corresponding module of
Ducc333https://gitlab.mpcdf.mpg.de/mtr/ducc.

We defined a MUSER-I configuration file for RASCIL, giving the longitude,
latitude, and altitude values of MUSER-I. In addition, we carefully configured
the antenna diameter, polarization mode, coordinates, and other related
parameters. In the simulated observations, the observation frequency is 1 GHz,
and the polarization mode is set to Stocks-I. The sky model used for the
simulation consists of 49 point sources (see Figure 5 (a)), and the flux of
point sources is 10 Jy/beam. The distance between the leftmost and rightmost
point sources is about 32’, which is consistent with the solar diameter. Based
on the results of the previous section, we selected 12:00 (CST) on December
11, 2020, for the simulation observations as well. The right ascension of
observation is 258.6503(degrees), and the declination of observation is
-23.0267(degrees). Meanwhile, the observation hour angle is 6.2296(radians).
The observation uv coverage is in Figure5(b). The imaging cell size is 6.81
arcsecond, the image size is $512\times 512$. Therefore the FOV of the image
is 0.968 degrees.

During the processing of simulated observational data, at first we adopted the
classical imaging method of radio interferometry (see Equation2), i.e.,
gridding the received visibilities, and then obtaining the dirty image by
Fourier transform. The non-coplanar baseline effects are not corrected in this
imaging process. Then the w-projection (Cornwell 2008; Muscat 2014) algorithm
is used to imaging, meanwhile to correct the non-coplanar baseline effects. In
w-projection, the number of $w$-layers is obtained by using the methods of
Humphreys and Cornwell (2011).

Figure 5: (a) sky model of 2020-12-11-12:00:00(CST).(b)UV coverage of MUSER-I
.

Figure 6: (a) dirty image with non-coplanar baseline effects.(b)dirty image
without non-coplanar baseline effects(w-projection). (c)the differences image
of dirty image with/without non-coplanar baseline effects.(d) the sources
centroid locations in dirty images.

After the simulation observation and imaging processing, we present the
corresponding results in Figure 6. Figure 6(a) and Figure 6(b) shows the dirty
images with/without non-coplanar baseline effects, respectively. Figure 6(c)
shows the result of subtracting Figure 6(a) and Figure 6(b).

In dirty images obtained by both methods, the point source at the phase center
(right ascension:258.6504, declination:-23.0264, unit:degrees) has a flux of
10.93 Jy/beam, but the point source at the location (258.6504, -22.7673) has a
flux of 10.43 Jy/beam in Figure6(a) and 10.81 Jy/beam in Figure6(b). Point
sources farther from the phase center, such as at the location (258.9313,
-22.7673) have a flux of 8.95 Jy/beam in Figure6(a) and 9.82 Jy/beam in
Figure6(b). It can be noticed that non-coplanar baseline effects would have
different effects on fluxes of point sources at different locations.
Therefore, in the dirty image without w term correction, we further
investigated the relationship between the fluxes of 49 point sources and their
distances from the phase center. In Figure7(a) the horizontal axis denotes the
distances from the locations of the point sources to the phase center, and the
vertical axis indicates that the flux attenuation of the point source varies
with the distance from the phase center. At a distance of 18.98 arc-minutes
from the phase center, the ratio is about 80$\%$.

Figure 7: (a) diagram of point source flux attenuation without w-term
correction.(b)statistical result of centroid location offset in dirty image
without w-term correction.

The point sources at the edge are more severely affected by the non-coplanar
baseline effects, and the locations of the point sources are shifted. By
calculating the centroid locations of the point sources in the dirty images,
we obtained the centroid locations offset distribution of the point sources in
both imaging methods. In Figure 6(d), the red markers are the centroid
locations of point sources in the dirty image without the non-coplanar
baseline effects correction, and the green markers are the centroid locations
of point sources in the dirty image, which is obtained by the w-projection
algorithm. It can be inferred that the point sources have offset in centroid
locations when non-coplanar baseline effects are not corrected in the imaging
process.

Therefore, we statistic the centroid locations offset data for all point
sources in Figure6(a), and the results are shown in Table2. We could see that
the centroid locations offset in right ascension is 0.494 arc seconds, and in
declination is 9.848 arc seconds at the location (258.9313, -23.2858).
Furthermore, Figure7(b) shows the relationship between the point source
location and the mean value of the point source centroid location offset. At a
distance 18.98 arc minutes from the phase center, the offset mean value is
about 7 arc seconds. It can be concluded that the farther the point sources
are from the phase center, the larger the centroid position offsets due to the
non-coplanar baseline effect.

Table 2: Offset of point sources centroid locations in the dirty image without
non-coplanar baseline effects correction (unit:arc-second)

RA$\backslash$ DEC(degrees) | -23.28578874 | -23.20799939 | -23.08696216 | -23.02644368 | -22.9659252 | -22.844888 | -22.76734809
---|---|---|---|---|---|---|---
258.9313419 | (0.494, 9.848) | (0.413, 6.573) | (0.851, 6.887) | (1.990, 4.108) | (0.846, 7.669) | (0.257, 9.146) | (-1.429, 9.196)
258.8472499 | (0.190, 4.327) | (-1.169, 3.986) | (0.124, 1.583) | (0.704, 2.468) | (0.286, 1.136) | (0.884, 4.402) | (0.321, 6.373)
258.7159846 | (-0.170, 1.803) | (0.608, -0.088) | (-0.005,-3.049) | (0.383, -1.652) | (0.462, -2.110) | (1.223, -0.275) | (0.759, 2.826)
258.6503521 | (0.025, 3.560) | (0.719, 2.469) | (-0.768, -0.367) | (0.066, -0.030) | (0.323, 1.459) | (-0.808, 1.815) | (0.075, 5.320)
258.5847195 | (-1.093, 6.825) | (-1.362, 3.791) | (-0.695, 1.950) | (-0.226, 1.361) | (-0.067, 2.821) | (-0.459, 4.698) | (-0.173, 7.392)
258.4534543 | (-0.513, 6.872) | (-1.224, 3.482) | (-0.790, 3.267) | (-0.772, 1.814) | (-0.240, 3.119) | (0.907, 4.749) | (-0.794, 8.811)
258.3693622 | (0.834, 7.954) | (-0.076, 6.122) | (-0.487, 3,295) | (-1.322, 5.746) | (-0.398, 3.530) | (-0.008, 7.637) | (-0.330, 7.976)

From the above imaging results, it can be concluded that MUSER is affected by
the non-coplanar baseline effect during imaging. This effect is more
noticeable when the $w$ value is larger. Moreover, since the non-coplanar
baseline effect is proportional to the FOV size, the edge part of the solar
disk is more affected than the central part (phase center).

##  4 Conclusions and Future work

In this study, we investigate the non-coplanar baseline effects on MUSER,
which has an almost 3000-meter long baseline. Experiments on simulation
observation have shown that the maximum $w$ values of MUSER are from 3156.74
(number of wavelengths) to 9412.12 (number of wavelengths), which would affect
the imaging quality of observation. We demonstrate the non-coplanar baseline
effects in the imaging process by two imaging algorithms, including the
classical Fourier transform algorithm and w-projection. It can be concluded
that the non-coplanar baseline effects should be considered in MUSER’s imaging
processing for almost all the year if we want to get a high dynamic and finely
solar image.

For radioheliograph imaging, a long baseline can improve the imaging
resolution. However, this study presents that although the field of view is
small for solar observations, non-coplanar effects can significantly impact
solar observation imaging, especially for measurements related to precise
positions.

Furthermore, to obtain higher quality observation images for MUSER,
investigating of imaging algorithms for correcting non-coplanar baseline
effects should be taken into account. There are two issues to be considered.
The first is to compare and choose proper methods for non-coplanar baseline
effects correction. Besides the w-projection mentioned above, there are
additional imaging algorithms that are used to correct the non-coplanar
baseline effects. As these algorithms have different features, it is necessary
to find one that is the best suitable for the MUSER data processing pipeline.
The second is to meet the requirements of real-time imaging. MUSER-I would
output 1.92 GB of observation data per minute and MUSER-II would reach 3.6 GB
per minute (Wang et al. 2015). How to correct the non-coplanar baseline
effects for these data in real-time is also a problem of MUSER’s data
processing pipeline. In the future, the focus will be on the running
efficiency of various imaging algorithms, and apparently, parallel computing
based on the multi-core CPU or GPU will be a good choice.

###### Acknowledgements.

This work is supported by the National SKA Program of China (2020SKA0110300),
the Joint Research Fund in Astronomy (U1831204, U1931141) under cooperative
agreement between the National Natural Science Foundation of China (NSFC) and
the Chinese Academy of Sciences (CAS), the Funds for International Cooperation
and Exchange of the National Natural Science Foundation of China
(11961141001), the National Natural Science Foundation of China (No.11903009).

## References

  * Ryle (1962) Ryle M.,1962 ,Nature,194,517
  * Bastian (2003) Tim S. Bastian, 2003, Innovative Telescopes and Instrumentation for Solar Astrophysics, proceedings of SPIE, Vol.4583, doi:10.1117/12.460293
  * Bhatnagar (1999) Sanjay Bhatnagar,1999,Notes from summer school on low frequency astronomy at NCRA, Pune
  * Bhatnagar et al. (2008) S. Bhatnagar , T. J. Cornwell,K. Golap , J. M. Uson,2008, A&A, 487,419
  * Cornwell (2008) Tim J. Cornwell, 2008, IEEE Journal of Selected Topics in Signal Processing, 2,793
  * Cornwell et al. (2005) Tim J. Cornwell , K. Golap , S. Bhatnagar ,2005, Astronomical Society of the Pacific Conference Series, 347
  * Gary (2003) D. Gary, 2003,Journal of The Korean Astronomical Society,36,135
  * Grechnev et al. (2003) V.V. Grechnev , S.V. Lesovoi , G. Ya. Smolkov et al.,2003, Solar Physics,216,239
  * Humphreys and Cornwell (2011) B. Humphreys , T. J. Cornwell,2011,Square Kilometre Array Memo,132
  * Kerdraon and Delouis (1997) Kerdraon Alain , Delouis Jean-Marc,1997, The Nançay Radioheliograph, Coronal physics from radio and space observations, 192
  * Kerdraon and Klein (2011) Kerdraon Alain , Klein Ludwig Karl,2011, General Assembly and Scientific Symposium, 2011 XXXth URSI
  * Mei et al. (2017) Mei Ying , Wang Feng , Wang Wei , Chen Linjie et al.,2017,PASP, 130,014503
  * Muscat (2014) Muscat Daniel,2014, arXiv preprint arXiv:1403.4209
  * Nakajima et al. (1994) H. Nakajima , M. Nishio , S. Enome, K. Shibasaki , T. Takano , Y. Hanaoka et al.,1994, Proc.IEEE,82,702
  * Offringa et al. (2014) A.R. Offringa , B. McKiley et al.,2014,MNRAS, 444,606
  * Perley (1989) R. A. Perley,1989, Astronomical Society of the Pacific Conference Series,6
  * Perley (1999) R. A. Perley, 1999,Synthesis imaging in radio astronomy II,180,383
  * Ramesh et al. (1998) R. Ramesh , K. R. Subramanian , M. S. Sundararajan , C. V. Sastry,1998, Sol.Phys,181,439
  * Sault (1999) R. J. Sault et al., 1999, A&A, 139,387
  * Taylor et al. (1999) Greg B.Taylor , Chris Luke Carilli, Richard A Perley,1999, Synthesis Imaging in Radio Astronomy II,180
  * Thompson et al. (2008) A. R. Thompson, J. M. Moran, G.W.Swenson Jr, 2008, John Wiley & Sons
  * Wang et al. (2015) Wang Feng, Mei Ying, Wang Wei , Ji, Kaifan,2015,PASP,127,950
  * Wang and Yan (2019) Wang Wei, Yan Yihua, 2019,Advances in Astronomy,2019
  * Wei et al. (2016) Wei Shoulin , Wang Feng , Deng Hui,Liu Cuiyin , DaiWei, Liang Bo ,Mei Ying,Shi Congming, Liu Yingbo ,Wu Jingping,2016,PASP,129,024001
  * Yan et al. (2015) Yan Yihua ,Chen Linjie,Yu Sijie,2015,International Astronomical Union. Proceedings of the International Astronomical Union,11,427
  * Yan et al. (2011) Yan Yihua ,Zhang Jian , Chen Zhijun , Wang Wei, Liu Fei, Geng Lihong, 2011,General Assembly and Scientific Symposium, 2011 XXXth URSI
  * Yan et al. (2004) Yan Yihua, Zhang Jian,Huang Guangli,2004, Radio Science Conference, 2004. Proceedings. 2004 Asia-Pacific
  * Yan et al. (2021) Yan Yihua ,Chen Zhijun,Wang Wei,Liu Fei, Geng Lihong,Chen Linjie , Tan Chengming ,Chen Xingyao et al., 2021,Frontiers in Astronomy and Space Sciences


A. Ivascenko
<EMAIL_ADDRESS>
# Semi-analytical model of cosmic ray electron transport

A. Ivascenko  F. Spanier  Universität Würzburg, ITPA, Lehrstuhl für
Astronomie, Emil-Fischer-Str. 31, 97074 Würzburg, Germany

(12 November 2010; 25 January 2011; 14 February 2011)

###### Abstract

We present a numerical extension to the analytical propagation model
introduced in Hein and Spanier (2008) to describe the leptonic population in
the galactic disc. The model is used to derive a possible identification of
the components that contribute to the leptonic cosmic ray spectrum, as
measured by PAMELA, Fermi and HESS, with an emphasis on secondary
$e^{+}-e^{-}$ production in collisions of cosmic ray particles with ambient
interstellar medium (ISM). We find that besides secondaries, an additional
source symmetric in $e^{+}$ and $e^{-}$ production is needed to explain both
the PAMELA anomaly and the Fermi bump, assuming a power-law primary electron
spectrum. Our model also allows us to derive constraints for some properties
of the ISM.

Recently the leptonic component of the cosmic ray spectrum has gained new
attention. New observations from ATIC, PAMELA and Fermi show a deviation from
a power-law in the form of an excess in both the electron and positron
spectra. Annihilating dark matter (Allahverdi et al., 2009) and nearby pulsars
(Büsching et al., 2008; Grimani, 2009; Blasi and Amato, 2011) (among other
hypotheses) have been proposed as possible sources of the excess leptons.
Regardless of the source, a propagation model is needed to connect the energy
spectrum measured on Earth with the injection spectra.

We present our numerical cosmic ray transport model in application to the high
energy electron transport in the ISM. Spatial and momentum diffusion, particle
escape, acceleration via Fermi I and continuous energy losses were taken into
account and their effects on the steady-state energy spectrum analyzed. In
solving the transport equation we employed quasi-linear transport theory, the
diffusion approximation and a separation of the spatial and momentum problem
to obtain the leaky-box-equation, which was then solved numerically. The
spatial problem was solved analytically in cylindrical and prolate spheroidal
coordinates.

The transport model was employed to calculate the spectrum of secondary
electrons in our galaxy. We assume the leptons from pioArXivn decay to be the
dominating component of the secondary spectrum. The lepton spectrum from pion
collisions of highly relativistic cosmic ray protons with thermal protons in
the ISM was calculated using the parametrization of pion production in p-p-
collisions from Kelner et al. (2006) and used as the injection spectrum for
the transport model. With realistic simulation parameters the resulting
positron flux in the vicinity of the solar system lies remarkably close to the
low-energy PAMELA data. Assuming a generic power-law primary $e^{-}$ spectrum
and an additional $e^{+}-e^{-}$ source allows us to fit the datasets from
PAMELA, Fermi and HESS very nicely and to put constraints on several transport
model parameters and the corresponding ISM properties.

##  1 The data

The leptonic cosmic ray spectrum gained new attention when the results of the
balloon experiment ATIC were published (Chang et al., 2008). The authors
claimed an anomalous excess in the spectrum just below 1 TeV. This started
wild speculations about the nature of the source for these leptons, which
should be in the relative vicinity of the solar system.

Those speculations were fueled further by new measurements from the PAMELA
satellite published in Adriani et al. (2009) claiming a rise in the positron
spectrum above 10 GeV. The lepton spectrum from the Fermi satellite (Abdo et
al., 2009) could not confirm the ATIC-peak quantitatively but showed a
deviation from a straight power-law in the same energy region.

The HESS atmospheric Cerenkov telescope extended the spectrum measurement
beyond 1 TeV showing a cut-off or a steeper power-law but could also not
confirm the ATIC excess (Aharonian et al., 2009).

Recently the Fermi lepton spectrum was extended to lower energies (7 GeV) in
Ackermann et al. (2010). Also the PAMELA results were updated with new
experimental data and different evaluation techniques in Adriani et al.
(2010).

To interpret these measurements, it is necessary to connect the injection
spectra at the (possible) sources with the observations in the solar system.
To do that we need a transport model that includes all the relevant processes
in the ISM (spatial and momentum diffusion, cooling, particle escape, etc.).

##  2 Transport model

The derivation of the transport equation follows closely Lerche and
Schlickeiser (1988) and Hein and Spanier (2008). Using Quasi Linear Theory and
the diffusion approximation the Vlasov equation is transformed into the
steady-state diffusion-convection equation for the phase space density
$f(\mathbf{x},p)$, which can be written in the form

$\mathcal{L}_{x}f+\mathcal{L}_{p}f+S(\textbf{x},p)=0$ (1)

with the spatial operator

$\mathcal{L}_{x}(\textbf{x},p)\equiv\nabla\left[\kappa(\textbf{x},p)\nabla\right]$
(2)

and the momentum operator

$\mathcal{L}_{p}(\textbf{x},p)\equiv\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left(p^{2}D(\textbf{x},p)\frac{\partial}{\partial
p}-p^{2}\dot{p}(\textbf{x},p)\right)-\frac{1}{\tau(\textbf{x},p)}\textrm{.}$
(3)

This equation describes the diffusion (coefficients $\kappa$ and $D$ in space
and momentum, respectively), acceleration and cooling (rates $\dot{p}$) and
escape (timescale $\tau$) of an injected source distribution $S$.

Following Lerche and Schlickeiser (1988) the equation can be solved by
separating the operators and the source term in their spatial and momentum
dependencies. Then the resulting steady state particle distribution can be
written as an infinite sum:

$f(\textbf{x},r)=\sum_{i}A_{i}(\textbf{x})R_{i}(p)$ (4)

The spatial coefficients $A_{i}(\textbf{x})$ can be acquired by solving the
diffusion equation in the appropriate geometry

$\frac{\partial}{\partial
u}A_{i}\exp(-\omega_{i}^{2}u)=\nabla\left(\kappa\nabla
A_{i}\exp(-\omega_{i}^{2}u)\right)\textrm{,}$ (5)

where $\omega_{i}^{2}$ take the role of inverse escape times for the momentum
modes and $u$ is a convolution variable.

In this work we used the analytical solution for cylindrical geometries as
provided by Wang and Schlickeiser (1987) since this offers the best reflection
of the symmetries of a spiral galaxy. Analytical solutions in spherical and
prolate spheroidal coordinates are also available (Schlickeiser et al., 1987;
Hein and Spanier, 2008) and can be used to describe particle transport in
elliptical galaxies and galaxy clusters.

The momentum modes $R_{i}(p)$ are solutions of the ODE, commonly referred to
as the leaky-box equation:

$\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[\frac{a_{2}p^{4}}{\kappa(p)}\frac{\partial R_{i}}{\partial
p}-\frac{a_{1}p^{3}}{\kappa(p)}R_{i}+\dot{p}R_{i}\right]-\omega_{i}^{2}\kappa(p)R_{i}=-Q(p)$
(6)

The parameters $a_{1}$ and $a_{2}$ represent the absolute scales for Fermi I
and II processes and depend mainly on the Alfven and shock speed,
respectively. As Lerche and Schlickeiser (1988) have shown, an analytical
solution for the equation can be found, if the cooling rate $\dot{p}$ scales
linearly with the particle energy, which is a good assumption for high energy
protons and heavier nuclei. Electrons above a few GeV, however, are cooled
primarily by synchrotron radiation and IC scattering, both scaling
quadratically in particle energy. Therefore the leaky-box equation was solved
using numeric relaxation methods.

To summarize, the model describes a particle distribution that is injected by
a homogeneous population of time-independent sources in the galactic disc
(injection region), diffuses throughout the galactic disc and halo
(confinement region), cooling by bremsstrahlung, adiabatic deceleration,
synchrotron radiation and inverse Compton scattering, and leaves the
confinement region at some point. The matter, photon and magnetic field
distributions are assumed to be homogeneous over the whole confinement region.

A model with a uniform source distribution is not a good choice to simulate
high energy leptons from discrete sources (SNR, PWN…), since the short energy
loss timescales make the source distribution critical for the resulting
spectra, as discussed in Blasi and Amato (2011). Conversely, production of
secondaries from p-p collisions in the galactic halo can be approximated very
well as a homogeneous source, since there are no dramatic variations in the CR
and ISM densities.

A similar analysis by Moskalenko and Strong (1998) used an analogue model as
in this paper, but obviously did not include the new data above 1 GeV.
Stephens (2001) has basically the same drawbacks, in addition this work has a
more sophisticated secondary production model. Compared with the recent works
in this area like Blasi and Amato (2011), we concentrate more on the low
energy end of the current PAMELA data and the constraints it can put on
propagation parameters if interpreted as secondaries.

##  3 Secondary leptons

Fortunately there is a lepton production process that can be approximated with
a homogeneous spatial distribution: the secondary leptons from the decay of
products of proton collisions. To get a galaxy-averaged flux of secondaries we
consider collisions of high energy cosmic ray protons with ambient thermal
matter, assuming that both populations are homogeneously distributed with
average densities.

The pion production cross-section and the electron source function for a
single p-p collision used here are parametrizations of results from the SIBYLL
event generator from Kelner et al. (2006). Folding the source function
$F_{e}\left(\frac{E_{e}}{E_{p}},E_{e}\right)$ with the cosmic ray spectrum
$J_{p}(E_{p})$ and the inelastic collision cross-section
$\sigma_{\mathrm{inel}}(E_{p})$ yields the electron flux:

$\frac{dN_{e}}{dE}=4\pi
n_{H}\int\sigma_{\mathrm{inel}}(E_{p})J_{p}(E_{p})F_{e}\left(\frac{E_{e}}{E_{p}},E_{e}\right)\frac{dE_{p}}{E_{p}}$
(7)

Since the source function is sharply peaked, the CR spectrum can be cut off
exponentially at the first knee ($\approx 1000$ TeV) with negligible
deviations to the lepton spectrum below 10 TeV:

$J_{p}(E_{p})=0.252\left(\frac{E_{p}}{\mathrm{TeV}}\right)^{-2.677}e^{\frac{-E_{p}}{10^{3}\
\mathrm{TeV}}}\mathrm{m}^{-2}\mathrm{sr}^{-1}\mathrm{s}^{-1}\mathrm{TeV}^{-1}$

This parametrisation of the CR spectrum was obtained by fitting the CR data
between $10^{10}$ and $10^{14}$ eV in Cronin et al. (1997).

The average thermal matter density $n_{H}$ is treated as a free parameter. The
resulting electron flux with a slightly harder power-law (spectral index
$2.62$) is used as the injection function in the transport model. Since the
parametrizations in Kelner et al. (2006) only consider p-p interactions, the
actual flux of secondaries is underestimated by the contribution from heavier
nuclei.

##  4 Parameter constraints

Figure 1: Parameter studies for the secondary positron spectrum compared to
the PAMELA $e^{+}$ fraction multiplied by a fit of the Fermi $e^{+}+e^{-}$
data. Combined systematic errors of PAMELA and Fermi are shown by the grey
area. The thick black line with $n_{H}=1.5\cdot 10^{3}\,\mathrm{m}^{-3}$,
$B=3\,\mathrm{\mu G}$ and $t=t_{\mathrm{esc}}$ yields the best fit of the
positron spectrum, if combined with an extra high energy source as shown in
Fig. 2. The thin and dashed black lines show the slope change with changing
escape losses.

The principal shape of the resulting steady-state lepton spectrum is shown in
Fig. 1. The cooling, diffusion and escape processes introduce breaks in the
injected power-law, so that multiple energy intervals with different spectral
indices are formed, each of them defined by a different process. At the first
break the injected power-law is steepened by $1/3$ as the diffusion/escape
time scale becomes shorter than the time scale of bremsstrahlung and adiabatic
cooling (linear in energy). The second break appears when the time scale for
synchrotron and IC losses becomes shorter that the diffusion time scale. These
quadratic losses steepen the spectrum by 1. The absolute flux value is
determined mainly by the ISM density, since the injected secondary spectrum
has a linear $n_{H}$ dependence (Eq. (7)), and to a lesser extent by escape
and cooling losses, which transport particles out of the simulation region.

So the model has three major parameters: the linear energy loss time scale
that reflects the thermal gas density $n_{H}$, the quadratic loss time scale
that incorporates the energy densities of the EM field (composed of a constant
CMB density of $\approx 0.3$ eV/cm3 and a variable magnetic field $B$) and the
escape time scale $t_{\mathrm{esc}}$ that corresponds to the size of the
confinement region.

Measurements of the positron spectrum, in particular the newest results from
PAMELA, allow us to constraint those parameters. A rather conservative
statement can be made, if we assume that the secondaries have a negligible
contribution to the positron spectrum. As shown in Fig. 1, assuming an average
ISM density of $\approx 10^{3}$ m-3, the confinement region has to be less
than the assumed $30$ kpc or the diffusion process has to be more efficient,
so that the particles leave the galaxy 10 times faster (black dotted line). A
higher ISM density would require an even more efficient particle escape, so as
not to overshoot the measured positron spectrum. Since the spectral shape
wouldn’t play any role, the magnetic field can be almost arbitrary.

*

Figure 2: A fit of Fermi (Ackermann et al., 2010) and HESS (Aharonian et al.,
2009) lepton data (solid red line) and PAMELA (Adriani et al., 2010) positron
data (solid green line). Grey areas show the systematic errors. The components
that contribute to the spectra are secondary positrons (dashed green), a
generic power-law primary spectrum (dashed red) and an additional symmetric
$e^{+}-e^{-}$ source (solid and dashed black).

On the other hand, the PAMELA data can be fitted very well with a
superposition of two power-laws with a spectral index of about $3.6$ for the
low energy part, corresponding remarkably well with the secondary spectrum
steepened by synchrotron and IC losses. Assuming the low energy positrons to
be secondaries puts much stricter constraints on the parameters, since the
second break in the spectrum has to appear at $\approx 1$ GeV to match the
data. The break energy is determined by the equilibrium of the
diffusive/escape losses and the radiative losses. The thick black line in Fig.
1 is the best fit with $n_{H}=1.5\cdot 10^{3}$ m-3, $B=3\,\mu G$ and
$t=t_{\mathrm{esc}}$. The radiative losses in the galactic halo are dominated
by IC scattering on the CMB photons, so that lowering the magnetic field has
no effect on the break energy. A slight flattening at the low energy end of
the spectrum is needed for the best fit of the new PAMELA data. This
corresponds to a halo radius of $30$ kpc and a spatial diffusion coefficient
$\kappa(E=1\,\mathrm{GeV})=4.5\cdot 10^{24}$ m2 s-1, yielding an average
confinement time $t_{\mathrm{esc}}=2\cdot 10^{17}$ s. The absolute flux fixes
the third parameter, the ISM density, to $1.5\cdot 10^{3}$ m-3. Taking into
consideration that about $99\%$ of the confinement region represent the
galactic halo with an average gas density of about $10^{3}$ m-3 (McKee and
Ostriker, 1977), the parameters seem very realistic.

##  5 Lepton excess

In Fig. 2 the latest Fermi lepton measurement, in which the spectrum has been
extended down to $7$ GeV, is shown. A generic fit of it was used to calculate
a positron flux from the PAMELA positron fraction measurement. This data is
used because it was taken simultaneously in the same solar cycle making
interpretation easier. The HESS high-energy data was shifted down by $15\%$ to
better coincide with the Fermi data. This is still within the systematic error
margins of HESS, not to mention the systematic error of Fermi being of the
same order of magnitude.

Once the secondary background has been subtracted from the PAMELA positron
flux, the excess can be fitted very well with a simple power-law with a
spectral index of $2.3$ (Fig. 2, black dashed line). Remarkably, if we assume
the source of these “extra” positrons to be symmetric in $e^{+}-e^{-}$, we can
also fit the Fermi bump leaving a primary lepton component of the shape
$E^{-3.2}\cdot\exp(-E/1.46$ TeV) (red dashed line). The (shifted) high-energy
HESS data provides a means to fix the cut-off energy of the “extra” leptons to
about 920 GeV.

With the help of our semi-analytic transport model we were able to identify
the low-energy PAMELA data as secondary positrons from collisions of CR
protons with ISM gas using very realistic values for the simulation
parameters. An indirect implication from that is that lepton transport is
dominated by radiative losses (synchrotron and IC) above 1 GeV.

Assuming a generic primary lepton spectrum of the form
$E^{-3.2}\cdot\mathrm{exp}(-E/1.46$ TeV), the PAMELA anomaly, the Fermi bump
and the high energy HESS data can all be fitted by a single component, that is
symmetric in $e^{+}-e^{-}$ with a power-law spectral index of $2.3$ and a cut-
off energy of 920 GeV. This could either be a lepton population injected into
the ISM by nearby source with $E^{-2.3}$, so that the propagation time scale
is shorter than the radiative loss time scale, or a more distant source with a
much harder spectrum $E^{-1.3}$ which is then steepened by radiative losses. A
supernova remnant would be a good candidate in the first case, pulsars of
different ages would suite both cases.

###### Acknowledgements.

We acknowledge the support by the Deutsche Forschungsgesellschaft through the
Graduate School “Astroplasmaphysik” and through grant SP 1124-3.

Edited by: R. Vainio
Reviewed by: C. Grimani and another anonymous referee

## References

  * Abdo et al. (2009) Abdo et al. (Fermi LAT Collaboration): Measurement of the Cosmic Ray $e^{+}+e^{-}$ Spectrum from 20 GeV to 1 TeV with the Fermi Large Area Telescope, Phys. Rev. Lett., 102, 181101, doi:10.1103/PhysRevLett.102.181101, 2009.
  * Ackermann et al. (2010) Ackermann et al. (Fermi LAT Collaboration): Fermi LAT observations of cosmic-ray electrons from 7 GeV to 1 TeV, Phys. Rev. D, 82, 092004, doi:10.1103/PhysRevD.82.092004, 2010.
  * Adriani et al. (2009) Adriani et al. (PAMELA Collaboration): An anomalous positron abundance in cosmic rays with energies 1.5–100 GeV, Nature, 458, 607–609, doi:10.1038/nature07942, 2009.
  * Adriani et al. (2010) Adriani et al. (PAMELA Collaboration): A statistical procedure for the identification of positrons in the PAMELA experiment, Astropart. Phys., 34, 1–11, doi:10.1016/j.astropartphys.2010.04.007, 2010.
  * Aharonian et al. (2009) Aharonian et al. (H.E.S.S. Collaboration): Probing the ATIC peak in the cosmic-ray electron spectrum with H.E.S.S., Astron. Astrophys., 508, 561–564, doi:10.1051/0004-6361/200913323, 2009.
  * Allahverdi et al. (2009) Allahverdi, R., Dutta, B., Richardson-McDaniel, K., and Santoso, Y.: A Supersymmetric B-L Dark Matter Model and the Observed Anomalies in the Cosmic Rays, Phys. Rev. D, 79, 075005, doi:10.1103/PhysRevD.79.075005, 2009.
  * Blasi and Amato (2011) Blasi, P. and Amato, E.: Positrons from pulsar winds, Astroph. Space Sci. Proc., 623–641, doi:10.1007/978-3-642-17251-950, 2011.
  * Büsching et al. (2008) Büsching, I., de Jager, O. C., Potgieter, M. S., and Venter, C.: A Cosmic-Ray Positron Anisotropy due to Two Middle-Aged, Nearby Pulsars?, Astrophys. J., 678, L39–L42, doi:10.1086/588465, 2008.
  * Chang et al. (2008) Chang et al. (ATIC Collaboration): An excess of cosmic ray electrons at energies of 300-800GeV, Nature, 456, 362–365, doi:10.1038/nature07477, 2008.

  * Cronin et al. (1997) Cronin, J. W., Gaisser, T. K., and Swordy, S. P.: Cosmic rays at the energy frontier., Scient. American, 276, 32–37, 1997.
  * Grimani (2009) Grimani, C.: Experimental clues on pulsar energy losses and the role of circumpulsar supernova fallback disks, Classical and Quantum Gravity, 26, 235 009, doi:10.1088/0264-9381/26/23/235009, 2009.
  * Hein and Spanier (2008) Hein, T. and Spanier, F.: Analytical view of diffusive and convective cosmic ray transport in elliptical galaxies, Astron. Astrophys., 481, 1–16, doi:10.1051/0004-6361:20078645, 2008.
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., and Bugayov, V. V.: Energy spectra of gamma rays, electrons, and neutrinos produced at proton-proton interactions in the very high energy regime, Phys. Rev. D, 74, 034 018+, doi:10.1103/PhysRevD.74.034018, 2006.
  * Lerche and Schlickeiser (1988) Lerche, I. and Schlickeiser, R.: On the energy spectrum of galactic primary cosmic rays - Results from the transport equation, Astroph. Space Science, 145, 319–354, 1988.
  * McKee and Ostriker (1977) McKee, C. F. and Ostriker, J. P.: A theory of the interstellar medium - Three components regulated by supernova explosions in an inhomogeneous substrate, Astrophys. J., 218, 148–169, doi:10.1086/155667, 1977\.
  * Moskalenko and Strong (1998) Moskalenko, I. V. and Strong, A. W.: Production and Propagation of Cosmic-Ray Positrons and Electrons, Astrophys. J., 493, 694, doi:10.1086/305152, 1998.
  * Schlickeiser et al. (1987) Schlickeiser, R., Sievers, A., and Thiemann, H.: The diffuse radio emission from the Coma cluster, Astron. Astrophys., 182, 21–35, 1987.
  * Stephens (2001) Stephens, S. A.: Origin of cosmic ray electrons, Adv. Space Res., 27, 687–692, doi:10.1016/S0273-1177(01)00110-7, 2001.
  * Wang and Schlickeiser (1987) Wang, Y.-M. and Schlickeiser, R.: Smearing of a beaming pattern by an isotropic cloud - an analysis with applications to nonpulsing X-ray sources and cosmic rays, Astrophys. J., Part 1, 313, 200–217, doi:10.1086/164962, 1987.