# Automated detection and mapping of crystal tilt using thermal diffuse
scattering in transmission electron microscopy

Mauricio Cattaneo  Knut Müller-Caspary  Juri Barthel  Katherine E. MacArthur
Nicolas Gauquelin  Marta Lipinska-Chwalek  Johan Verbeeck  Leslie J. Allen
Rafal E. Dunin-Borkowski  Ernst Ruska Centre (ER-C), Forschungszentrum Jülich
GmbH, 52425 Jülich, Germany  Department of Chemistry and Centre for
NanoScience, Ludwig-Maximilians-University Munich, Butenandtstr. 11, 81377
Munich, Germany  Electron Microscopy for Materials Science (EMAT), University
of Antwerp, 2020 Antwerp, Belgium  School of Physics, University of Melbourne,
Parkville, Victoria 3010, Australia  Central Facility for Electron Microscopy,
RWTH Aachen University, 52074 Aachen, Germany

###### Abstract

Quantitative interpretation of transmission electron microscopy (TEM) data of
crystalline specimens often requires the accurate knowledge of the local
crystal orientation. A method is presented which exploits momentum-resolved
scanning TEM (STEM) data to determine the local mistilt from a major zone
axis. It is based on a geometric analysis of Kikuchi bands within a single
diffraction pattern, yielding the centre of the Laue circle. Whereas the
approach is not limited to convergent illumination, it is here developed using
unit-cell averaged diffraction patterns corresponding to high-resolution STEM
settings. In simulation studies, an accuracy of approximately 0.1 mrad is
found. The method is implemented in automated software and applied to
crystallographic tilt and in-plane rotation mapping in two experimental cases.
In particular, orientation maps of high-Mn steel and an epitaxially grown
La${}_{\text{0.7}}$Sr${}_{\text{0.3}}$MnO${}_{\text{3}}$-SrTiO${}_{\text{3}}$
interface are presented. The results confirm the estimates of the simulation
study and indicate that tilt mapping can be performed consistently over a wide
field of view with diameters well above 100 nm at unit cell real space
sampling.

††journal: Ultramicroscopy

##  1 Introduction

Scanning transmission electron microscopy (STEM) has established itself as a
method for structural characterization as well as chemical composition mapping
with atomic-scale resolution. A diversity of imaging techniques has been
developed, such as annular dark field (ADF) imaging providing $Z$-contrast
[1], bright field (BF) and annular bright field (ABF) imaging of light
elements [2], elemental mapping using electron energy-loss spectroscopy (EELS)
[3] and energy-dispersive x-ray spectroscopy (EDX) [4].

Common to these methods is a signal dependence on the crystallographic
orientation of the specimen with respect to the incident electron beam. Even
relatively small mistilts away from a crystal zone axis can affect the
dynamical scattering of the electron wave significantly [5, 6, 7].
Quantitative conclusions can often be drawn only via the comparison of
experimental data with accompanying simulations, which require accurate tilt
information as input to model scattering dynamics realistically. In the case
of EDX spectroscopy, for example, a deliberate sample tilt is applied in order
to suppress electron channelling [8, 9]. Furthermore, structure retrieval by,
e.g., inverse multislice [10, 11, 12, 13] can benefit either from the
knowledge of local tilt as initial input, or as a sanity check of the
reconstruction result.

Furthermore, research partly focuses on the local tilt distribution itself.
Examples include materials where local variations of the orientation due to,
e.g., stacking faults, interfaces, dislocations and grain boundaries are
related to materials properties such as elasticity or ferroelectricity [14,
15, 16]. A technique to measure and map the local sample tilt would thus be
beneficial for such applications, in particular if the data providing tilt
information can be acquired simultaneously with, e.g., spectroscopic or
Z-contrast data. In this regard, the advent of momentum resolution in STEM
using ultra-fast cameras [17, 18, 19, 20, 21, 22] has dramatically expanded
the versatility and flexibility of STEM, as complete diffraction patterns can
be acquired at each scan position up to a certain spatial frequency. In that
respect, four-dimensional STEM (4D-STEM) data sets have already been used to
determine specimen tilt and thickness in an automated way by means of neural
networks analyzing the elastic scattering signal [23]. Moreover, a correlative
approach comparing mainly the elastic scattering signal in experimental and
simulated position-averaged convergent beam electron diffraction (PACBED)
patterns has also shown a high sensitivity to sample mistilt [24]. Because
these approaches are crystal and zone axis specific, they require prior
knowledge. In the method presented here, the amount of needed prior
information is significantly reduced while maintaining a high evaluation speed
and accuracy.

At least for specimens investigated at sufficiently large thicknesses and
temperatures, typically above 10 nm at 300 K, thermal diffuse scattering (TDS)
constitutes a considerable fraction of the dark-field intensity. Besides the
smooth diffuse background, Kikuchi bands can be observed in diffraction
patterns of crystals at high scattering angles [25]. Since Kikuchi bands are
aligned with reciprocal lattice planes, their exploitation to measure crystal
tilt is widespread in electron diffraction. For example, crystallographic
orientation mapping via Kikuchi bands has a long tradition in scanning
electron microscopy (SEM) using electron back-scatter diffraction (EBSD) [26,
27, 28]. The kernel of the existing Kikuchi band detection schemes in EBSD
usually involves a Hough transform to locate the Kikuchi band crossings.
However, Hough-based methods can become imprecise in cases where rectangular-
shaped Kikuchi bands do not constitute the dominant features in diffraction
patterns. A strong zero-order Laue zone in diffraction patterns recorded in
transmission significantly complicates the application of Hough transforms,
especially at high beam energies in the (S)TEM.

In this work, we present an alternative approach based on Kikuchi band
recognition from azimuthal intensity profiles in an annular dark-field region
of a diffraction pattern. While requiring only a minimum of input, such as the
expected number of strong Kikuchi bands, the method performs very robustly in
typical high-resolution STEM scenarios in, or close to, zone-axis orientation
and allows us to quantify mistilt with well below 1 mrad accuracy and
precision. Whereas, methodologically, each diffraction pattern is treated
independently, the application of this scheme also provides the ability to map
local tilts in large 4D STEM data sets. The procedure was implemented in a
computer program to run fully automatically and computationally efficiently,
such that the detection of the centre of the Laue circle takes on the order of
several milliseconds per diffraction pattern on current personal computer
hardware.

In the following, the method to measure the tilt between a crystal axis and
the incident beam direction based on single diffraction patterns is introduced
conceptually. Major steps of the analysis include the enhancement of Kikuchi
band contrast in the dark-field region using a renormalization scheme, the
identification of Kikuchi band crossings from azimuthal profiles, as well as
iteratively eliminating systematic errors arising from the initially unknown
centre of the Laue circle. Frozen phonon multislice simulation studies are
performed in order to assess accuracy and precision of the method in
dependence on electron dose and mistilt. The method is then applied to two
experimental cases demonstrating tilt mapping (1) across a large field of view
spanning more than 100 nm in a bent cross-section containing an
La${}_{\text{0.7}}$Sr${}_{\text{0.3}}$MnO${}_{\text{3}}$/SrTiO3 interface, and
(2) to a smaller scale of approximately 10 nm across a twin boundary in high-
Mn steel. For both examples, 4D-STEM data sets were recorded with atomic
column resolution. Before the tilt analysis, the diffraction patterns were
averaged over scan areas corresponding to the size of one projected unit cell.
Mistilt maps are thus generated from PACBED patterns with unit cell sampling.
Due to the robustness and accuracy of the method, small changes of sample tilt
are resolved over a large field of view as well as with high sampling close to
interfaces and grain boundaries.

##  2 Methods

###  2.1 Conceptual outline and problem statement

The measurement of sample mistilt from a zone axis orientation based on
Kikuchi bands, as proposed here, consists of locating the crossing point of
the bands by geometrical construction and determining its distance from the
center of the bright-field disc. In principle, this analysis could be
performed exploiting the full information on tilt present in a two-dimensional
diffraction pattern, e.g., by cross-correlation against an expected band
pattern. However, the intensity distribution and details like the shape and
width of Kikuchi bands depend on the crystal structure, sample thickness,
zone-axis orientation, and further experimental parameters. This means a full
two-dimensional analysis of the Kikuchi-band pattern obtained in transmission
electron diffraction which would require extensive information as input to an
automated routine. Since we are aiming for a fast method which requires only a
minimum of prior knowledge, the evaluation approach is simplified by
determining points in the diffraction plane through which Kikuchi bands pass.

Figure 1: (a) Diffraction pattern of a steel sample recorded with 300 keV
electrons and at 7.2 mrad mistilt from the [110] zone axis. Kikuchi bands
intersect in point $\bf{P}$ whereas the incident beam direction corresponds to
point $\bf{C}$. The logarithm of intensity is plotted to enhance the
visibility of Kikuchi bands. (b) Re-normalized signal of the pattern in (a)
within an annular region that is centered around the point $\bf{Q}=\bf{C}$.
(c) Diffraction pattern of a SrTiO3 sample recorded with 300-keV electrons at
a large mistilt of about 75 mrad from the [100] zone axis. (d) Signal of the
pattern in (c) within an annular region around $\bf{Q}$ in the detector
center. Arrows indicate two traces of one Kikuchi band in the annular mask.

In a preliminary step, an annular mask is applied to the diffraction pattern
so as to separate the data including thermal diffuse scattering (TDS) and thus
Kikuchi bands at medium to high angles from the bright field and low-angle
diffraction. An example with a relatively small mistilt of 7.2 mrad is shown
in Fig. 1 (a,b) for the full and masked experimental diffraction pattern,
respectively. In the following, we denote the locations of the undiffracted
beam by $\mathbf{C}$, the centre of the Laue circle by $\mathbf{P}$ and the
centre of the annular mask by $\mathbf{Q}$. As indicated by the arrows in Fig.
1 (b), Kikuchi bands passing through the mask leave two traces at
approximately opposite sides of the annulus with azimuths differing by
approximately $\pi$. Figure 1 (c,d) shows an example with a comparably large
mistilt of 75 mrad, and a ring mask that is neither centered around
$\mathbf{C}$ nor $\mathbf{P}$. In such cases, the azimuth difference
corresponding to a single Kikuchi band can deviate significantly from $\pi$.
Moreover, the centering of the mask away from $\mathbf{C}$ leads to a strongly
varying mean scattered intensity in dependence of the azimuth $\varphi$, and
the traces of two bands can become very close and might even overlap in some
cases. In addition, there can be a much higher intensity close to the inner
rim of the mask.

These circumstances complicate an automatic and robust detection of Kikuchi
bands by determining the azimuthal positions of their traces in the annular
mask. Furthermore, the signal itself decreases along the radial direction,
i.e. towards larger scattering angles where it approaches a Rutherford type of
decay of the scattering cross-section [29]. This would inherently increase the
weight of scattering close to the inner rim of the annulus if not corrected
for during the evaluation. In that respect, it is noteworthy that the
azimuthal resolution increases with the distance from the mask center. In the
following, solutions to these problems are described which provide an
enhancement of the Kikuchi-band signal, facilitate locating band positions
within the annular region, and finally enable the determination of the
position of the common band crossing, that is, the centre of the Laue circle.

###  2.2 Enhancing Kikuchi band contrast

Since Kikuchi bands are not a dominant feature in the bright field and also in
the low-angle dark-field region, a short camera length is beneficial for
accessing large scattering angles such that the Bragg reflections within the
zero-order Laue zone in the center cover less than a quarter of the recorded
diffraction pattern, as shown by the examples in Fig. 1. For a given camera
length and diffraction projection on the pixelated detector, we assume that
both Kikuchi-band crossing $\mathbf{P}$ and bright field $\mathbf{C}$ are
inside the inner radius of the annular mask. Consequently, the geometrical
approach developed in the following is restricted to tilt variations limited
by the inner radius of the annular mask at the present stage. However, since
this already covers several tens of milliradians tilt variation, it is
applicable to the majority of practically relevant cases. The outer radius of
the mask is made as large as possible reaching out to the edge of the detector
area.

