# A Categorization of Complexity Classes for Information Retrieval and
Synthesis Using Natural Logic

Greg Coppola
coppola.ai
Research. Develop. Meme

###### Abstract

Given the emergent reasoning abilities of large language models, _information
retrieval_ is becoming more complex. Rather than just _retrieve a document_ ,
modern information retrieval systems adverstise that they can _synthesize an
answer_ based on potentially many different documents, conflicting data
sources, and using _reasoning_. But, different kinds of questions have
different answers, and different answers have different complexities. In this
paper, we introduce a novel framework for analyzing the complexity of a
question answer based on the _natural deduction calculus_ as presented in
[Prawitz, 1965]. Our framework is novel both in that no one to our knowledge
has used this logic as a basis for complexity classes, and also in that no
other existing complexity classes to these have been delineated using any
analagous methods either. We identify three decidable fragments in particular
called the _forward_ , _query_ and _planning_ fragments, and we compare this
to what would be needed to do proofs for the _complete_ first-order calculus,
for which theorem-proving is long known to be undecidable.

###### Contents

  1. 1 Introduction
  2. 2 Background
    1. 2.1 First-Order Theorems
    2. 2.2 The Church-Turing Thesis
    3. 2.3 The Natural Deduction Calculus
  3. 3 The Forward Fragment
    1. 3.1 Fragment Definition
      1. 3.1.1 Horn Clauses
      2. 3.1.2 Disjunctive Normal Form
      3. 3.1.3 Conjoined Conclusions
    2. 3.2 Datalog’s Safety Restriction
    3. 3.3 Analysis of the Fragment
  4. 4 The Query Fragment
    1. 4.1 Motivation
    2. 4.2 Existential Quantifiers as Queries
    3. 4.3 Complexity of the Full Fragment
    4. 4.4 Useful Best-Effort Fragments
  5. 5 The Planning Fragment
    1. 5.1 Example of Reasoning Under Uncertainty
    2. 5.2 Disjunctive Normal Form
    3. 5.3 Two-Player Games
    4. 5.4 Relevant Empirical Results
  6. 6 Discussion
    1. 6.1 Review of Studied Fragments
    2. 6.2 The Remaining Fragments
  7. 7 Appendix: The Natural Deduction Calculus

##  1 Introduction

##### Information Retrieval Now Requires Theorem-Proving

Given the reasoning abilities of large language models, _information
retrieval_ is becoming more complex. Rather than just retrieve a document,
modern information retrieval systems are asked to _synthesize an answer_ based
on multiple data sources. One major information company researches both how to
find a “needle in a haystack,” as well as how to do “geometry” as responses to
queries.

##### A Framework Based on Natural Deduction

In this paper, we seek a _general framework_ for understanding the _complexity
classes_ of different kinds of _answers_ that a user can ask for. Our analysis
is based on the rules for the _natural deduction calculus_ as presented in
[Prawitz, 1965]. In [Coppola, 2024], we showed how to implement a _logical
graphical model_ , in which each node in a _Bayesian Network_ is a sentence in
a logical language, and each has a probability attached. We asserted that a
_forward pass_ in our network corresponded to inference using a _specific
subset_ of the rules in [Prawitz, 1965]’s calculus. We now expand on this.

##### Analyzing the Reasoning Abilities of Transformers

A major theme recently has become the realization that transformers are
reasoning to some extent, and to investigation as to the limits of their
reasoning abilities, e.g. [Almeekam et al., 2023, Kambhampati, 2023]. We take
the perspective that, if the end result is _reasoning_ , then the optimal way
to study this is through formal logic, especially the _first-order logic_ ,
which is taken for granted to be sufficient to describe science and
mathematics from the perspective of the philosophy of science (see, e.g.
[Pelletier, 2000]). Our insight is that natural and useful classes and models
of computation can be derived from the _natural_ calculus, corresponding to
the difference in modern systems between doing a _forward pass_ through a
grapical model, _querying_ , and _planning_.

##  2 Background

###  2.1 First-Order Theorems

A theorem is a pair ${\bf d}=\left({\bf A},{\bf C}\right)$, where ${\bf A}$ is
set of assumptions and ${\bf C}$ is a set of conclusions, such that each ${\bf
A}$ and ${\bf C}$ is a subset of the possible sentences in a logical language
$\ell$ of interest. A proof of ${\bf d}$ in the calculus $\Gamma_{\ell}$ is a
sequence of deduction steps $\gamma_{\bf d}=[d_{1},...,d_{n}]$, that derive
${\bf d}$. Given a theorem ${\bf d}$ and a sequence $\gamma_{\bf d}$, we can
trivially verify whether $\gamma_{\bf d}$ constitutes a valid proof of ${\bf
d}$. In such a case, we say that $\left({\bf A},{\bf C}\right)$ is provable in
$\Gamma_{\ell}$. We say that a set ${\bf C}$ is true given ${\bf A}$, if every
model satisfying ${\bf A}$, that honors the meaning of the logical connectives
($\land$, $\lor$, $\rightarrow$, $\forall$, $\exists$ and $\bot$), must also
satisfy ${\bf C}$. The first-order calculus is so useful because it is
consistent, meaning that everything provable is true and complete, meaning
that everything true is provable [Gödel, 1931, Gentzen, 1934].

###  2.2 The Church-Turing Thesis

While the task of verifying whether $\gamma_{\bf d}$ is a valid proof of ${\bf
d}$ is trivial, the task of deciding whether ${\bf d}$ has a proof is related
to the halting problem, and is undecidable in general [Turing, 1937, Church,
1936]. That is, there is no universal program that can take an arbitrary
theorem $\left({\bf A},{\bf C}\right)$ and say whether it has a proof. If we
remove the quantifiers $\forall$ and $\exists$, we are left with the
propositional calculus. Proving a theorem in this calculus corresponds to
deciding boolean satisfiability, and this is decidable, but is NP-hard in
general [Cook, 1971], which is to say $\Omega(2^{N})$ where $N$ is the number
of boolean variables. The Church-Turing thesis is that anything that is
_computable_ is computable by a machine isomorphic to a Turing machine. Since
a Turing machine’s program is isomorphic to a theorem, we propose to use the
theory of _theorem-proving_ to identify natural computing complexity classes.

###  2.3 The Natural Deduction Calculus

[Prawitz, 1965] develops the natural deduction calculus of [Gentzen, 1934].
This is a _complete_ and _consistent_ calculus for first-order logic, in which
there are _twelve_ rules, one _Introduction_ and one _Elimination_ rule for
each of the six logical connectives $\land$, $\lor$, $\rightarrow$, $\forall$,
$\exists$ and $\bot$. In [Coppola, 2024], we show how to implement a subset of
these inference rules, which we will call the _Forward Fragment_ in Section 3.
It was always wondered why the _natural deduction_ calculus was considered
_natural_ , but intuitively it did seem that there was something _natural_
about it. Our work in [Coppola, 2024] suggests a reason for this: the logic is
_natural_ because it corresponds to how we would make a graphical model. That
is, inferences in a graphical model structure correspond to certain of the
rules in [Prawitz, 1965]’s calculus. The rules of the natural deduction
calculus are reviewed in detail in the Appendeix in Section 7.

##  3 The Forward Fragment

In this section, we investigate the _Forward Fragment_ of the deduction
calculus. This is the fragment implemented in [Coppola, 2024], and corresponds
to a single forwards pass through a graphical network, in which the nodes are
indexed by sentences in the first-order language.

###  3.1 Fragment Definition

####  3.1.1 Horn Clauses

Consider the first-order language with a set of predicate symbols, function
symbols, and constants. Suppose that that we have $n$ different functions
$P_{i}$ of at most $k$ variables, and suppose that $C$ is a function of $k$
open variables. Then, the _Forward_ fragment, implemented in [Coppola, 2024],
involves a system of _quantified Horn clauses_ of the form:

$\forall x_{1},\ldots,\forall
x_{k},\Bigg{[}\bigwedge_{i=1}^{n}P_{i}(x_{1},\ldots,x_{k})\Bigg{]}\rightarrow
C(x_{1},\ldots,x_{k})$ (1)

We emphasize that every variable mentioned $x_{1},...,x_{k}$ must appear in
$C$. This is called the _safety_ restriction of _Datalog_ , which is very
crucial to the distinction in complexity classes, and is further discussed in
Section 3.2. An example of an inference that uses the form is the following:

$\forall x_{1},x_{2},\Bigg{[}like(x_{1},x_{2})\wedge
like(x_{2},x1)\Bigg{]}\rightarrow friends(x_{1},x_{2})$ (2)

This says that if any $x_{1}$ and $x_{2}$ both like _each other_ , then they
are _friends_.

####  3.1.2 Disjunctive Normal Form

The presence of multiple ways to derive the same predicate amounts to a
_disjunction_ over the possible ways to reach that conclusion, and thus $m$
statements of the form 1 can be written as a statement with a disjoined
premise as:

$\forall x_{1},\ldots,\forall
x_{k},\Bigg{[}\bigvee_{i=1}^{m}\bigwedge_{i=1}^{n}P_{i,j}(x_{1},\ldots,x_{k})\Bigg{]}\rightarrow
C(x_{1},\ldots,x_{k})$ (3)

Thus, any sentence in _disjunctive normal form_ , apart from the restriction
on quantification, can be used as a premise. Since any statement in the
propositional calculus can be written in disjunctive normal form [Andrews,
1986], any statement obeying the restriction on quantification can be a
premise to a deduction rule, so this is a very powerful framework.

####  3.1.3 Conjoined Conclusions

Finally, if we have a statements of the form 3, note that in some cases, the
premise can be the same, and so a database of different statements like 3 can
be re-written as a single statement like:

$\forall x_{1},\ldots,\forall
x_{k},\Bigg{[}\bigvee_{i=1}^{n}\bigwedge_{i=1}^{n}P_{i}(x_{1},\ldots,x_{k})\Bigg{]}\rightarrow\Bigg{[}\bigwedge_{i=1}^{n}C(x_{1},\ldots,x_{k})\Bigg{]}$
(4)

That is, a set of statements of the form 1 and 4 are equivalent in the proofs
that they allow.

###  3.2 Datalog’s Safety Restriction

In this _Forward_ fragment, we are incorporating an important restriction on
the forms of clauses, which is that we can only mention in the
$P_{i}(x_{1},...,x_{k})$ variables that have also been mentioned in
$C(x_{1},...,x_{k})$. That is, in this fragment, we cannot introduce new
universally quantified variables used in $P_{i}$ but not in $C$, and means we
cannot say that we are _truly_ implementing _$\forall$ -Elimination_, for
example, because the _$\forall$ -Elimination_ rule in [Prawitz, 1965]’s system
does _not_ include this restriction. A Datalog rule is considered to be safe
if every variable that appears in the head of the rule (the conclusion part)
also appears in a positive literal in the body of the rule (the condition
part) [Abiteboul et al., 1995].

###  3.3 Analysis of the Fragment

This fragment allows all statements that use only the following rules:
_$\land$ -Elimination_, _$\land$ -Introduction_, _$\lor$ -Introduction_,
_$\forall$ -Elimination_ and _$\rightarrow$ -Elimination_, but subject to the
_safety_ restriction on quantification as discussed above. Because this
fragment follows the safety restriction of Datalog, it amounts to _Horn
Satisfiability_ , which is a well-known efficient fragment, where the time
taken to do inference is _linear_ in the number of variables _total_ in the
theory, which is really very efficient, when compared to the other fragments,
as we will see. [Coppola, 2024] implements this fragment, and shows how we can
not only do logical inference, but even assign probabilities if this fragment
is stuck to. That is, we give probabilities to a similar fragment as that for
which [Pereira and Shieber, 1987] investigated theorem-provability. That is,
we can distinguish between theories that are more or less _likely_. And, we
can determine a _generative_ probability for the data set, so that it can be
compressed [Sutskever, 2023].

##  4 The Query Fragment

###  4.1 Motivation

Suppose that instead of 1, suppose we have an integer $K>k$, and:

$\forall x_{1},\ldots,\forall
x_{K},\Bigg{[}\bigwedge_{i=1}^{n}P_{i}(x_{1},\ldots,x_{k},x_{k+1},...,x_{K})\Bigg{]}\rightarrow
C(x_{1},\ldots,x_{k})$ (5)

