# Statistical analysis of intermittency and its association with proton
heating in the near Sun environment.

Nikos Sioulas Earth, Planetary, and Space Sciences, University of California,
Los Angeles
Los Angeles, CA 90095, USA  Marco Velli Earth, Planetary, and Space Sciences,
University of California, Los Angeles
Los Angeles, CA 90095, USA  Rohit Chhiber Department of Physics and Astronomy,
Bartol Research Institute, University of Delaware
Newark, DE 19716, USA  NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  Loukas Vlahos Deprtment of Physics, Aristotle
University of Thessaloniki
GR-52124 Thessaloniki, Greece  William H. Matthaeus Department of Physics and
Astronomy, Bartol Research Institute, University of Delaware
Newark, DE 19716, USA  Riddhi Bandyopadhyay Department of Astrophysical
Sciences, Princeton University,
Princeton, NJ 08544, USA  Manuel E. Cuesta Department of Physics and
Astronomy, Bartol Research Institute, University of Delaware
Newark, DE 19716, USA  Chen Shi Earth, Planetary, and Space Sciences,
University of California, Los Angeles
Los Angeles, CA 90095, USA  Trevor A. Bowen Space Sciences Laboratory,
University of California,
Berkeley, CA 94720-7450, USA  Ramiz A. Qudsi Department of Physics and
Astronomy, Bartol Research Institute, University of Delaware,
Newark, DE 19716, USA  Michael L. Stevens Harvard-Smithsonian Center for
Astrophysics,
Cambridge, MA 02138, USA,  Stuart D. Bale  Space Sciences Laboratory,
University of California,
Berkeley, CA 94720-7450, USA  Physics Department, University of California,
Berkeley, CA 94720-7300, USA

###### Abstract

We use data from the first six encounters of Parker Solar Probe and employ the
Partial Variance of Increments ($PVI$) method to study the statistical
properties of coherent structures in the inner heliosphere with the aim of
exploring physical connections between magnetic field intermittency and
observable consequences such as plasma heating and turbulence dissipation. Our
results support proton heating localized in the vicinity of, and strongly
correlated with, magnetic structures characterized by $PVI\geq 1$. We show
that on average, such events constitute $\approx 19\%$ of the dataset, though
variations may occur depending on the plasma parameters. We show that the
waiting time distribution ($WT$) of identified events is consistent across all
six encounters following a power-law scaling at lower $WTs$. This result
indicates that coherent structures are not evenly distributed in the solar
wind but rather tend to be tightly correlated and form clusters. We observe
that the strongest magnetic discontinuities, $PVI\geq 6$, usually associated
with reconnection exhausts, are sites where magnetic energy is locally
dissipated in proton heating and are associated with the most abrupt changes
in proton temperature. However, due to the scarcity of such events, their
relative contribution to energy dissipation is minor. Taking clustering
effects into consideration, we show that smaller scale, more frequent
structures with PVI between, $1\lesssim PVI\lesssim 6$, play the major role in
magnetic energy dissipation. The number density of such events is strongly
associated with the global solar wind temperature, with denser intervals being
associated with higher $T_{p}$.

Parker Solar Probe, proton heating, solar wind, turbulence, intermittency

##  1 Introduction

The solar wind is a strongly magnetized, weakly collisional stream of charged
particles expanding at supersonic speeds from the outermost layer of the sun,
the corona (Parker, 1958). For coronal temperatures of order $\sim 100\ eV$
and adiabatic cooling, the temperature at 1 AU would be expected to decrease
to below $\sim 1\ eV$. Observations of the solar wind, however, show values in
the order of $\sim 10\ eV$ at 1 AU (Gazis et al., 1994; Richardson et al.,
1995) and ion temperature decays as a function of radial distance as $T\sim
r^{-\gamma}$, with the exponent, attaining values in the range,
$0.5\lesssim\gamma\lesssim 1$ (Richardson et al., 1995; Stansby et al., 2018),
a much slower rate than spherically symmetric adiabatic expansion models would
predict (i.e. $T\sim r^{-4/3})$. A complete description of the dynamics of the
solar wind must therefore include non-adiabatic heating processes that
contribute to its internal energy. Most theoretical investigations suggest
collisionless dissipation mechanisms that can be broadly classified into two
categories: wave-particle interactions and current sheets/reconnection (Velli
et al., 1989; Matthaeus et al., 1999; Cranmer & van Ballegooijen, 2005;
Cranmer et al., 2007; Lionello et al., 2014; Isliker et al., 2019; Sioulas et
al., 2020a, b; Perez et al., 2021; González et al., 2021). Despite all these
efforts, the non-adiabatic expansion of the solar wind remains one of the
major outstanding problems in the field of space physics.

The turbulent cascade in space plasma systems is a multi-scale process during
which the energy entering the system at large scales is removed from the
turbulent cascade at the scale of the ion gyroradius (kinetic scales) due to
non-linear interactions among the fluctuations (Biskamp, 2003; Coleman, 1968).
This is not the complete story though as, both simulations (Matthaeus &
Montgomery, 1980; Meneguzzi et al., 1981) and observations (Marsch & Tu, 1990;
Chhiber et al., 2020) indicate that turbulence also produces intermittency. In
the solar wind, fluctuations that extend over several decades in the frequency
domain, from the inverse of the solar rotation period down to the electron
gyrofrequency, abound. In this sense, the solar wind is regarded as the
exquisite turbulent plasma laboratory, as it offers the chance to study
coherent structures, such as current sheets and vortices, being dynamically
generated as the end result of the turbulent cascade (Veltri, 1999; Greco et
al., 2010).

Coherent structures and spatial intermittency arise from the non-self-
similarity of the nonlinear cascade, i.e. the probability distribution
functions of fluctuations of a field $\phi$,
$\Delta\phi=\phi(\ell+\Delta\ell)-\phi(\ell),$ at a given length , $\ell$,
display larger and larger tails with respect to a Gaussian distribution as the
lengths become smaller and smaller (Castaing et al., 1990). This behavior,
often described as multifractal, is associated with the occurrence at small
scales of coherent structures, superposed on a background of random
fluctuations. Indeed such structures that persist in time longer than the
surrounding stochastic fluctuations are characterized by non-Gaussian
statistics and constitute only a minor fraction of the entire dataset (Osman
et al., 2012; Bruno, 2019). They nevertheless strongly influence the
dissipation, heating, transport, and acceleration of charged particles
(Matthaeus & Velli, 2011; Karimabadi et al., 2013; Tessein et al., 2013;
Bandyopadhyay et al., 2020a). In recent years, several numerical simulations
(Biskamp & Müller, 2000; Parashar et al., 2009; Lehe et al., 2009; Servidio et
al., 2012; Jain et al., 2021; Sioulas et al., 2022) have suggested localized,
intermittent heating in structures such as reconnection sites and current
sheets.

Additionally, observational studies have indicated a statistical link between
coherent magnetic field structures and elevated ion temperatures, both in the
near-Earth solar wind (Leamon et al., 2000; Osman et al., 2012; Yordanova et
al., 2021), as well as in the near-Sun (Qudsi et al., 2020) environment. As
Parker Solar Probe (PSP) (Fox et al., 2016) moves closer to the solar wind
sources, it offers an unprecedented opportunity to study intermittency and the
physics of dissipation in the near sun solar wind environment. In this
article, we make use of the PVI method and take advantage of the high-
resolution magnetic field and particle PSP data from the first six encounters
to perform a statistical study of intermittency. Furthermore, the relationship
between solar wind turbulence and the concept of intermittent heating in the
nascent solar wind environment is investigated.

The outline of this work is as follows. Sections 2, 3 define the PVI and LET
methods respectively, and provide some background; Section 4 presents the
selection of PSP data and their processing. The results of our analysis are
presented in 5; in subsection 5.1 the statistics of discontinuities in the
magnetic field obtained by the PVI method are analyzed, while in subsection
5.2, the effects of the identified discontinuities in the proton temperature
of the solar wind are investigated. A summary of the results along with the
conclusions is finally given in Section 6.

##  2 Partial Variance of Increments (PVI)

When confronted with a dataset that samples a turbulent space plasma system,
an important task is to find the subset of the data that corresponds to the
underlying coherent structures. In recent years, a plethora of methods have
been suggested for the identification of intermittent structures and
discontinuities in the magnetic field. Some of those include the Phase
Coherence Index method (Hada et al., 2003) and the wavelet-based Local
Intermittency Measure (LIM) (Bruno et al., 1999). In this paper, we employ a
simple and well-studied method that has been effectively used in the past for
the study of intermittent turbulence and the identification of coherent
structures, both in simulations (Servidio et al., 2011), and observations
(Chasapis et al., 2018; Hou et al., 2021), namely, the Partial Variance of
Increments (PVI). The advantage of the PVI method is that it provides an easy-
to-implement tool that measures the sharpness of a signal relative to the
neighborhood of that point. For lag $\tau$, the normalized $PVI$ at time t is
defined as:

$PVI(t,\tau)\ =\
\frac{|\Delta\textbf{B}(t,\tau)|}{\sqrt{\langle|\Delta\textbf{B}(t,\tau)|^{2}\rangle}},$
(1)

where, $|\Delta\textbf{B}(t,\tau)|\ =\ |\textbf{B}(t\ +\ \tau)\ -\
\textbf{B}(t)|$ is the magnitude of the magnetic field vector increments, and
$\langle...\rangle$ denotes average over a window that is a multiple of the
estimated correlation time. PVI is a threshold method, so to proceed with the
analysis, one can impose a threshold, ${\theta}$ on the PVI and select
portions in which $PVI>{\theta}$. One could easily see that
${\langle}PVI^{4}{\rangle}$ is related to the kurtosis of the magnetic field,
and thus a possible way to calibrate the PVI method is to exclude all the
values above ${\theta}$ and compute the kurtosis of the remaining data. The
process can be repeated several times, each time lowering the value of
${\theta}$ until the kurtosis of the remaining signal is equal to the value
expected for Gaussian random signals. Greco et al. (2018) have shown that
increments with $PVI>3$ lay on the “heavy tails” observed in the distribution
of increments and can thus be associated with Non-Gaussian structures. By
increasing the threshold value $\theta$, one can thus identify the most
intense magnetic field discontinuities like current sheets and reconnection
sites. Finally, note that the method is insensitive to the mechanism that
generates the coherent structures. This means that the PVI can be implemented
for the identification of any form of sharp gradients in the vector magnetic
field. A more comprehensive review of PVI, as well as a comparison with
aforementioned methods, appropriate for identifying discontinuities, can be
found in Greco et al. (2018).

##  3 Local Energy Transfer (LET)