Before an azimuthal scan is extracted from the area selected by the annular
mask, the mask is partitioned in concentric rings of widths between one and
two detector pixels. For each thin annular partition, the mean value of
intensity is subtracted and the remaining signal variation is normalized with
respect to the standard deviation. Each partition is then of equal strength in
signal variation. For the subsequent integration of the signal along the
radial direction, this renormalization avoids stronger weight of inner parts
of the annular mask compared to outer parts due to the decay of the electron
scattering cross-section towards larger angles. The result of renormalization
within the annular mask is shown in Figs. 1 (b) and (d), where the Kikuchi
bands show only weak intensity variation along the radial direction, but a
strong variation occurs azimuthally.

###  2.3 Azimuthal scan and Kikuchi band detection

A one-dimensional azimuthal intensity profile depending on $\varphi$ is
extracted from the annular mask by integrating the renormalized data along the
radial direction. Sufficient azimuthal sampling is given, for example, by the
number of pixels crossed by the inner circle limiting the annular mask. For
each azimuthal sample the intensity is effectively integrated over a larger
azimuthal range (1) in order to suppress the fine structure of Kikuchi bands
which could possibly impede subsequent peak finding and (2) to further improve
signal-to-noise ratio. The size of these azimuthal sectors should be close to
the size of the Kikuchi bands within the annular mask. In that respect,
sectors spanning five to ten degrees were found to be good choices in all
cases analyzed so far.

Figure 2 (a) shows the azimuthal scan that has been directly extracted from
the pattern in Fig. 1 (b) in blue, using sectors with a width of six degrees.
The three major and three minor Kikuchi bands in this pattern add up to twelve
peaks in total. The azimuthal positions of the peaks are measured by detecting
the local maxima, and the mean radius of the annular mask is used as radial
position.

Figure 2: Azimuthal scan (a) as extracted from the [110] steel diffraction
pattern in Fig. 1(b) and (b) as extracted from the [100] SrTiO3 pattern in
Fig. 1(d). The azimuthal sampling was done in steps of one degree and binning
into sectors spanning six degrees (blue curves). Kikuchi bands appear as
narrow peaks in the profile. A strong background modulation is present because
parts of the azimuthal scan are closer to the bright field than others. This
background modulation is removed in the red curves by means of a high-pass
filter. Black circles indicate the result of local peak detection.

The azimuthal profile in Fig. 2 (b) has been extracted directly from the
example with larger mistilt in Fig. 1(d). It shows a strong one-fold, low-
frequency modulation with an amplitude similar to that of the set of peaks
caused by Kikuchi bands which impedes robust automatic peak detection and
stems from the large deviation of $\mathbf{Q}$ from $\mathbf{C}$. Such a low-
frequency modulation is also present in the blue curve of Fig. 2(a), although
with a much lower amplitude. A higher intensity occurs for azimuths where the
mask is closer to the bright field and covers lower scattering angles. Besides
difficulties in automated peak detection, this effect causes a significant
shift of peaks located at strong slopes of the background. In order to prevent
systematic errors arising from this effect, a high-pass filter is applied to
the azimuthal scan which effectively removes the low-frequency background.
This high-pass filter sets the lowest-order Fourier coefficients to zero
resulting in the red profiles in Fig. 2, which are now dominated by the peaks
caused by Kikuchi bands in terms of signal variation.

The black circles in Fig. 2(a) indicate the positions of the twelve strongest
local maxima of the red curve expected for six Kikuchi bands in the case of a
steel sample investigated close to [110] orientation. Likewise, the eight
highest local maxima occurring in the red curve of Fig. 2(b) mark the
locations expected from four strong Kikuchi bands for [100] SrTiO3. In the
example of Fig. 2 (a), peaks have comparable heights, stand out significantly
from the background and can hence be detected straight forwardly by
established methods. However, the situation is different for the high-pass
filtered azimuthal profile in Fig. 2 (b) which shows a large number of local
maxima at azimuths not related to the centers of Kikuchi bands. A
classification of a local maximum as a Kikuchi band or as a peak caused by
noise related to, e.g., low counting statistics, is not obvious even for a
human observer at least in some cases. For example, setting a threshold based
on peak height, width, area or any other quantity regarded as peak strength,
may easily reject actual Kikuchi-bands while accepting noise-related peaks.
Therefore, the detection routine is allowed to find more than the expected
number of local maxima and registers them in a sorted list by means of a
parameter of strength, e.g. by the peak height. This approach allows for the
identification of falsely detected peaks in the next step, based on a
criterion that leads to the expected geometrical situation of Kikuchi bands,
namely that they are crossing in one point of the diffraction pattern.
Finally, note the regular spacing of peaks in Fig. 2 a due to the nearly
coinciding centres $\mathbf{P}$ and $\mathbf{Q}$ in Fig.1 b, as well as the
distorted spacings of maxima on the $\varphi$ axis in Fig. 2 b where the mask
centering deviates significantly from the centre of the Laue circle according
to Fig. 1 d.

###  2.4 Determination of the Kikuchi band intersection

For $n$ Kikuchi bands in a diffraction pattern the annular mask is ideally
crossed twice by each band, generating $2n$ peaks in the azimuthal scan. In
such cases, a straightforward approach to associate peak pairs to Kikuchi
bands, is to connect one peak with the $n^{\mathrm{th}}$ peak in the azimuthal
sequence by a straight line. Thus $n$ lines are constructed which have in
total

$\eta=\frac{n(n-1)}{2}$ (1)

crossings, of which all coincide at or close to the centre of the Laue circle
in the ideal situation, as shown by the blue markers in Fig. 3. The procedure
to calculate the crossing point for two pairs of different azimuths is
detailed in A.

This procedure relies entirely on the correct detection of all peaks and, in
particular, of the complete set of pairs of peaks related to the Kikuchi
bands. The association of peaks to pairs based on the azimuthal sequence fails
to produce a set of crossing points forming a compact cloud, if Kikuchi bands
are only partially detected. For example, this can happen if only one of the
two intersections of a band with the annular mask is found, if two bands
overlap and produce only one instead of two peaks, or when a false recognition
occurs due to detecting peaks arising from noise. An example for a false
relation of points to Kikuchi bands is shown by the red markers in Fig. 3,
where only one of the two intersections of weak bands with the annulus were
included. This leads to a false sequential association of pairs and thus to a
wide scatter of crossing points quite far away from the common crossing of
Kikuchi bands. The residual scatter of the crossing points for a correct
selection of peaks is smaller and expected to be in the range of the width of
Kikuchi bands.

Figure 3: (a) PACBED pattern of [110] fcc steel with larger mistilt and (b)
azimuthal scan demonstrating an erroneous (red) and a correct (blue)
association of detected peaks to pairs, which are assumed to represent the two
traces of Kikuchi bands in the annular mask. The pairs deduced from the list
of peaks marked by red triangles lead to the wrong result because the second
trace of two Kikuchi bands is not included. The resulting set of crossing
points shows a large scatter and is far away from the crossing point of the
Kikuchi bands. The list of peaks marked by blue squares leads to a set of
crossing points with low scatter at the crossing point of the Kikuchi bands.

In most cases, the set of the $2n$ strongest peaks found in the azimuthal scan
provides also the best solution, especially when there are $n$ Kikuchi bands
that are clearly visible and sample mistilt is low. In order to stabilize the
solution for more challenging scenarios, as the one shown in Fig. 3, more than
$2n$ peaks are taken into consideration. In particular, the list of peaks
found in the azimuthal scan is sorted with respect to the peak height, and
extended to include also weaker ones which are only potentially related to
Kikuchi bands. Several sets of $2n$ peaks are created from this extended list
by randomly exchanging mostly weaker peaks. From these sets, the one leading
to the lowest scatter of crossing points is chosen as the solution to identify
the location of the Kikuchi-band crossing.

###  2.5 Iterative error correction

On startup of the evaluation procedure, the position $\bf P$ of the centre of
the Laue circle in a diffraction pattern is unknown. Especially in the case of
large sample mistilt, i.e. for a large distance between points $\bf P$ and
$\bf C$ in Fig. 1, the annular mask applied in the previous steps will usually
be centered around a point $\bf Q$ in the diffraction pattern that may be
rather remote from the initially unknown point $\bf P$. Consequences are that
(1) Kikuchi bands are more difficult to resolve in parts of the mask close to
$\bf P$ and (2) some bands cause peaks with an asymmetric shape in the
azimuthal profile. The first issue can lead to a serious lack of detected
pairs, whereas the second one causes a small systematic shift of the result
towards the center of the mask. The two problems are illustrated schematically
in Fig. 4 with two bands and a significant offset of their crossing from the
center position of the annular mask. Errors due to these issues can be quite
significant, i.e. larger than the scatter of crossing points found in the
previous step.

Figure 4: Illustration of systematic errors in Kikuchi band detection due to
an off-centered annular mask or significant sample mistilt. (a) Situation
showing two bands crossing under a small angle close to the edge of the
annular mask, which makes it difficult to resolve their intersections with the
annulus in region “S1”, whereas two well separated band sections occur in the
opposite region “S2”. (b) Situation with two perpendicular Kikuchi bands,
where the band labelled “KB 1” is shifted off the mask center along band “KB
2”. Symmetric peaks are registered in azimuthal sectors for “KB 2”, which
crosses the annular mask along the radial direction. In contrast, the peaks
due to “KB 1” are asymmetric, which leads to a systematic shift of the deduced
crossing point.

The two problems described above are minimum if the annular mask is centered
at the crossing of Kikuchi bands. Therefore, an iterative scheme is applied in
which the annular mask is centered to the result of the previous evaluation
and a new evaluation is performed. In most cases only one or two iterations
are needed until the correction of the mask position is smaller than the
scatter of crossing points.

###  2.6 Software implementation

The evaluation procedure described above has been implemented in a computer
program using the C++ language. With commercial hardware, each iteration of
finding a Kikuchi intersection for a single PACBED takes on the order of ten
milliseconds using a single CPU thread. Parallel programming via OpenMP is
applied to further speed up the analysis making use of the fact that the
evaluation of multiple diffraction patterns can be performed independently. A
set of 1000 diffraction patterns is evaluated by four threads in less than 25
seconds.

The program can take one or multiple diffraction patterns as input in form of
two-dimensional arrays of intensity values together with parameters describing
the size of the array. The inner radius limiting the annular mask and the
number of expected Kikuchi bands $n$ are the required input ensuring stability
of the procedure. Defining the outer radius of the annular mask is optional
but limiting the width of the annulus has shown improved performance in some
cases. Output is the position of the Kikuchi band intersection $P_{N}$ in the
input array as determined after $N$ iterations and an error estimate for the
position corresponding to the root-mean-square deviation of all crossing
points of detected Kikuchi bands. The mistilt is then calculated as the
distance of the output position $P_{N}$ to the known position $C$ of the
incident probe direction considering scaling parameters which describe the
diffraction projection in physical units, e.g. in milliradians per pixel.

##  3 Simulation studies

In order to assess the accuracy and precision of the mistilt measurement,
diffraction patterns were simulated using the multislice method [30] as
implemented in the Dr. Probe software [31] for a large range of mistilts.
Thermal diffuse scattering was simulated according to the quantum excitation
of phonons (QEP) model [32] assuming an Einstein model of uncorrelated thermal
vibration. As a test case, diffraction patterns were calculated for 40 nm
thick fcc Fe viewed along the $[110]$ zone axis based on the model of
Straumanis & Kim [33]. An aberration-free incident electron probe of 200 keV
kinetic energy and 25 mrad convergence semi-angle was assumed. PACBED patterns
were obtained by summation of patterns resulting from more than 100 probe
positions distributed uniformly over the projected unit cell. The diffraction
patterns extend out to scattering angles of about 250 mrad due to a fine
sampling of the supercell of $5\times 7$ unit cells with $800\times 800$
points perpendicular to the incident beam direction. The grid $x$-axis was
along the $[\bar{1}10]$ direction and its $y$-axis along $[001]$. This setup
also realises a step size of approximately 1 mrad in the diffraction patterns,
which is on the order of magnitude that can be achieved with fast pixelated
direct electron detectors and convenient camera lengths in experiments. The
potentials were projected along $[110]$ within volume slices that are thin
enough to contain just one atomic layer. For each layer 100 positional
configurations were calculated by including random displacements with mean
squared amplitudes of $\langle u^{2}\rangle=0.007\,\text{\AA}^{2}$[34]. By
random selection of positional configurations from the pre-calculated set for
each slice during the propagation of the probe function through a thick
crystal in the multislice scheme, the scattering for each probe position is
essentially calculated with an individual positional configuration. The
summation of patterns from many probe positions form a PACBED pattern which
thus also includes a summation over positional configurations achieving
convergence for the signal related to thermal diffuse scattering. A pattern
calculated in this way shows three strong Kikuchi bands along the main
reciprocal lattice plains and weaker bands in between.