This statement is _unsafe_ in the terms of Datalog (see Section 3.2), because
the variables $x_{k+1},...,x_{K}$ are _not mentioned_ in the conclusion, and
thus _unbound_. Such a statement can be useful, for example, in saying that if
$x_{1}$ and $x_{2}$ both want the same thing $x_{3}$, then $x_{1}$ and $x_{2}$
are _in competition_ :

$\forall x_{1},x_{2},x_{3},\Bigg{[}want(x_{1},x_{3})\wedge
want(x_{2},x_{3})\Bigg{]}\rightarrow compete(x_{1},x_{2})$ (6)

Thus, the restriction that we have put in Section 3.2 is something that we
want to relax.

###  4.2 Existential Quantifiers as Queries

We can re-write 5 as:

$\forall x_{1},\ldots,\forall x_{k},\Bigg{[}\exists
x_{k+1},...,x_{K},\bigwedge_{i=1}^{n}P_{i}(x_{1},\ldots,x_{k},x_{k+1},...,x_{K})\Bigg{]}\rightarrow
C(x_{1},\ldots,x_{k})$ (7)

Each existentially quantified variable $\exists x_{k+i}$ corresponds to a
_query_ as to whether any element in the domain can be found that satisfies
the premise.

###  4.3 Complexity of the Full Fragment

If our goal is to do _full deterministic_ theorem-proving, the worst-case
bound is very bad, for reasons that we will now investigate.

##### Quantifiers in Same Sentence

Suppose that the $x_{i}$ are universally quantified and the $y_{j}$ are
existentially quantified, but we will still write $\exists$ for clarity. Now,
consider an implication with a single existential:

$\left[\exists y_{1},P(x_{1},y_{1})\right]\rightarrow C(x_{1})$ (8)

If the domain that $\exists y_{1}$ ranges over is finite, which we assume that
it is, then we can _implement_ $\exists$ using _$\vee$ -Introduction_, with a
disjunction over $D$ elements in the domain:

$\left[\bigvee_{i=1}^{D}(x_{1},c_{i})\right]\rightarrow\exists
y_{1}P(x_{1},y_{1})$ (9)

This is computable, but quickly becomes a problem from an efficiency
perspective. For example, to implement a premise with two exisentially
quantified variables in the premise like:

$\left[\exists y_{1},y_{2},P(x_{1},y_{1},y_{2})\right]\rightarrow C(x_{1})$
(10)

This would require $\Theta(D^{2})$ simple rules to implement. In general, if
we begin with a graph of size $G$, and a domain of size $D$, then to implement
$N$ $\exists$ would require work $\Theta(GD^{N})$.

###  4.4 Useful Best-Effort Fragments

We must stress that the exponential blow-up in computational cost described
above is only in the _worst-case_ assuming _full theorem-proving_. We discuss
two efficient but useful fragments of the query calculus that we believe can
be helpful in practice.

##### Shallow Queries

Suppose instead, we can say that, when eliminating a $\exists$, we are
interested only in facts that are _already in_ an existing _statically_ stored
database, without requiring any _dynamic_ theorem-proving. In other words, we
can make a single database query to see if there is any $x_{3}$ such that
$want(x_{1},x_{3})$ and $want(x_{2},x_{3})$. This can be done in $O(1)$ time
relative to the complexity of the graph, possibly even using a database index.

##### Probabilistic Ranking

In traditional, non-probabilistic theorem-proving, if we want to prove
$\exists xP(x)$, we have no a priori way to _order_ the $x$ in terms of which
would be _most likely_ make $P(x)$ true. Such a ranking could cut down on the
search time. We could also use a heuristic like _try only the most probable
candidate_ , which would eliminate the exponential blow up in the search
space, by only allowing $1$ try per $\exists$.

##### A* Search

Otherwise, if we have a _ranking_ of different candidates, and we really want
to search exhaustively, or as exhaustively as possible, we speculate we can
view this as analogous to searching through a _maze_ and use _A* search_ as in
[Lehnert et al., 2024].

##  5 The Planning Fragment

The natural rule of _$\lor$ -Elimination_ we also call _reasoning by cases_ ,
and now we will investigate why that is, and see the relation to _planning_
under _uncertain conditions_.

###  5.1 Example of Reasoning Under Uncertainty

##### Set-Up

Let us consider the example of a party that is planning an _enjoyable
excursion_ , where the goal at the end is to be _satisfied_ with their
experience. This particular party ${\bf p}$ has their own individual
preferences. And, they are considering a trip to a particular beach town ${\bf
t}$. Suppose the party will be _satisfied_ with their excursion if the
excitement level $l$ with excursion $e$ is at least a $7$ out of $10$:

$\forall e,\exists l,\Bigg{[}excited({\bf p},e,l)\wedge l\geq
7\Bigg{]}\rightarrow satisfied({\bf p},e)$ (11)

Here, we assume that we can use $\exists$. If the party visits the beach town
${\bf t}$ they can go to the beach:

$\Bigg{[}visit({\bf p},e,{\bf t})\Bigg{]}\rightarrow visit({\bf
p},e,beach({\bf t}))$ (12)

If they visit the beach, and it is sunny, the happiness level will be a $10$.

$\Bigg{[}visit({\bf p},e,beach({\bf t}))\wedge sunny(e)\Bigg{]}\rightarrow
excited({\bf p},e,10)$ (13)

But, if at the beach, it is not sunny, the party will be _very_ unhappy:

$\Bigg{[}visit({\bf p},e,beach({\bf t}))\wedge\neg sunny(e)\Bigg{]}\rightarrow
excited({\bf p},e,1)$ (14)

Now, there is a favorite _restaurant_ of the party’s in t, which they enjoy
going to when it rains:

$\Bigg{[}visit({\bf p},e,restaurant({\bf t}))\wedge\neg
sunny(e)\Bigg{]}\rightarrow excited({\bf p},e,7)$ (15)

Going to this restaurant is not as fun as the beach, but it is still fun on a
rainy day. However, going to the restaurant if it is _not_ sunny does not make
the party happy, as they will wonder why they are not at the beach:

$\Bigg{[}visit({\bf p},e,restaurant({\bf t}))\wedge
sunny(e)\Bigg{]}\rightarrow excited({\bf p},e,3)$ (16)

##### Observation

We are now in a situation where the party actually _can_ guarantee that they
will be satisfied with their excursion, becuase they can go to the town ${\bf
t}$, and if it is sunny, they can go to the beach, and if it rains, they can
go to their restaurant, and they will be happy either way. The fact that they
will be happy _either way_ requires the ability to reason by _cases_ , and
evaluate their position _in either case_.

###  5.2 Disjunctive Normal Form

In order to handle the motivating case, and others as well, it suffices to
allow, instead of 3, that a conclusion can now itself contain a disjunction as
in:

$\forall x_{1},\ldots,\forall
x_{k},\Bigg{[}\bigvee_{i=1}^{n}\bigwedge_{i=1}^{n}P_{i}(x_{1},\ldots,x_{k})\Bigg{]}\rightarrow\bigvee_{i=1}^{n}C(x_{1},\ldots,x_{k})$
(17)

Of course, it is less information to learn that $A\lor B$ than to learn $A$,
because the latter allows us to conclude $A$, and the former does not, and
this is why we must _resaon by cases_.

###  5.3 Two-Player Games

In general, interaction between the agent and the environment can be viewed as
a two-player game, where the agent plays themselves, and _reality_ plays as
the opponent. In an _antagoistic_ two-player game, we assume that the
_opponent_ always does their best to make the life of the _protagonist_ as bad
as possible. In the case of the _protagonist_ against _reality_ , _reality_
might not always try to make things as bad as possible for the _protagonist_ ,
but rather reality will react according to a _probability distribution_. In
either case, this is similar to a game of _chess_ or _go_ [Silver et al.,
2016].

###  5.4 Relevant Empirical Results

Because of the $\Omega(2^{N})$ explosion in complexity of solving a boolean
satisfiability problem [Cook, 1971], we expect that a _large language model_
of finite size will not be able to handle problems corresponding to boolean
satisfiability, if the input grows large enough. This has been investigated
and shown in [Almeekam et al., 2023, Kambhampati, 2023].

##  6 Discussion

###  6.1 Review of Studied Fragments

In the _forward_ fragment of Section 3, we investigate a system that contains
the rules:

  * •

forward fragment

    * –

_$\land$ -Introduction_

    * –

_$\land$ -Elimination_

    * –

_$\lor$ -Introduction_

    * –

_$\rightarrow$ -Elimination_

    * –

_$\forall$ -Elimination_ (limited by the _safety_ requirement on
quantification)

In the _query_ fragment of Section 4, we weaken the restriction on
quantification and add _$\exists$ -Introduction_:

  * •

query fragment

    * –

all those from the _direct fragment_ , plus

    * –

_$\forall$ -Elimination_ (without limitation)

    * –

_$\exists$ -Introduction_

Adding in _reasoning by cases_ creates the _planning fragment_ of Section 5:

  * •

planning fragment

    * –

all those from the _query fragment_ , plus

    * –

_$\vee$ -Elimination_

###  6.2 The Remaining Fragments

The remaining fragments are the ones that [Prawitz, 1965] called _improper_ :

  * •

improper rules

    * –

_$\lor$ -Elimination_

    * –

_$\rightarrow$ -Introduction_

    * –

_$\forall$ -Introduction_

    * –

_$\exists$ -Elimination_

    * –

_$\bot$ -Introduction_

In [Coppola, 2024] we implemented the _direct_ fragment, and showed how to
assign probabilities to conclusions. We said that one feature of the research
is that it shows how to _unify logical and probabilistic reasoning_ , by
providing a _unified model_ of the two. This is because, while we can show
what happens in a _forward pass_ of a _direct_ fragment, we can also show how
this fragment that is implemented related to a _complete_ and _consistent_
calculus for the first-order logic. However, it is important to understand the
reason that we can _not_ implement an entire logic as a _forward pass_ in a
graphical network. To do proofs in the unbound fragments, one must implement
_partial_ strategies for these rules that [Prawitz, 1965] called _improper_.
Such a list must always be incomplete, because full theorem proving in the
first-order calculus is undecidable in general [Church, 1936, Turing, 1937].

##  7 Appendix: The Natural Deduction Calculus

##### Introduction

This is a review of [Prawitz, 1965]’s formulation of the _natural deduction
calculus_. [Prawitz, 1965] distinguishes between deduction rules that are
_proper_ and _improper_. We instead call them _simple_ and _complex_ , to
avoid the value judgment. We will point out for each rule which category
[Prawitz, 1965] gave it. The _complex_ ones are harder to implement.

##### The Logical Rules

This is a _complete_ and _consistent_ calculus for first-order logic, in which
there are _twelve_ rules, one _Introduction_ and one _Elimination_ rule for
each of the six logical connectives $\land$, $\lor$, $\rightarrow$, $\forall$,
$\exists$ and $\bot$.

#####  $\land$-Introduction

The rule of $\land$-Introduction says that if we have proved $A$ and $B$, we
can conclude $A\land B$:

$\left\\{\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{A,B\right\\}\right)\right\\}\rightarrow\left({\bf
A}_{0},{\bf C}_{0}\cup\left\\{A\wedge B\right\\}\right)$ (18)

This inference rule is _simple_.

#####  $\land$-Elimination

The rule of $\land$-Elimination says that if we have proved $A\land B$, we can
conclude both $A$ and $B$:

$\left\\{\left({\bf A}_{0},{\bf C}_{0}\cup\left\\{A\wedge
B\right\\}\right)\right\\}\rightarrow\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{A,B\right\\}\right)$ (19)

This inference rule is simple. Note that $\land$ is the only logical symbol
for which both rules are simple.

#####  $\lor$-Introduction

The rule of $\lor$-Introduction says that if we have either $A$ or $B$, we can
conclude $A\lor B$:

$\left\\{\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{A\right\\}\right)\right\\}\rightarrow\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{A\right\\}\cup\left\\{A\lor B\right\\}\right)$ (20)

This inference rule is simple. This corresponds to a disjunction gate in the
QBBN, and can be learned to make a statistical model [Coppola, 2024].

#####  $\lor$-Elimination