As already noted in Sec. 1, coherent structures constitute a minor fraction of
the dataset. At the same time coherent structures correspond to sites where
significant heating and dissipation takes place, contributing a
disproportionate amount to the internal energy of the solar wind. Once the
subset of the data-set corresponding to coherent structures has been
identified, a method that can allow us to quantify their contribution to
energy dissipation and heating of the solar wind is necessary. The Kolmogorov-
Yaglom law, extended to isotropic MHD (Politano & Pouquet, 1998a, b) (PP98)
has been used to evaluate the average energy transfer rate per unit mass, at
scale $\Delta t$. Taking into account the conservation law of the respective
inviscid invariants, the scaling law can be directly derived from the
dynamical MHD equations for the third-order moment. Assuming homogeneity,
scale separation, isotropy, and time-stationarity PP98 yields the linear
scaling of the mixed third-order moment of the field fluctuations

Figure 1: (a) PDFs of $PVI$ for lag, $\tau\ =0.837s$ (b) Fraction of PVI
events exceeding a given threshold $\theta$ for $E_{1}-E_{6}$.

$Y(\Delta t)=\langle\Delta
v_{r}(|\Delta\textbf{v}|^{2}+|\Delta\textbf{b}|^{2})-\Delta
b_{r}(\Delta\textbf{v}\cdot\Delta\textbf{b})\rangle\\\
=-\frac{4}{3}\langle\epsilon\rangle V_{SW}\Delta t,$ (2)

where, $\langle...\rangle$ indicates spatial averaging, and
$\langle\epsilon\rangle$ is the mean energy transfer rate. The mixed third-
order moment is computed using the scale-dependent increments of the plasma
velocity field v, and the magnetic field given in velocity units,
$\textbf{b}=\frac{\textbf{B}}{\sqrt{\mu_{0}m_{p}n_{p}}}$, where $\mu_{0}$ is
the magnetic permeability of vacuum, $m_{p}$ is the proton mass, and $n_{p}$
is the number density of protons. The subscript r indicates the longitudinal
component, i.e., measured in the sampling direction. Note, that the validity
of the Taylor’s hypothesis (Taylor, 1938), has been assumed for converting
spatial to temporal scales through

$\ell=V_{SW}\Delta t,$ (3)

The linear relation in Eq. 2 has been confirmed under various plasma
conditions revealing the turbulent dynamics and providing an estimate of the
energy transfer rate of space plasmas throughout the heliosphere (Sorriso-
Valvo et al., 2018, 2019; Bandyopadhyay et al., 2018, 2020b, 2021; Hernández
et al., 2021). In recent years, it was suggested (Sorriso-Valvo et al., 2018,
2019) that by omitting the averaging process in Eq. 2, a proxy of the local
energy transfer rates (LET) at a given scale $\Delta t$ can be obtained as

$\small{\epsilon^{\pm}(t,\Delta t)=-\frac{3}{4}\frac{\Delta
v_{r}(|\Delta\textbf{v}|^{2}+|\Delta\textbf{b}|^{2})-\Delta
b_{r}(\Delta\textbf{v}\cdot\Delta\textbf{b})}{V_{SW}\Delta t}}.$ (4)

The LET is composed of two additive terms. The first term,
$\epsilon^{e}=-3/(4V_{SW}\Delta t)[\Delta
v_{r}(|\Delta\textbf{v}|^{2}+|\Delta\textbf{b}|^{2})$], is associated with the
magnetic and kinetic energy advected by the velocity fluctuations. The second
term, $\epsilon^{c}=$-$3/(4V_{SW}\Delta t)[\Delta
b_{r}(\Delta\textbf{v}\cdot\Delta\textbf{b})]$, is associated with the cross-
helicity coupled to the longitudinal magnetic fluctuations. Note that several
important scaling contributions are neglected in this procedure of estimating
LET. Such contributions in Eq. 2, are suppressed by averaging over a large
sample. Consequently, this approach provides only a local proxy, indicating
the contribution to the local energy transfer rate (Kuzzay et al., 2019)

Figure 2: PDFs of WT between (a) $PVI>3$, (b) $PVI>6$ events for lag
$\tau=0.837~{}{\rm seconds}$.

##  4 DATA AND ANALYSIS PROCEDURES

We analyze data from first six encounters ($E_{1}-E_{6}$) of PSP from 2018 to
2020, covering heliocentric distances $0.1\lesssim R\lesssim 0.25$ au. We use
magnetic field data from the FIELDS fluxgate magnetometers Bale et al. (2016).
To estimate the PVI time-series at time-lags, $\tau\ =\ 0.837~{}{\rm seconds}$
magnetic field data have been resampled to a cadence of $0.837~{}{\rm
seconds}$ using linear interpolation. As outlined in Sec 2, in order to
compute the variance, a moving average over a window that is a multiple of the
correlation time is required. The correlation time can be estimated using the
e-folding technique by considering the time it takes for the autocorrelation
function to drop to $e^{-1}$ of its maximum value (Matthaeus & Goldstein,
1982; Krishna Jagarlamudi et al., 2019). For encounters $E_{1}-E_{6}$, the
correlation time was estimated to be between $500\lesssim t\lesssim
2000~{}{\rm seconds}$. Accordingly, we perform the ensemble averaging over a
window of 8 hours for all six encounters. Several different averaging windows,
ranging $2-12~{}{\rm hours}$ have been implemented, with qualitatively similar
results on our analysis. Plasma data from Solar Probe Cup (SPC), part of the
Solar Wind Electron, Alpha and Proton (SWEAP) suite (Kasper et al., 2016) have
also been analyzed to obtain the bulk velocity and radial proton
temperature/thermal speed measurements at $\sim 0.837~{}{\rm seconds}$
resolution. The radial temperature and bulk velocity time-series have been
pre-processed to eliminate spurious spikes and outliers using the Hampel
filter (Davies & Gather, 1993). An important effect that should be taken into
account when analyzing solar wind particle data is the anisotropy in the
parallel ($T_{||}$) and perpendicular ($T_{\perp}$) temperatures with respect
to the background magnetic field (Huang et al., 2020; Hellinger et al., 2011).

Figure 3: An example indicating clustering of intermittent structures
associated with increased proton temperature. The shaded areas point to the
location of $PVI\geq 3$ events. From top to bottom, the magnitude of the
magnetic field $|B|$, the radial component of the magnetic field ($B_{R}$),
the tangential and normal components of the magnetic field ($B_{T}$) and
($B_{N}$), in green and red respectively, the PVI time-series for lag,
$\tau=0.837~{}{\rm seconds}$ the temperature $T$, the radial component of the
proton bulk velocity $V_{R}$, the tangential and normal components of the
proton bulk velocity ($V_{T}$), and proton number density ($n$) are shown.

SPC was designed to sample the radial temperature. Consequently, to ensure
that the temperature we were measuring was indeed radial temperature all
intervals for which the angle between $B_{r}\hat{r}$ and B was greater than
thirty degrees, $\theta_{B_{r}\hat{r},\textbf{B}}\geq 30$ have been discarded.
Finally, intervals for which the solar wind was outside of the field of view
of SPC, have also been manually discarded.

##  5 Results

###  5.1 Statistical properties of intermittent structures.

One of the major questions one has to address when studying intermittency is
the nature of the physical processes that initiate the coherent structure
production at the origin. Accordingly, in order to gain insight into the
statistics of intermittent magnetic structures in the solar wind, we follow
the process described in Sec. 4 to estimate the PVI time-series for time-lag,
${\tau}\ =\ 0.837~{}{\rm seconds}$. In Fig. 1a we show the Probability Density
Functions (PDFs) of PVI values. The most probable value is close to $\sim\
0.3$, indicating that the majority of the detected events can be characterized
as non-intermittent. In Fig. 1b, the fraction of the entire data-set occupied
by PVI events exceeding the threshold $PVI>\theta$, is shown for
$E_{1}-E_{6}$. For reasons that will become obvious in Sec 5.2, the fraction
of events characterized by a PVI index greater than unity $f_{PVI\geq 1}$ is
also shown. The fraction is consistent for $E_{1}-E_{4}$ attaining a value of
$f_{PVI\geq 1}\approx 19.1\%$, but decreases to $f_{PVI\geq 1}\approx 18\%$
for $E_{5},E_{6}$. Beyond this point, the fraction of the dataset occupied by
coherent structures characterized by higher values of PVI decreases rapidly.

Figure 4: Binned average of radial proton temperature plotted against PVI
(blue) along with the number of points per bin (red) for the first six
encounters of PSP. Notice a rough upward trend in mean temperature at higher
PVI bins.

Another method that can provide an insight into the statistics of the solar
wind coherent structures is the Waiting Time (WT) distribution analysis. In
the case of the PVI time-series, we define the waiting time as the time passed
between the end and the start of two subsequent events for which the value of
PVI stays above a threshold $\theta$. Note, that PVI events have a finite
duration. Therefore, subsequent times for which the PVI time-series stay above
the threshold are considered as part of the same event. Also note that in
order to maintain an adequate sample size, we have imposed a restriction on
the minimum counts per bin. Consequently, bins with fewer than ten counts have
been discarded. The waiting time interval is itself a new random variable, the
distribution of which, is independent of the original random variable’s
distribution. A simple inspection of the distribution shape can then reveal
whether or not the underlying mechanism can be classified as a random
Poissonian-type of process, or it possesses “memory” indicating strong
correlation and clustering. In the first case, the distribution is better
described by an exponential, while in the latter distribution, scales like a
power-law (Greco et al., 2009, 2010). Figs. 2a,b, show the PDF’s of WT between
intermittent PVI events with lag $\tau\ =0.837~{}{\rm seconds}$ and threshold
${\theta}_{1}=3$, ${\theta}_{2}=6$, respectively, for the first six
encounters. At lower WT’s, the best fit analysis indicates that WT
distributions are better described by a power-law. The index of the power-law
fit attains values in the range $a\ \in[-0.8,\ -0.52]$, progressively getting
softer as the threshold value, $\theta$, increases. In contrast, for events
further apart in time, the distribution is better described by an exponential.
The change between power-law and exponential in the distribution can be
interpreted as the breaking-point between $intercluster$ and $intracluster$
waiting times Greco et al. (2010). This change seems to coincide with an ill-
posed, due to the power-law nature of the distribution, mean value of WT,
$\langle WT\rangle$ Chhiber et al. (2020). The WT distribution analysis was
repeated by estimating the PVI time-series using a different time-lag, $\tau\
=\ 8.37,\ 83.7~{}{\rm seconds}$, still sampling though, overtimes with
$0.837~{}{\rm seconds}$ cadence. The resulting waiting time distributions (not
shown here) are once again remarkably similar between all six encounters and
thus resemble the ones reported in Chhiber et al. (2020) for $E_{1}$. This
provides a strong indication that high PVI valued coherent structures, are not
evenly distributed within the solar-wind, but rather tend to be strongly
correlated and form clusters, creating alternating regions of very low and
very high magnetic field activity respectively (see also (Dudok de Wit et al.,
2020; Chhiber et al., 2020; Bale et al., 2021; Yordanova et al., 2021)). One
such example, out of the thousands of clustering events we were able to
recover is presented in Fig. 3. To emphasize the clustering of coherent
structures, blue, vertical lines have been added to indicate the location of
events characterized by a magnetic PVI index, $PVI\geq 3$. Notice the elevated
proton temperature at regions where coherent structures abound, an effect that
is discussed in detail in Sec. 5.2.