Test mistilts were distributed uniformly over a range up to 40 mrad tilt
magnitude. Patterns obtained directly from the simulations correspond to
infinite dose conditions. Scenarios representing a finite dose were simulated
by adding noise following Poisson statistics. The total dose assumed in the
simulation of one PACBED pattern translates to incident probe currents where
we assume here a dwell time of 1 ms per probe position and 100 probe positions
contributing to each simulated PACBED.

Figure 5: Histograms of the overall radial error of mistilt measurements with
simulated patterns for different probe currents. (a) Detection error for a
direct evaluation, without iterative refinement. A systematic error of
approximately 0.1 mrad is observed. (b) Improved accuracy after iterative
refinement showing no significant systematic error.

The histograms in Fig. 5 show the radial component of the mistilt detection
error (a) without and (b) with iterating the tilt measurement using the result
of one cycle as a start for the next. Without iteration a systematic error of
approximately -0.2 mrad occurs, which might already be tolerable for many
applications. The measurement becomes perfectly accurate if the iterative
scheme is applied, since the distribution is symmetric around zero. This
behaviour of the automated evaluation procedure is independent of noise added
due to a finite electron dose. Also the width of the distributions plotted in
Fig. 5 shows no strong dependence on electron counting noise but is reduced by
the iterative scheme. A precision of about $\pm 0.13$ mrad is deduced from the
range containing 68% of the test cases, corresponding to the 1-$\sigma$ level
of a normal distribution. A slight decrease of the precision is found for the
lowest probe current of 10 pA assumed in the simulations. Reduction of the
incident probe current by another order of magnitude leads to a failure of the
evaluation in over 50% of the cases under the applied conditions of the test,
since the low signal of the Kikuchi bands could not be distinguished from
noise anymore. It should be noted though that the probe current given in this
test is not a fundamental limit because the signal-to-noise ratio of the
Kikuchi bands also depends on the detection dwell-time, the size and
efficiency of detector pixels, as well as the number of probe positions
included in a PACBED. This means that the automated tilt measurement can also
be applied in cases with lower probe currents or in low-dose experiments by
averaging over a larger area of the scan field and by applying pattern
binning.

##  4 Experimental results

###  4.1 Data acquisition

Experimental data has been acquired with an FEI Titan 80-300 [35] and an FEI
Titan G2 80-200 [36] electron microscope at accelerating voltages of 300 and
200 kV, respectively. The probe convergence semi-angle was limited to 25 mrad
by an aperture and probe aberrations were minimized within this angular range
using a double-hexapole corrector providing below 0.1 nm spatial resolution.

Momentum-resolved STEM data of a
La${}_{\text{0.7}}$Sr${}_{\text{0.3}}$MnO${}_{\text{3}}$ epitaxial layer grown
on an SrTiO3 substrate was recorded with a Merlin for EM Medipix 3 single chip
detector (Quantum Detectors, Oxfordshire, UK) with a $256\times 256$ pixel
array [18, 37]. The TEM sample was prepared by focussed ion beam (FIB) milling
in cross-section geometry. A second experimental case study of high-Mn steel
has been performed using an electron microscope pixel array detector (EMPAD)
(Thermo Fisher Scientific) with $128\times 128$ pixels [21]. The TEM sample
was prepared by a lift-out method in a dual-beam focused ion beam (FIB) FEI
Helios NanoLab 460F1 [38].

In order to demonstrate the applicability of the method in practice, we
present two case studies in the following. Firstly, it is shown that specimen
bending occurring in thin foils of a perovskite heterostructure can be
quantified at large fields of view (FOV). This example also shows that
atomically resolved 4D-STEM with $10^{6}$ scan pixels is feasible despite a
slow frame rate in the kHz range, i.e. for a total scan time on the order of
1000 seconds. Secondly, we consider a perfect twin boundary in a high-Mn steel
specimen. Since this interface is purely crystallographic, it exhibits no
compositional or strain gradients but an in-plane rotation between the
adjacent grains. Observations made in this example provide information about
the spatial resolution of the method.

###  4.2 Bending of specimens at large fields of view

Figure 6(a) shows the ADF signal obtained from 4D-STEM data at a scan raster
of $400\times 3500$, which corresponds to a field of view of 17 nm $\times$
146 nm. The left-hand side contains an epitaxial
La${}_{\text{0.7}}$Sr${}_{\text{0.3}}$MnO${}_{\text{3}}$ layer (bright), grown
on an SrTiO3 substrate on the right hand side [39]. Subsets of the data are
shown in Fig. 6(b) for five substrate regions A-E labelled in (a), being
approximately 30 nm apart from each other.

Figure 6: Momentum-resolved STEM acquisition of epitaxial
La${}_{\text{0.7}}$Sr${}_{\text{0.3}}$MnO${}_{\text{3}}$ (bright layer, left)
on SrTiO${}_{\text{3}}$ (right). (a) Virtual detector ADF image of a $17\times
146\,$nm scan. (b) Five selected areas from (a) with real space scan images
and corresponding PACBED patterns. Green and yellow dots in the PACBEDs mark
detected Kikuchi bands and deduced band crossings, respectively. The r.m.s.
scatter $\sigma$ of crossing points is given. (c) Map of the measured out-of-
plane mistilt. (d) Map of estimated error of the tilt measurement.

Although atomic resolution is present in all ADF images, their contrast varies
significantly and is maximum in regions C and D, which correspond to an
approximate zone axis condition. A selection of PACBED patterns from different
regions of the scan is displayed in Fig. 6 (b). In particular, a variation of
crystal tilt is identified therein by the apparent vertical shift of the
crossing of the Kikuchi bands when comparing the patterns in sequence from
region A to region E. In the example patterns shown, green dots mark the
Kikuchi bands detected by the automated mistilt measurement. The measurement
routine provides an internal error estimate $\sigma$ which corresponds to the
root-mean-squared (r.m.s.) deviation of line crossings found in the analysis,
e.g., as shown in Fig. 3. Error estimates in this experimental study are well
below one milliradian, whereas the mistilt magnitude from $[100]$ zone axis
reaches up to $9$ mrad in some regions.

From the scan a set of $40\times 350$ PACBED patterns was generated, where
each PACBED is calculated from an area corresponding to one projected pseudo-
cubic unit of the perovskite. The application of the iterative evaluation
procedure took 20 minutes in total on a consumer-grade computer. The result of
the systematic mistilt analysis for the whole region of Fig. 6 (a) is shown in
Fig. 6 (c). A colour scale is used to display the modulus and direction of
mistilt in terms of brightness and hue, respectively. The change from blue to
green along the vertical scanning direction corresponds to a change of the
vertical tilt component. This dominant variation of tilt is approximately a
linear function of the scan position and corresponds to a bending of the
specimen around a horizontal axis by approximately $\pm 9\,$mrad, as plotted
in Fig. 7.

Although the global trend of the tilt variation is linear, significant local
deviations from this trend are observed in the epitaxial layer. These areas
also appear with higher brightness in the ADF image of Fig. 6 (a), which is
indicative for enhanced static disorder due to defects. Additional diffuse
intensity in the diffraction pattern due to Huang scattering [40] associated
with disorder hampers Kikuchi band detection by lowering their contrast. A map
of the measurement error estimated by the automated procedure is plotted in
Fig. 6 (d). A good precision of 0.3 to 0.6 mrad is obtained in the substrate
while the majority of measurements in the epitaxial layer has an error below 2
mrad. Significantly larger errors are found in some areas of the epitaxial
layer, especially at the bottom of the map, where also the measured tilts
scatter strongly. This indicates that the algorithm could not converge in
these cases within a predefined limit of 30 iterations, most likely due to a
low contrast of the Kikuchi-band pattern.

Figure 7: Line profiles of sample mistilt along the arrow in Fig. 6 in the
epitaxial layer (red) and the substrate (blue). (a) Mistilt $x$-component,
perpendicular to the interface, and (b) $y$-component, parallel to the
interface, show a different variation with sample position. Whereas there is
almost no change of the $x$-component, the $y$ component exhibits a linear
trend, indicating sample bending of 20 mrad over 140 nm. Oscillations of the
sample mistilt around the linear trend are observed in the epitaxial layer for
both components.

It is instructive to interpret the obtained tilt map in terms of strain-
induced tilt variations caused by the epitaxial layer grown with a lattice
mismatch of 0.8 %. For this purpose, we calculated a line profile of the
mistilt along the direction indicated by the green arrow in Fig. 6(a)
separately for the substrate and epitaxial layer. These line profiles are
plotted in Fig. 7 for (a) the component of tilt perpendicular to the interface
and (b) the component of tilt parallel to the interface. The blue profile was
measured from substrate areas and shows a constant bending of the lamella
along the interface with a tilt gradient of 143 $\mu$rad/nm and a standard
deviation of 1 mrad, as indicated by error bars. In contrast to that, the red
curves measured on the layer material show a wave-like deviation from the
linear trend of up to 7 mrad, which are likely caused by thin foil relaxation
and defects. Larger error bars of the layer measurement reflect both, the
slightly worse precision of tilt detection and an increased tilt variation in
the layer material also along the direction perpendicular to the interface.

###  4.3 Tilt measurement in high-Mn steel

The second case study is a tilt mapping across a twin boundary in austenitic
steel. The results in Fig. 8 show the mapping of a perfect $\Sigma$3 twin
boundary viewed along the common [110] direction of the two grains with the
twin boundary in a {111} plane and oriented horizontally in Fig. 8(a). From
the methodological point of view, this structure is interesting since it
provides an interface of a perfect twin which is well understood, however,
without adding complexity due to compositional gradients or strain. Besides
mapping the tilt induced by bending of the thinned TEM sample, as in the
previous case, this second case provides additional insight regarding the
attainable spatial resolution of the tilt measurement, because the
crystallographic orientation changes abruptly at the twin boundary. The change
of orientation is apparent in the PACBED patterns as an in-plane mirror with
respect to the [111] interface plane, which is marked by the dashed line in
Fig. 8(b).

Figure 8: Mistilt mapping across a perfect twin boundary in high-Mg steel. (a)
Virtual ADF image calculated from the 4D-STEM data set. Two grains with
different in-plane orientation are separated by a horizontal grain boundary.
(b) Example PACBED patterns of three selected areas from (a) with positions A
and C taken on either side of the boundary and B taken from an area that
includes the boundary. (c) Out-of-plane orientation map showing a small sample
bending of approximately 2 mrad over the entire field of view.

Figure 8(c) shows the orientation map in a scan of 14 nm $\times$ 14 nm. The
sample appears mostly flat with tilt variations of 1 to 2 mrad. Although these
orientation changes are small, they are smooth across the twin boundary and
measured with very good precision of $\sigma<0.1$ mrad. Above this level of
uncertainty, no abrupt change of the out-of-plane crystallographic orientation
was found at the grain boundary. This indicates that the two grains are indeed
perfect crystallization twins.

Pattern B displays an obvious problem for the procedure of tilt measurement at
grain boundaries with significant orientation changes. There are partial bands
on either side of the {111} plane, i.e. no counterpart is found on the
opposite side to support a pair-wise relation of peaks in the azimuthal scan
to a Kikuchi band. This causes the measurement routine to fail, which is
reflected by the white area at the twin boundary in Fig. 8(c). The width of
the white area corresponds to that of one unit cell on either side of the twin
boundary. Already at a distance of two unit cells from the boundary, the
routine works well with good accuracy and precision, i.e. the mistilt at these
positions doesn’t deviate much from the results obtained deeper in the grains,
and the error estimates are likewise small. Accordingly, we deduce that the
result in one grain is not influenced by the difference in orientation of the
other grain at a distance of two unit cells, which could be interpreted as a
spatial resolution of better than 0.8 nm in this case.