$\lor$-Elimination says that if we have proved $A\lor B$, and we have proven
that from $A$ we can conclude $C$, and we have proven that from $A$ we can
conclude $C$, then we can conclude $C$.

$\left\\{\begin{aligned} &\left({\bf A}_{0}\cup\left\\{\right\\},{\bf
C}_{0}\cup\left\\{A\lor B\right\\}\right)\\\ &\left({\bf
A}_{0}\cup\left\\{A\right\\},{\bf C}_{0}\cup\left\\{C\right\\}\right)\\\
&\left({\bf A}_{0}\cup\left\\{B\right\\},{\bf
C}_{0}\cup\left\\{C\right\\}\right)\end{aligned}\right\\}\rightarrow\left({\bf
A}_{0},{\bf C}_{0}\cup\left\\{C\right\\}\right)$ (21)

This is complex, and actually involves three sets of assumptions, ${\bf
A}_{0}$, ${\bf A}_{0}\cup\left\\{A\right\\}$ and ${\bf
A}_{0}\cup\left\\{B\right\\}$. Although this deduction is complex, we can show
in a sense that this rule is less complex than the others, in that we can
always decide proofs in a logic that only uses the simple rules and this rule,
but it becomes NP-hard, see Section 5.

#####  $\rightarrow$-Introduction

The rule of $\rightarrow$-Introduction says that if we know we can prove $B$
from $A$, then we can prove $A\rightarrow B$: $\rightarrow$-Introduction is
complex and says:

$\left\\{\left({\bf A}_{0}\cup\left\\{A\right\\},{\bf
C}_{0}\cup\left\\{B\right\\}\right)\right\\}\rightarrow\left({\bf
A}_{0}\cup\left\\{\right\\},{\bf C}_{0}\cup\left\\{A\rightarrow
B\right\\}\right)$ (22)

This is complex because it involves a change of assumptions from ${\bf
A}_{0}\cup\left\\{A\right\\}$ to ${\bf A}_{0}$. This is a rare rule to use in
practice.

#####  $\rightarrow$-Elimination

The rule of $\rightarrow$-Elimination is simple. It corresponds to a single
factor forward inference in the QBBN. Starting with the simple
$\rightarrow$-Elimination, we have:

$\left\\{\left({\bf A}_{0},{\bf C}_{0}\cup\left\\{A,A\rightarrow
B\right\\}\right)\right\\}\rightarrow\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{A,A\rightarrow B\right\\}\cup\left\\{B\right\\}\right)$ (23)

This is simple. This deduction is the bedrock of logic, i.e., modus ponens.

#####  $\forall$-Introduction

The rule of $\forall$-Introduction asserts that if we can prove $A$ without
specifically referencing any particular instance of $x$, then $A$ can be
generalized to $\forall xA$:

$\left\\{\left({\bf A}_{0},{\bf
C}_{0}\cup\\{A\\}\right)\right\\}\rightarrow\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{\forall x.A\right\\}\right)$ (24)

This rule is considered complex due to the introduction of a universal
quantifier, expanding the scope of $A$ to all instances of $x$. The
assumptions remain unchanged.

#####  $\forall$-Elimination

The rule of $\forall$-Elimination allows for the instantiation of a
universally quantified statement to a specific instance. Given $\forall xA$,
we can deduce $A[t/x]$ for any term $t$:

$\left\\{\left({\bf A}_{0},{\bf C}_{0}\cup\left\\{\forall
x.A\right\\}\right)\right\\}\rightarrow\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{A[t/x]\right\\}\right)$ (25)

This rule is simple.

#####  $\exists$-Introduction

The rule of $\exists$-Introduction posits that if $A[t/x]$ is proven for some
term $t$, then $\exists xA$ is also proven:

$\left\\{\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{A[t/x]\right\\}\right)\right\\}\rightarrow\left({\bf
A}_{0},{\bf C}_{0}\cup\left\\{\exists x.A\right\\}\right)$ (26)

This rule is simple.

#####  $\exists$-Elimination

The rule of $\exists$-Elimination involves deducing a conclusion $B$ from an
existential premise $\exists x.A$, requiring a temporary assumption $A[c/x]$
for a new constant $c$:

$\left\\{\left({\bf A}_{0}\cup\left\\{\exists x.A\right\\},{\bf
C}_{0}\right),\left({\bf A}_{0}\cup\left\\{A[c/x]\right\\},{\bf
C}_{0}\right)\right\\}\rightarrow\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{B\right\\}\right)$ (27)

This rule is complex because it involves the introduction of a new assumption
for a hypothetical instance $c$ that satisfies $A$, and the set of assumptions
changes during the process.

#####  $\bot$-Introduction

The rule of $\bot$-Introduction allows for the introduction of a contradiction
($\bot$) when both a statement $A$ and its negation $\neg A$ are proven:

$\left\\{\left({\bf A}_{0},{\bf C}_{0}\cup\left\\{A,\neg
A\right\\}\right)\right\\}\rightarrow\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{\bot\right\\}\right)$ (28)

This rule is complex as it derives a contradiction, indicating a fundamental
inconsistency within the assumptions.

#####  $\bot$-Elimination

Applying the rule of $\bot$-Elimination, any statement $B$ can be concluded
from a contradiction ($\bot$):

$\left\\{\left({\bf A}_{0},{\bf
C}_{0}\cup\left\\{\bot\right\\}\right)\right\\}\rightarrow\left({\bf
A}_{0},{\bf C}_{0}\cup\left\\{B\right\\}\right)$ (29)

This rule is simple.

## References

  * [Abiteboul et al., 1995] Abiteboul, S., Hull, R., and Vianu, V. (1995).  Foundations of Databases.  Addison-Wesley.
  * [Almeekam et al., 2023] Almeekam, K., Marquez, M., Sreedharan, S., and Kambhampati, S. (2023).  On the planning abilities of large language models - a critical investigation.  In Thirty-seventh Conference on Neural Information Processing Systems.  URL https://openreview.net/forum?id=X6dEqXIsEW.
  * [Andrews, 1986] Andrews, P. B. (1986).  An Introduction to Mathematical Logic and Type Theory: To Truth Through Proof.  Academic Press, Inc., Orlando, San Diego, New York, Austin, Boston, London, Sydney, Tokyo, Toronto.
  * [Church, 1936] Church, A. (1936).  An unsolvable problem of elementary number theory.  American Journal of Mathematics, 58(2):345–363.
  * [Cook, 1971] Cook, S. A. (1971).  The complexity of theorem-proving procedures.  Proceedings of the third annual ACM symposium on Theory of computing.
  * [Coppola, 2024] Coppola, G. (2024).  The quantified boolean bayesian network: Theory and experiments with a logical graphical model.
  * [Gentzen, 1934] Gentzen, G. (1934).  Untersuchungen über das logische schließen.  Mathematische Zeitschrift, 39:176–210, 405–431.
  * [Gödel, 1931] Gödel, K. (1931).  Über formal unentscheidbare sätze der principia mathematica und verwandter systeme i.  Monatshefte für Mathematik, 38(1):173–198.  English title: On Formally Undecidable Propositions of Principia Mathematica and Related Systems I.
  * [Kambhampati, 2023] Kambhampati, S. (2023).  Can llms really reason and plan?  Communications of the Association for Computing Machinery (CACM) Blog.  URL https://cacm.acm.org/blogs/blog-cacm/276268-can-llms-really-reason-and-plan/fulltext.
  * [Lehnert et al., 2024] Lehnert, L., Sukhbaatar, S., Mcvay, P., Rabbat, M., and Tian, Y. (2024).  Beyond a*: Better planning with transformers via search dynamics bootstrapping.
  * [Pelletier, 2000] Pelletier, F. J. (2000).  A history of natural deduction and elementary logic textbooks.  Logical consequence: Rival approaches, 1:105–138.
  * [Pereira and Shieber, 1987] Pereira, F. C. N. and Shieber, S. M. (1987).  Prolog and Natural-Language Analysis.  CSLI Lecture Notes 10. Center for the Study of Language and Information - CSLI.
  * [Prawitz, 1965] Prawitz, D. (1965).  Natural Deduction: A Proof-Theoretical Study.  Stockholm Studies in Philosophy 3. Almqvist & Wiksell, Stockholm; Göteborg; Uppsala.  Acta Universitatis Stockholmiensis.
  * [Silver et al., 2016] Silver, D., Huang, A., Maddison, C. J., Guez, A., Sifre, L., Van Den Driessche, G., Schrittwieser, J., Antonoglou, I., Panneershelvam, V., Lanctot, M., et al. (2016).  Mastering the game of go with deep neural networks and tree search.  Nature, 529(7587):484–489.
  * [Sutskever, 2023] Sutskever, I. (2023).  An observation on generalization.  Simons Institute, YouTube.  Accessed: 2024-01-29.
  * [Turing, 1937] Turing, A. (1937).  On computable numbers, with an application to the entscheidungsproblem.  Proceedings of the London Mathematical Society, s2-42(1):230–265.


# Numerical Modeling of Time Dependent Diffusive Shock Acceleration

S. Aerdker,11footnotetext: Corresponding author. L. Merten  J. Becker Tjus  D.
Walter  F. Effenberger  H. Fichtner

###### Abstract

Motivated by cosmic ray (CR) re-acceleration at a potential Galactic Wind
Termination Shock (GWTS), we present a numerical model for time-dependent
Diffusive Shock Acceleration (DSA). We use the stochastic differential
equation solver (DiffusionSDE) of the cosmic ray propagation framework
CRPropa3.2 with two modifications: An importance sampling module is introduced
to improve statistics at high energies by keeping the simulation time short.
An adaptive time step is implemented in the DiffusionSDE module. This ensures
to efficiently meet constraints on the time and diffusion step, which is
crucial to obtain the correct shock spectra. The time evolution of the
spectrum at a one-dimensional planar shock is verified against the solution
obtained by the grid-based solver VLUGR3 for both energy-independent and
energy-dependent diffusion. We show that the injection of pre-accelerated
particles can lead to a broken power law spectrum in momentum if the incoming
spectrum of CRs is harder than the reaccelerated spectrum. If the injected
spectrum is steeper, the shock spectrum dominates at all energies. We finally
apply the developed model to the GWTS by considering a spherically symmetric
shock, a spiral Galactic magnetic field, and anisotropic diffusion. The time-
dependent spectrum at the shock is modeled as a basis for further studies.

##  1 Introduction

The cosmic-ray energy spectrum and composition is by now well-studied from GeV
up to ZeV energies. The origin of these highly energetic, charged particles is
still largely unknown. Particularly unresolved is the question of the energy
range between $10^{15}\,\mathrm{eV}$ and $3\times 10^{18}\,\mathrm{eV}$. This
energy range is defined by two kinks in the spectrum, the "knee" and the
"ankle". Up to $\sim 10^{15}\,\mathrm{eV}$ it is believed that CRs are
accelerated in the Galaxy, most likely at Supernova Remnants (SNRs) (see e.g.
[1]). Above the ankle ($\sim 10^{18}\,\mathrm{eV}$) [1] CRs have gyroradii
larger than the Galaxy and are clearly of extra-galactic origin. The breaks in
the spectrum may indicate changes in the contributions to the spectrum and a
transition from Galactic to extra-galactic origin of CRs. The break at
$\approx 200-300\,\mathrm{GV}$ especially at lower energies ([1] and
references therein) are also attributed to a change in the transport
properties.

In general, the observed power-law spectrum indicates that CRs are accelerated
by stochastic processes. There are two possible scenarios of how to accelerate
particles stochastically, both developed early-on by Fermi. Fermi Second Order
Acceleration [2] was suggested first as acceleration on moving magnetized
clouds in the Galaxy. Fermi First Order Acceleration [3, 4, 5] is more
efficient: CRs are scattered on both sides of a shock and may cross the shock
front repeatedly. Each time they cross the shock, they are accelerated
depending on their energy, $\Delta E\propto E$. The energy gain can be derived
by Lorentz boosts between the upstream and downstream reference frames. The
slope of the stationary shock spectrum is purely dependent on the shock’s
compression ratio.