Figure 5: Average $T_{p}$ conditioned on the temporal separation from PVI
events that exceed a PVI threshold. The black dashed line indicates the mean
proton temperature ($\tilde{T}_{p}$) for each encounter.  Figure 6: Superposed
average of $T_{p}$ conditioned on the temporal separation from PVI events that
exceed a PVI threshold for $E_{1}-E_{6}$

###  5.2 Intermittent heating of the young solar wind.

The purpose of this section is to investigate the idea that local coherent
structures contribute to the heating of the solar wind. For this reason, we
follow the method described in Sec. 4 and carry out several diagnostics on the
correlation between the PVI and proton temperature ($T_{p}$) time-series. As a
first step, we use binned statistics to interpret $T_{p}$ as a function of
PVI. In Fig. 4, we use 600 PVI bins, to present the average proton temperature
per bin (i.e., $\langle T_{p}(\theta_{i}\leq PVI\leq\theta_{i+1})\rangle$,
where $\theta_{i}$ is the PVI threshold) plotted against the center of the
bin. Uncertainty bars are also shown as vertical lines indicating the standard
error of the sample (Gurland & Tripathi, 1971). The uncertainty is thus
estimated as $\sigma_{i}/\sqrt{n}$, where $\sigma_{i}$ is the standard
deviation of the samples inside the bin. The number of samples included in
each PVI bin is also shown as a dotted red line. For all encounters, a
statistically significant, positive correlation between $T_{p}$ and PVI can be
observed for $PVI\leq 3$. Beyond this point, the trend is eroded for $E_{5}$.
For the remaining encounters, the trend extends, with good statistical
significance to $PVI\approx 6$. For greater values of PVI, $6\leq PVI\leq 10$,
though on average the temperature increases with increasing PVI, there is also
a high degree of variability that could most probably be attributed to the low
number of recorded events, as indicated by the red dotted line. For most
orbits, the local temperature peaks at $PVI\sim 10$. Between the lowest and
highest PVI values, the difference in mean temperature is of the order of
$\sim 10^{5}K$, but, can range as high as $\sim 2\cdot 10^{5}K$ (e.g.,
$E_{1},\ E_{2},\ and\ E_{4}$). Note, however, the high degree of variability
in $\tilde{T_{p}}$ for bins in the $PVI\geq 6$ range. Although bins of
relatively high average temperature are the rule, several bins of low
temperature, comparable to one of the bins of $PVI\approx 2$, can be observed.

Figure 7:  Lagged cross correlation between $T_{p}^{{}^{\prime}}$ and
$PVI^{{}^{\prime}}$ for different thresholds in the PVI timeseries. The
correlation coefficient at zero lag is also shown. Note that the correlation
functions have been shifted in the vertical direction for clarity.

To further elucidate the relationship between $T_{p}$, and magnetic field
discontinuities, we estimate averages of $T_{p}$ constrained on the temporal
separation from PVI events that lay in a given PVI bin. This can be expressed
as (Osman et al., 2012; Tessein et al., 2013; Sorriso-Valvo et al., 2018):

$\tilde{T}_{p}(\Delta t,\ \theta_{i},\ \theta_{i+1})\ =\\\ \langle
T_{p}(t_{PVI}\ +\Delta t)|\theta_{i}\leq PVI\leq\theta_{i+1}\rangle.$ (5)

Here, $\Delta t$ is the temporal lag relative to the location of the main PVI
event taking place at time $t_{PVI}$, and $\theta\ =\ [0,\ 1,\ 2,\ 4,\ 6,\ 8,\
10]$. The conditional average of temperature at different lags is presented in
Fig. 5 for six different PVI bins. The dip in the mean temperature at lag
equal to $t\ =0s$, for $PVI\leq 1$ (red line), suggests that no significant
heating of the solar wind occurs at times where the magnetic field is
relatively smooth. This observation is reinforced by the fact, that in such
regions, the plasma temperature is lower than the mean temperature of the
respective encounter, indicated by the black dashed line. On the other hand,
by increasing the threshold value $\theta$, we can make several observations
regarding the nature of intermittent dissipation in the solar wind. For
encounters $E_{1}-E_{4},\&\ E_{6}$, and $\theta\geq 1$, we can observe a
global maximum in mean $T_{p}$ close to zero lag, indicating a higher
probability of increased proton temperature in the vicinity of coherent
structures. The peak is then followed by a gradual decrease in mean $T_{p}$ as
we move further away from the main discontinuity. The fact that $T_{p}$ stays
elevated near the main event can, most probably, be attributed to the
clustering of coherent structures noted ins Sec. 5.1. The rate of decrease is
distinct for each bin, with the steepest gradients in $T_{p}$ observed around
the sharpest discontinuities. This is better illustrated in Fig. 6, where a
superposed average of $E_{1}-E_{6}$ is presented.

Figure 8:  PDF’s of $T_{p}$ conditioned on the percentage fraction of coherent
structures with $PVI\geq 1$, $f_{PVI\leq 1}$ for the first six encounters.
Vertical lines indicate the median temperature for each of the PDFs. As the
percentage of high PVI valued events increases, we can see a significant
increase in the probability density for higher proton temperatures. Note that
PDFs don’t change when we remove $PVI\geq 6$ events (black, dashed line).

Figure 9:  Fraction of coherent structures of PVI index greater than $\theta$,
$f_{PVI\geq\theta}$, where (a) $\theta$ = 4, (b) $\theta$ = 5, (c) $\theta$= 6
as a function of the fraction of coherent structures of PVI index greater than
one $f_{PVI\geq 1}$ and mean proton temperature of the 5-minute interval
$\tilde{T_{p}}$. The data points were binned according to $f_{PVI\geq\theta}$
and $f_{PVI\geq 1}$, and the mean value inside each bin was calculated, which
is reflected in the color. The numbers in the plot indicate the number of
5-minute intervals inside each bin. Bins with 0 counts are shown in dark grey.

On the contrary, $E_{5}$ exhibits a distinct behavior. The mean proton
temperature is $T_{p}\ \sim 3\cdot 10^{5}K$ marking the lowest observed proton
temperature between all encounters, and the average $T_{p}$ peaks at lower
thresholds of PVI. These results motivated us to perform a more thorough
investigation of $E_{5}$. A possible explanation for the observed discrepancy
during E5 is the crossing of PSP through the Heliospheric Current Sheet
($HCS$). Several statistical studies using superposed epoch analysis suggest
that the HCS is characterized by a local minimum in proton temperature, and a
local maximum in the proton density (Borrini et al., 1981; Suess et al., 2009;
Liou & Wu, 2021; Shi et al., 2022). As a result of the extended intervals
characterized by an increased proton density in the vicinity of the HCS,
transient mechanisms such as adiabatic proton heating can obscure the effects
of non-adiabatic, intermittent heating. Moreover, recent observations (Chen et
al., 2021), support a poorly developed turbulence in the proximity of the HCS,
indicated by the relatively lower amplitude of the turbulent fluctuations.
These results are corroborated by the fact that $E_{5}$, is characterized by a
relatively low coherent structure density. More specifically, the precise
value for the fraction of the data-set occupied by events of $PVI\geq 1$, was
estimated at $f_{PVI\geq 1}\approx 18\%$, namely, $6.1\%$ lower than the mean
value of $\langle f_{PVI\geq 1}\rangle=19.1\%$ for the rest of the encounters.
The poorly developed turbulence in the vicinity of the HCS could act as a
contributing factor to the ineffective heating of the protons in that region.
Note, however that other features of the solar wind, including the Alfvénic
nature of the fluctuations (e.g. through cross-helicity) (Stansby et al.,
2019); solar wind speed (Burlaga & Ogilvie, 1970; Lopez & Freeman, 1986; Shi
et al., 2021); plasma beta (Vech et al., 2017; Huang et al., 2020); the
presence of waves (Howes et al., 2012; He et al., 2015; Mozer et al., 2021)
have been shown to be closely correlated with proton temperature. Taking all
the aforementioned effects into account, we can understand that the
intermittent heating in the solar wind, could be obscured by other transient
heating mechanisms that are strongly dependent on the nature of the solar wind
under study, and could explain the discrepancy for $E_{5}$.

The results presented up to this point provide strong evidence for the
intermittent heating of the solar wind. However, for a more complete
understanding, it would be convenient to quantify the relative contribution of
different types of coherent structures, as indicated by the PVI index, to the
solar wind’s internal energy. To do this, a careful analysis that allows for
clustering effects, as well as transient effects (e.g. HCS crossing) to be
subtracted is required. One such analysis involves the study of lagged cross-
correlation between $T_{p}$ and the PVI time-series, as it offers the chance
to infer whether or not the two quantities are changing at the same time.
Here, the lag refers to how far the series is offset and its sign determines
which series is shifted, in this case, PVI being the lagged quantity. The goal
is to correlate the fluctuating parts of PVI and $T_{p}$. We thus have to
subtract the mean

$PVI^{{}^{\prime}}\ =\ PVI(t)\ -\ \langle PVI\rangle,$ (6)
$T_{p}^{{}^{\prime}}\ =\ T_{p}(t)\ -\ \langle T_{p}\rangle,$ (7)

where $\langle...\rangle$ denotes the time average over the entire data set.
We can therefore estimate the lagged cross-correlation as

$Corr(T_{p}^{{}^{\prime}},PVI^{{}^{\prime}})\ =\ \frac{\langle
PVI^{{}^{\prime}}(t\ +\ \tau)\cdot
T_{p}^{{}^{\prime}}(t)\rangle}{\sqrt{\sigma_{PVI}\cdot\sigma_{T_{p}}}},$ (8)

where $\sigma$ is the variance of the subscripted quantity. In Fig. 7, the
lagged cross-correlation between PVI and temperature, for different $PVI$ bins
is presented. Note that each curve has been shifted vertically by $0.2\cdot
i$, where, $i\ =\ 0,1,..,5$ for better visualization.