##  5 Discussion

A geometric approach to measure crystal orientation from diffraction patterns
has been presented. By choosing Kikuchi bands as the source of the signal for
orientation mistilt measurement, the evaluation procedure requires that
Kikuchi bands are separated from other much stronger features in electron
diffraction patterns, such as the bright field and Bragg reflections. A
sufficiently large sample thickness is required, in most cases more than 10
nm, to generate significant Kikuchi bands at room temperature. The purpose of
the method is to measure mistilt in orientations close to a high symmetry zone
axis as frequently used in high-resolution STEM. In principle it is beneficial
for the geometric analysis that the diffraction pattern covers a large angular
range which can be achieved by choosing a short camera length. In our
experiments the method produced reliable results in cases where the bright-
field diameter is smaller than half the edge length of the detector. In
addition, shorter camera lengths allow for larger mistilts to be analyzed and
better precision. A compromise should be found when extending the angular
range of detection because the signal drops towards larger angles (decay of
scattering factors) and non-linear projection distortions can complicate
angular calibration. Even on simulated patterns, as well as under optimized
experimental conditions we found a limitation of precision of about 0.1 mrad
root-mean-square for the mistilt which is mainly due to the finite width of
the Kikuchi bands.

It is instructive to study whether Kikuchi band detection is facilitated by
heating specimens. Simulations show the disappearance of higher-order Laue-
zone (HOLZ) rings with increasing temperature, which may otherwise lead to
systematic errors or failures of the procedure. In particular, this would be
the case when at least one HOLZ ring is partially but not completely included
in the azimuthal scan. The diffuse scattering signal increases with increasing
sample temperature on the one hand. However, the simulations do not show a
likewise increased contrast of Kikuchi bands on the other hand. Apart from the
disappearance of the HOLZ signal, heating the sample does not seem to offer
significant advantages.

At lower beam energies and large mistilt Kikuchi bands can appear curved in
the diffraction pattern, since their edges follow conic sections in the strict
sense. As a consequence, the prediction of the crossing point of the bands
with the present method could become slightly inaccurate. In this case, a
final step evaluating the full two-dimensional course of the Kikuchi bands
could be necessary, e.g., assuming hyperbolic shapes. However, in the cases
presented here, sufficiently high beam energies have been used such that this
error is negligible.

A significant step in the pattern analysis is reducing a two-dimensional
problem, the intersection of bands, to a one-dimensional problem, measuring
peak positions in an azimuthal intensity profile. This is beneficial for the
speed and complexity of the numerical evaluation. However, a problem may arise
when associating pairs of peaks to bands. The smallest root-mean-square
distance of band intersections is used as the main criterion for finding the
correct association. There are rare cases, even under good experimental
conditions, where this criterion yields a deceptively good but actually wrong
result. An example is depicted in Fig. 9. Note that the computer only operates
based on the 1D scan and doesn’t perceive the 2D pattern like a human
observer. These initially wrong results are, however, healed by the refinement
step using another azimuthal scan extracted from a different area of the 2D
pattern. In all cases observed so far, this procedure converged to the correct
result.

Figure 9: Simulated PACBED pattern of Fe [110] fcc steel and example of a
false detection of the Kikuchi band intersection with a very small root-mean-
square deviation. This is a rare case, where by chance, the depicted wrong
association of pairs of peaks to bands leads to the smallest deviation of
common intersections. Given the underlying 2-dimensional pattern, the error is
obvious, but not detectable in the evaluation of a single azimuthal scan.
Using the iterative scheme, this false detection is nevertheless eliminated.

Applying the method to 4D-STEM data allows for mapping local mistilt. The two
experimental case studies presented in the results section both show
systematic variations in local orientation. Each pixel in the maps corresponds
to a scanned area from which diffraction patterns have been averaged and then
analyzed independently. The locality of tilt information, however, needs to be
interpreted with care due to the significant thickness- and orientation-
dependent broadening of the electron probe in the specimen. An empirical
assessment of the spatial resolution for mistilt variation is found in the
mapping of Fig. 8 for a perfect twin-boundary in steel. Already at a distance
of one unit cell from the twin boundary to either side, the measured tilts are
consistent with the rest of the grain, which suggests that the mistilt
information is generated from a projected area with a radius of approximately
one unit cell in the present case. For patterns generated from areas closer to
the grain boundary, e.g. pattern B in Fig. 8(b), the method fails for obvious
reasons. In this case all Kikuchi bands except one terminate at the pole and
do not generate the expected signal on the opposite side of the azimuthal
scan. In the case of high-resolution STEM of steel, a continuous Kikuchi band
along a direction is therefore generated as soon as there is at least one
atomic column from the scan position in both directions to form a crystal
plane. For the most intense Kikuchi bands, the neighboring columns have
distances of one or two unit cells. The spatial resolution for local tilt
mapping is therefore determined by the spacing of atomic columns in the
crystallographic planes corresponding to the Kikuchi bands used in the
analysis. Since all orientation mappings presented in this study were
generated from PACBED patterns of at least unit cell areas, tilt information
in one data point can be seen as independent from its neighbors.

##  6 Summary and Conclusion

To measure small tilts from the zone axis geometry, an approach based on the
signal from Kikuchi bands in diffraction patterns was chosen. This approach
requires a minimum of a-priori information and can thus be used with little
preliminary work. The need to record Kikuchi bands does not preclude
simultaneous evaluation of bright field or elastic diffraction signal as long
as both are recorded with sufficient area on the detector. This means that,
for example, with a momentum-resolved STEM measurement (4D-STEM), data
relevant to material science as well as supporting data important for later
modelling, such as tilt, can be recorded simultaneously.

The evaluation procedure has been implemented numerically and provides fast
and reliable results with accuracy on the order of 0.1 mrad mistilt. However,
the application of the method to typical (S)TEM diffraction patterns is
limited to small mistilts of less than 5° and requires a crystal thickness
sufficient to form Kikuchi bands. As a further development it could be
investigated whether the need for minimal inputs can be completely avoided,
e.g. by training a neural network with multiple patterns, to obtain a fully
automatic method.

Two examples have demonstrated that the method can be used to map sample
mistilts over large areas and interfaces. The resulting information allows for
the measurement of orientation changes across interfaces and to determine
grain boundary parameters. This is important in modelling, e.g., when
diffraction effects need to be separated from that of chemical gradients. The
high speed of the evaluation allows fast feedback for an online sample-tilt
optimization. This would be a step towards automating high resolution
transmission electron microscopy in finding zone-axis orientation or applying
deliberate mistilts for, e.g., quantitative EDX.

##  7 Acknowledgements

M.C. acknowledges funding from the Jülich Melbourne Postgraduate Academy
(JUMPA). K.M.-C. acknowledges funding from the Deutsche Forschungsgemeinschaft
(DFG) under grant number EXC 2089/1 - 390776260 (Germany’s Excellence
Strategy) and the Helmholtz association under contract VH-NG-1317. J.V. and
K.M.-C. received funding from the FWO under contract G042920N. M.L.-C
acknowledges financial support of the DFG within the Collaborative Research
Center (SFB) 761 “Steel – ab initio”. The collection of EMPAD data took place
within a cooperation with Thermo Fisher Scientific which is gratefully
acknowledged.

##  Appendix A Calculation of line crossings from positions in an azimuthal
scan

Given is a list of $2n$ unequal azimuthal angles $\varphi_{i}$ corresponding
to positions

${\bf X}_{i}={\bf Q}+r\begin{pmatrix}\cos\varphi_{i}\\\
\sin\varphi_{i}\end{pmatrix}$ (2)

on a circle of radius $r$ with center $\bf Q$ in a two-dimensional plane. The
list is sorted, such that $\varphi_{i}<\varphi_{i+1}$ for all $i\geq 1$ and
$i<2n.$ $n$ lines are constructed by connecting pairs of positions ${\bf
X}_{i}$ and ${\bf X}_{i+n}$ for $i\leq n.$ A line that contains two points
${\bf X}_{i}$ and ${\bf X}_{j}$ is given by the formula

${\bf L}_{ij}(t)={\bf X}_{i}+t({\bf X}_{j}-{\bf X}_{i}),t\in\mathbb{R}.$ (3)

The crossing point of two different lines ${\bf L}_{ij}(t_{1})$ and ${\bf
L}_{kl}(t_{2})$ with $k\neq i,$ $j=i+n,$ and $l=k+n$ is found by solving the
system of two equations given by

${\bf 0}={\bf L}_{ij}(t_{1})-{\bf L}_{kl}(t_{2})$ (4)

for the parameters $t_{1}$ and $t_{2}$ and calculating the respective position
via Eq. (3). The solution is found as

$\begin{pmatrix}t_{1}\\\ t_{2}\end{pmatrix}=A^{-1}{\bf B},$ (5)

with the $2\times 2$ matrix $A$ given by

$A=\left({\bf X}_{j}-{\bf X}_{i}\ {\bf X}_{k}-{\bf X}_{l}\right)$ (6)

and the right side vector

${\bf B}={\bf X}_{k}-{\bf X}_{i}.$ (7)