There are some limitations to DSA, one of which is that particles must be able
to reach the shock from downstream to upstream again, see, e.g. [6]. In order
to cross the shock, CRs must have been accelerated to energies that give them
a gyroradius larger than the shock width, which is also known as the injection
problem. Also, CRs have to be contained in the acceleration region to
experience ample acceleration. This implies a diffusive particle transport
with diffusion low enough that CRs are not escaping too quickly. An upper
bound of the maximum energy produced by a source can be estimated by the
Hillas criterion and depends on the particles rigidity and the size of the
acceleration region [7]. There are other limiting factors like the lifetime of
the accelerator or loss processes like synchrotron radiation.

In the transition region from Galactic to extragalactic CRs the Galactic Wind
Termination Shock (GWTS) might contribute to the sources of CRs via re-
acceleration [8, 9]. A possible Galactic wind can be driven by, e.g.,
radiation or the CRs themselves (e.g. [10, 11, 12, 13, 14]). It may form a
shock when the supersonic wind slows down due to interaction with the
Intergalactic Medium (IGM).

CRs with their origin in the Galactic disk are advected outwards and are
accelerated to higher energies at the GWTS. Either they leave the Galaxy or a
fraction of the high-energy CRs may be able to propagate back to the Galaxy
against the Galactic wind. The idea of particle acceleration at the GWTS was
already discussed e.g. by [15], [8], [9], and most recently by [16]. With the
assumption that about two solar masses per year are advected outwards, the
GWTS cannot be supported by the wind longer than $100\,\mathrm{Myr}$ [17]. The
finite lifetime of the shock may also impact the shock spectrum and its
contribution to the observed cosmic-ray spectrum on Earth.

Not only the origin of CR but also their transport properties change in the
transition region. CRs can be treated as diffusion below a threshold energy.
Particles undergo a random walk due to scattering on magnetic field
turbulence. CRs with energies above $\sim 10^{18}\,\mathrm{eV}$ engage less in
interaction with the turbulent Galactic magnetic field and thus, cannot be
described properly by the diffusion approximation on these scales [18]. Where
exactly the transitions between diffusive and ballistic transport as well as
Galactic and extra-galactic CRs lie is still unclear.

Methods for the simulations of CR transport have to reflect that transition:
For high-energy extra-galactic CRs the equation of motion is integrated. In
this ballistic regime, individual particle trajectories in arbitrary magnetic
field configurations can be simulated precisely. However, for Galactic CRs
with energies below $10^{16}\,\mathrm{eV}$ this approach becomes
computationally costly as most of the computation is spend on resolving the
particles’ gyration. Therefore, for Galactic CRs the transport equation is
commonly used to describe diffusive transport in space and momentum.

The simplified, assuming isotropic (in momentum space) CR distributions,
Fokker-Planck equation, often called Parker transport equation (e.g. [19])

$\frac{\partial n}{\partial t}+\vec{u}\cdot\nabla
n=\nabla\cdot\left(\hat{\kappa}\nabla
n\right)+\frac{1}{p^{2}}\frac{\partial}{\partial p}\left(p^{2}D\frac{\partial
n}{\partial p}\right)+\frac{1}{3}\left(\nabla\cdot\vec{u}\right)\frac{\partial
n}{\partial\ln p}+S(\vec{x},p,t),$ (1.1)

describes the time evolution of the cosmic-ray density $n$ in space $\vec{x}$
and momentum $p=|\vec{p}|$. Individual terms describe advection of CR with a
background flow $\vec{u}$, spatial diffusion described by the diffusion tensor
$\hat{\kappa}$, momentum diffusion described by the coefficient $D$, adiabatic
energy changes and cosmic-ray sources $S$.

In the ensemble-averaged picture, acceleration at shocks results from
adiabatic heating. The Lorentz boosts between the particle frame and the frame
in which the background flow is at rest are implicitly done in the derivation
of the transport equation for accelerating background flows (see [20]).

The spatial diffusion tensor is often defined in a local coordinate frame,
describing diffusion parallel and perpendicular to the local magnetic field.
The momentum diffusion coefficient describes acceleration of CR in the
presence of magnetic field turbulence (Second-order Fermi acceleration). The
knowledge of both, spatial and momentum diffusion tensors, is essential for
the complete description of CR transport.

There is not yet a closed description of the diffusion tensor. Depending on
the assumed magnetic field, coherent background and turbulent component, and
the CRs properties analytical or numerical approximations might be available
[21, 22]. On the other, hand it can also be used as a free parameter and
fitted to match, e.g., the observed primary and secondary ratios in a model of
Galactic cosmic-ray transport. Since all of these approaches come with their
own caveats, we will use normalized diffusion coefficients
$\kappa_{\parallel}=1$ throughout this work.

Considering the GWTS mentioned at the beginning, re-accelerated CRs are
expected to be in the transition from diffusive to ballistic transport. Close
to the shock due to strong turbulence induced by the CRs themselves the
diffusive approach is valid for higher energies than in the interstellar or
IGM [4, 23]. Particles escaping the acceleration region, however, cannot
necessarily be described by the diffusion approximation. In order to assess
the contribution of the GWTS to the cosmic-ray spectrum they must be
propagated back to the Galaxy using either the ballistic or diffusive approach
depending on their energy. Thus, the spectrum observed on Earth is further
modulated by the energy-dependence and spatial variability of the diffusive
transport.

The CR propagation framework CRPropa3.2 [24] offers the possibility to
simulate acceleration and propagation in both regimes. Also, the structure of
CRPropa3.2, which is described in more depth in Section 2, makes it relatively
easy to define arbitrary magnetic field and shock configurations in up to
three dimensions. Together with the GWTS scenario, this motivates an in-depth
study of DSA using the diffusion approach of CRPropa3.2 [25]. In test
scenarios we explore different questions related to the re-acceleration of CR
at the GWTS: How does the spectrum change when a finite shock lifetime is
considered? What are the effects on the spectrum when CRs are already pre-
accelerated to a power-law? To what extent does a finite acceleration region
affect the maximal energy that can be reached?

First, we give a short overview of CRPropa3.2 and its diffusion approach which
is based on Stochastic Differential Equations (SDEs) in Section 2. A new
module is presented that was specifically implemented to enhance statistics
when simulating DSA with CRPropa3.2.

In Section 3 we show that DSA simulations with CRPropa3.2 result in the
expected spectral slope for acceleration at one-dimensional planar shocks. We
validate the obtained stationary spectra with predictions from theory and
other ensemble-averaged approaches to DSA. The time evolution of the spectrum
at the shock is compared to simulations that integrate the transport equation
using the finite difference code VLUGR3 [26, 27, 28]. We clarify constraints
for simulating DSA based on SDEs. Based on those findings, we consider a
finite shock lifetime, energy-dependent diffusion and the injection of pre-
accelerated spectra. The effects of each modification on the spectrum at the
shock or the downstream boundary are analysed separately.

In Section 4 we take anisotropic diffusion into account and consider oblique
shocks. Finally, we present the time-dependent spectrum at a simple model for
the GWTS: A spherical symmetric shock and spiral magnetic field.

##  2 Methods

In this Section we briefly introduce SDEs and how they relate to the transport
equation (1.1). We explain how they are used in CRPropa3.2 to solve the
transport equation. For the simulation of DSA a new module,
CandidateSplitting, was implemented to enhance statistics at high energies.

###  2.1 Stochastic Differential Equations

One way to simulate DSA is the solution of the transport equation (1.1), a
partial differential equation in reduced phase-space. The solution directly
results in the distribution function of CRs. We followed this approach using
the well-established VLUGR3 [26, 27] code to cross check our results. Another
way is to make use of the connection between Fokker-Planck type equations and
SDEs. In general, Fokker-Planck type equations describe the transition
probability from one state to a set of other states. SDEs describe dynamical
systems which are subject to noise.

In the following we discuss briefly the set of SDEs implemented in CRPropa3.2.
For an in-depth explanation the reader is referred to [25] or e.g. [29].

In the absence of second-order Fermi acceleration processes, i.e. $D=0$, and
spatially constant diffusion the system of SDEs equivalent to the transport
equation (1.1) is given by

$\displaystyle\mathrm{d}\vec{x}$
$\displaystyle=\vec{u}(x)\,\mathrm{d}t+\hat{B}\,\mathrm{d}\vec{\omega}_{t},$
(2.1) $\displaystyle\mathrm{d}p$
$\displaystyle=-\frac{p}{3}\nabla\cdot\vec{u}\,\mathrm{d}t,$ (2.2)

where $\mathrm{d}\vec{\omega}_{t}=\sqrt{dt}\,\vec{\eta}$ is a Wiener process
with the Gaussian distributed random numbers $\eta_{\mathrm{i}}$. The tensor
$\hat{B}$ corresponds to spatial diffusion. The diffusion tensor becomes
diagonal when drift terms described by off-diagonal elements are neglected.
With this assumption $\hat{B}$ is given by
$B_{\mathrm{ij}}=\delta_{\mathrm{ij}}\sqrt{2\kappa_{\mathrm{ij}}}$. Equation
(2.1) describes the spatial displacement due to advection and stochastic
fluctuations. Equation (2.2) describes the adiabatic energy change. Without
diffusion in momentum space it is not an SDE but an ordinary differential
equation.

SDEs can also be written in integral form by making use of the Itô integral
[30]. In this form they can be approximated numerically (see e.g. [31, 32]).
The SDE approach has the advantage that it is easy to implement and to extend
to arbitrary geometries. We describe the first-order Euler-Maruyama scheme
implemented in CRPropa3.2 for numerical solution of the SDE (Eq. (2.1)) in the
following Section.

###  2.2 Simulation with CRPropa3.2

CRPropa3.2 has a modular structure with the Candidate as central element. When
the EoM is used, the candidate simply corresponds to the propagated particle.
When the diffusive transport is used, the candidate corresponds to a phase-
space element that is propagated. The phase-space element is also called
pseudo-particle in the following. The candidate module holds all information
about the (pseudo-)particle that can be used and altered by other modules in
the ModuleList that are successively called each simulation step. The user can
add different kinds of modules to the ModuleList: Propagation, Acceleration,
Interaction, Boundaries and Observers. Other modules handle the simulation
environment (MagneticField, AdvectionField), the injection of candidates
(Source) and the output. The modular structure makes it easy to set up
simulations for various scenarios and to add new modules.

To model DSA, we use the DiffusionSDE module [25]. Pseudo-particles are
propagated parallel and perpendicular to a magnetic field, depending on the
diffusion tensor

$\hat{\kappa}=\kappa_{\parallel}\begin{pmatrix}\epsilon&0&0\\\ 0&\epsilon&0\\\
0&0&1\end{pmatrix}\left(\frac{E}{E_{0}}\right)^{\alpha},$ (2.3)

which is constant in space and can depend on energy. The diffusion
coefficients for the normal and binormal direction of the magnetic background
field line can be assumed to be the same for large curvature radii of the
magnetic field lines
$\kappa_{\mathrm{n}}=\kappa_{\mathrm{b}}\equiv\kappa_{\perp}=\epsilon\kappa_{\parallel}$.
In Section 3.4 energy-dependent diffusion is considered. Other than that,
diffusion is considered to be energy-independent. The SDE 2.1 is integrated
using the Euler-Maruyama Scheme (e.g. [31])

$\vec{x}_{n+1}-\vec{x}_{n}=\vec{u}\Delta
t+\sqrt{2\hat{\kappa}}\Delta\vec{\omega}_{n},$ (2.4)

with time step $\Delta t$ and random variables
$\Delta\vec{\omega}_{n}=\vec{\omega}_{t_{n+1}}-\vec{\omega}_{t_{n}}$ drawn
from a normal distribution with an expected value of zero and variance of
$\Delta t$. SDE methods have no constraints on the used time step. However,
the choice of time step can be crucial to obtain correct results. See [33] for
a discussion. The diffusive step
$\sqrt{2\hat{\kappa}}\cdot\Delta\vec{\omega}_{n}$ is calculated in the
orthonormal base of the magnetic field, while the advective step
$\vec{u}\Delta t$ is calculated in the lab frame.

To ensure one-dimensional propagation in Section 3 the magnetic field is set
constant in x-direction and diffusion is allowed only parallel to the magnetic
field, i.e. $\epsilon=0$. In Section 4 perpendicular diffusion is taken into
account as well as a more complex magnetic field.

Acceleration is handled by the AdiabaticCooling module according to the
momentum equation (Eq. (2.2)) based on the divergence of the shock profile
specified in the AdvectionField class. The advection field and magnetic fields
used for simulations are specified in the following sections.