For $PVI\leq 1$, we can observe a minimum of the correlation coefficient at
zero lag providing an additional confirmation that no temperature enhancements
occur in such regions. On the other hand, for higher thresholds, there is a
trend, with the correlation coefficient attaining lower values as we move
further away from the discontinuity. In general, the maximum value of the
coefficient tends to peak at higher values with the increase of PVI. In many
cases, though, the three highest PVI bins peak at comparable values. For
$E_{5}$, even though the correlation coefficient stays relatively low, peaks
still emerge at the location of the principal discontinuity, at zero lag. We
can thus understand that the strongest discontinuities in the magnetic field
are related to the most abrupt changes in temperature, and thus play a prime
role in determining the local plasma dynamics. Finally, note the emergence of
local peaks at some distance from the main discontinuity, becoming more
pronounced as we move to higher PVI thresholds. We can attribute this
observation to the combination of two factors. The first factor is the
relatively low frequency of high PVI valued events. Due to the scarcity of
such events, around which we perform the conditional averaging, the number of
samples gets progressively lower as the PVI threshold increases. For example,
for PVI index, $PVI\geq 8$, the number of recorded events for $E_{1}-E_{6}$
ranges between 140 to 277. Consequently, even though there are enough events
at the largest PVI bins to safely determine the mean and variance in the
correlation coefficient, single strong events of PVI index, $PVI\geq 6$, shown
to be highly correlated to local temperature changes, located at a given
distance from the main discontinuity might be controlling the correlation. The
second factor is the clustering of coherent structures reported in Sec.5.1.
For the clustering effects, it would be important to also take into account
events of PVI index, $PVI\geq 6$, with the number of recorded events for
$E_{1}-E_{6}$ ranging between 456 to 820. As shown in Fig. 2b, the WT
distribution for $PVI\geq 6$ is described by a power-law at lower waiting
times indicating that such events tend to be highly correlated and form
clusters. Consequently, considering clustering effects, we can understand that
intense events taking place close to the main event can strongly affect the
average value of the mean $T_{p}$ at certain lags, resulting in the local
peaks observed in Fig. 7.

Finally, we move on to study the macroscopic effects of intermittent heating.
For this reason, we divide the data into $5\ minute$ sub-intervals, and for
each sub-interval, the fraction of the data that correspond to events of PVI
index greater than unity, as well as, the mean temperature of the sub-interval
are estimated. This allows us to estimate the PDFs of the mean proton
temperature, $T_{p}$, conditioned on the percentage fraction of $PVI\geq 1$
events within the sub-interval, $p(T_{p}|f_{PVI\geq 1}\leq\theta\%)$, where,
$\theta\ =\ 1,\ 10,\ 15,\ 20$. As shown Fig. 8, the probability density for a
higher mean temperature is increased, as the fraction of coherent structures
in the interval increases. Additionally, the median value of proton
temperature within each bin is estimated, shown as a vertical line of the same
color to the respective PDF. It is clear from Fig. 8 that for denser in
coherent structure intervals, the median of $T_{p}$ moves to higher values.
Notice that even when we exclude all events characterized by $PVI\geq 6$,
usually associated with reconnection exhausts, from our analysis, both the
PDFs and the median $T_{p}$ practically remain intact. Taking into account the
steep temperature gradients associated with events of high PVI index, shown in
Fig. 8, the results of Fig. 8 may at first seem contradictory. However, the
negligible contribution of events characterized by a high PVI index can be
explained in terms of the relatively low frequency of occurrence of such
events, as shown in Fig. 1b. This indicates, that as a result of their
scarcity, coherent structures characterized by a PVI index that exceeds a
given threshold, dissipate a negligible amount to the internal energy of the
solar wind. To better illustrate the relative contribution, of coherent
structures at a given threshold $\theta$, $PVI\geq\theta$, we present in Fig.
9a,b,c the fraction of coherent structures of PVI index greater than
$\theta=4,\ 5,\ 6$ respectively as a function of the fraction of coherent
structures of PVI index greater than unity, $f_{PVI\geq 1}$ and mean proton
temperature of the 5-minute interval $\tilde{T_{p}}$. In these figures, the
color indicates the mean proton temperature $T_{p}$ of the 5-minute interval,
while, the numbers in the plot indicate the number of 5-minute intervals
inside each bin. As shown in Fig. 9c, the vast majority of the intervals are
characterized by a very low number density, $f_{PVI\geq 6}$, as well as, a
clear increase in $T_{p}$ with increasing $f_{PVI\geq 1}$. On the other hand,
for $\theta\geq 6$ even though intervals with the highest observed temperature
are characterized by a relatively high number of $f_{PVI\geq 6}$, not a clear
trend with increasing $f_{PVI\geq 6}$ is observed. At the same time, such
intervals constitute only a small, almost negligible fraction of the entire
dataset. On the contrary, as shown Fig. 9a,b the number of intervals with a
high coherent structure number density and elevated $\tilde{T}_{p}$ steadily
increases as we lower the threshold to $\theta=4,5$, respectively.

Figure 10:  Sum of absolute values of LET, associated with PVI events of a
given PVI threshold, $\sum|\epsilon^{\pm}_{PVI\geq\theta}|$. The legend
indicates the estimates of the percentage fraction of total LET for coherent
structures at a PVI threshold,
$100\cdot\frac{\sum|\epsilon^{\pm}_{PVI\geq\theta}|}{\sum|\epsilon^{\pm}_{PVI\geq
0}|}$.

Finally, to further corroborate the previous results, and quantify the
relative contribution of coherent structures we take advantage of the LET
method (Sec. 3). In Fig. 10, we estimate the sum of local energy transfer
associated with coherent structures characterized by a PVI index that exceeds
a given threshold. Note, that LET is a signed quantity. Consequently, for the
purposes of this study, the absolute value of LET is considered (Sorriso-Valvo
et al., 2018). In agreement with (Osman et al., 2012), the strongest PVI
events, $PVI\geq 5$, contribute $\sim 10.5\ \%$ of the internal energy.
However, as the PVI threshold is increased to $PVI\geq 6$, the contribution of
coherent structures decreases to $\sim 4.8\ \%$. These results indicate that
in addition to the amount of energy dissipated per identified event, future
studies should also consider the frequency of occurrence of such events when
considering intermittent dissipation in the solar wind. Consequently, even
though, intense events of $PVI\geq 6$, usually linked to reconnection
exhausts, are the ones that determine the local plasma dynamics, due to the
scarcity of such high PVI valued coherent structures, it is rather the
accumulative effect of many small scale dissipation events that determine the
global temperature of the solar wind. More specifically, our results, suggest
that in the near sun solar wind environment the number density of coherent
structures characterized by a PVI index in the range $1\lesssim PVI\lesssim
6$, play the major role in magnetic energy dissipation and thus constitute an
important factor in determining the global temperature of the solar wind.

##  6 Summary and conclusions

In this work, we have analyzed magnetic field and particle data from the first
six encounters of the PSP mission. Our goal was to study the statistics of
intermittency and further elucidate the nature of turbulent dissipation in the
neighborhood of the solar wind sources. As a first step, an effort was made to
understand the nature of the mechanism that is responsible for the generation
of intermittency and coherent structures in the solar wind. We have shown that
coherent structures, corresponding to $PVI\geq 1$, constitute only $\approx
19\%$ of the dataset. As a follow-up to the study of Chhiber et al. (2020), we
studied the waiting time distributions by applying thresholds on the PVI time-
series. We have confirmed that intermittent magnetic field structures are not
evenly distributed in the solar wind but rather tend to strongly cluster,
forming regions characterized by a high magnetic field variability followed by
intervals for which the magnetic field is relatively smooth. This observation
is also reinforced by the power-law nature of the waiting time distributions
at low waiting times, indicating the presence of an intracluster population.
The power law is followed by an exponential at longer waiting times,
suggesting a second intercluster population of coherent structures in the
solar wind’s magnetic field. Additionally, the power-law scaling of the WT
distributions is indicative of clusters that do not have a typical size or
distance, except for the limiting size given by the exponential cutoff. We
moved on to investigate the correlation of identified coherent structures to
the proton temperature as a proxy to intermittent dissipation in the solar
wind. In Sec. 5.2, we have provided strong evidence that supports the theory
of plasma heating by coherent magnetic structures generated by the turbulent
cascade. More specifically, in agreement with previous studies (Servidio et
al., 2012; Qudsi et al., 2020; Yordanova et al., 2021), the present analysis
suggests that the strongest discontinuities, $PVI\geq 6$ in the magnetic
field, are linked to extreme dissipation events. In the past, both theoretical
(e.g., Servidio et al. (2011)) and observational studies (e.g., (Hou et al.,
2021)), have demonstrated that such events are most likely to be associated
with magnetic reconnection sites. A visual inspection of tens of high PVI
valued events recovered during $E_{1}-E_{6}$ supports this theory, as in most
cases, such events are related to reversals of the magnetic field and zero-
crossings of at least one of the magnetic field components, as well as
significant energy dissipation and proton heating. We have demonstrated in our
analysis that such high PVI value events can considerably alter the local
plasma dynamics by increasing $T_{p}$ in a very short amount of time. However,
as indicated by Fig. 8, and, Fig. 9, due to the relatively small fraction of
these structures in our dataset, their contribution in energy dissipation is
practically insignificant. Consequently, it is rather the accumulative effect
of many small-scale coherent structures of PVI index in the $1\leq PVI\leq 6$
range that determines $T_{p}$ on a macroscopic scale. This means $T_{p}$ is
strongly related to the number density of coherent structures in our dataset,
with high-density intervals being associated with higher $T_{p}$. In
accordance with our study, Hou et al. (2021), have recently analyzed
measurements from the Magnetospheric Multiscale Spacecraft (MMS) mission to
show that despite having a significant role in local dissipation, magnetic
reconnection events are not the major contributor to the magnetosheath’s
internal energy. This study provides an additional confirmation that the
number density of intense coherent structures such as reconnection sites and
current sheets is an important parameter that should be taken into account in
numerical investigations that study the heating of astrophysical and space
plasmas.

Finally, it would be important to compare our results to those of previous
studies investigating intermittent heating at different heliospheric
distances. This comparison can allow us to infer the relative importance of
intermittent heating on different phases of the development of the turbulent
cascade. Such an example is Fig. 3 of Osman et al. (2012), who studied the
effects of intermittent dissipation on $T_{p}$ at 1AU. The caveat here is,
that Osman et al. (2012) have considered values of $PVI\geq 2.4$ to interpret
$T_{p}$ as function of $f_{PVI\geq 2.4}$. However, even though not clearly
outlined in Osman et al. (2012), events of $1\leq PVI\leq 2.4$ seem to also
contribute to the heating of the solar wind protons (e.g., see Fig. 4).
Nevertheless out of this comparison, it is clear that the effects of
intermittent heating on $T_{p}$ close to the solar wind sources are not as
pronounced as in the outer parts of the heliosphere. A possible explanation
for the decreased efficiency of intermittent heating in the young solar wind
is that PSP mostly samples intervals for which the angle between the
background magnetic field and the solar wind flow $\Theta_{VB}$ is parallel.

In the past, intermittency in the inertial range of MHD turbulence has been
shown to be highly anisotropic. For parallel intervals, the statistical
signature of the magnetic field fluctuations is that of a non-Gaussian
globally scale-invariant process, in contrast to multi-exponent statistics
observed when the local magnetic field is perpendicular to the flow direction
(Horbury et al., 2008; Osman et al., 2014). This result can be interpreted as
a decreased level of intermittency for parallel intervals and consequently a
decreased efficiency of the intermittent heating mechanism. We will soon
revisit the topic by studying the evolution of intermittent heating as a
function of helispheric distance, as well as $\Theta_{VB}$ (Sioulas et al. in
preparation).