## References

  * [1] S. J. Pennycook, L. A. Boatner, Chemically sensitive structure-imaging with a scanning transmission electron microscope, Nature 336 (1988) 565–567.  doi:10.1038/336565a0.
  * [2] S. D. Findlay, N. Shibata, H. Sawada, E. Okunishi, Y. Kondo, Y. Ikuhara, Dynamics of annular bright field imaging in scanning transmission electron microscopy, Ultramicroscopy 110 (2010) 903–923.  doi:10.1016/j.ultramic.2010.04.004.
  * [3] M. Bosman, V. J. Keast, J. L. Garcia-Munoz, A. J. D’Alfonso, S. D. Findlay, L. J. Allen, Two-dimensional mapping of chemical information at atomic resolution, Physical Review Letters 99 (8) (2007) 086102.  doi:10.1103/PhysRevLett.99.086102.
  * [4] A. J. D’Alfonso, B. Freitag, D. Klenov, L. J. Allen, Atomic-resolution chemical mapping using energy-dispersive X-ray spectroscopy, Physical Review B 81 (10) (2010) 100101.  doi:10.1103/PhysRevB.81.100101.
  * [5] S. E. Maccagnano-Zacher, K. A. Mkhoyan, E. J. Kirkland, J. Silcox, Effects of tilt on high-resolution ADF-STEM imaging, Ultramicroscopy 108 (8) (2008) 718–726.  doi:10.1016/j.ultramic.2007.11.003.
  * [6] D. Zhou, K. Müller-Caspary, W. Sigle, F. F. Krause, A. Rosenauer, P. A. van Aken, Sample tilt effects on atom column position determination in ABF–STEM imaging, Ultramicroscopy 160 (2016) 110–117.  doi:10.1016/j.ultramic.2015.10.008.
  * [7] K. E. MacArthur, A. B. Yankovich, A. Béché, M. Luysberg, H. G. Brown, S. D. Findlay, M. Heggen, L. J. Allen, Optimizing experimental conditions for accurate quantitative energy-dispersive X-ray analysis of interfaces at the atomic scale, Microscopy and Microanalysis 27 (2021) 528–542.  doi:10.1017/S1431927621000246.
  * [8] K. E. MacArthur, H. G. Brown, S. D. Findlay, L. J. Allen, Probing the effect of electron channelling on atomic resolution energy dispersive X-ray quantification, Ultramicroscopy 182 (2017) 264–275.  doi:10.1016/j.ultramic.2017.07.020.
  * [9] N. R. Lugg, G. Kothleitner, N. Shibata, Y. Ikuhara, On the quantitativeness of EDS STEM, Ultramicroscopy 151 (2015) 150–159.  doi:10.1016/j.ultramic.2014.11.029.
  * [10] Z. Chen, Y. Jiang, Y.-T. Shao, M. E. Holtz, M. Odstrčil, M. Guizar-Sicairos, I. Hanke, S. Ganschow, D. G. Schlom, D. A. Muller, Electron ptychography achieves atomic-resolution limits set by lattice vibrations, Science 372 (6544) (2021) 826–831.  arXiv:https://science.sciencemag.org/content/372/6544/826.full.pdf, doi:10.1126/science.abg2533.
  * [11] H. Sha, Y. Ma, G. Cao, J. Cui, W. Yang, Q. Li, R. Yu, Sub-nanometer-scale mapping of crystal orientation and depth-dependent structure of dislocation cores in SrTiO3, Nature Communications 14 (1) (2023) 162.  doi:10.1038/s41467-023-35877-7.
  * [12] B. Diederichs, Z. Herdegen, A. Strauch, F. Filbir, K. Müller-Caspary, Exact inversion of partially coherent dynamical electron scattering for picometric structure retrieval, Nature Communications 15 (1) (2024) 101.  doi:10.1038/s41467-023-44268-x.
  * [13] W. Van den Broek, C. T. Koch, Method for retrieval of the three-dimensional object potential by inversion of dynamical electron scattering, Physical Review Letters 109 (2012) 245502.  doi:10.1103/PhysRevLett.109.245502.
  * [14] Z. Liao, N. Gauquelin, R. J. Green, K. Müller-Caspary, I. Lobato, L. Li, S. V. Aert, J. Verbeeck, M. Huijben, M. N. Grisolia, V. Rouco, R. El Hage, J. E. Villegas, A. Mercy, M. Bibes, P. Ghosez, G. A. Sawatzky, G. Rijnders, G. Koster, Metal-insulator-transition engineering by modulation tilt-control in perovskite nickelates for room temperature optical switching, Proceedings of the National Academy of Sciences 115 (38) (2018) 9515–9520.  doi:10.1073/pnas.1807457115.
  * [15] M. Calcagnotto, D. Ponge, E. Demir, D. Raabe, Orientation gradients and geometrically necessary dislocations in ultrafine grained dual-phase steels studied by 2d and 3d EBSD, Materials Science and Engineering: A 527 (10) (2010) 2738–2746.  doi:https://doi.org/10.1016/j.msea.2010.01.004.
  * [16] A. Strauch, B. März, T. Denneulin, M. Cattaneo, A. Rosenauer, K. Müller-Caspary, Systematic errors of electric field measurements in ferroelectrics by unit cell averaged momentum transfers in STEM, Microscopy and Microanalysis 29 (2023) 499–511.  doi:10.1093/micmic/ozad016.
  * [17] K. Müller, H. Ryll, I. Ordavo, S. Ihle, L. Strüder, K. Volz, J. Zweck, H. Soltau, A. Rosenauer, Scanning transmission electron microscopy strain measurement from millisecond frames of a direct electron charge coupled device, Applied Physics Letters 101 (2012) 212110.  doi:10.1063/1.4767655.
  * [18] R. Plackett, I. Horswell, E. N. Gimenez, J. Marchal, D. Omar, N. Tartoni, Merlin: a fast versatile readout system for Medipix3, Journal of Instrumentation 8 (2013) C01038.  doi:10.1088/1748-0221/8/01/C01038.
  * [19] K. Müller-Caspary, A. Oelsner, P. Potapov, Two-dimensional strain mapping in semiconductors by nano-beam electron diffraction employing a delay-line detector, Applied Physics Letters 107 (2015) 072110.  doi:10.1063/1.4927837.
  * [20] H. Ryll, M. Simson, R. Hartmann, P. Holl, M. Huth, S. Ihle, Y. Kondo, P. Kotula, A. Liebel, K. Müller-Caspary, A. Rosenauer, R. Sagawa, J. Schmidt, H. Soltau, L. Strüder, A pnCCD-based, fast direct single electron imaging camera for TEM and STEM, Journal of Instrumentation 11 (2016) P04006.  doi:10.1088/1748-0221/11/04/P04006.
  * [21] M. W. Tate, P. Purohit, D. Chamberlain, K. X. Nguyen, R. Hovden, C. S. Chang, P. Deb, E. Turgut, J. T. Heron, D. G. Schlom, D. C. Ralph, G. D. Fuchs, K. S. Shanks, H. T. Philipp, D. A. Muller, S. M. Gruner, High dynamic range pixel array detector for scanning transmission electron microscopy, Microscopy Microanalysis 22 (2016) 237–249.  doi:10.1017/S1431927615015664.
  * [22] D. Jannis, C. Hofer, C. Gao, X. Xie, A. Béché, T. J. Pennycook, J. Verbeeck, Event driven 4D STEM acquisition with a Timepix3 detector: Microsecond dwell time and faster scans for high precision and low dose applications, Ultramicroscopy 233 (2022) 113423.  doi:10.1016/j.ultramic.2021.113423.
  * [23] W. Xu, J. M. LeBeau, A deep convolutional neural network to analyze position averaged convergent beam electron diffraction patterns, Ultramicroscopy 188 (2018) 59–69.  doi:10.1016/j.ultramic.2018.03.004.
  * [24] J. M. LeBeau, S. D. Findlay, L. J. Allen, S. Stemmer, Position averaged convergent beam electron diffraction: Theory and applications, Ultramicroscopy 110 (2010) 118–125.  doi:10.1016/j.ultramic.2009.10.001.
  * [25] S. Nishikawa, S. Kikuchi, Diffraction of cathode rays by mica, Nature 121 (3061) (1928) 1019–1020.  doi:10.1038/1211019a0.
  * [26] J. A. Venables, C. J. Harland, Electron back-scattering patterns — A new technique for obtaining crystallographic information in the scanning electron microscope, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 27 (5) (1973) 1193–1200.  doi:10.1080/14786437308225827.
  * [27] P. W. Trimby, Orientation mapping of nanostructured materials using transmission Kikuchi diffraction in the scanning electron microscope, Ultramicroscopy 120 (2012) 16–24.  doi:10.1016/j.ultramic.2012.06.004.
  * [28] G. L. Burton, S. Wright, A. Stokes, D. R. Diercks, A. Clarke, B. P. Gorman, Orientation mapping with kikuchi patterns generated from a focused stem probe and indexing with commercially available EDAX software, Ultramicroscopy 209 (2020) 112882.  doi:10.1016/j.ultramic.2019.112882.
  * [29] L. Reimer, H. Kohl, Transmission Electron Microscopy, Springer, New York, NY, 2008.  doi:10.1007/978-0-387-40093-8.
  * [30] J. M. Cowley, A. F. Moodie, The scattering of electrons by atoms and crystals. I. A new theoretical approach, Acta Crystallographica 10 (10) (1957) 609–619.  doi:10.1107/S0365110X57002194.
  * [31] J. Barthel, Dr. Probe: A software for high-resolution STEM image simulation, Ultramicroscopy 193 (2018) 1–11.  doi:10.1016/j.ultramic.2018.06.003.
  * [32] B. D. Forbes, A. V. Martin, S. D. Findlay, A. J. D’Alfonso, L. J. Allen, Quantum mechanical model for phonon excitation in electron diffraction and imaging using a Born-Oppenheimer approximation, Physical Review B 82 (2010) 104103.  doi:10.1103/PhysRevB.82.104103.
  * [33] M. E. Straumanis, D. C. Kim, Lattice constants, thermal expansion coefficients, densities and perfection of structure of pure iron and of iron loaded with hydrogen, International Journal of Materials Research 60 (4) (1969) 272–277.  doi:10.1515/ijmr-1969-600404.
  * [34] L.-M.Peng, G. Ren, S. L. Dudarev, M. J. Whelan, Debye–Waller Factors and Absorptive Scattering Factors of Elemental Crystals, Acta Crystallographica Section A 52 (3) (1996) 456–470.  doi:10.1107/S010876739600089X.
  * [35] Ernst Ruska-Centre for Microscopy and Spectroscopy with Electrons (ER-C), Fei titan 80-300 stem, Journal of large-scale research facilities 2 (2016) A42.  doi:10.17815/jlsrf-2-67.
  * [36] Ernst Ruska-Centre for Microscopy and Spectroscopy with Electrons (ER-C), Fei titan g2 80-200 crewley, Journal of large-scale research facilities 2 (2016) A43.  doi:10.17815/jlsrf-2-68.
  * [37] G. W. Paterson, R. J. Lamb, R. Ballabriga, D. Maneuski, V. O’Shea, D. McGrouther, Sub-100 nanosecond temporally resolved imaging with the Medipix3 direct electron detector, Ultramicroscopy 210 (2020) 112917.  doi:10.1016/j.ultramic.2019.112917.
  * [38] M. Kruth, D. Meertens, K. Tillmann, FEI helios nanolab 460f1 fib-sem, Journal of large-scale research facilities JLSRF 2 (03 2016).  doi:10.17815/jlsrf-2-105.
  * [39] M. Keunecke, F. Lyzwa, D. Schwarzbach, V. Roddatis, N. Gauquelin, K. Müller-Caspary, J. Verbeeck, S. J. Callori, F. Klose, M. Jungbauer, V. Moshnyaga, High-$t_{\mathrm{C}}$ interfacial ferromagnetism in srmno3/lamno3 superlattices, Advanced Functional Materials 30 (18) (2020) 1808270.  doi:10.1002/adfm.201808270.
  * [40] K. Huang, X-ray reflexions from dilute solid solutions, Proc. Roy. Soc. Lond. A 190 (1020) (1947) 102–117.  doi:10.1098/rspa.1947.0064.


1]Institut für Theoretische Physik und Astrophysik
Universität Würzburg
Am Hubland
97074 Würzburg

M. Weidinger
<EMAIL_ADDRESS>
13 January 2010

# Modelling the steady state spectral energy distribution of the BL-Lac Object
PKS 2155-30.4 using a selfconsistent SSC model