Free-escape boundaries $L_{\pm}$ can be specified in the upstream and
downstream region. Observer modules detect pseudo-particles when they cross
those boundaries or the TimeEvolutionObserver is used to detect pseudo-
particles’ positions and energies at given times $t_{i}$ during the simulation
in the acceleration region. All units are normalized so that
$x_{0}/(v_{0}\cdot t_{0})=1$.

###  2.3 Analysis

The energy spectrum $J(E)=\mathrm{d}N/\mathrm{d}E$ of particles in the
acceleration region is approximated by the histogram $\Delta N/\Delta E$, with
$\Delta N$ being the number of pseudo-particles in each energy bin $\Delta E$.
Since we expect a power-law spectrum energy is binned with equal distance in
logarithmic space. The error for each bin is then given by $\Delta
J=J/\sqrt{\Delta N}$. The spectrum can be weighted by $E^{2}$ to highlight the
spectral slope $s=-2$ predicted for acceleration at strong shocks.

###  2.4 Candidate Splitting

Diffusive Shock Acceleration produces a power-law spectrum in energy and the
spectral slope depends only on the shock compression ratio. When mono-
energetic candidates are injected at the shock and accelerated until they
leave the acceleration region, the number of candidates decreases with
increasing energy. Especially for less efficient acceleration at shocks with a
low compression ratio, statistics at high energies may be so bad that it is
difficult to evaluate whether the stationary solution is already reached or to
determine the resulting spectral slope.

The number of injected candidates can be increased in order to get sufficient
statistics at high energies, which comes at large computational costs. Another
way to overcome this problem is to split candidates in $n_{\mathrm{split}}$
copies once they cross energy boundaries in log-space. For DSA the optimal
splitting number depends on the expected spectral slope to compensate for the
loss of candidates at higher energies. Each splitted candidate is assigned
with a weight $w=1/n_{\mathrm{split}}$ during the simulation. In the later
analysis candidates are weighted accordingly to obtain the correct spectra. In
order to determine the error of the spectrum for each energy bin all
candidates can be used, which reduces the uncertainty at high energies.

In Fig. 1 resulting shock spectra simulated with and without candidate
splitting are compared at the same simulation time ($t=400\,t_{0}$). The
injected number of candidates is $N_{0}=10^{6}$ for the simulation without and
$N_{0}=10^{4}$ with candidate splitting. For the latter, the actual number of
simulated candidates increases during the simulation. Still, computational
time goes down by a factor of $25$222CPU Time (Apple M1 Pro) for
$N_{0}=10^{4}$ candidates with CandidateSplitting: $194\pm 7.92\,\mathrm{s}$.
CPU Time for $N_{0}=10^{6}$ candidates without CandidateSplitting: $5212\pm
75\,\mathrm{s}$. Both averaged over 8 runs. when candidate splitting is used,
since fewer candidates are simulated in low-energy bins.

Figure 1: Left: Spectra obtained with and without using the CandidateSplitting
module at simulation time $t=400\;t_{0}$ evaluated at the shock, $x=[0,2]$.
The shock compression ratio is $q=4$ and both simulations reproduce the
expected spectral slope $s=-2$. The simulation with candidate splitting
reaches higher energies up to $10^{7}\;E_{0}$ with smaller errors. Right:
Weighted spectra (top) and relative error $(J_{-2}-J)/J$ where $J_{-2}$ is the
predicted spectrum (bottom).

Splitting of candidates at energy bins in logarithmic space is one way of
importance sampling used also in other approaches to DSA [34]. Another way is
the introduction of artificial drift terms pushing particles to the shock [33,
35]. Such drift terms need to be extracted in the analysis. The candidate
splitting method is, however, easier to apply to more complex configurations
of the magnetic field and advection field.

Candidate splitting is implemented as an independent module in CRPropa3.2. It
is possible to define individual energy bins and number of splits. For use in
diffusive shock acceleration only the minimal and maximal energy as well as
the expected spectral index need to be specified.

##  3 One-Dimensional Diffusive Shock Acceleration

In the ensemble averaged picture, energy gain at shocks is described by
Equation (2.2). To apply this adiabatic description a finite shock width is
considered since the advective speed $u(x)$ must be continuously
differentiable. A smooth shock profile can also be motivated if interactions
between the particles and background plasma are taken into account [36]. The
one-dimensional shock profile is implemented

$\displaystyle
u(x)=\frac{u_{1}+u_{2}}{2}-\frac{u_{1}-u_{2}}{2}\tanh\left(\frac{x}{L_{\mathrm{sh}}}\right),$
(3.1)

with upstream and downstream velocity, $u_{1}$ and $u_{2}=u_{1}/q$, and shock
width $L_{\mathrm{sh}}$. In the following, all units are normalized so that

$x=\frac{\tilde{x}}{x_{0}},u=\frac{\tilde{u}}{u_{0}},t=\frac{\tilde{t}}{t_{0}},\kappa=\frac{\tilde{\kappa}}{\kappa_{0}},E=\frac{\tilde{E}}{E_{0}},$
(3.2)

with $x_{0}/v_{0}=t_{0}$. The velocity profile $u(x)$ is also used in other
studies of linear and non-linear DSA [37, 34, 28]. Figure 2 illustrates $u(x)$
for a shock compression ratio $q=u_{1}/u_{2}=4$ and different shock widths.
The narrower the shock width, the better the approximation of the ideal shock.
However, from the numerical perspective the region where $\partial u/\partial
x\neq 0$ must be large enough compared to the integration step length for
candidates to experience acceleration. Otherwise energy gains may be
underestimated (see Section 3.2 for a detailed analysis).

Figure 2: One-dimensional advection field with shock at $x=0$ and a shock
compression ratio $q=4$. The smaller the shock width compared to the
integration step length, the better the ideal shock solution is approximated.

###  3.1 Stationary Solution

Simulating with CRPropa 3.2, pseudo-particles are injected at $t=0$ at the
shock and propagated until a maximum simulation time or until they leave the
acceleration region downstream through free-escape boundaries.

In order to obtain a stationary solution to compare the results with shock
spectra predicted from theory, $s=2-3q/(q-1)$ (e.g. [38]), continuous
injection upstream must be assumed. The inherent structure of CRPropa 3.2 does
not allow for continuous injection during the simulation itself. To construct
the stationary solution at the shock afterwards CRPropa’s
TimeEvolutionObserver module is used. Positions and energies of pseudo-
particles are stored at times $t_{0}+n\Delta T$. With the assumption of
continuous injection, the solution at time $t$ is constructed by summing over
all such time snaps to approximate the time-integrated solution. This approach
is valid given that the solution does not change too much during one time
interval $\Delta T$. The time intervals $\Delta T$ can be chosen in
logarithmic space so that the time resolution decreases over simulation time,
when number density and spectrum undergo smaller changes. The resulting time
evolution of the spectrum at the shock and number density are shown in Fig. 3
compared to the solutions obtained by the integration of the Fokker-Planck
equation (1.1) using the same parameters. For the integration of the partial
differential equations VLUGR3 [26, 27] is used. The spectra are weighted by
energy so that the predicted slope lies horizontally in the plot. The
resulting spectra and number densities are in good agreement. Differences may
result from binning in energy and space for the CRPropa3.2 simulation and from
resolution in energy and time for the VLUGR3 simulation.

Figure 3:  Time evolution of spectrum at the shock, $x=[0,2]$ (left) and
number density in the acceleration region (right). Particles with energy
$E_{0}$ are continuously injected at the shock from time $t=0$. Upstream speed
is $u_{1}=1$ with a compression of $q=4$ at the shock. The diffusion
coefficient is constant in space and energy, $\kappa=1$. Free escape
boundaries are at $L_{\pm}=\pm 100$. Spectrum and number density resulting
from simulation with CRPropa3.2 (dots) are compared to the solutions obtained
by integrating the transport equation with VLUGR3 (lines).

The time evolution can be interpreted in the following ways: The shock is
active from time $t=0$ on and accelerated particles entering the shock region
with energy $E_{0}$. The longer the shock is active the higher the maximum
energy the particles can reach and the cut-off of the spectrum moves to higher
and higher energies. Having an infinite acceleration region, the stationary
solution of a planar one-dimensional shock with compression $q=4$ is a power-
law with slope $s=-2$. With increasing simulation time, the approximation of
the stationary solution gets better and the solution relaxes to the power-law
for ever higher energies. In Fig. 3 the stationary solution at $t=200$ is
reached for energies $E<10^{3}E_{0}$. Theoretically there is no limit for the
maximal energy, due to the infinitely large shock front all particles will
eventually cross the shock front again with $t\rightarrow\infty$.

In order to evaluate the approximation of the stationary solution the spectrum
is fitted at $t=500$ (not shown in Fig. 3) for energies below $10^{4}E_{0}$.
Up to that energy candidates are splitted by the CandidateSplitting module.
The resulting spectral slope of $-2.054\pm 0.003$ matches the expected value
$s=-2$ nicely.

The time evolution of the number density shows how particles gradually migrate
into the downstream region due to the advection field. At $t=200$ the
stationary solution is reached close to the shock, $x<30$, at later times also
further into the downstream region. The free-escape boundary at $x=100$ does
not influence the number density profiles at the shock. In principle, free
escape boundaries are only needed if the acceleration region has a finite
size. In Fig. 17 we show how the resulting spectrum is affected by free-escape
boundaries. Only a fraction of particles makes it into the upstream region
$x<0$ against the advection flow. Acceleration is efficient for a constant,
low diffusion coefficient as assumed here, however, particles are unlikely to
escape upstream.

Analytical solutions exist only for the special case of
$v^{2}/\kappa=\mathrm{const.}$ across the shock [39, 38], resulting in a drop
of spatial diffusion by $1/16$ for a strong shock. Spatial dependence of the
diffusion coefficient is not considered here, since the resulting drift term
which adds to the advective step in Equation (2.1) is currently not
implemented in CRPropa3.2. For SDEs the strong additional drift term that is
expected at the shock may also lead to numerical difficulties [34]. The
integration of the transport equation using VLUGR3 is verified against
existing analytical solutions in earlier works [28].

###  3.2 Constraints

In order to reproduce the predicted spectrum the simulation set-up must
fulfill constraints on the integration step length and shock width. For SDEs
Achterberg and Krülls [37] and later Achterberg and Schure [34] present a
thorough analysis of the choice of shock width and step length and the
resulting spectral slope. The simulations resulting in Fig. 3 take their
findings into account. In the following we show that DSA simulations with
CRPropa3.2 are subject to the same constraints. These constraints are not
purely numerical but can also be motivated from physics.

Considering a diffusion coefficient constant in space and energy, the
advective step length, diffusive step length and shock width must fulfill the
inequality

$u\Delta t<L_{\mathrm{sh}}\lesssim\sqrt{\kappa\Delta t},$ (3.3)

to obtain the correct spectral index [37]. The first inequality ensures that
pseudo-particles experience the gradient of the advection field and therefore
gain energy. On the other hand, the diffusive step length $\sqrt{\kappa\Delta
t}$ \- a measure for the stochastic step - must be larger than the shock width
to increase the likelihood of pseudo-particles to cross the shock front
multiple times. When diffusion is too low, particles will not make it back to
the shock against the advection flow. If the inequality is not fulfilled, the
modelled energy gain may be less efficient leading to steeper spectra, or
depending on the ratio of advective step to the shock width, even more
efficient leading to flatter spectra than predicted.

Obeying Equation (3.3) alone is not sufficient to reproduce the ideal shock
spectra. If the shock width and with that the integration time step are chosen
too large, pseudo-particles experience a smoothly changing advection field
instead of a discrete shock. For diffusion independent on space and energy
Krülls and Achterberg [37] calculated the expected spectral slopes for such
smooth velocity gradients. Also, when a finite acceleration region is
considered, diffusion must not only be large enough so that particles are able
to diffusive back to the shock but also low enough so that they are contained
in the acceleration region and do not escape too quickly.

#### Advective Step