As noted in Sec. 1, the advantage of the PVI method to capture all kinds of
discontinuities can also be its biggest drawback. Namely, two equally valued
PVI events may correspond to different types of magnetic field
discontinuities, that have contrasting effects on the proton temperature. For
example, different types of coherent structures may instigate heating
mechanisms that preferentially operate close but not inside magnetic field
discontinuities. This could strongly affect the results of this study, as
heating may not coincide in time with strong discontinuities in the magnetic
field. Thus to further categorize such discontinuities in the near sun solar
wind environment, several additional magnetic field and plasma parameters
should be taken into account (e.g., Sorriso-Valvo et al. (2018)). We hope that
future studies will revisit the topic by performing a more thorough analysis,
considering such parameters.

We would like to thank the anonymous referee for their insightful suggestions
and comments, which helped improve and clarify this manuscript. This research
was funded in part by the FIELDS experiment on the Parker Solar Probe
spacecraft, designed and developed under NASA contract NNN06AA01C; the NASA
Parker Solar Probe Observatory Scientist grant NNX15AF34G and the HERMES DRIVE
NASA Science Center grant No. 80NSSC20K0604.

## References

  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., 204, 49, doi: 10.1007/s11214-016-0244-5
  * Bale et al. (2021) Bale, S. D., Horbury, T. S., Velli, M., et al. 2021, ApJ, 923, 174, doi: 10.3847/1538-4357/ac2d8c
  * Bandyopadhyay et al. (2018) Bandyopadhyay, R., Chasapis, A., Chhiber, R., et al. 2018, ApJ, 866, 106, doi: 10.3847/1538-4357/aade04
  * Bandyopadhyay et al. (2020a) Bandyopadhyay, R., Matthaeus, W. H., Parashar, T. N., et al. 2020a, The Astrophysical Journal Supplement Series, 246, 61, doi: 10.3847/1538-4365/ab6220
  * Bandyopadhyay et al. (2020b) Bandyopadhyay, R., Goldstein, M. L., Maruca, B. A., et al. 2020b, ASTROPHYSICAL JOURNAL SUPPLEMENT SERIES, doi: 10.3847/1538-4365/ab5dae
  * Bandyopadhyay et al. (2021) Bandyopadhyay, R., Chasapis, A., Gershman, D. J., et al. 2021, MNRAS, 500, L6, doi: 10.1093/mnrasl/slaa171
  * Biskamp (2003) Biskamp, D. 2003, Magnetohydrodynamic Turbulence
  * Biskamp & Müller (2000) Biskamp, D., & Müller, W.-C. 2000, Physics of Plasmas, 7, 4889, doi: 10.1063/1.1322562
  * Borrini et al. (1981) Borrini, G., Gosling, J. T., Bame, S. J., Feldman, W. C., & Wilcox, J. M. 1981, J. Geophys. Res., 86, 4565, doi: 10.1029/JA086iA06p04565
  * Bruno (2019) Bruno, R. 2019, Earth and Space Science, 6, 656, doi: https://doi.org/10.1029/2018EA000535
  * Bruno et al. (1999) Bruno, R., Bavassano, B., Bianchini, L., et al. 1999, in ESA Special Publication, Vol. 9, Magnetic Fields and Solar Processes, ed. A. Wilson & et al., 1147
  * Burlaga & Ogilvie (1970) Burlaga, L. F., & Ogilvie, K. W. 1970, ApJ, 159, 659, doi: 10.1086/150340
  * Castaing et al. (1990) Castaing, B., Gagne, Y., & Hopfinger, E. J. 1990, Physica D Nonlinear Phenomena, 46, 177, doi: 10.1016/0167-2789(90)90035-N
  * Chasapis et al. (2018) Chasapis, A., Yang, Y., Matthaeus, W. H., et al. 2018, ApJ, 862, 32, doi: 10.3847/1538-4357/aac775
  * Chen et al. (2021) Chen, C., Chandran, B., Woodham, L., et al. 2021, Astronomy & Astrophysics, 650, doi: 10.1051/0004-6361/202039872
  * Chhiber et al. (2020) Chhiber, R., Goldstein, M. L., Maruca, B. A., et al. 2020, The Astrophysical Journal Supplement Series, 246, 31, doi: 10.3847/1538-4365/ab53d2
  * Coleman (1968) Coleman, Paul J., J. 1968, ApJ, 153, 371, doi: 10.1086/149674
  * Cranmer & van Ballegooijen (2005) Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265, doi: 10.1086/426507
  * Cranmer et al. (2007) Cranmer, S. R., van Ballegooijen, A. A., & Edgar, R. J. 2007, ApJS, 171, 520, doi: 10.1086/518001
  * Davies & Gather (1993) Davies, L., & Gather, U. 1993, Journal of the American Statistical Association, 88, 782, doi: 10.1080/01621459.1993.10476339
  * Dudok de Wit et al. (2020) Dudok de Wit, T., Krasnoselskikh, V. V., Bale, S. D., et al. 2020, ApJS, 246, 39, doi: 10.3847/1538-4365/ab5853
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7, doi: 10.1007/s11214-015-0211-6
  * Gazis et al. (1994) Gazis, P. R., Barnes, A., Mihalov, J. D., & Lazarus, A. J. 1994, J. Geophys. Res., 99, 6561, doi: 10.1029/93JA03144
  * González et al. (2021) González, C. A., Tenerani, A., Matteini, L., Hellinger, P., & Velli, M. 2021, The Astrophysical Journal Letters, 914, L36, doi: 10.3847/2041-8213/ac097b
  * Greco et al. (2018) Greco, A., Matthaeus, W. H., Perri, S., et al. 2018, Space Sci. Rev., 214, 1, doi: 10.1007/s11214-017-0435-8
  * Greco et al. (2009) Greco, A., Matthaeus, W. H., Servidio, S., & Dmitruk, P. 2009, Phys. Rev. E, 80, 046401, doi: 10.1103/PhysRevE.80.046401
  * Greco et al. (2010) Greco, A., Perri, S., & Zimbardo, G. 2010, Journal of Geophysical Research (Space Physics), 115, A02203, doi: 10.1029/2009JA014690
  * Gurland & Tripathi (1971) Gurland, J., & Tripathi, R. C. 1971, The American Statistician, 25, 30.  http://www.jstor.org/stable/2682923
  * Hada et al. (2003) Hada, T., Koga, D., & Yamamoto, E. 2003, Space Sci. Rev., 107, 463, doi: 10.1023/A:1025506124402
  * He et al. (2015) He, J., Tu, C., Marsch, E., et al. 2015, ApJ, 813, L30, doi: 10.1088/2041-8205/813/2/L30
  * Hellinger et al. (2011) Hellinger, P., Matteini, L., Štverák, Š., Trávníček, P. M., & Marsch, E. 2011, Journal of Geophysical Research (Space Physics), 116, A09105, doi: 10.1029/2011JA016674
  * Hernández et al. (2021) Hernández, C. S., Sorriso-Valvo, L., Bandyopadhyay, R., et al. 2021, ApJ, 922, L11, doi: 10.3847/2041-8213/ac36d1
  * Horbury et al. (2008) Horbury, T. S., Forman, M., & Oughton, S. 2008, Phys. Rev. Lett., 101, 175005, doi: 10.1103/PhysRevLett.101.175005
  * Hou et al. (2021) Hou, C., He, J., Zhu, X., & Wang, Y. 2021, ApJ, 908, 237, doi: 10.3847/1538-4357/abd6f3
  * Howes et al. (2012) Howes, G. G., Bale, S. D., Klein, K. G., et al. 2012, ApJ, 753, L19, doi: 10.1088/2041-8205/753/1/L19
  * Huang et al. (2020) Huang, J., Kasper, J. C., Vech, D., et al. 2020, The Astrophysical Journal Supplement Series, 246, 70, doi: 10.3847/1538-4365/ab74e0
  * Huang et al. (2020) Huang, J., Kasper, J. C., Vech, D., et al. 2020, ApJS, 246, 70, doi: 10.3847/1538-4365/ab74e0
  * Isliker et al. (2019) Isliker, H., Archontis, V., & Vlahos, L. 2019, The Astrophysical Journal, 882, 57, doi: 10.3847/1538-4357/ab30c6
  * Jain et al. (2021) Jain, N., Büchner, J., Comişel, H., & Motschmann, U. 2021, ApJ, 919, 103, doi: 10.3847/1538-4357/ac106c
  * Karimabadi et al. (2013) Karimabadi, H., Roytershteyn, V., Wan, M., et al. 2013, Physics of Plasmas, 20, 012303, doi: 10.1063/1.4773205
  * Kasper et al. (2016) Kasper, J. C., Abiad, R., Austin, G., et al. 2016, Space Sci. Rev., 204, 131, doi: 10.1007/s11214-015-0206-3
  * Krishna Jagarlamudi et al. (2019) Krishna Jagarlamudi, V., Dudok de Wit, T., Krasnoselskikh, V., & Maksimovic, M. 2019, ApJ, 871, 68, doi: 10.3847/1538-4357/aaef2e
  * Kuzzay et al. (2019) Kuzzay, D., Alexandrova, O., & Matteini, L. 2019, Phys. Rev. E, 99, 053202, doi: 10.1103/PhysRevE.99.053202
  * Leamon et al. (2000) Leamon, R. J., Matthaeus, W. H., Smith, C. W., et al. 2000, ApJ, 537, 1054, doi: 10.1086/309059
  * Lehe et al. (2009) Lehe, R., Parrish, I. J., & Quataert, E. 2009, The Astrophysical Journal, 707, 404, doi: 10.1088/0004-637x/707/1/404
  * Lionello et al. (2014) Lionello, R., Velli, M., Downs, C., Linker, J. A., & Mikić, Z. 2014, ApJ, 796, 111, doi: 10.1088/0004-637X/796/2/111
  * Liou & Wu (2021) Liou, K., & Wu, C.-C. 2021, ApJ, 920, 39, doi: 10.3847/1538-4357/ac1586
  * Lopez & Freeman (1986) Lopez, R. E., & Freeman, J. W. 1986, J. Geophys. Res., 91, 1701, doi: 10.1029/JA091iA02p01701
  * Marsch & Tu (1990) Marsch, E., & Tu, C.-Y. 1990, Journal of Geophysical Research: Space Physics, 95, 8211, doi: https://doi.org/10.1029/JA095iA06p08211
  * Matthaeus & Goldstein (1982) Matthaeus, W. H., & Goldstein, M. L. 1982, J. Geophys. Res., 87, 10347, doi: 10.1029/JA087iA12p10347
  * Matthaeus & Montgomery (1980) Matthaeus, W. H., & Montgomery, D. 1980, Annals of the New York Academy of Sciences, 357, 203, doi: 10.1111/j.1749-6632.1980.tb29687.x
  * Matthaeus & Velli (2011) Matthaeus, W. H., & Velli, M. 2011, Space Science Reviews, 160, 145, doi: 10.1007/s11214-011-9793-9
  * Matthaeus et al. (1999) Matthaeus, W. H., Zank, G. P., Oughton, S., Mullan, D. J., & Dmitruk, P. 1999, ApJ, 523, L93, doi: 10.1086/312259
  * Meneguzzi et al. (1981) Meneguzzi, M., Frisch, U., & Pouquet, A. 1981, Phys. Rev. Lett., 47, 1060, doi: 10.1103/PhysRevLett.47.1060
  * Mozer et al. (2021) Mozer, F. S., Vasko, I. Y., & Verniero, J. L. 2021, ApJ, 919, L2, doi: 10.3847/2041-8213/ac2259
  * Osman et al. (2014) Osman, K., Kiyani, K. H., Chapman, S. C., & Hnat, B. 2014, in AGU Fall Meeting Abstracts, Vol. 2014, SH52A–03
  * Osman et al. (2012) Osman, K. T., Matthaeus, W. H., Wan, M., & Rappazzo, A. F. 2012, Phys. Rev. Lett., 108, 261102, doi: 10.1103/PhysRevLett.108.261102
  * Parashar et al. (2009) Parashar, T. N., Shay, M. A., Cassak, P. A., & Matthaeus, W. H. 2009, Physics of Plasmas, 16, 032310, doi: 10.1063/1.3094062
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664, doi: 10.1086/146579
  * Perez et al. (2021) Perez, J. C., Chandran, B. D. G., Klein, K. G., & Martinović, M. M. 2021, Journal of Plasma Physics, 87, 905870218, doi: 10.1017/S0022377821000167
  * Politano & Pouquet (1998a) Politano, H., & Pouquet, A. 1998a, Geophysical Research Letters, 25, 273
  * Politano & Pouquet (1998b) —. 1998b, Phys. Rev. E, 57, R21, doi: 10.1103/PhysRevE.57.R21
  * Qudsi et al. (2020) Qudsi, R. A., Maruca, B. A., Matthaeus, W. H., et al. 2020, ApJS, 246, 46, doi: 10.3847/1538-4365/ab5c19
  * Richardson et al. (1995) Richardson, J. D., Paularena, K. I., Lazarus, A. J., & Belcher, J. W. 1995, Geophys. Res. Lett., 22, 325, doi: 10.1029/94GL03273
  * Richardson et al. (1995) Richardson, J. D., Paularena, K. I., Lazarus, A. J., & Belcher, J. W. 1995, Geophysical Research Letters, 22, 325, doi: https://doi.org/10.1029/94GL03273
  * Servidio et al. (2012) Servidio, S., Valentini, F., Califano, F., & Veltri, P. 2012, Phys. Rev. Lett., 108, 045001, doi: 10.1103/PhysRevLett.108.045001
  * Servidio et al. (2011) Servidio, S., Dmitruk, P., Greco, A., et al. 2011, Nonlinear Processes in Geophysics, 18, 675, doi: 10.5194/npg-18-675-2011
  * Shi et al. (2022) Shi, C., Velli, M., Tenerani, A., Réville, V., & Rappazzo, F. 2022, arXiv preprint arXiv:2201.02894
  * Shi et al. (2021) Shi, C., Velli, M., Panasenco, O., et al. 2021, A&A, 650, A21, doi: 10.1051/0004-6361/202039818
  * Sioulas et al. (2020a) Sioulas, N., Isliker, H., & Vlahos, L. 2020a, ApJ, 895, L14, doi: 10.3847/2041-8213/ab9092
  * Sioulas et al. (2022) —. 2022, A&A, 657, A8, doi: 10.1051/0004-6361/202141361
  * Sioulas et al. (2020b) Sioulas, N., Isliker, H., Vlahos, L., Koumtzis, A., & Pisokas, T. 2020b, MNRAS, 491, 3860, doi: 10.1093/mnras/stz3259
  * Sorriso-Valvo et al. (2018) Sorriso-Valvo, L., Carbone, F., Perri, S., et al. 2018, Solar Physics, 293, doi: 10.1007/s11207-017-1229-6
  * Sorriso-Valvo et al. (2019) Sorriso-Valvo, L., De Vita, G., Fraternale, F., et al. 2019, Frontiers in Physics, 7, doi: 10.3389/fphy.2019.00108
  * Stansby et al. (2019) Stansby, D., Horbury, T. S., & Matteini, L. 2019, MNRAS, 482, 1706, doi: 10.1093/mnras/sty2814
  * Stansby et al. (2018) Stansby, D., Salem, C., Matteini, L., & Horbury, T. 2018, Sol. Phys., 293, 155, doi: 10.1007/s11207-018-1377-3
  * Suess et al. (2009) Suess, S. T., Ko, Y. K., von Steiger, R., & Moore, R. L. 2009, Journal of Geophysical Research (Space Physics), 114, A04103, doi: 10.1029/2008JA013704
  * Taylor (1938) Taylor, G. I. 1938, Proceedings of the Royal Society of London Series A, 164, 476, doi: 10.1098/rspa.1938.0032
  * Tessein et al. (2013) Tessein, J. A., Matthaeus, W. H., Wan, M., et al. 2013, The Astrophysical Journal, 776, L8, doi: 10.1088/2041-8205/776/1/l8
  * Vech et al. (2017) Vech, D., Klein, K. G., & Kasper, J. C. 2017, ApJ, 850, L11, doi: 10.3847/2041-8213/aa9887
  * Velli et al. (1989) Velli, M., Grappin, R., & Mangeney, A. 1989, Phys. Rev. Lett., 63, 1807, doi: 10.1103/PhysRevLett.63.1807
  * Veltri (1999) Veltri, P. 1999, Plasma Physics and Controlled Fusion, 41, A787, doi: 10.1088/0741-3335/41/3a/071
  * Yordanova et al. (2021) Yordanova, E., Vörös, Z., Sorriso-Valvo, L., Dimmock, A. P., & Kilpua, E. 2021, ApJ, 921, 65, doi: 10.3847/1538-4357/ac1942