M. Weidinger  M. Rüger  F. Spanier  [

(6 April 2009; 19 October 2009; 23 November 2009)

###### Abstract

In this paper we present a fully selfconsistent SSC model with particle
acceleration due to shock and stochastic acceleration (Fermi-I and Fermi-II-
Processes respectively) to model the quiescent spectral energy distribution
(SED) observed from PKS 2155. The simultaneous August/September 2008
multiwavelength data of H.E.S.S., Fermi, RXTE/SWIFT and ATOM give new
constraints to the high-energy peak in the SED concerning its curvature. We
find that, in our model, a monoenergetic injection of electrons at
$\gamma_{0}=910$ into the model region, which are accelerated by Fermi-I- and
Fermi-II-processes while suffering synchrotron and inverse Compton losses,
finally leads to the observed SED of PKS 2155-30.4 shown in H. E. S. S. and
Fermi-LAT collaborations (2009). In contrast to other SSC models our
parameters arise from the jet’s microphysics and the spectrum is evolving
selfconsistently from diffusion and acceleration. The $\gamma_{0}$-factor can
be interpreted as two counterstreaming plasmas due to the motion of the blob
at a bulk factor of $\Gamma=58$ and opposed moving upstream electrons at
moderate Lorentz factors with an average of $\gamma_{u}\approx 8$.

[INTRODUCTION] Among the class of active galactic nuclei (AGN), blazars are
showing a spectral energy distribution (SED) that is strongly dominated by
nonthermal emission across a wide range of wavelengths, from radio waves to
gamma rays, and rapid, large-amplitude variability. Presumably, these
characteristics are due to a highly relativistic jet which covers a small
angle to the line-of-sight, emitting the observable Doppler-boosted
synchrotron and inverse Compton radiation.
In SSC models the characteristic double humped spectra of blazars are
explained by electrons in the jet emitting synchrotron radiation while being
accelerated in a magnetic field, which gives the first peak in the SED. These
high energy electrons upscatter the very same synchrotron photons to TeV
energies due to the inverse Compton effect, resulting in a second peak in the
SED. Another approach to the explanation of the double humped structure are
proton initiated electromagnetic cascades (e.g., Mannheim, 1993) or external
Compton models.
The key issue is to understand which physical mechanisms are leading to such
SEDs. In particular that means to explain the ultrarelativistic electron
spectra within the jet, which are believed to be responsible for the gamma
radiation.
The high peaked BL Lac objects (HBLs) as a subclass of blazars show a peak in
their SED in the X-ray regime, suggesting that an inverse Compton peak should
occur at correspondingly high gamma-ray energies. In fact, a large fraction of
the known nearby HBLs have already been discovered with Cerenkov telescopes,
such as H.E.S.S., MAGIC, and VERITAS. Since 2008 the Fermi satellite measures
at these high gamma-ray energies. The energy range of the Fermi data is
slightly different from the H.E.S.S. and VERITAS-Telescopes which gives new
constraints to the SEDs.
The first Fermi data published is from PKS 2155-30.4, a HBL at redshift
$z=0.117$ (luminosity distance: $d_{\text{L}}=1.67\cdot
10^{27}\leavevmode\nobreak\ \text{cm}$) (H. E. S. S. and Fermi-LAT
collaborations, 2009).

We present a selfconsistent SSC model that is not only able to model the SED
of PKS 2155-30.4 shown in H. E. S. S. and Fermi-LAT collaborations (2009) but
also to partly explain the “ad-hoc” injected particle spectra of many SSC
models. Therefore we introduce and solve the kinetic equation describing the
synchrotron-self-Compton emission numerically in two different zones within
the jet (see section 1). We use the exact Klein-Nishina cross section which is
important at the relevant very high gamma-energies to describe the inverse
Compton radiation and energy losses of the electrons. The emphasis lies on the
accurate treatment of the two possible particle acceleration mechanisms
(Fermi-I- and Fermi-II) which are able to produce high energy electrons as
well as on the selfconsistent treatment of the radiation processes.

##  1 THE MODEL

###  1.1 Model geometry

We extend the well-established SSC model by Fermi-I and Fermi-II acceleration
mechanisms to a selfconsistent SSC model with two zones in a nested setup.
Both regions (the acceleration- and the radiation zone) forming the blob are
assumed to be spherical and homogeneous containing isotropically distributed
non-thermal electrons and a randomly oriented magnetic field. The acceleration
zone is assumed to be spatially significantly smaller than the surrounding
radiation zone. Furthermore every electron leaving the accelerationzone enters
the radiation zone. These assumptions are common place in SSC models (e.g.,
Kirk et al., 1998).
To derive the kinetic equations describing the time evolution of
$n_{e}(\gamma)$, $N_{e}(\gamma)$ ($n_{e}$ in the acceleration zone, $N_{e}$ in
the radiation zone) as the differential electron densities we use the one
dimensional diffusion approximation (eq. (1)) of the relativistic Vlasov
equation (e.g. Schlickeiser, 2002), which is applicable due to the assumptions
made above.

$\displaystyle\frac{\partial}{\partial t}f(p,t)$
$\displaystyle=\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[F\left(p,f,\frac{\partial}{\partial
p}f\right)\right]+S(p,t)\leavevmode\nobreak\ ,$ (1)

where $f(p,t)$ is a particle distribution function, with the particle’s
momentum value $p$. $F$ describes the contributing processes, such as
synchrotron radiation or acceleration, in momentum space. Catastrophic
particle gains and losses are considered via $S(p,t)$.
Making use of the relativistic approximation $E\approx pc=\gamma mc$ and the
relation $n(p,t)=4\pi p^{2}f(p,t)$ one can derive the kinetic equations
governing the model.

###  1.2 Kinetic equations

####  1.2.1 Acceleration zone

While the blob propagates through the jet, electrons are continuously injected
into the acceleration zone when considering the blob’s rest frame, leading to
an injection function

$\displaystyle Q_{\text{inj}}(\gamma,t)$
$\displaystyle:=Q_{0}\delta(\gamma-\gamma_{0})\vartheta{(t-t_{0})}\leavevmode\nobreak\
\leavevmode\nobreak\ ,\leavevmode\nobreak\ t_{0}=0$ (2)

which we assume monoenergetic and time independent. These low- to mid-energy
electrons are accelerated systematically and stochastically due to Fermi-I and
Fermi-II processes while suffering synchrotron and inverse Compton losses.
Energy losses due to inverse Compton scattering are calculated using the full
Klein-Nishina cross section, see eq. (1.2.2). This leads to
$P_{\text{IC}}(\gamma)$ given in eq. (8) with the corresponding radiation
field $n_{\text{PH}}$ of the acceleration zone. Due to the non equilibrium of
magnetic and radiative energy in the acceleration zone the energy losses via
inverse Compton scattering can become quite significant and must not be
neglected, also the Thomson limit is not appropriate here. The synchrotron
losses are calculated using eq. (3) from Ginzburg and Syrovatskii (1969) for
isotropic particle distributions

$\displaystyle P_{s}(\gamma)$
$\displaystyle=\frac{1}{6\pi}\frac{\sigma_{\text{T}}B^{2}}{mc}\gamma^{2}=\beta_{s}\gamma^{2}$
(3)

with the Thomson cross section $\sigma_{\text{T}}$. According to Schlickeiser
(1984) particle acceleration via parallel shockfronts and stochastic
acceleration caused by scattering at Alfvén waves leads to

$\displaystyle F\left(p,f,\frac{\partial}{\partial p}f\right)$
$\displaystyle=p^{4}\frac{v_{A}^{2}}{9K_{||}}\frac{\partial f}{\partial
p}+p^{3}\frac{v_{S}^{2}}{4K_{||}}f$ (4)

for the function $F$. With the parallel spatial diffusion coefficient
$K_{||}$, which is momentum independent for hard spheres and the
characteristic speeds $v_{A}$ for the Alfvén mediated stochastic acceleration
and $v_{S}$ for parallel shockfronts. Substituting $p\rightarrow\gamma$ in eq.
(1) and eq. (4) according to the relativistic approximation mentioned above,
one will finally find eq. (1.2.1); the kinetic equation of the acceleration
zone.

$\displaystyle\frac{\partial n_{e}(\gamma,t)}{\partial t}=$
$\displaystyle\frac{\partial}{\partial\gamma}\left[(\beta_{s}\gamma^{2}+P_{\text{IC}}(\gamma)-t_{\text{acc}}^{-1}\gamma)\cdot
n_{e}(\gamma,t)\right]+$
$\displaystyle+\frac{\partial}{\partial\gamma}\left[[(a+2)t_{\text{acc}}]^{-1}\gamma^{2}\frac{\partial
n_{e}(\gamma,t)}{\partial\gamma}\right]+$
$\displaystyle+Q_{0}(\gamma-\gamma_{0})-t_{\text{esc}}^{-1}n_{e}(\gamma,t)\leavevmode\nobreak\
,$ (5)

where the characteristic acceleration timescale $t_{\text{acc}}$ is given by

$\displaystyle t_{\text{acc}}$
$\displaystyle=\left(\frac{v_{s}^{2}}{4K_{||}}+2\frac{v_{A}^{2}}{9K_{||}}\right)^{-1}\leavevmode\nobreak\
.$ (6)

Eq. (6) for $t_{\text{acc}}$ is a direct consequence of the derivation of eq.
(1.2.1) out of eq. (1) using eqn. (4) and (3). The expression in eq. (6)
includes the analytical timescale for non-relativistic shock acceleration.
According to Bednarz and Ostrowski (1996) and especially to Ellison et al.
(1990) the acceleration timescale for parallel relativistic shock waves
decreases approximately by a factor of $3$. We did not take into account this
behavior for it is unclear how the analytical expression looks like in that
case. Secondly we are omitting the energy dependency of $t_{\text{acc}}$ using
hard spheres for the plasma instabilities anyway. This issue is irrelevant for
the modelling (for we are setting numerical values for $t_{\text{acc}}$) and
the type of energy spectrum (a powerlaw) produced by Fermi-I acceleration is
identical in the non-relativistic and relativistic case (Sokolov et al.,
2004). But it has to be kept in mind for the interpretation of the results.
The parameter $a\approx v_{s}^{2}/v_{A}^{2}$ determines the ratio of shock to
stochastic acceleration. $t_{\text{esc}}=\eta R_{\text{acc}}/c$ is the
characteristic timescale for electrons escaping from the acceleration region,
where $\eta$ is an empirical factor set to $\eta=10$ and $R_{\text{acc}}$ the
radius of the acceleration sphere. All escaping electrons enter the radiation
zone downstream the jet. The seperation in two zones can firstly explain the
injected electron spectra and secondly takes account of a much more confined
shock region for Fermi-I acceleration will probably not occur in the whole
blob region when considering physical sources.
Our model can be compared with the model presented by Katarzyński et al.
(2006). The kinetic equation (eq. 3 in their paper) is almost similar to the
kinetic equation in the acceleration zone eq. 1.2.1. One major difference to
our model is their sole use of stochastic acceleration. In fact their model is
the limit of our model for $a\to 0$. Additionally they limit themselves to
radiation in the acceleration zone, which is useful when not taking into
account shock acceleration. Besides that there are number of minor differences
regarding the exact treatment of inverse Compton losses and the derivation of
escape rates.
Due to the small spatial extent the acceleration zone does not contribute to
the SED directly, i.e. $n_{\text{ph}}(\nu)$ is only calculated in order to
determine the inverse Compton loss rate for the electrons in the acceleration
zone.

####  1.2.2 Radiation zone

The electrons are not accelerated here. Thus the kinetic equation takes the
simple form

$\displaystyle\frac{\partial N_{e}(\gamma,t)}{\partial t}=$
$\displaystyle\frac{\partial}{\partial\gamma}\left[\left(\beta_{s}\gamma^{2}+P_{\text{IC}}(\gamma)\right)\cdot
N_{e}(\gamma,t)\right]-$
$\displaystyle-\frac{N_{e}(\gamma,t)}{t_{\text{rad,esc}}}+\left(\frac{R_{\text{acc}}}{R_{\text{rad}}}\right)^{3}\frac{n_{e}(\gamma,t)}{t_{\text{esc}}}\leavevmode\nobreak\
\text{.}$ (7)

Electrons in the radiation zone suffer synchrotron
($\propto\beta_{s}\gamma^{2}$) and inverse Compton losses (eq. (8)), other
energy losses are irrelevant in jetsystems of such low electron density (see
e.g. Böttcher and Chiang, 2002).

$\displaystyle P_{\text{IC}}(\gamma)$
$\displaystyle=m^{3}c^{7}h\int_{0}^{\alpha_{max}}{d\alpha\alpha\int_{0}^{\infty}{d\alpha_{1}N_{\text{ph}}(\alpha_{1})\frac{dN(\gamma,\alpha_{1})}{dtd\alpha}}}$
(8)

The integrals in eq. (8) are solved numerically using the full Klein-Nishina
cross section for a single electron given in eq. (1.2.2). The photon energies
are rewritten in terms of the electrons rest mass, i.e. $h\nu=\alpha mc^{2}$
for the scattered photons and $h\nu=\alpha_{1}mc^{2}$ for the target photons.
The integration bounds of the outer integral in eq. (8) are a direct
consequence of the kinematics. The non-trivial dependency of
$P_{\text{IC}}(\gamma)$ from $N_{\text{ph}}$ and thus of $N_{e}$ from eq. (8)
makes the numerical treatment of the kinetic equations inevitable leading to a
time resolved model. The loss rates for the electron distribution of PKS 2155
in the steady state are shown in Fig. 1, which indicates that in our case the
inverse Compton losses would be slightly overestimated in the often used
Thomson limit for all Lorentzfactors $\gamma$ because of the dependency on the
photon field $P_{\text{thom}}\propto\int{d\nu\nu N_{\text{ph}}}\gamma^{2}$ and
its special shape due to the modified injection at $\gamma_{0}=910$. This
would not be the case for low energetic injected electrons and the resulting
powerlaw-like photon distribution. For high Lorentzfactors $\gamma$ however
the deviation of the Thompson limit for the inverse Compton scattering to the
real Klein-Nishina treatment becomes more significant in each case.

Figure 1: Lossrates due to the inverse Compton effect for an electron of
Lorentzfactor $\gamma$ in the model photon field of PKS 2155 (Thompson limit,
red and used Klein-Nishina treatment, black) compared to the synchrotron
losses, blue. Other losses like adiabatic cooling, pair production are
irrelevant at typical SSC configurations.

Again electrons escaping the blob are taken into account via
$t_{\text{esc,rad}}=\eta R_{\text{rad}}/c$ with the empirical factor $\eta$
set to $\eta=10$. Both electronic PDEs are connected via the catastrophic
particle loss/gain-term. The factor
$\left(R_{\text{acc}}/R_{\text{rad}}\right)^{3}$ ensures particle
conservation.

To determine the time-dependent spectral energy distribution of blazars we
solve the differential equation for the differential photon number density,
obtained from the radiative transfer equation, including the corresponding
terms with respect to the SSC model,

$\displaystyle\frac{\partial N_{\text{ph}}(\nu,t)}{\partial t}$
$\displaystyle=R_{s}-c\alpha_{\nu}N_{\text{ph}}(\nu,t)+R_{c}-\frac{N_{\text{ph}}(\nu,t)}{t_{\text{ph,esc}}}\leavevmode\nobreak\
\text{.}$ (9)

To describe the synchrotron photon production rate $R_{s}$ in a convenient way
we use the well known Melrose approximation (see e.g., Pohl, 2002)

$\displaystyle R_{s}$ $\displaystyle=1.8\frac{\sqrt{3}e^{3}B_{\bot}}{h\nu
mc^{2}}\int{d\gamma
N_{e}(\gamma,t)\left(\frac{\nu}{\nu_{c}(\gamma)}\right)^{\frac{1}{3}}e^{-\frac{\nu}{\nu_{c}(\gamma)}}}\leavevmode\nobreak\
.$ (10)

with the characteristic frequency eq. (11) for an electron with a Lorentz
factor of $\gamma$.

$\displaystyle\nu_{c}(\gamma)$ $\displaystyle=\frac{3\gamma^{2}eB_{\bot}}{4\pi
mc}$ (11)

In optically thick regimes the emitted synchrotron radiation is absorbed by
the emitting electrons itself. This is described by the synchrotron self
absorption coefficient,

$\displaystyle\alpha_{\nu}$
$\displaystyle=\frac{1}{12}\frac{c}{\nu^{2}eB}P_{\nu}(\gamma_{c})\frac{N_{e}(\gamma_{0},t)}{\gamma_{0}^{2}}\text{.}$
(12)

(with $\gamma_{c}=f(\nu_{c})^{-1}$). Here we made use of the monochromatic
approximation (e.g., Felten and Morrison, 1966) for the synchrotron power:

$\displaystyle P_{\nu}(\nu,\gamma)$
$\displaystyle=\frac{\sqrt{3}e^{2}B_{\bot}}{mc^{2}}\frac{\nu}{\nu_{c}{\gamma}}\int_{0}^{\infty}{d\nu^{\prime}K_{\frac{5}{3}}(\nu^{\prime})}$
(13)

In SSC models the second hump in the SED of a blazar is due to inverse Compton
scattered photons by the synchrotron radiation emitting electrons themselves.
Here the full Klein-Nishina cross section from Blumenthal and Gould (1970) is
used to calculate the inverse Compton photon production rate.

$\displaystyle R_{c}=$ $\displaystyle\int d\gamma\,N_{e}(\gamma)\;\times$
$\displaystyle\times\int
d\alpha_{1}\left[N_{\text{ph}}(\alpha_{1})\frac{dN(\gamma,\alpha_{1})}{dtd\alpha}-N_{\text{ph}}(\alpha)\frac{dN(\gamma,\alpha)}{dtd\alpha_{1}}\right]$
(14)

To fully exploit the Klein-Nishina cross section, eq. (1.2.2), we used the
approximate inverse Compton spectrum of a single electron scattering off a
unit density photon field (e.g., Jones, 1968; Jauch and Rohrlich, 1976):

$\displaystyle\frac{dN(\gamma,\alpha_{1})}{dtd\alpha}=$
$\displaystyle\frac{2\pi r_{0}^{2}c}{\alpha_{1}\gamma^{2}}\bigg{[}2q\ln
q+(1+2q)(1-q)+$ $\displaystyle+\frac{1}{2}\frac{(4\alpha_{1}\gamma
q)^{2}}{(1+4\alpha_{1}\gamma q)}(1-q)\bigg{]}\leavevmode\nobreak\ ,$ (15)

with the electron’s Lorentz radius $r_{0}=e^{2}/(mc^{2})$, the scattering
parameter $q={\alpha}/(4\alpha_{1}\gamma^{2}(1-\alpha/\gamma))$ and $0\approx
1/(4\gamma^{2})<q\leq 1$. Due to momentum and energy conservation this
equation is valid for $\alpha_{1}<\alpha\leq
4\alpha_{1}\gamma^{2}/(1+4\alpha_{1}\gamma)$.
The last catastrophic term in eq. (9) describes photons escaping from the
emitting region, where

$\displaystyle t_{\text{ph, esc}}=\frac{3R_{\text{rad}}}{4c},$ (16)

is the approximate escape time, with $R_{\text{rad}}$ the radius of the
emitting blob. The escape time is chosen to be the light crossing time of the
photons.
The photon lossrate due to the pair production of electrons and positrons is
not taken into account for two reasons. Firstly it is insignificant compared
to the dominating synchrotron and inverse Compton processes. This is a
consequence of the relatively low density (Böttcher and Chiang, 2002).
Secondly it would violate the selfconsistency of our model for positrons are
not treated, hence violating energy conservation.

To compute the SEDs in our model we must shift the frame of reference from the
blob to the observer. For a sphere of radius $R$ the observed flux at distance
$r$ is

$\displaystyle F_{\nu}^{obs}(r)$ $\displaystyle=\pi
I_{\nu}^{obs}\frac{R^{2}}{r^{2}}\leavevmode\nobreak\ \text{.}$ (17)

With the Lorentz boosted intensity $I_{\nu}^{obs}=\delta^{3}I_{\nu}^{blob}$
due to the bulk motion with a doppler factor $\delta$ of the blob and the
Lorentz transformed, red shifted frequency $\nu^{obs}=\delta/(1+z)\nu$. Where
$I_{\nu}^{blob}$ is calculated from the photon unit density

$\displaystyle I_{\nu}^{blob}$ $\displaystyle=\frac{h\nu
c}{4\pi}N_{\text{ph}}(\nu)$ (18)

for homogenous spheres.

##  2 NUMERICS

In our model we numerically solve the kinetic equations forward in time in
order to obtain a model SED. The downstream motion of the electrons induces
the sequence of solving the acceleration zone’s equation before the kinetic
equation of the radiation zone in each time step. The simple Euler scheme was
found adequate to do the time integration.
In the acceleration zone we had to combine the Crank-Nicholson scheme (Press,
2002) with Godunov’s method to provide both correct treatment of the
characteristics and stability for the derivation in $\gamma$. In the radiation
zone the characteristic flows, due to the absence of acceleration, only in one
direction making the Crank-Nicholson scheme sufficient.
With our carefully tested code it is possible to calculate the dynamics of
SEDs in a range of 20 orders of magnitude. The implemented code complies
particle conservation in each zone alone and both together as well as the
conservation of the total energy (i.e. of the electrons and the photons) over
typical simulation times with a maximum error of O(5%). For negligible
stochastic acceleration, i.e. $a\rightarrow\infty$, and without a radiation
field, i.e. no inverse Compton losses, the steady state solution for the
kinetic equation yields

$\displaystyle
n_{\text{e,steady}}(\gamma)=C\frac{1}{\gamma^{2}}\left(\frac{1}{\gamma}-\beta_{\text{s}}t_{\text{acc}}\right)^{\frac{t_{\text{acc}}-t_{\text{esc}}}{t_{\text{esc}}}}$

with $\gamma_{\text{max}}=\left(t_{\text{acc}}\beta_{\text{s}}\right)^{-1}$
and the constant $C$ determined by the injection function $Q$. The implemented
numeric model converges against this solution for sufficient simulation time.
Additionally it was tested against the steady state analytical solution with
Fermi-II processes given in Schlickeiser (1984) with no significant
deviations. Setting $t_{\text{esc,rad}}\rightarrow\infty$ and neglecting
inverse Compton scattering the spectral index of the powerlaw part of the
electron distribution in the radiation zone is (analytically) reduced by one
compared to the one in the acceleration zone, which also was confirmed by the
implemented code. The inverse Compton scattering rate was confirmed against
the approximate analytical results (low energetic Thompson regime and extreme
Klein-Nishina limit) before implementing. Concerning the photon distributions
we validated the expected spectral indices in the steady state solution for
the different frequency regimes, which together with the energy conservation
between electrons and photons approves the integrity of the model.
A detailed description of the used numeric techniques as well as the
implemented model also in context with the variability of the sources will be
given in a paper yet to be published.

##  3 RESULTS

The recent Fermi data give new constraints on the gamma-ray peak of the HBL
PKS 2155-30.4 concerning its curvature. This is leading to a deep dip between
the optical/X-ray and the gamma-ray peak. We are able to model the SED of PKS
2155-30.4 with our model by setting

$\displaystyle\gamma_{0}$ $\displaystyle=910$ (19)

for the monoenergetic injection into the acceleration zone. This is rather
unusual but required to model the SED of PKS 2155. Such moderate but not small
Lorentz factors can be explained e.g. by two counterstreaming plasmas. If the
upstream electrons would be at rest, the bulk doppler factor of $\delta=116$
would automatically lead to $\gamma_{0}=\Gamma\approx 58$. Assuming
speculatively that the upstream electrons moving in the opposite direction of
the blob with a mean velocity of $v_{u}$ hence a upstream Lorentz-Factor
$\gamma_{u}$ the $\gamma_{0}$ factor in the blob’s rest frame must be
calculated according to the relativistic superposition:

$\displaystyle\gamma_{0}$
$\displaystyle=\sqrt{1-\left(\frac{\sqrt{\Gamma^{2}-1}\Gamma+\sqrt{\gamma_{u}^{2}-1}\gamma_{u}}{\Gamma^{2}+\gamma_{u}^{2}-1}\right)^{2}}^{-1}$
(20)

Solving eq. (20) for our setup we find $\gamma_{u}\approx 8$ for the upsteam
electrons which are streaming towards the blob. The numerically solved steady
state electron density in the acceleration zone is shown in Fig. 2. We also
show the time development for a “switched on” injection, i.e.
$\left.n_{\text{e}}(\gamma)\right|_{t<0}=0\leavevmode\nobreak\
\forall_{\gamma}$ and $Q=Q_{0}\delta(\gamma-\gamma_{0})\vartheta(t)$, until
the steady state is reached.

Figure 2: Steady state electron distribution and its time development in the
acceleration zone modelling the SED of PKS 2155-30.4 shown in Fig. 3 as
arising from the injection function $Q_{0}$. The corresponding intrinsic times
are $t=1000\leavevmode\nobreak\ \text{s}$ (dashed black),
$t=5000\leavevmode\nobreak\ \text{s}$ (dashed blue), $t=1\cdot
10^{4}\leavevmode\nobreak\ \text{s}$ (dashed red), $t=2\cdot
10^{4}\leavevmode\nobreak\ \text{s}$ (dashed green). The steady state with its
rising and falling powerlaw and the exponential cutoff at $\gamma\approx
10^{5}$ is reached at about $t=10^{5}\leavevmode\nobreak\ \text{s}$.

In Fig. 2 it can clearly be seen tvhat accelerating electrons using Fermi-I
and Fermi-II processes leads to powerlaw electron distributions with an
exponential cut-off often used as the ad-hoc injection function in onezone-SSC
models (Chiang and Böttcher, 2002) (right side of $\gamma_{0}$ in Fig. 2),
thus explaining them using the diffusion theory derived from plasma physics.
By injecting electrons with eq. (19) and significant stochastic acceleration
(i.e. $a=O(1)$) we are also able to produce rising electron spectra before
decreasing in a power-law and an exponential cut-off, like introduced in
Böttcher and Chiang (2002) (left side of $\gamma_{0}$ in Fig. 2). The Fermi-II
processes are responsible for the rising power-law and exponential cut-off,
whereas the ratio of $t_{\text{acc}}/t_{\text{esc}}$ determines the spectral
index of the power-law at $\gamma>\gamma_{0}$. It can clearly be seen from
Fig. 2 that the convergence against the steady state solution for the electron
density begins relatively rapid while slowing down eventually. The simulation
time when the steady state is reached corresponds to the escape time of the
electrons in the acceleration zone. When concerning variability and time
resolved lightcurves of blazars this is an advantage of the twozone model
because the rising part in such lightcurves corresponds partially to the
escape time of the acceleration zone (while the falling part is connected to
the response time of the system, $t_{\text{rad,esc}}$).

Table 1: Chosen parameters for the model SED shown in Fig. 4 to fit the data
(H. E. S. S. and Fermi-LAT collaborations, 2009) of PKS 2155-30.4.
$Q_{0}(\text{cm}^{-3})$ | $B(\text{G})$ | $R_{\text{acc}}(\text{cm})$ | $R_{\text{rad}}(\text{cm})$ | $t_{\text{acc}}/t_{\text{esc}}$ | $a$ | $\Gamma$
---|---|---|---|---|---|---
$5.25\cdot 10^{4}$ | $0.29$ | $3.0\cdot 10^{13}$ | $6.3\cdot 10^{14}$ | $1.55$ | $1$ | $58$
|  |  |  |  |  |

An acceleration zone electron density as shown by the solid black curve in
Fig. 2, leads to the desired broken power-law electron spectrum in the
radiation zone which finally is able to model the SED of PKS 2155-30.4 (see
Fig. 3 and Fig. 4). We used the parameters in Table 1 for the model SED in
Fig. 4 (black, solid line). The black dashed curve in Fig. 4 corresponds to a
fit assuming a black body for the thermal contribution of the host galaxy thus
the ATOM optical data is not to be taken into account for the SSC modelling.
The curvature and deep dip in the model SED is a direct consequence of the
rising part in the electron density of the acceleration zone. Thus it can be
modeled by varying the ratio $a$ of shock to stochastic acceleration. All the
parameters in Table 1 are consistent with the limits given via other
observations and statistics, e.g. determination of $\Gamma$ using superluminal
motion of Quasar jets.
The recent Fermi, H.E.S.S. and ATOM data (H. E. S. S. and Fermi-LAT
collaborations, 2009) have been averaged over a period of 14 days and show a
lowstate of the HBL PKS 2155-30.4. This is confirmed by the Aharonian et al.
(2005) data of H.E.S.S. a few years ago which show the same flux level as the
recent data. We used the EBL studies described in Primack et al. (2005) to do
the EBL deabsorption for the H.E.S.S. datapoints, a correction of the Fermi
data is not necessary.

Figure 3: Time evolution due to the switched on injection function until the
steady state of the SED of PKS2155 is reached (see also Fig. 4). The intrinsic
times are $t=1\cdot 10^{4}\leavevmode\nobreak\ \text{s}$ (dashed black),
$t=2\cdot 10^{4}\leavevmode\nobreak\ \text{s}$ (dashed blue), $t=5\cdot
10^{4}\leavevmode\nobreak\ \text{s}$ (dashed red), $t=1\cdot
10^{5}\leavevmode\nobreak\ \text{s}$ (dashed green), the complete steady state
is reached at about $t=2\cdot 10^{6}\leavevmode\nobreak\ \text{s}$, which
correlates to the response time of the system due to $t_{\text{esc,rad}}$.

The time development of the SED due to a switched on injection of electrons
into the acceleration zone at time $t_{0}=0$ is shown in Fig. 3. It can
clearly be seen that the final state of the model SSC correlates with the
response time of the radiation zone $t_{\text{esc,rad}}$ and that the
convergence again begins fast and slows down rapidly at higher simulation
times. [CONCLUSIONS] Our model is able to explain the injection function of
many onezone-SSC models as shock and stochastic acceleration of electrons
upstream the jet entering the blob while continuously suffering synchrotron
losses. By introducing Fermi-II acceleration we get rid of the sharp cut-off
introduced in Kardashev (1962) or Kirk et al. (1998) which probably does not
occur in physical sources. Additionally we are able to model relatively
complex electron densities with increasing and decreasing parts through the
stochastic acceleration of electrons, only by varying the monoenergetic
injection to higher $\gamma_{0}$. In contrast to the ad-hoc injection of some
onezone-SSC models such Lorentz factors have a physically reasonable, but
highly speculative, explanation as upstream previously accelerated but already
partially cooled electrons. These electrons are averagely moving in the
opposite direction of the blob with a mean Lorentz factor of
$\gamma_{u}\approx 8$ resulting, together with the motion of the blob, in
$\gamma_{0}\approx 900$ for the monoenergetic injection function used in the
acceleration zone of our model.

As recent data points out, these complex electron distributions are necessary
to model the new constraints concerning the gamma-ray peak of blazar’s SEDs if
one does not simply shift the synchrotron-peak to achieve the inverse Compton
spectrum (e.g., Kataoka et al., 2000). The curvature of the peak, and thus the
deep dip between the two humps, is a direct consequence of the rising part in
the responsible electron distribution within the blob. This constraint rules
out many SSC models, which are not able to produce such electron spectra.
With our model we are able to form the curvature of the gamma-ray peak and the
dip by varying the influence of the Fermi-II processes. The shape and position
of the synchrotron peak in the model SED is dominated by $t_{\text{acc}}$ and
$R_{\text{acc}},R_{\text{rad}}$ as well as $B$. For the parameters concerning
the acceleration arise from plasmaphysics considerations we gain insight into
the jets microphysics while modelling observed SEDs. We have also shown that
in such environments the Thomson approximation for the inverse Compton effect
can not always be applied, especially when considering time resolution and
hence non equilibria of the energy distribution in the blob.

Here we only introduced steady state solutions of our model, but due to the
spatially relatively small acceleration region, which is at least an oder of
magnitude smaller than the emitting region, this twozone-SSC model is able to
selfconsistently model the rising part in the lightcurves of flaring blazars
which are connected to the behavior in the acceleration zone, especially the
energy transport from low to high energies. This, together with the
consequences of the model geometry on the observable SEDs and lightcurves of
blazars like in Sokolov et al. (2004), will be subject of a following paper.

Figure 4: Lowstate of PKS 2155-30.4 with the simultaneous data of ATOM, SWIFT,
RXTE, Fermi and H.E.S.S of the August/September 2008 campaign from H. E. S. S.
and Fermi-LAT collaborations (2009) (red triangles and circles). The 2003
H.E.S.S. data (blue circles) is also shown, proofing the lowstate of
PKS2155-30.4. The VHE data have been deabsorbed using Primack et al. (2005).
The dashed black curve shows a thermal fit for the contribution of the host-
galaxy. Our model SSC fit, arising from the steady state electron distribution
in the radiation zone is shown in the solid black curve, a moderate energy
injection at $\gamma_{0}\approx 910$ into the acceleration zone together with
stochastic and systematic acceleration is needed to meet the curvature of the
VHE peak given via the Fermi data.

## References

  * Aharonian et al. (2005) Aharonian et al.: Multi-wavelength observations of PKS 2155-304 with HESS, 442, 895–907, 10.1051/0004-6361:20053353, 2005.
  * Bednarz and Ostrowski (1996) Bednarz, J. and Ostrowski, M.: The acceleration time-scale for first-order Fermi acceleration in relativistic shock waves, 283, 447–456, 1996.
  * Blumenthal and Gould (1970) Blumenthal, G. R. and Gould, R. J.: Bremsstrahlung, Synchrotron Radiation, and Compton Scattering of High-Energy Electrons Traversing Dilute Gases, Reviews of Modern Physics, 42, 237–271, 10.1103/RevModPhys.42.237, 1970\.
  * Böttcher and Chiang (2002) Böttcher, M. and Chiang, J.: X-Ray Spectral Variability Signatures of Flares in BL Lacertae Objects, Astrophys. J., 581, 127–142, 10.1086/344155, 2002.
  * Chiang and Böttcher (2002) Chiang, J. and Böttcher, M.: Synchrotron and Synchrotron Self-Compton Spectral Signatures and Blazar Emission Models, Astrophys. J., 564, 92–96, 10.1086/324294, 2002.
  * Ellison et al. (1990) Ellison, D. C., Reynolds, S. P., and Jones, F. C.: First-order Fermi particle acceleration by relativistic shocks, Astrophys. J., 360, 702–714, 10.1086/169156, 1990.
  * Felten and Morrison (1966) Felten, J. E. and Morrison, P.: Omnidirectional Inverse Compton and Synchrotron Radiation from Cosmic Distributions of Fast Electrons and Thermal Photons, Astrophys. J., 146, 686–708, 1966.
  * Ginzburg and Syrovatskii (1969) Ginzburg, V. L. and Syrovatskii, S. I.: Developments in the Theory of Synchrotron Radiation and its Reabsorption, Ann. Rev. Astron. Astrophys., 7, 375–+, 10.1146/annurev.aa.07.090169.002111, 1969.
  * H. E. S. S. and Fermi-LAT collaborations (2009) H. E. S. S. and Fermi-LAT collaborations: Simultaneous observations of PKS 2155-304 with H.E.S.S., Fermi, RXTE and ATOM: spectral energy distributions and variability in a low state, ArXiv e-prints, 2009.
  * Jauch and Rohrlich (1976) Jauch, J. M. and Rohrlich, F.: The theory of photons and electrons. The relativistic quantum field theory of charged particles with spin one-half, Texts and Monographs in Physics, New York: Springer, 1976, 2nd ed., 1976.
  * Jones (1968) Jones, F. C.: Calculated Spectrum of Inverse-Compton-Scattered Photons, Physical Review, 167, 1159–1169, 10.1103/PhysRev.167.1159, 1968.
  * Kardashev (1962) Kardashev, N. S.: Nonstationariness of spectra of young sources of nonthermal radio emission, Soviet Astronomy - AJ, 39, 393–409, 1962.
  * Kataoka et al. (2000) Kataoka, J., Takahashi, T., Makino, F., Inoue, S., Madejski, G. M., Tashiro, M., Urry, C. M., and Kubo, H.: Variability Pattern and the Spectral Evolution of the BL Lacertae Object PKS 2155-304, Astrophys. J., 528, 243–253, 10.1086/308154, 2000.
  * Katarzyński et al. (2006) Katarzyński, K., Ghisellini, G., Mastichiadis, A., Tavecchio, F., and Maraschi, L.: Stochastic particle acceleration and synchrotron self-Compton radiation in TeV blazars, 453, 47–56, 10.1051/0004-6361:20054176, 2006.
  * Kirk et al. (1998) Kirk, J. G., Rieger, F. M., and Mastichiadis, A.: Particle acceleration and synchrotron emission in blazar jets, 333, 452–458, 1998.
  * Mannheim (1993) Mannheim, K.: The proton blazar, 269, 67–76, 1993.
  * Pohl (2002) Pohl, M.: Einfuehrung in die Hochenergieastrophysik, Shaker Verlag, Aachen, 2002\.
  * Press (2002) Press, W. H.: Numerical recipes in C++ : the art of scientific computing, Numerical recipes in C++ : the art of scientific computing by William H. Press. xxviii, 1,002 p. : ill. ; 26 cm. Includes bibliographical references and index. ISBN : 0521750334, 2002.
  * Primack et al. (2005) Primack, J. R., Bullock, J. S., and Somerville, R. S.: Observational Gamma-ray Cosmology, in: High Energy Gamma-Ray Astronomy, edited by F. A. Aharonian, H. J. Völk, & D. Horns, vol. 745 of American Institute of Physics Conference Series, pp. 23–33, 10.1063/1.1878394, 2005.
  * Schlickeiser (1984) Schlickeiser, R.: An explanation of abrupt cutoffs in the optical-infrared spectra of non-thermal sources. A new pile-up mechanism for relativistic electron spectra, 136, 227–236, 1984.
  * Schlickeiser (2002) Schlickeiser, R.: Cosmic ray astrophysics, Astronomy and Astrophysics Library; Physics and Astronomy Online Library. Berlin: Springer. ISBN 3-540-66465-3, 2002, XV + 519 pp., 2002.
  * Sokolov et al. (2004) Sokolov, A., Marscher, A. P., and McHardy, I. M.: Synchrotron Self-Compton Model for Rapid Nonthermal Flares in Blazars with Frequency-dependent Time Lags, Astrophys. J., 613, 725–746, 10.1086/423165, 2004.