Acceleration at shocks of different widths $L_{\mathrm{sh}}$ is simulated with
decreasing time steps $\delta T$. The normalized diffusion coefficient is
$\kappa=1$ and the upstream velocity $u_{1}=1$. The shocks compression ratio
is $q=2$ leading to a spectral slope $s=-4$ in order to compare with the
predictions from [37]. The slope is determined by linear fits in log-space up
to the maximum energy of the CandidateSplitting module. Time is chosen such
that the cut-off of the spectrum has minimal impact on the slope approximating
the stationary solution. The results are shown in Fig. 4.

The expected spectral slopes depending on the shock widths are indicated by
lines. With increasing resolution the simulated spectrum well approximates the
predicted stationary solution. For sufficiently small shock widths, here
$L_{\mathrm{sh}}=0.008$, a spectral slope close to that of an ideal shock is
obtained. The larger the shock width, the less efficient the acceleration
which leads to steeper spectra. Depending on the shock width, with decreasing
time step, the simulation approaches the predicted stationary solution from
steeper ($L_{\mathrm{sh}}=0.128$) or flatter ($L_{\mathrm{sh}}=0.008$)
spectra. Thus, the compression ratio of the shock is over- or underestimated.

We show a more detailed figure with entangled information about the shock
width and time step in the Appendix in an analogous way to [37].

#### Diffusive Step

A similar analysis can be done for the dependence of the diffusion coefficient
in relation to the shock width. The larger the diffusive step, the better
becomes the approximation of an ideal shock. If the diffusive step however is
too large, pseudo-particles may miss the shock region and acceleration is
underestimated as a consequence.

For a constant shock width of $L_{\mathrm{sh}}=0.004$ the diffusion
coefficient was varied over several simulation runs. The resulting spectral
slope was fitted in the energy range $[E_{0},10^{3}E_{0}]$ given by the
CandidateSplitting. The results are shown in Fig. 4 in comparison to the
results from [34]. With large diffusive step length the spectrum may also be
affected by free-escape boundaries. For the simulations boundaries were chosen
such that they do not influence the time evolution of the spectrum for the
simulated diffusion coefficients. Likely, that is why we do not find the
spectral slope to decrease again for $\epsilon<0.04$. In the Appendix we show
how the free-escape boundaries influence the spectrum at the shock depending
on the diffusion coefficient.

Figure 4: Left: Resulting spectral slope depending on the time step $\Delta t$
in relation to the shock width $L_{\mathrm{sh}}$. The upstream speed
$u_{1}=1$. The compression ratio is $q=2$, thus for an ideal shock we expect
$s=-4$. With larger shock widths, the shock produces steeper spectra as
calculated by [37]. Slopes of the predicted stationary solutions are indicated
by the colored lines for the respective shock width. Right: Resulting spectral
slope depending on the ratio of upstream velocity and shock width to diffusion
coefficient. Data from [34] is shown for comparison. Here $u_{1}=1$, $q=4$ and
$L_{\mathrm{sh}}=0.004$ with $\Delta t=0.001$ is used while $\kappa$ is
varied. The free escape boundaries are $L_{\pm}=10$. With higher diffusion
coefficient the spectral slope approaches the predicted value $s=-2$. We
already show the results for a transversal shock wave and taking perpendicular
diffusion into account, which is discussed in Section 4.1.

We conclude that for the simulation of ideal shocks, the shock width must be
small compared to the advective field and the time step should be chosen such
that the advective step is at most as large as one fourth of the shock width.
The diffusive step must be larger than the shock width, but not too large to
still contain particles in the acceleration region. Those constraints may be
difficult to fulfill when modeling a physical scenario and may require very
small time steps close to the shock.

###  3.3 Finite Lifetime

We model two scenarios that lead to a shock spectrum that is not stationary:
Either the source of particles reaching the shock is only active for a short
period of time or the shock itself has a finite lifetime. Considering the
GWTS, the CR flux that reaches the shock from the Galactic disk can be assumed
to be stationary. The GWTS itself however cannot be maintained for an infinite
time since too much mass would be advected out of the Galaxy.

To model a burst-like particle source, candidates are injected at $t=0$ and
propagated until they leave the acceleration region through the free-escape
boundaries $L_{\pm}$. The time evolution depending on the source duration can
be derived similar to the stationary spectrum in Section 3.1. Energy and
position of candidates are stored in time intervals $\Delta T$ and the
resulting spectrum is approximated by summing of those time snapshots from
$[t-\Delta T_{\mathrm{src}},t]$ with $\Delta T_{\mathrm{src}}$ being the
source lifetime. For more details we point to [9]. Figure 5 shows the time
dependent spectrum at the shock for two different source duration $\Delta
T_{\mathrm{src}}$.

Figure 5: Candidates are injected at $t=0$ and simulated until they leave the
acceleration region through free-escape boundaries $L_{\pm}=\pm 100\;x_{0}$. A
finite lifetime $\Delta T_{\mathrm{scr}}=20$ (left) and $\Delta
T_{\mathrm{scr}}=80$ (right) of the particle source is assumed. The spectrum
at the shock $x=[0,2]$ is shown.

Considering a finite lifetime of the shock, we already know the time dependent
spectrum at the shock from Section 3. In order to assess the contribution to
the overall cosmic-ray spectrum, not only the spectrum at the shock, but also
the spectrum of particles escaping downstream, or even upstream is of
interest. Considering a free-escape boundary at $L_{+}=100$, Fig. 6 shows the
time-dependent spectrum of particles escaping for two different lifetimes of
the shock $\Delta T_{\mathrm{sh}}$. With a constant, low diffusion coefficient
$\kappa=1$ acceleration is efficient but particles do not escape upstream. In
the following Section we introduce energy-dependent diffusion: With increasing
energy, diffusion is higher and particles are more likely to escape upstream.

Figure 6: Candidates are injected at $t=0$ and simulated until they leave the
acceleration region through free-escape boundaries $L_{\pm}=\pm 100\;x_{0}$. A
finite lifetime $\Delta T_{\mathrm{scr}}=100$ (left) and $\Delta
T_{\mathrm{scr}}=400$ (right) of the shock is assumed. The spectrum of
particles escaping through the downstream boundary is shown.

###  3.4 Energy-Dependent Diffusion

So far, a diffusion coefficient constant in energy was considered. A more
realistic description of diffusive motion takes energy-dependence into
account. Assuming that the diffusion coefficient is proportional to
$(E/E_{0})^{\alpha}$ and $\alpha>0$, the diffusive step becomes larger with
increasing energy. At first glance it is then easy to satisfy the inequality
in Equation (3.3). However, the analysis in Section 3.2 revealed that the
diffusive step should not become too large and free-escape boundaries should
be dropped or must be set far away from the shock, otherwise high-energy
particles escape too quickly.

In order to keep the diffusive time step within a reasonable range, an energy-
dependent adaptive time step is implemented in the DiffusionSDE module. Within
a range of specified time steps
$[\mathrm{d}t_{\mathrm{min}},\mathrm{d}t_{\mathrm{max}}]$, it is chosen such
that the advective step is smaller than $1/4\,L_{\mathrm{sh}}$ and the
diffusive time step smaller than $100\,L_{\mathrm{sh}}$. This also leads to
better performance in case of the energy independent diffusion coefficient
($\alpha=0$), since the time step in the downstream region can be chosen
higher than in the upstream region by a factor of the shock compression ratio
$q$.

Figure 7 shows the time evolution of the spectrum at the shock and number
density in the acceleration region. Again the compression ratio is $q=4$ and
the predicted spectral slope $s=-2$. The diffusion coefficient is
$\kappa(E)=\kappa_{E_{0}}\left(E/E_{0}\right)^{\alpha}$ with $\alpha=1$ and
$\kappa_{E_{0}}=1$. The results obtain with the diffusion approach applied in
CRPropa3.2 are compared to the simulations using VLUGR3. The same times as in
Fig. 3 are used.

With $\alpha>0$ the diffusion coefficient $\kappa(E)$ is equal or greater than
in the energy-independent case shown in Section 3.1. With that the
acceleration time scale gets larger. The cut-off of the energy spectrum at the
shock is at lower energies compared to Fig. 3.

With higher diffusion at high energies, more particles make it in the upstream
region against the advection flow. Since high energies are only reached at
later times, the upstream number density also evolves over time in contrast to
Fig. 3 for energy-independent diffusion. Thus, the likelihood of escaping
upstream gets higher with energy, thus acceleration time, and increasing
$\alpha$.

Figure 7: Time evolution of the spectrum (right) at the shock and the number
density (left) close to the shock with energy-dependent diffusion coefficient
$\kappa=\kappa_{E_{0}}\left(E/E_{0}\right)$ compared to the results obtained
by integrating the transport equation using VLUGR3.

The cut-off energy for a planar shock is still only determined by the free-
escape boundaries and acceleration time. Here no free escape boundaries are
used, so that the finite acceleration region does not influence the spectrum
at the shock. In physical scenarios, high energy particles are likely to
escape the acceleration region.

In Fig. 8a the spectra at the shock for different values of $\alpha$ are
compared at time $t=200$. The stronger the energy dependence, the higher the
diffusion at high energies and the longer it takes to accelerate particles to
the same energy.

In case of energy-dependent diffusion, the upstream spectrum differs from the
downstream spectrum. In Fig. 8b we show the time evolution of the upstream
spectrum at $x=-5$ for $\alpha=1$.

Figure 8: Left: Spectrum at the shock at time $t=200$ with energy-dependent
diffusion $\kappa=\kappa_{0}\left(E/E_{0}\right)^{\alpha}$ and different
values of $\alpha$. Energy spectrum obtained with CRPropa3.2 (dots) is
compared against those resulting from VLUGR3 (lines). Right: Upstream
spectrum, $x=[-6,-5]$, with energy-dependent diffusion ($\alpha=1$). For
energy-dependent diffusion the upstream spectrum differs from the downstream
spectrum. Again CRPropa3.2 results are compared to VLUGR3.

###  3.5 Injection of Pre-accelerated Spectra

So far mono-energetic pseudo-particles of energy $E_{0}$ have been injected at
the shock. In a physical scenario of particle acceleration at the GWTS,
particles are assumed to arrive with an energy distribution of a power law
[8]. Depending on the acceleration process a range of spectral slopes is
possible.

We injected different spectra at the shock and investigated the impact on the
resulting re-accelerated spectra. According to [38] steeper spectra than the
one produced by the shock converge to the shock spectrum. When, on the other
hand, the shock produces a steeper spectrum than the one continuously
injected, the injected spectrum does not change. This results in broken power-
law spectra when the injected spectrum has a finite cut-off energy as shown in
Fig. 10.

Figure 9: Time evolution of the spectrum at the shock $x=[0,2]$ with injected
spectra of slope $s=-1$ (left) and $s=-3$ (right). The injected spectra have a
hard cut-off at $10^{2}E_{0}$. For $s=-1$ up to $10^{2}E_{0}$ the spectral
slope remains the same, for higher energies the spectrum evolves to the one
produced by the shock, $s=-2$. For $s=-3$ the spectrum converges to the shock
spectrum for all energies.

In Fig. 9a a harder spectrum than produced by the shock is continuously
injected up to the maximum energy $10^{2}\,E_{0}$. Up to that energy, the
injected spectrum prevails. At the shock, particles are accelerated
independent on their energy and the time evolution of the spectrum at the
shock becomes visible for $E>10^{2}\,E_{0}$. In Fig. 9b a flatter spectrum
than produced by the shock is continuously injected up to the maximum energy
$10^{2}\,E_{0}$. The spectrum converges to the one produced by the shock over
time, since all particles independent on their energy experience acceleration.
Thus, at the shock particles are re-accelerated to higher energies.

In case of harder spectra, broken power-law spectra may emerge. In case of
softer spectra, the spectrum converges to the shock spectrum. If the injected
spectrum is the same as produced by the shock it remains the same. This is
summarized in Fig. 10: The stationary spectrum up to $10^{4}E_{0}$ at the
shock with different injected spectra from $10^{-1}E_{0}-10^{2}E_{0}$ is
shown.

Figure 10:  Spectra at $t=400$ at the shock for different injected spectra
from $10^{-1}E_{0}$ up to $10^{2}E_{0}$. For $s<-2$ the spectrum converges to
the shock spectrum. For $s>-2$ the spectrum breaks at $10^{2}E_{0}$.

##  4 Anisotropic Diffusion, Oblique Shocks and Spherical Symmetry