# Deductive Verification via the Debug Adapter Protocol

Gidon Ernst Johannes Blau Ludwig-Maximilians-Universität
Munich, Germany<EMAIL_ADDRESS>University of Melbourne
Melbourne, Australia Toby Murray University of Melbourne
Melbourne, Australia<EMAIL_ADDRESS>
###### Abstract

We propose a conceptual integration of deductive program verification into
existing user interfaces for software debugging. This integration is well-
represented in the “Debug Adapter Protocol”, a widely-used and generic
technology to integrate debugging of programs into development environments.
Commands like step-forward and step-in are backed by steps of a symbolic
structural operational semantics, and the different paths through a program
are readily represented by multiple running threads of the debug target inside
the user interface. Thus, existing IDEs can be leveraged for deductive
verification debugging with relatively little effort. We have implemented this
scheme for SecC, an auto-active program verifier for C, and discuss its
integration into Visual Studio Code.

##  1 Introduction

Deductive verification of programs with respect to strong requirements relies
on human proof engineering effort. The user has to provide the primary
correctness specifications (e.g. procedure contracts), as well as auxiliary
annotations (e.g. loop invariants), key lemmas, and other proof hints. This is
much facilitated by modern Integrated Development Environments (IDEs) for
formal methods tools and by advances in verification technology. Over the
recent years, the term “push-button” has been coined, suggesting perhaps that
proof automation is nowadays good enough to not burden the user with internal
details. However, proofs for code that is already correctly annotated are
fundamentally different from the typical trial an error to find these. The
ability to dig into the causes of a verification failure is not just nice-to-
have—it is crucial have access to as much information as possible.

Related Work: To that end, state-of-the-art IDEs for formal development offer
different features: Dafny [14], for example, highlights those annotations
which cannot be proven, and the Boogie Verification Debugger [13] gives
structured access to concrete counterexamples. In VeriFast [12], one can
inspect the symbolic state and a tree representation of the paths explored.
Why3 [5] shows the generated verification conditions in a nice and structured
way and offers interactive as well as automatic proof steps. In contrast,
general purpose interactive theorem provers like Isabelle/PIDE [20], KIV [7],
Rodin [19], and PVS [15] (to name a few with a sophisticated user interface),
tend to expose proof internals in detail. The latter paper [15] nicely
compares some popular formal IDEs and their features.

The approaches mentioned are rather different from the experience of
traditional, concrete debugging of programs in IDEs like Eclipse, IntelliJ, or
Visual Studio Code, where the main features are breakpoints, single stepping,
and inspection of data at runtime. Recently, loose coupling between IDEs and
language-specific toolchains has become popular, based on the Language Server
Protocol (LSP)111https://microsoft.github.io/language-server-protocol and the
Debug Adapter Protocol (DAP).222https://microsoft.github.io/debug-adapter-
protocol LSP is used in several of the above mentioned formal IDEs. For VDM a
debugger for concrete model executions has been developed [18] using the DAP.
The KeY Symbolic Execution Debugger [10] is a feature-rich tool for Java
verification built on the Eclipse platform.

Contribution: In this paper, we propose and describe an integration of
deductive program verification into general purpose IDEs on top of the DAP,
which can be used to interactively navigate to specific parts of a proof. The
first contribution is an abstract characterization of how that integration
works conceptually (Sec. 2). As the second contribution, we briefly describe
an ongoing implementation effort of this scheme for SecC, an autoactive
verifier for correctness and security of C programs (Sec. 3) inside Visual
Studio Code. Our conclusion is that implementing debugging for an existing
verification tool that is based on symbolic execution is relatively easy and
straight-forward when using the DAP.

##  2 Conceptual Model of a Debug Server

Figure 1: Schematics of the Debug View inside an IDE.

In this section we develop a conceptual model of a debug server. The server
exposes certain operations to the IDE and maintains the state of running
processes that are being debugged (the debug targets). Fig. 1 depicts the
integration between the IDE (as the client, on the left) and a debug server
(on the right) in terms of the graphical front-end of the IDE and the messages
exchanged between the two components. Typically, the debugging perspective of
an IDE shows the program’s source code and the breakpoints within. Moreover,
there is a view that shows a current state of the program, in terms of runtime
values of variables (top-left), possibly organized according to the structure
of lexical scopes.

In Fig. 1 we anticipate a debug target that is executed symbolically, such
that the program variables x and y are assigned symbolic expression over
logical variables. The example program, a function abs to compute the absolute
value of parameter x, has been halted at a breakpoint in line 6,
consequentially, variable y is set to $-x_{\\_}0$ wrt. logical variable
$x_{\\_}0$ capturing the initial value of x. The synthetic variable path shows
the current path condition, a formula that captures the constraints for
reaching the respective code location, here the negated test of the
conditional. Below, the two possible branches through function abs are
represented as execution thread with identifiers 00 and 01, that are subsumed
by a parent thread 0. The views in the window on the left are populated from
data that is requested by the IDE from the server.