In the following, we extend our analysis to anisotropic diffusion, oblique
planar shocks and spherical symmetry. The diffusion coefficient is held
constant in space and energy. Note that since the diffusion coefficient is
defined in the local magnetic field coordinates, diffusion changes with the
direction of the magnetic field even when $\kappa$ is held constant. In
Section 4.3 a spherical shock profile with a spiral magnetic field is modeled.

###  4.1 Perpendicular Shock Wave

We consider a perpendicular shock wave with the magnetic field perpendicular
to the advection flow. If only parallel diffusion is considered, particles
would not be able to cross the shock and, therefore, do not experience
acceleration. The advective field is still given by Equation (3.1), the
magnetic field now points in y-direction, $\vec{B}=B_{0}\vec{e}_{\mathrm{y}}$.

The perpendicular diffusion coefficient now simply takes the role of parallel
diffusion. We can find the same dependency between perpendicular diffusion
$\kappa_{\perp}=\epsilon\kappa_{\parallel}$ and shock width as shown in Fig. 4
in Section 3.2. The parallel part of the diffusion tensor does thus not
influence the acceleration process.

###  4.2 Oblique Shock Wave

First, only the advection field is considered to be oblique and the magnetic
field remains parallel to the shock normal. Only the parallel component to the
shock normal of the advection field, $u_{\mathrm{x}}$, is shocked, so that it
is now given by

$\displaystyle u_{\mathrm{x}}$
$\displaystyle=\frac{u_{\mathrm{x},1}+u_{\mathrm{x},2}}{2}-\frac{u_{\mathrm{x},1}-u_{\mathrm{x},2}}{2}\tanh\left(\frac{x}{L_{\mathrm{sh}}}\right),$
(4.1) $\displaystyle u_{\mathrm{y}}$ $\displaystyle=\mathrm{const.}$ (4.2)

The $y$-component is varied from $u_{\mathrm{y}}=0$ (one-dimensional case) to
$u_{\mathrm{y}}=5u_{\mathrm{x,1}}$, while $u_{\mathrm{x,1}}=1$ like. The angle
$\theta$ between shock normal and advection flow varies from $\theta_{1}=0$ to
$\theta_{1}\approx 0.44\pi$. From theory, the spectrum still depends only on
the shock compression ratio
$q=u_{1}\mathrm{cos}\theta_{1}/u_{2}\mathrm{cos}\theta_{2}$ until the shock
vanishes for $\theta\rightarrow\pi/2$. Considering the constraints on the time
step and shock width, particles now experience varying "effective" shock
widths when they cross the shock with different angles to the shock normal.
Thus, it becomes difficult to obtain the expected shock spectra for large
angle $\theta$. The effective shock width can be written as
$L_{\mathrm{sh,eff}}=L_{\mathrm{sh}}/\mathrm{cos}\theta$. With a shock width
of $L_{\mathrm{sh}}=0.001$ and adaptive time step $\Delta
t\in\left[L_{\mathrm{sh}}/4,L_{\mathrm{sh}}/16\right]$ the predicted spectral
slope $s=-2$ is obtained. The time evolution of the spectrum at the shock is
independent on the angle between the shock and advection field and the amount
of perpendicular diffusion to the magnetic field. Candidates are now injected
upstream $x_{\mathrm{src}}=-1$, so that they are not trapped inside the shock
for large angles of $\theta$.

Now the magnetic field is considered to be parallel to the advection field.
Both the advection field and magnetic field break at the shock. Here,
curvature drifts due to the bent magnetic field lines at the shock are
neglected. Currently drift terms are not included in the implementation of the
diffusion coefficient in CRPropa3.2. Still, we expect to see acceleration at
the shock up to the predicted spectrum. Figure 11 shows that with parallel
diffusion only for large angles between magnetic field and shock normal, the
spectral slope at the shock is steeper than predicted as acceleration is less
efficient. However, introducing diffusion perpendicular to the magnetic field,
particles are again able to cross the shock even if the magnetic field is
almost perpendicular to the shock normal. A spectral slope $s=-2$ is obtained
even for almost perpendicular shocks.

Figure 11: Resulting spectra at the shock for different angles between
advection field and shock normal at $t=400$. The magnetic field is considered
parallel to the advection field. With parallel diffusion only (left) and
increasing angle between magnetic field and shock normal, efficient
acceleration becomes difficult and the predicted spectral slope $s=-2$ is only
obtained for small angles. Introducing a small amount of perpendicular
diffusion (right) with $k_{\perp}=0.2$ leads to the predicted spectral slope.

###  4.3 Spherical symmetry

We construct a simplified model of the GWTS, namely a spherical symmetric
shock with spiral magnetic field that takes the findings of the previous
sections into account. With large shock radii astrophysical shocks can in
principle be approximated by planar shocks with all implications described
previously. However, for sufficiently small shock radii e.g. [38] showed, that
the spectral slope in spherical geometry is influenced by the energy- and
spatial dependence of the diffusion coefficient. In contrast to the one-
dimensional planar shock, for spherical shocks a limit for the maximum energy
is expected, since the shock front is not infinitely large.

In the following, the advection field is assumed to be constant upstream,
$r<r_{\mathrm{sh}}$, and to decrease with $r^{-2}$ downstream,
$r>r_{\mathrm{sh}}$, modeled by

$\displaystyle\vec{u}(r)=u_{1}\left\\{1+\left[\left(\frac{R_{\mathrm{sh}}}{2r}\right)^{2}\frac{1}{1+\mathrm{e}^{-\frac{r-R_{\mathrm{sh}}}{L_{\mathrm{sh}}}}}\right]\right\\}\vec{e}_{\mathrm{r}}.$
(4.3)

Analogous to Equation (3.1), a finite shock width $L_{\mathrm{sh}}$ is
considered. Similar profiles were also used in other studies of spherical DSA
[40] and applied to the GWTS [15, 9]. The advection field is shown in Fig. 12
for different shock widths.

Figure 12: Radial shock profile. The advective field is constant upstream
$u(r)=u_{1}$, drops by $1/q$ at the shock and decreases downstream with
$r^{-2}$ for $r>R_{\mathrm{sh}}$.

For the magnetic field, analogous to [9], an Archimedean spiral is considered.
It is given by

$\displaystyle\vec{B}=\pm
B_{0}\left[\left(\frac{r_{0}}{r}\right)^{2}\vec{e}_{\mathrm{r}}-\frac{\Omega
r_{0}^{2}\sin{\theta}}{rv_{\mathrm{w}}}\vec{e}_{\phi}\right]$ (4.4)

with the constant wind velocity $v_{\mathrm{w}}$ and $\Omega
r_{0}\sin{\theta}$ being the rotational velocity at a reference radius $r_{0}$
at latitude $\theta$. The magnetic field lines for different values of
$\Omega$ are illustrated in Fig. 13. For $\Omega=0$ the magnetic field is
parallel to the shock normal and the advection field. For $\Omega\neq 0$ the
magnetic field is not parallel to the wind profile and the break of the
magnetic field at the shock is neglected. For high values of $\Omega$ the
magnetic field is almost perpendicular to the shock normal and resembles the
transversal shock wave investigated in Section 4.1.

In Fig. 14 spectrum and number density are shown for $\Omega=10$. Here,
perpendicular diffusion was also taken into account to achieve sufficient
acceleration. In both simulations the stationary shock spectrum matches the
predicted spectral slope $s=-2$ for a strong shock. Due to the decelerating
flow for $r>R_{\mathrm{sh}}$ the number density profile differs from the
planar scenario.

Figure 13: Face on view of the Archimedean Spiral magnetic field for
increasing values of $\Omega=[0.1,1.0,10]$ from left to right. For small
angular velocities the magnetic field is parallel to the shock normal, for
high angular velocity magnetic field lines become perpendicular.  Figure 14:
Time evolution of the number density (right) and spectrum at the shock,
$r=[r_{\mathrm{sh}},r_{\mathrm{sh}}+2]$ (left). Free escape boundary is
downstream $R_{+}=200$. Both parallel and perpendicular diffusion are
considered $\kappa_{\parallel}=1$, $\kappa_{\perp}=0.8$ and the magnetic field
is almost perpendicular to the shock normal, $\Omega=10$.

Figure 15 shows the time evolution of the spectrum at the shock and the number
density in the acceleration region for $\Omega=1$. Different values for the
perpendicular diffusion are considered. Introducing diffusive motion
perpendicular to the magnetic field, leads to higher diffusion in total. The
time evolution of the spectrum at the shock is slower with higher
perpendicular diffusion.

Figure 15:  Time evolution of the spectrum at the shock,
$r=[r_{\mathrm{sh}},r_{\mathrm{sh}}+2]$ for different values of perpendicular
diffusion, $\epsilon=[0,0.2,0.8]$ from left to right. Free escape boundary is
downstream $R_{+}=200$. The lower the perpendicular diffusion, the slower the
time evolution at the shock.

##  5 Summary and Outlook

In this paper we present a thorough analysis of modeling DSA with the
DiffusionSDE module of CRPropa as a basis for future studies of the GWTS.
CRPropa3.2 provides solvers for both transport regimes, and thus offers the
opportunity to model acceleration at the GWTS and propagation back to the
Galaxy in a coherent framework. Taking into account the constraints discussed
in Section 3.2, time-dependent acceleration at shocks in three dimensions can
be simulated with CRPropa3.2. Arbitrary magnetic field and advection field
configurations can be used. Diffusion can be constant or energy-dependent as
well as anisotropic. The new CandidateSplitting module speeds up the
simulation significantly.

We validated the resulting time-dependent spectra at a one-dimensional planar
shock against those obtained from integrating the Fokker-Planck equation with
the finite difference solver VLUGR3 and presented first applications
introducing energy-dependent and anisotropic diffusion as well as spherical
shock geometry.

With the Acceleration module, CRPropa3.2 also offers the possibility to model
DSA in the ballistic regime. This was first introduced for second-order Fermi
acceleration [41] and later extended to DSA [24]. Each time a particle crosses
the shock front, the module calculates the energy gain performing Lorentz
transformations from the particle frame to the rest frame of the scatter
centers, the upstream or downstream background flow. However, this approach is
less flexible than the simulation based on the transport equation or the
corresponding SDEs.

In Section 3 we show that considering a one-dimensional planar parallel shock,
simulations lead to the expected spectral slope at the shock. We validated the
correct time evolution of the spectrum at the shock by comparing with
simulations based on the integration of the transport equation (see [28]). In
Section 3.2 we summarize constraints on the shock width, diffusion and time
step first described by [37] and [34]. We find the same dependence between
spectral slope, shock width and advective step, as well as diffusive step.
Likewise we found that depending on the shock width for too small time steps,
the compression ratio of the shock is over- or underestimated (see Fig. 4 and
Fig. 16). Presumably, this is due to the pseudo-particles only meeting the
shock region by chance: The smaller the shock width the stronger the velocity
gradient they encounter. This leads to harder spectra at least for those
pseudo-particles experiencing acceleration. Section 3.2 explains in detail how
to set up DSA simulations using SDEs in general and CRPropa3.2 specifically.

Based on the findings in Section 3.2 for DSA an alternative adaptive time step
was implemented in the DiffusionSDE module. This ensures that within a user
specified range the largest possible time step based on the inequalities (Eq.
(3.3)) for advective and diffusive step is used. Another extension to
CRPropa3.2 is the CandidateSplitting module. High energies are only reached by
a small fraction of candidates. For better statistics at high energies,
candidates crossing energy boundaries are split into $n$ copies depending on
the expected spectral slope. This way of importance sampling significantly
reduces computation time.

We approached more physical scenarios by considering energy-dependent
diffusion, pre-accelerated spectra and shocks with finite lifetime. Depending
on the life-time of the shock, the spectrum seen by a downstream observer
differs from the stationary spectrum. Energy or spatial dependence of the
diffusion coefficient would further modulate the downstream spectrum over
time. The effects on the shock spectrum discussed in sections 3.3 to 3.5 were
deliberately considered separately and are the basis for future combined
studies, potentially adding spatial dependent diffusion.