The simple, semi-formal model below omits inessential details of realistic
languages (e.g., memory, see [6, 3, 12, 9] here for details). The key point is
that the ideas presented here translate to any formal system that can be
described by a symbolic structural operational semantics [16].

### Representation of States and Execution Steps.

The state associated with a debug target primarily consists of a
_configuration_ $c$, that describes the progress of the execution.
Configurations bundle up programs $p$, defined over the set of program
variables $X$, with _symbolic assignment_ $s\colon X\to E$ from program
variables to logical expressions $E$ and a _path condition_ $\varphi$ as a
logical formula that describes the branches taken so far. Here, we capture
execution using three kinds of configurations:

$\displaystyle\text{Configurations }c$
$\displaystyle\Coloneqq\big{(}s,\varphi,p_{\\_}1;\ldots;p_{\\_}k\big{)}\mid\big{(}c_{\\_}1\parallel\cdots\parallel
c_{\\_}n\big{)}\mid\big{(}\varphi\land\lnot\psi,\lightning\big{)}$

_Sequential compositions_ execute programs $p_{\\_}1;\ldots;p_{\\_}k$ from
symbolic assignment $s$ and path condition $\varphi$, where for $k=0$ this
execution branch terminates. _Parallel compositions_ of configurations capture
multiple branches that arise e.g. from conditionals. _Proof obligations_ ,
marked by $\lightning$, verify that condition $\psi$ holds for the path
condition $\varphi$, where $\varphi\land\lnot\psi\ \mathbf{unsat}$ denotes
that the obligation discharges successfully. In practice, such proof
obligations would be annotated with additional contextual information.

Symbolic execution proceeds by unwinding a small-step relation
$c\stackrel{{\scriptstyle\sigma}}{{\longrightarrow}}c^{\prime}$ between
configurations according to a _schedule_ $\sigma$, which reflects the user’s
choice where to step next. Recursively, the first entry in the schedule, for
$1\leq i\leq n$, resolves which part of a parallel configuration performs the
next step:

$\displaystyle\big{(}c_{\\_}1\parallel\cdots c_{\\_}i\cdots\parallel
c_{\\_}n\big{)}$ $\displaystyle\stackrel{{\scriptstyle\langle
i\rangle\cdot\sigma}}{{\longrightarrow}}\big{(}c_{\\_}1\parallel\cdots
c^{\prime}_{\\_}i\cdots\parallel c_{\\_}n\big{)}\qquad\text{ if }\qquad
c_{\\_}i\stackrel{{\scriptstyle\sigma}}{{\longrightarrow}}c^{\prime}_{\\_}i$

In the base case, execution of sequential composition steps with the empty
schedule $\sigma=\langle\rangle$. Nondeterminism is captured by producing
multiple branches, for example:

$\displaystyle\big{(}s,\varphi,\mathbf{assume}\ e;\ldots\big{)}$
$\displaystyle\stackrel{{\scriptstyle\langle\rangle}}{{\longrightarrow}}\big{(}s,\varphi\land
s(e),\ldots\big{)}$ $\displaystyle\big{(}s,\varphi,\mathbf{assert}\
e;\ldots\big{)}$
$\displaystyle\stackrel{{\scriptstyle\langle\rangle}}{{\longrightarrow}}\big{(}s,\varphi\land
s(e),\ldots\big{)}\parallel\big{(}\varphi\land\lnot s(e),\lightning\big{)}$
$\displaystyle\big{(}s,\varphi,x\coloneqq e;\ldots\big{)}$
$\displaystyle\stackrel{{\scriptstyle\langle\rangle}}{{\longrightarrow}}\big{(}s[x\mapsto
s(e)],\varphi,\ldots\big{)}$ $\displaystyle\big{(}s,\varphi,\mathbf{if}\ e\
\mathbf{then}\ p_{\\_}1\ \mathbf{else}\ p_{\\_}2;\ldots\big{)}$
$\displaystyle\stackrel{{\scriptstyle\langle\rangle}}{{\longrightarrow}}\big{(}s,\varphi\land
s(e),p_{\\_}1;\ldots\big{)}\parallel\big{(}s,\varphi\land\lnot
s(e),p_{\\_}2;\ldots\big{)}$

where $s(e)$ denotes evaluation of program expression $e$ in symbolic state
$s$. Of course, configurations with an unsatisfiable path condition or empty
sequential compositions can be soundly dropped from their surrounding parallel
context (which our implementation does eagerly), e.g., the proof obligation
from an $\mathbf{assert}$ when $s(e)$ follows from $\varphi$, or either of the
branches of an $\mathbf{if}$ in case it is unreachable. Loops can be unwound
interactively (not discussed here) or summarized by invariants $I$ as shown

$\displaystyle\big{(}s,\varphi,\mathbf{while}\ e\ \mathbf{do}\
p;\ldots\big{)}$
$\displaystyle\stackrel{{\scriptstyle\langle\rangle}}{{\longrightarrow}}\big{(}s,\varphi\land\lnot
s(I);\lightning\big{)}\parallel\big{(}\hat{s},\varphi\land\hat{s}(e)\land\hat{s}(I),p;\mathbf{assert}\
I\big{)}\parallel\big{(}\hat{s},\varphi\land\lnot\hat{s}(e)\land\hat{s}(I);\ldots\big{)}$

where the latter produces three successor configurations, 1) to prove
invariant $I$ initially, 2) to preserve $I$ over an arbitrary iteration, where
$\hat{s}=s[x_{\\_}1\mapsto x^{\prime}_{\\_}1,\ldots,x_{\\_}n\mapsto
x^{\prime}_{\\_}n]$ introduces fresh logical variables
$x_{\\_}1^{\prime},\ldots,x^{\prime}_{\\_}n$ for the program variables
modified in loop body $p$, and 3) to continue with the code after the loop.

### Initialization.

The initial symbolic configurations for an entire translation unit is
exemplified below. A C file has a list of global variables $g_{\\_}i$ that are
to be initialized in sequence, which we represent as a program $g$ such that
$g\equiv\\{g_{\\_}1\coloneqq e_{\\_}1;\ldots;g_{\\_}n\coloneqq e_{\\_}n\\}$. A
top-level procedure declaration of the form
$f(y_{\\_}1,\ldots,y_{\\_}m)~{}\\{q\\}$, with precondition/postcondition pair
$P,Q$, parameters $y_{\\_}j$, and implementation $q$ can be mapped onto a
configuration $c_{\\_}f=\big{(}s,\mathit{true},\mathbf{assume}\
P;q;\mathbf{assert}\ Q\big{)}$ where $s_{\\_}f=[g_{\\_}i\mapsto
x_{\\_}i,y_{\\_}j\mapsto y_{\\_}j,\ldots]$ initializes all these variables in
the respective scopes to fresh logical ones. The verification of the main
procedure may additionally assume that the globals were just initialized,
which can be represented as
$c_{\\_}\texttt{main}=\big{(}s_{\\_}\texttt{main},\mathit{true},g;q\big{)}$
where the body $q$ of main is prefixed by the sequence $g$ of global
initializers.

### Dispatching Requests in the Debug Server.

We outline how to realize the operations that implement the main requests
issued by the IDE in reference to a top-level configuration $C$. In addition,
we track a set $B$ of breakpoints, which are program locations, and we denote
the location of a program by $\mathit{loc}(p)$. In the following, we denote by
$\sigma(C)=\big{(}s,\varphi,\ldots\big{)}$ the sub-configuration triggered by
$\sigma$, which is necessarily a sequential one in our simple model.

  * •

The request $\texttt{Launch}(\mathit{file})$ for program stored in
$\mathit{file}$ produces the initial configuration as a big parallel
composition $C\coloneqq\big{(}c_{\\_}{f_{\\_}1}\parallel\cdots\parallel
c_{\\_}{f_{\\_}n}\parallel c_{\\_}\texttt{main}\big{)}$ where the respective
sub-configurations are constructed wrt. the procedure declarations and globals
as outlined above.

  * •

The request GetThreads returns the currently running “threads”, which are
those parts of parallel compositions that can be stepped. This information can
be represented in terms of possible schedules for the next step, i.e., the set
GetThreads returns $\\{\sigma\mid\exists\ c^{\prime}.\
C\stackrel{{\scriptstyle\sigma}}{{\longrightarrow}}c^{\prime}\\}$. As
suggested in Fig. 1, the hierarchical structure of these identifiers can be
exploited for more complex operations, such as stepping all branches that
share a common prefix schedule.

  * •

The request $\texttt{GetVariables}(\sigma)$ returns the variable assignment
$s$ stored in the sequential configuration $\sigma(C)$ as remarked above. In
our implementation for SecC, we use this request to add synthetic variables,
such as $\texttt{path}\mapsto\varphi$ where $\varphi$ is the path constraint
of configuration $\sigma(C)$ for that thread, and we additionally include a
representation of the symbolic heap and information about the attacker level
for security proofs.

  * •

The request $\texttt{SetBreakpoints}(b)$ simply sets $B\coloneqq b$ to the set
of specified breakpoints. The protocol knows not only source breakpoints but
also function breakpoints (triggered by calls), and exception breakpoints
(when an exception is thrown), which we have not used so far.

  * •

The requests $\texttt{Next}(\sigma)$ and $\texttt{StepIn}(\sigma)$ execute a
single transition of a given thread $\sigma$ according to the rules for
$C\stackrel{{\scriptstyle\sigma}}{{\longrightarrow}}c^{\prime}$, where
$c^{\prime}$ is taken as the next state, $C\coloneqq c^{\prime}$. The
difference between these two commands in a concrete execution is that the
first proceeds over function calls in one atomic step, whereas the second
jumps into functions. This behavior can be mirrored in a modular deductive
verifier, where Next dispatches such a call using a given function contract,
whereas StepIn inlines the call and disregards such contracts, similarly for
proving loops with invariants or unfolding a finite number of iterations Our
implementation supports Next only so far, but e.g. KIV and KeY support both
interactions in their respective GUIs.

  * •

The request $\texttt{Continue}(\sigma)$ executes multiple transition of a
given thread $\sigma$ until the corresponding configuration
$\sigma(C)=\big{(}s,\varphi,p;\ldots\big{)}$ with $\mathit{loc}(p)\in B$,
i.e., the program execution has reached a breakpoint, or until
$\sigma(C)=\big{(}s,\varphi\big{)}$ is final with no residual program.

  * •

The request $\texttt{StepBack}(\sigma)$ undoes the latest corresponding
transition (or sequence of transitions) of a particular sub-configuration.
This can be realized by keeping a history $c_{\\_}0,c_{\\_}1,\ldots$ of
previous top-level configurations, which in practice is facilitated by the
fact that often tools are implemented in functional languages and do not use
destructive modification of states.

  * •