In Section 4 we took anisotropic diffusion into account and showed that planar
oblique shocks produce decent results, when diffusion parallel to the shock
normal is high enough. A prototype for simulating the GWTS was presented in
Section 4.3: A radial shock profile and a Archimedean spiral magnetic field.
The time-dependent spectrum obtained at the GWTS as well as at the upstream
and downstream escape boundaries can be used in further studies, investigating
the portion of particles that is able to propagate back to the Galaxy (as done
for a stationary shock spectrum by [9]). Only high-energy CRs are able to
propagate against the outflowing wind, so that energy-dependent diffusion must
be considered. We already saw in Section 3.4 that energy-dependent diffusion
leads to a higher fraction of particles being able to diffusive back into the
upstream region against the wind.

For high-energy particles leaving the acceleration region the diffusion
approximation might not be valid anymore. At this point we stress again the
advantage of CRPropa3.2: Follow up simulations can be done using both
diffusive and ballistic propagation depending on the candidates energy within
the same framework.

Not only energy-dependent diffusion but also a spatial dependence might impact
the acceleration process as well as the cosmic-ray transport in the Galactic
magnetic field. A spatially varying diffusion tensor, however, induces an
additional drift term to Equation (2.1) which is not yet implemented in
CRPropa3.2. Also, this term adds to the left-hand side of the inequality in
Equation (3.3) discussed in Section 3.2. This makes it more difficult to
choose a valid time step for a given shock width. Achterberg and Schure [34]
even found that a second-order scheme might be necessary for handling the
strong drift terms that would occur at the shock when the diffusion
coefficient drops over the shock.

The GWTS is most likely not a perfect sphere but like the termination shock of
the heliosphere a spheroid. Also, it does not necessarily enclose the complete
Galaxy. With a more variable wind profile CRs may even propagate back to the
Galaxy more easily. A more realistic wind profile presumably changes the time
evolution of the spectrum produced by the shock itself. When a finite shock
lifetime is considered this may have an impact on the contribution to the
shock spectrum at Earth.

With the radial wind profile and Archimedean spiral magnetic field we modeled
a quasi three-dimensional shock configuration. For the Galaxy, there exist
more realistic magnetic field models, like [42, 43] or [44]. The magnetic
field structure certainly has a great impact on the possibility of CRs
travelling back to the Galaxy. Merten et al. [9] showed that the length of the
magnetic field lines and the amount of diffusion parallel to the magnetic
field lines is critical for a contribution of CRs re-accelerated at the GWTS
to the spectrum. The magnetic field also has an impact on the arrival
direction of CRs. With CRPropa3.2 the magnetic field can easily be exchanged
and the impact on the spectrum and arrival direction on Earth can be studied
for various magnetic fields.

##  Appendix A Advective Step

Figure 16:  Spectral slope of the stationary solution at the shock depending
on the time step (indicated by color) and the ratio of upstream advective step
and shock width $n$. The expected spectral index for finite shock widths is
indicated by horizontal lines (data from [37]).

##  Appendix B Diffusive Step

Figure 17: Slope of the stationary spectra at the shock depending on the
diffusion coefficient $\kappa$ for different free-escape boundaries $L_{\pm}$.
With the diffusion coefficient too large compared to the acceleration region,
acceleration is less efficient and the spectral slope decreases.

## Acknowledgments

We acknowledge support from the DFG within the Collaborative Research Center
SFB1491 "Cosmic Interacting Matters - From Source to Signal". Additional
support from DFG projects EF98/4-1 and FI706/26-1 is acknowledged.

## References

  * [1] J. Becker Tjus and L. Merten, _Closing in on the origin of Galactic cosmic rays using multimessenger information_ , _Phys. Rep._ 872 (2020) 1 [2002.00964].
  * [2] E. Fermi, _On the Origin of the Cosmic Radiation_ , _Phys. Rev._ 75 (1949) 1169.
  * [3] W.I. Axford, E. Leer and J.F. McKenzie, _The structure of cosmic ray shocks_ , _A &A_ 111 (1982) 317.
  * [4] A.R. Bell, _The acceleration of cosmic rays in shock fronts - I._ , _MNRAS_ 182 (1978) 147.
  * [5] A.R. Bell, _The acceleration of cosmic rays in shock fronts - II._ , _MNRAS_ 182 (1978) 443.
  * [6] R. Kallenbach, M. Hilchenbach, S.V. Chalov, J.A. Le Roux and K. Bamert, _On the “injection problem” at the solar wind termination shock_ , _A &A_ 439 (2005) 1.
  * [7] A.M. Hillas, _The Origin of Ultra-High-Energy Cosmic Rays_ , _ARA &A_ 22 (1984) 425.
  * [8] S. Thoudam, J.P. Rachen, A. van Vliet, A. Achterberg, S. Buitink, H. Falcke et al., _Cosmic-ray energy spectrum and composition up to the ankle: the case for a second Galactic component_ , _A &A_ 595 (2016) A33 [1605.03111].
  * [9] L. Merten, C. Bustard, E.G. Zweibel and J. Becker Tjus, _The Propagation of Cosmic Rays from the Galactic Wind Termination Shock: Back to the Galaxy?_ , _ApJ_ 859 (2018) 63 [1803.08376].
  * [10] D. Breitschwerdt, J.F. McKenzie and H.J. Voelk, _Galactic winds. I. Cosmic ray and wave-driven winds from the galaxy._ , _A &A_ 245 (1991) 79.
  * [11] J.E. Everett, E.G. Zweibel, R.A. Benjamin, D. McCammon, L. Rocks and I. Gallagher, John S., _The Milky Way’s Kiloparsec-Scale Wind: A Hybrid Cosmic-Ray and Thermally Driven Outflow_ , _ApJ_ 674 (2008) 258 [0710.3712].
  * [12] M. Uhlig, C. Pfrommer, M. Sharma, B.B. Nath, T.A. Enßlin and V. Springel, _Galactic winds driven by cosmic ray streaming_ , _MNRAS_ 423 (2012) 2374 [1203.1038].
  * [13] C. Bustard, E.G. Zweibel and E. D’Onghia, _A Versatile Family of Galactic Wind Models_ , _ApJ_ 819 (2016) 29 [1509.07130].
  * [14] S.A. Mao and E.C. Ostriker, _Galactic Disk Winds Driven by Cosmic Ray Pressure_ , _ApJ_ 854 (2018) 89 [1801.06544].
  * [15] J.R. Jokipii and G. Morfill, _Ultra-High Energy Cosmic Rays in a Galactic Wind and Its Termination Shock_ , _ApJ_ 312 (1987) 170.
  * [16] P. Mukhopadhyay, E. Peretti, N. Globus, P. Simeon and R. Blandford, _Reacceleration of Galactic Cosmic Rays Beyond the Knee at the Termination Shock of a Cosmic-Ray-Driven Galactic Wind_ , _arXiv e-prints_ (2023) arXiv:2301.08902 [2301.08902].
  * [17] C. Bustard, E.G. Zweibel and C. Cotter, _Cosmic Ray Acceleration by a Versatile Family of Galactic Wind Termination Shocks_ , _ApJ_ 835 (2017) 72 [1610.06565].
  * [18] A. Kääpä, K.-H. Kampert and J. Becker Tjus, _Flux predictions in the transition region incorporating the effects from propagation of cosmic rays in the Galactic magnetic field_ , _EPJ Web Conf._ 283 (2023) 03006.
  * [19] R. Schlickeiser, _Cosmic ray astrophysics_ , Springer Science & Business Media (2002).
  * [20] J.G. Kirk, R. Schlickeiser and P. Schneider, _Cosmic-Ray Transport in Accelerating Flows_ , _ApJ_ 328 (1988) 269.
  * [21] P. Reichherzer, J. Becker Tjus, E.G. Zweibel, L. Merten and M.J. Pueschel, _Turbulence-level dependence of cosmic ray parallel diffusion_ , _MNRAS_ 498 (2020) 5051 [1910.07528].
  * [22] P. Reichherzer, L. Merten, J. Dörner, J. Becker Tjus, M.J. Pueschel and E.G. Zweibel, _Regimes of cosmic-ray diffusion in Galactic turbulence_ , _SN Applied Sciences_ 4 (2022) 15 [2104.13093].
  * [23] E. Amato and P. Blasi, _Non-linear particle acceleration at non-relativistic shock waves in the presence of self-generated turbulence_ , _MNRAS_ 371 (2006) 1251 [astro-ph/0606592].
  * [24] R. Alves Batista, J. Becker Tjus, J. Dörner, A. Dundovic, B. Eichmann, A. Frie et al., _CRPropa 3.2 - an advanced framework for high-energy particle propagation in extragalactic and galactic spaces_ , _JCAP_ 2022 (2022) 035 [2208.00107].
  * [25] L. Merten, J. Becker Tjus, H. Fichtner, B. Eichmann and G. Sigl, _CRPropa 3.1 - a low energy extension based on stochastic differential equations_ , _JCAP_ 2017 (2017) 046 [1704.07484].
  * [26] J.G. Blom and J.G. Verwer, _VLUGR3: A vectorizable adaptive grid solver for PDEs in 3D - Part I: Algorithmic aspects and applications_ , _App. Num. Math._ 16 (1994) 129.
  * [27] J.G. Blom and J.G. Verwer, _VLUGR3: A vectorizable adaptive-grid solver for PDEs in 3D - Part II: Code description_ , _ACM Transactions on Mathematical Software (TOMS)_ 22 (1996) 329.
  * [28] D. Walter, F. Effenberger, H. Fichtner and Y. Litvinenko, _A nonlinear model of diffusive particle acceleration at a planar shock_ , _Physics of Plasmas_ 29 (2022) 072302.
  * [29] A. Kopp, I. Büsching, R.D. Strauss and M.S. Potgieter, _A stochastic differential equation code for multidimensional Fokker-Planck type problems_ , _Computer Physics Communications_ 183 (2012) 530.
  * [30] K. Itô, _On stochastic differential equations_ , no. 4, American Mathematical Soc. (1951).
  * [31] P.E. Kloeden and E. Platen, _Stochastic differential equations_ , vol. 23, Springer (1992).
  * [32] C. Gardiner, _Stochastic methods_ , vol. 4, Springer Berlin (2009).
  * [33] R.D.T. Strauss and F. Effenberger, _A Hitch-hiker’s Guide to Stochastic Differential Equations. Solution Methods for Energetic Particle Transport in Space Physics and Astrophysics_ , _SSRv_ 212 (2017) 151 [1703.06192].
  * [34] A. Achterberg and K.M. Schure, _A more accurate numerical scheme for diffusive shock acceleration_ , _MNRAS_ 411 (2011) 2628 [1103.3049].
  * [35] M. Zhang, _Calculation of Diffusive Shock Acceleration of Charged Particles by Skew Brownian Motion_ , _ApJ_ 541 (2000) 428.
  * [36] J.A. le Roux, M.S. Potgieter and V.S. Ptuskin, _A transport model for the diffusive shock acceleration and modulation of anomalous cosmic rays in the heliosphere_ , _J. Geophys. Res._ 101 (1996) 4791.
  * [37] W.M. Kruells and A. Achterberg, _Computation of cosmic-ray acceleration by Ito’s stochastic differential equations_ , _A &A_ 286 (1994) 314.
  * [38] L.O. Drury, _An introduction to the theory of diffusive shock acceleration of energetic particles in tenuous plasmas_ , _Rep. Prog. Phys_ 46 (1983) 973.
  * [39] I.N. Toptygin, _Acceleration of particles by shocks in a cosmic plasma._ , _SSRv_ 26 (1980) 157.
  * [40] V. Florinski and J.R. Jokipii, _Cosmic-Ray Spectra at Spherical Termination Shocks_ , _ApJ_ 591 (2003) 454.
  * [41] T. Winchen and S. Buitink, _Energy spectrum of fast second order Fermi accelerators as sources of ultra-high-energy cosmic rays_ , _Astropart. Phys._ 102 (2018) 25 [1612.03675].
  * [42] R. Jansson and G.R. Farrar, _A New Model of the Galactic Magnetic Field_ , _ApJ_ 757 (2012) 14 [1204.3662].
  * [43] J. Kleimann, T. Schorlepp, L. Merten and J. Becker Tjus, _Solenoidal Improvements for the JF12 Galactic Magnetic Field Model_ , _ApJ_ 877 (2019) 76 [1809.07528].
  * [44] P. Terral and K. Ferrière, _Constraints from Faraday rotation on the magnetic field structure in the Galactic halo_ , _A &A_ 600 (2017) A29 [1611.10222].