The request $\texttt{Evaluate}(\sigma,e)$ is issued to inspect the value of
arbitrary expressions $e$ within a state. The response consists of evaluation
$s(e)$ wrt. the logical variables from
$\sigma(C)=\big{(}s,\varphi,\ldots\big{)}$. Since the format of the result is
just a string, further information can be computed with the help of a solver
and included in the response, such as whether $\varphi\implies s(e)$ holds
when $e$ is boolean, or concrete values for $e$ and its free variables if
$\varphi\land s(e)\ \mathbf{sat}$.

##  3 Debug Server Implementation for SecC

SecC333https://covern.org/secc is an autoactive verifier programs for
functional correctness and security of C programs. It is built around Security
Concurrent Separation Logic (SecCSL) [6], which can express value-dependent
security properties of concurrent heap-manipulating programs. The tool is
currently used to verify a variety of small case studies. Internally, SecC is
based on a symbolic execution engine for Separation Logic [3], which is
similar to that of VeriFast (the latter is described nicely in [12]). Thus,
SecC lends itself to the approach outlined in Sec. 2.

SecC is implemented in the Scala programming
language,444https://bitbucket.org/covern/secc/ which runs on the Java Virtual
Machine so that we can rely on the mature library
lsp4j,555https://github.com/eclipse/lsp4j which fully abstracts the DAP (and
also the Language Server Protocol, LSP) in Java. Creating a debug server with
lsp4j simply amounts to implementing a particular Java interface whose
operations correspond to DAP requests, all protocol data structures are
available as Java classes, too. Within the client, here Visual Studio Code,
some additional effort has to be spent to register the respective language
extension, and to provide the necessary hooks that spin up the debug server
with an appropriate configuration.

Figure 2: Debugging an insecure concurrent program in VSCode with SecC.

The integration was developed as a VS Code extension in a Master thesis
project by the second author over the course of roughly six months, resulting
in about 1000 LoC of Scala for the server and 400 LoC of TypeScript for the
extension, albeit getting a working initial version for a toy language was a
matter of a few days. In addition to debugging, the VS Code extension provides
syntax highlighting, verify-on-save, a debug console that lets one inspect the
current state and evaluate expressions, and a graphical view of the symbolic
execution tree. Note that these extended features cannot be realized via the
DAP alone but require LSP functionality as well as the extension facilities
inside VS Code. The extension is currently available in binary form (file
secc-0.2.0.vsix on Bitbucket) and will be released as open source soon.

A screenshot of the debugging perspective is shown in Fig. 2, for an example
from [6, Sec. 2] that has a defect. The symbolic store $s$ appears under
“Vars” and the synthetic variables, including $\varphi$ as “Path” under
“State” on the left; in addition there is a list of symbolic heap chunks
describing the memory. Stepping line 14 with the controls shown at the top-
right subsequently leads to a verification failure. Informally, the code
writes a secret value to a public memory location OUTPUTREG. This can be
recognized from the data shown as follows: Value $d_{\\_}7$ stored in
rec->data (item 1 under Heap) is classified information (item 3:
$d_{\\_}7\mathbin{\dblcolon}\texttt{high}$ under “Path”), whereas the memory
location OUTPUTREG is public (item 2:
$\texttt{OUTPUTREG}_{\\_}2\mapsto[\texttt{low}]\ldots$ under “Heap”).

While we have not done a systematic evaluation or larger case study in this
new SecC IDE yet, the integration was useful for the second challenge of the
VerifyThis 2021 competition [8]. During the competition, it was helpful to
investigate the symbolic states while developing the correctness proof
interactively, for example to determine some subtle arithmetic constraints, or
to debug the unfolding/folding of memory predicates.

##  4 Discussion & Conclusion

We have shown a concept to embed symbolic execution engines into existing IDEs
for interactive debugging via the established Debug Adapter Protocol (Sec. 2)
Our implementation for the autoactive verifier SecC proved to be straight-
forward and low effort (Sec. 3).

The approach inherits as a limitation the exponential path explosion from the
nondeterministic execution when branches are not joined. This is a problem
with many conditionals in sequence, but we have not yet been impeded by this
limitation. At a conceptual level, it is not entirely clear how to remedy the
approach proposed here with ideas that defer splitting up branches to the SMT
solver as it is done in Boogie [2] and generally in Horn clause verifiers [4].
Our approach can nevertheless complement such ideas, e.g., by stepping
selectively only that thread corresponding to a particular procedure or branch
of interest to investigate precisely those proof obligations that fail with
the more efficient encoding of [2].

The vanilla formulation of Sec. 2 repeatedly traverses the tree of
configurations down to the currently executing leaf and rebuilds it on the way
back. Zippers [11] avoid this (probably perceived) inefficiency, and initial
experiments with changing the implementation suggests that Zippers lead to
quite elegant code, too. This idea has indeed been followed before [17].

Overall, we think that the proposed approach is general and flexible enough,
to be used to retrofit existing verification tools and languages with a
symbolic interactive debugger. By relying on existing infrastructure, such an
undertaking is well within the reach of short-term projects. By relying on
established interaction paradigms, the approach brings software development
practice and program verification a step closer together. For future work we
would like to investigate how to integrate concrete symbolic and concrete
debugging techniques, and we plan to conduct a larger case study inside the
SecC IDE to evaluate the benefits of the proposed approach in practice.

_Acknowledgement._ We thank the reviewers for their suggestions to improve the
presentation.

## References

  * [1]
  * [2] Mike Barnett & K Rustan M Leino (2005): _Weakest-precondition of unstructured programs_.  In: Proc. of Program Analysis for Roftware Tools and Engineering (PASTE), ACM, pp. 82–87, 10.1145/1108792.1108813.
  * [3] Josh Berdine, Cristiano Calcagno & Peter W O’Hearn (2005): _Symbolic execution with Separation Logic_.  In: Proc. of Asian Symposium on Programming Languages and Systems (APLAS), LNCS 3780, Springer, pp. 52–68, 10.1007/115754675.
  * [4] Nikolaj Bjørner, Arie Gurfinkel, Ken McMillan & Andrey Rybalchenko (2015): _Horn clause solvers for program verification_.  In: Fields of Logic and Computation II, LNCS 9300, Springer, pp. 24–51, 10.1007/978-3-319-23534-92.
  * [5] François Bobot, Jean-Christophe Filliâtre, Claude Marché, Guillaume Melquiond & Andrei Paskevich (2013): _The Why3 platform_.  Technical Report, LRI, CNRS & Univ. Paris-Sud & INRIA Saclay.  Available at https://hal.inria.fr/hal-00822856.
  * [6] G. Ernst & T. Murray (2019): _SecCSL: Security Concurrent Separation Logic_.  In: Proc. of Computer Aided Verification (CAV), LNCS 11562, Springer, pp. 208–230, 10.1007/978-3-030-25543-513.
  * [7] G. Ernst, J. Pfähler, G. Schellhorn, D. Haneberg & W. Reif (2015): _KIV—Overview and VerifyThis competition_.  Software Tools for Technology Transfer (STTT) 17(6), pp. 677–694, 10.1007/s10009-014-0308-3.
  * [8] Gidon Ernst, Marieke Huisman, Wojciech Mostowski & Mattias Ulbrich (2019): _VerifyThis–Verification competition with a human factor_.  In: Proc. of Tools and Algorithms for the Construction and Analysis of Systems (TACAS), LNCS 11429, Springer, pp. 176–195, 10.1007/978-3-030-17502-312.
  * [9] José Fragoso Santos, Petar Maksimović, Sacha-Élie Ayoun & Philippa Gardner (2020): _Gillian, Part I: A multi-language platform for symbolic execution_.  In: Proc. of Programming Language Design and Implementation (PLDI), ACM, pp. 927–942, 10.1145/3385412.3386014.
  * [10] Martin Hentschel, Reiner Hähnle & Richard Bubel (2016): _The interactive verification debugger: Effective understanding of interactive proof attempts_.  In: Proc. of Automated Software Engineering (ASE), ACM, pp. 846–851, 10.1145/2970276.2970292.
  * [11] Gérard Huet (1997): _The zipper_.  Journal of functional programming 7(5), pp. 549–554, 10.1017/S0956796897002864.
  * [12] Bart Jacobs, Jan Smans, Pieter Philippaerts, Frédéric Vogels, Willem Penninckx & Frank Piessens (2011): _VeriFast: A powerful, sound, predictable, fast verifier for C and Java_.  In: Proc. of NASA Formal Methods (NFM), Springer, pp. 41–55, 10.1007/978-3-642-20398-54.
  * [13] Claire Le Goues, K Rustan M Leino & Michał Moskal (2011): _The Boogie Verification Debugger_.  In: Proc. of Software Engineering and Formal Methods (SEFM), Springer, pp. 407–414, 10.1007/978-3-642-24690-628.
  * [14] K. Rustan M. Leino & Valentin Wüstholz (2014): _The Dafny Integrated Development Environment_.  In: Proc. of Formal Integrated Development Environment (F-IDE), EPTCS 149, pp. 3–15, 10.4204/EPTCS.149.2.
  * [15] Paolo Masci & César A. Muñoz (2019): _An Integrated Development Environment for the Prototype Verification System_.  In: Proc. of Formal Integrated Development Environment (F-IDE), EPTCS 310, pp. 35–49, 10.4204/EPTCS.310.5.
  * [16] Gordon D Plotkin (2004): _The origins of structural operational semantics_.  The Journal of Logic and Algebraic Programming (JLAP) 60, pp. 3–15, 10.1016/j.jlap.2004.03.009.
  * [17] Norman Ramsey & Joao Dias (2006): _An applicative control-flow graph based on Huet’s zipper_.  Electronic Notes in Theoretical Computer Science (ENTCS) 148(2), pp. 105–126, 10.1016/j.entcs.2005.11.042.
  * [18] Jonas Kjær Rask, Frederik Palludan Madsen, Nick Battle, Hugo Daniel Macedo & Peter Gorm Larsen (2021): _Visual Studio Code VDM Support_.  In: Proc. of Overture Workshop, pp. 35–49.  Available at https://arxiv.org/abs/2101.07261.
  * [19] Laurent Voisin & Jean-Raymond Abrial (2014): _The Rodin platform has turned ten_.  In: Proc. of Abstract State Machines, Alloy, B, TLA, VDM, and Z (ABZ), LNCS 8477, Springer, pp. 1–8, 10.1007/978-3-662-43652-31.
  * [20] Makarius Wenzel (2012): _Isabelle/jEdit–A Prover IDE within the PIDE framework_.  In: Proc. of Intelligent Computer Mathematics (AISC/MKM/Calculemus), LNCS 7362, Springer, pp. 468–471, 10.1007/978-3-642-31374-538.