off a switch” 7 Request permissions when needed 15 As a key to a door 3 End-
Users should be notified 11 As a bridge 2 Trust & reputation 9

##  Appendix E Questionnaire for the Survey Study With End-Users

[Answer options to close-ended questions were randomized where appropriate.]

  * •

On average, how many hours do you spend on your primary smartphone on a daily
basis? Please estimate your daily average usage in hours. [Slider 0–24]

  * •

For how many years have you been using a smartphone? (Numbers only)

  * •

Please name five smartphone apps that you have frequently used in the past
year? Use commas to separate the items. For example, app1, app2, app3, app4,
app5.

  * •

In this survey, we use the term “permissions” or “permission dialogs” to refer
to the following example dialogs in smartphone operating systems like Android
and iOS. [Example dialogs from Android and iOS documentation were included]

  * •

In one sentence, why do smartphone apps request permissions? There is no right
or wrong answer, please use your own words, as we are not looking for a
correct or a perfect answer.

  * •

In one sentence, how do you explain what permissions requested by apps do in
smartphones? There is no right or wrong answer, please use your own words, as
we are not looking for a correct or a perfect answer.

  * •

Please select the top five permission requests that you have frequently seen
in the past year on your smartphone? $-$ Photos or Gallery $-$ Location $-$
Camera $-$ Internet $-$ Contacts $-$ Microphone $-$ Sensors or Motion $-$
Calendars $-$ Bluetooth $-$ Biometrics (e.g., Touch ID or Face ID) $-$ Storage
or files $-$ WiFi or Network state

  * •

If a smartphone app uses permissions, it means that the app . . .

$-$ Can see the content of all the files on the device you are using

$-$ Is not a risk to infect your device with a computer virus

$-$ Will automatically prompt you to update your web browser software if it is
out of date

$-$ Can access certain resources on your device

$-$ I don’t know

  * •

Why do you choose to grant permissions requested by smartphone apps? (Select
all that apply)

$-$ I want to use a specific feature that requires permission

$-$ I think the app will not work otherwise

$-$ I trust the app developer

$-$ Because the app is popular

$-$ I will not be able to grant permissions later

$-$ I have nothing to hide

$-$ I want the permission screen to go away

$-$ Nothing bad will happen

$-$ The app developer already has the requested information about me

$-$ I don’t know

  * •

To what extent do you agree or disagree with the following statements? [Scale:
Strongly disagree, Somewhat disagree, Neither agree nor disagree, Somewhat
agree, Strongly agree]

$-$ There are certain permissions that I am sensitive about accepting.

$-$ I do not care about what permissions smartphone apps ask for.

$-$ There have been times when I have decided NOT to install a smartphone app
because of its permission requests.

$-$ Granting or not granting a permission request from a smartphone app is my
choice and responsibility.

$-$ I need clearer explanations describing each permission request I receive.

$-$ Current permissions requested by smartphone apps are broad and generic.

$-$ I trust smartphone apps with fewer permission requests.

$-$ If a smartphone app does not contain permission requests, it is safe to
use, and I feel comfortable using it.

$-$ I feel overwhelmed by the number of permission requests I get on my
smartphone.

$-$ I would grant permission requested by an app developed by someone or a
software company that I trust and has a good reputation.

$-$ There is no need for permission requests on smartphones.

  * •

Thinking about the five smartphone apps that you have frequently used in the
past year, to what extent do you agree or disagree with the following
statements? [Scale: Strongly disagree, Somewhat disagree, Neither agree nor
disagree, Somewhat agree, Strongly agree]

$-$ I am aware of the permissions used in the smartphone apps that I
frequently use. $-$ Permissions requested by the smartphone apps that I
frequently use are aligned with my expectations.

$-$ I understand why certain permissions are requested by smartphone apps that
I frequently use.

  * •

Have you ever removed a permission that you have granted?

$-$ Yes, I have removed a permission that I granted before. Please explain why
and how.

$-$ No, I have not ever removed a permission that I granted before. Please
explain why.

  * •

Do you think permissions requested by a smartphone app you use could harm you?

$-$ Yes, I can see harm in permission requests. Please explain why

$-$ No, I do not see any harm in permission requests. Please explain why.

  * •

In your opinion, to what extent does granting the following permissions on
your device have privacy implications for you? [Scale: None at all, A little,
A moderate amount, A lot, A great deal]

$-$ Photos or Gallery $-$ Location $-$ Camera $-$ Internet $-$ Contacts $-$
Microphone $-$ Sensors or Motion $-$ Calendars $-$ Bluetooth $-$ Biometrics
(e.g., Touch ID or Face ID) $-$ Storage or files $-$ WiFi or Network state

  * •

Please rate the following statements. [Scale: None at all, A little, A
moderate amount, A lot, A great deal]

$-$ How often do you check the permissions requested by an app before
installing it on your smartphone?

$-$ How often do you read the description and explanation of permission
dialogs?

$-$ How often do you contact the developer of a smartphone app to ask for more
information and clarifications about the app’s permissions?

$-$ How often do you change the settings of your smartphone to manage the
permissions of a smartphone app? For example, to remove an already given
permission or to grant an already removed permission.

  * •

To what extent do you agree or disagree with the following statements? [Scale:
Strongly disagree, Somewhat disagree, Neither agree nor disagree, Somewhat
agree, Strongly agree]

$-$ After installing an app on my smartphone, it would bother me if the app
started asking me to grant permissions.

$-$ If I installed an app on my smartphone, I would think twice before
granting any of the app permissions.

$-$ If I accepted app permissions, my activities on my smartphone would be
monitored or tracked at least part of the time.

$-$ If I accepted app permissions, I would be concerned that the app would
know more information about me.

$-$ If I accepted app permissions, I would be concerned that the app would
monitor my activities on my smartphone.

$-$ If I accepted app permissions, others would know about me more than I
would be comfortable with.

$-$ If I accepted app permissions, information about me that I consider to be
private or sensitive would be more easily accessible to others than I would
want.

$-$ If I accepted app permissions, information about me would be out there.
Also, if that information were used, my privacy would be invaded.

$-$ If I accepted app permissions, I would be concerned that the app might use
my personal information for other purposes without notifying me or getting my
authorization.

$-$ If I accepted app permissions, I would be concerned that the app might use
my information for other purposes.

$-$ If I accepted app permissions, I would be concerned that the app might
share my personal information with other entities without getting my
authorization.

##  Appendix F Additional Plots for the Survey Study With End-Users

Figure 4. End-User participants’ answers to “In your opinion, to what extent
does granting the following permissions on your device have privacy
implications for you?”

End-User participants’ answers to “In your opinion, to what extent does
granting the following permissions on your device have privacy implications
for you?” Each bar shows end-users answers to a Likert question.

Figure 5. End-User participants’ answers to questions about privacy
ramifications of permissions.

End-User participants’ answers to questions about privacy ramifications of
permissions. Each bar shows end-users answers to a Likert question.

Figure 6. End-User participants’ answers to questions rooted in the developer
study.

End-User participants’ answers to questions rooted in the developer study.
Each bar shows end-users answers to a Likert question.

Figure 7. End-User participants’ answers to “Please select the top five
permission requests that you have frequently seen in the past year on your
smartphone?”

End-User participants’ answers to “Please select the top five permission
requests that you have frequently seen in the past year on your smartphone?”
Each bar shows end-users answers to a Likert question.

Figure 8. End-User participants’ answers to “Why do you choose to grant
permissions requested by smartphone apps? (Select all that apply).”

End-User participants’ answers to “Why do you choose to grant permissions
requested by smartphone apps? (Select all that apply).”


# Camera PPG waveforms at the forehead

A.C. den Brinker, H. Der Sarkissian, J.H. Wülbern, B. Balmaekers,
M. Padalko, J. Sénégas, R. Springorum and C. Possanzini111A.C. den Brinker, H.
Der Sarkissian, B. Balmaekers and M. Padalko are with Philips Research,
Eindhoven (NL). J.H. Wülbern and J. Sénégas are with Philips Research, Hamburg
(De). R. Springorum and C. Possanzini are with MR R&D, Philips Healthcare,
Best (NL).

## Abstract

In order to obtain insights into the feasibility of replacing ECG-guided
triggering in magnetic resonance imaging (MRI) by a system based on video
photoplethysmography (PPG), PPG and ECG data were collected from volunteers in
an MRI scanner. PPG waveforms obtained using remote camera PPG directed at the
forehead are studied in qualitative and quantitative sense over a number of
volunteers. The data analysis considers variations in PPG waveforms across
volunteers, modelling of the waveforms in Fourier series, dependencies of
waveforms and features on the interbeat interval (IBI) and breath-holding, and
models for ECG-blind estimation of R-peak position. The main findings are that
the PPG waveform depends on the volunteer and that its shape changes with IBI
and does not depend on breath-holding in the given scenario. Low-order
harmonic models provide accurate approximations to the PPG waveform, where for
higher IBI the waveform shows more temporal details. Accurate predictions (20
ms std) of the delays between markers in ECG and PPG appear feasible from a
single PPG feature.

##  1 Introduction

Photoplethysmography (PPG) is a technique that measures the pulsation of blood
volume based on the subtle color changes of the human skin, which led to
various applications in the field of physiological and health monitoring
reviewed in [1]. Originally, it was introduced as a contact-measurement and is
used, e.g., in form of the finger tip cardiac sensor in MRI. More recently, a
camera-based variant of PPG has been introduced to remotely monitor the heart
rate as described in [2, 3].

Widespread use of cardiac MRI in standard clinical care is still hindered by
complicated planning and workflow. The application of ECG patches adds to this
complexity. Their application requires working time of skilled staff, and it
adds to the exam time which lowers patient throughput. Patients need to be
touched which increases risk of infection and many have to be shaved which
lowers patient comfort. Though ECG signal quality is mostly sufficiently
robust for cardiac triggering, sometimes electrodes have to be repositioned,
optimal positioning depends on the skill of the operator, and occasionally ECG
triggering is not possible at all, especially in obese patients.

Technically, ECG signal quality in the MR system is disturbed by several
effects: (i) the magneto-hydrodynamic effect, i.e., Lorentz forces on ions in
the flowing blood exposed to the stationary magnetic field cause an additional
electric field superimposed on the ECG as described in [4]; (ii) switched
magnetic gradient fields for spatial encoding of the MR signal and (iii)
switched control currents to (de-)tune radio frequency coils induce additional
voltages that add to the ECG (see [5]). Hence, the ECG signal measured in the
MR bore becomes distorted, which prevents the ECG to be diagnostic and lowers
the robustness of ECG triggering especially in high and ultrahigh field MR
systems.

Triggering based on markers in a remote PPG signal does not have these
disadvantages and has been proposed as an alternative [6, 7]. Replacement of
ECG triggers by a PPG-based trigger mechanism requires detailed knowledge of
the PPG signal, its characteristics and distribution of properties over the
pertinent population. This paper provides first insights into PPG waveforms
specifically targeting this application.

The outline of the paper is as follows. We start with describing the data
collection. We continue with describing the data analysis methods with
emphasis on personalized prototype construction, modeling and PPG features.
The results of the analysis are shown in Section 4, featuring intraperson
variations and R-peak predictability. We close with a discussion and
conclusions. The paper is supported by further online material.

##  2 Methods: data collection

The data were acquired on volunteers on MR systems located in test bays at the
Philips Innovation Center, BIU MR in Best (NL). These test-bays are equipped
with _3T Ingenia Elition X_ systems including a VitalEye product camera, which
provides camera-based respiratory-triggered MRI and a view into the bore. This
study has been approved by the Internal Committee Biomedical Experiments of
Philips Research, and informed consent has been obtained from each volunteer.

Volunteers were set up as for a conventional cardiac exam using the anterior
and posterior coil, ECG-electrodes, and in some volunteers the finger-tip PPU.
The cardiac exam included several periods where the participants were asked to
hold their breath.

The ECG-signal was manually calibrated after the volunteer was moved into the
bore. The software of the commercialized MR system was used to determine the
R-peaks of the ECG waveform and the peaks of the contact-PPG waveform.

The standard mirror was replaced by a system having main mirror and side
mirror for improved view of the subjects face.

The available data was collected in two sessions. In the first session, 9
volunteers participated. For a second session 20 volunteers were recruited.
One of the volunteers in the second session had an extremely noisy PPG; this
data has been excluded from further analysis. For the second session the age
range was 30-85 years, with weight 60-102 kg and height 162-196 cm. In total
we report on data from 28 participants.

##  3 Methods: signal analysis

###  3.1 Prototype creation

For prototype creation, PPG signal segments created by an automated skin
selection mechanism were used. Segments with 50% overlap in time were taken
and analyzed. The segments were cut into cardiac cycles. The signal in each a
cardiac cycle was resampled to a uniform grid (e.g., $N_{s}=100$ samples in a
cycle), and the median was taken over the different cycles where the number of
cycles is denoted as $N_{c}$. The median is called the prototype PPG waveform;
next to the median, the interquartile range (IQR) was determined to obtain a
measure for the variability. Two methods for dividing the segments into
cardiac cycles were used providing waveforms $p_{i}(k)$, where $i$ gives the
method ($i=1,2$) and $k$ defines the sample within the cardiac cycle
$k=1,\cdots,N_{s}$. The two methods are as follows.

Method 1: Extraction by ECG-based segmentation
In this extraction method, the ECG signal was used for segmentation. The
R-peaks of the ECG signal were identified and the PPG signal was resampled and
segmented according to this timing.

Method 2: Extraction using PPG signal only
In this extraction method, no external timing information was used. As marker
for time zero, the downgoing zero-crossing was used. This instance closely
corresponds to the steepest slope of the zero-mean PPG signal.

The prototypes of both methods are denoted as $p_{1}$ and $p_{2}$,
respectively. Cardiac cycles with a duration much shorter or longer than the
median IBI were discarded. The prototypes $p_{1}$, $p_{2}$ are compared in the
results section. To make such a comparsion, the variation in the measurements
need to be taken into account. The variation of the waveform in the various
cycles was tracked and was quantified by the interquartile range (IQR) at each
sample in the cycle. For clarity and simplicity, only the IQR of the $p_{1}$
prototype is shown in the graphs. In view of the additional processing needed
to determine segmentation in the second method, one would expect the second
method to provide more noisy results.

The cardiac cycles were also separated based on whether or not they belonged
to a breath-hold epoch. In this way, two PPG prototypes were created for each
participant with one reflecting the average PPG waveform during breath-hold
periods and the other in other periods. In order to limit the amount of work
associated with the required manual annotation, this was done for a limited
number of participants.

For further intraperson data analysis, the cardiac cycles were separated into
3 bins depending on the IBI. Per bin a prototype was constructed, yielding
three prototype per volunteer corresponding to a low, medium and high IBI.
Systematic changes in features of these prototypes (signal strength and
position of the maximum and zero-crossings) were studied.

###  3.2 Harmonic model

The discrete-time prototype was modeled in a Fourier series

$p(k)=\sum_{m=0}^{M}A_{m}\cos(2\pi k/Ns+\phi_{m})+e(k)$ (1)

where $M$ is the order of the harmonic model, $m$ denotes the harmonic
component ($m=0$ being the fundamental), $A_{m}$ the amplitude, $\phi_{m}$ the
phase, and $e(k)$ the non modeled part. In the results section the modeled
energy contained in this series is considered as a function of $M$, as well as
its dependence on the IBI.

###  3.3 Feature extraction

From the PPG prototype, a number of feature can be extracted. Although
numerous features can be defined, we consider only four. All of them are
related to specific instances in time. The first three are the M, F and D
feature being the position of the maximum M of the PPG signal, the negative-
going zero-crossing F, and the positive-going zero-crossing D. The zero-
crossing F closely corresponds to steepest decrease in PPG signal assumed to
be the highest incoming flux of the blood volume or its induced tissue effect.
D corresponds to a decay of blood flux (at the forehead) assumed to occur in
the diastolic phase. All three features are markers from the prototype itself
and indicated in Fig. 1. The fourth feature is called ZH and is the position
of the positive-going zero-crossing of the phase of the analytic signal $a$
associated with $p$. The feature ZH typically precedes the maximum M.

Figure 1: PPG prototype markers M, F and D.

##  4 Results

###  4.1 Prototypes: extraction method

In Fig. 2, the constructed PPG prototypes are shown for two volunteers, where
the results from both prototype construction methods are shown. It shows that
the ECG-blind method provides prototype waveforms agreeing well with those
obtained by using ECG timing information. Apparently, the method of signal
segmenting by zero-crossings does not compromise the prototype construction.
Qualitatively, the prototypes are biphasic responses, where the positive peak
is almost always sharper than the negative one and the downward slope is
steeper than the rising one. We also observe that different volunteers provide
different waveforms, e.g., the distances between zero-crossings and other
waveform details are volunteer dependent.

Figure 2: PPG prototype waveforms of Volunteers 2-16 and 2-18. For both
volunteers, two prototypes were constructed where one method is ECG agnostic,
and the other is not. The interquartile range (IQR) of the prototype
associated ECG-based segmentation is shown as the dashed lines.

One of the participants showed clearly deviant PPG behavior. The common PPG
prototype has a relatively narrow peak shortly after the R-peak maximum. This
can be understood from the fact that the R-peak is the onset of the heart
muscle contraction, which is followed by opening of the valves and a fresh
blood bolus being pushed through the vessels, finally to arrive at the
observation site where the PPG signal reflects the blood volume change.
Deviant behavior was defined as a PPG were the maximum of the prototype is
relatively broad (in comparison to the valley) and the maximum occurs at or
before the R-peak. It is assumed that some signal disturbance caused this
atypical behavior. The data from this participant have been excluded from the
next analyses.

###  4.2 Prototypes: breath-holding

For 6 participants of the second session, prototypes were created for the
breath-holding epochs and non-breath-holding epochs. The latter epochs thus
include preparation and recovery of breath-holding phases. In Fig. 3, the
constructed PPG prototypes are shown for two volunteers, where the prototypes
for the two conditions are given. The constructed prototypes have not been
smoothed nor normalized (in time or amplitude). One PPG prototype is cleaner
and of shorter duration, the other longer and more noisy (as inferred from the
IQR). For all considered participants, the prototype waveforms from these
different epochs agree closely similar to those shown in Fig. 3, i.e., with
curves falling inside the IQR range.

Figure 3: PPG prototype waveforms of Volunteers 2-01 and 2-07. For both
volunteers, two prototypes were constructed for the breath-holding and non-
breathholding (normal) epochs. The interquartile range associated with the
normal prototype is indicated by the dashed lines.

###  4.3 Prototypes: IBI dependence

For analysis of intraperson variation of the PPG waveforms, the prototypes
from ECG-based segmentation were used. Each cardiac cycle was labeled
according to its IBI. Three IBI bins were created: the 25% lowest and highest
IBI (denoted as Q1 and Q4) and the remaining mid range (Mid 50%). This was
done for each participant separately and for each bin, the corresponding PPG
prototype was created. In most cases, the median IBI of both outer bins was
less than 5% from median IBI in the central bin. No clear differences between
the prototypes were observed in these cases. For two volunteers having a large
IBI variation within the session, the prototypes are shown in Fig. 4 where
amplitude normalization was applied for clarity of the graphs. These plots
suggest that the PPG waveform changes with increasing IBI mainly by time
stretching of its rising part; the down-going part appears hardly affected.

Figure 4: PPG prototype waveforms of Volunteers 1-06 and 1-08. For both
volunteers, three prototypes were constructed based on IBI subdivision. The
interquartile range of the prototype associated with the central IBI bin is
shown as the dashed lines.

For a more detailed analysis, the prototypes were smoothed by modeling the
prototype as a harmonic series consisting of five components: fundamental and
four harmonics. The removal of higher harmonic results in smooth prototypes
without losing essential information as is shown in the next section.

From the smoothed prototypes, four features were extracted. These are the
positions of the markers M, F and D (relative to the R-peak) and the signal
amplitude at the maximum. These were established for all three bins, with the
bins denoted as D, N and I reflecting the decreased, normal and increased IBI
categories. Each participant has a median IBI in these three categories
denoted as $T_{D}(p)$, $T_{N}(p)$ and $T_{I}(p)$ with $p$ indicating the
participant. Considered were cases with small (S) and large (L) increases and
decreases in duration quantified by the ratios

$R_{D}(p)=T_{D}(p)/T_{N}(p)\mbox{ and }R_{I}(p)=T_{I}(p)/T_{N}(p).$ (2)

The D data was split into two sets: participants with $R_{D}$ less than the
median or not, splitting the D data in 13 and 14 participants, respectively.
The sets are denote as LD and SD corresponding to participants exhibiting a
large and small decrease of IBI. Similarly, the I data was split into two
sets: participants with $R_{I}$ greater than the median or not, splitting the
I data in 13 and 14 participants, respectively. The sets are called LI and SI
indicated a large and small increase of IBI. A scatter plot of the ratios
$R_{I}$ and $R_{D}$ as a function of pulse rate (derived from the N data) is
shown in Fig. 5 and clarifies the subdivision into the four categories.

Figure 5: Intraperson IBI variation categories LD, SD, SI and LI characterized
by the median heart rate of the interquartile range (vertical axis) and the
ratio of the outer quartiles and interquartile median IBIs (horizontal axis).
Blue asterisks: $R_{D}$, red asterisks: $R_{I}$.

Boxplots of the amplitude and the position change associated with four
categories of IBI change are given in Fig. 6. The plots indicate the RD
distance changes significantly and consistently with the IBI change. For the
other metrics (amplitude, RM and RF distance) it would require more data to
support the hypothesis of a change, and if so, the change would be small
(amplitude change of less than 1 dB, position shifts of several ms).

Figure 6: Boxplots of feature variation in the four intraperson IBI variation
categories LD, SD, SI and LI (see text). Top row: amplitude change and
position shift of the maximum. Bottom row: position shift of F and D marker.

###  4.4 Harmonic model

The prototype PPG waveforms were modeled in a Fourier series of different
order $M$, see (1). The fraction of unmodeled energy (in dB) is shown as a
function of order of the harmonic series over all 25 volunteers as a boxplot
in Fig. 7. Beyond order $M=2$, central 50% ranges of adjacent orders overlap
to a high degree. It shows that with $M=2$ almost all energy is captured by
the truncated Fourier series: for 50% of the cases the modeling error is less
than 20 dB.

Figure 7: Boxplot of unmodeled energy in the harmonic model for PPG prototype
waveforms as a function of the number of harmonics.

The modeling was also analyzed for dependence on IBI. This is done by creating
3 IBI bins and comparing the results of two outer bins, i.e., containing the 9
lowest and highest IBIs. In Fig. 8, the unmodeled energy of this subdivision
is shown suggesting that for low IBI (high heart rate) the modeling with 2
components gives better results than in case of high IBI. For orders $M=0$ and
$2$, the distributions overlap. In other words, the first harmonic component
is stronger for high IBI than for low.

Figure 8: Boxplot of unmodeled energy in the harmonic model ($M=0,1,2$) for
PPG prototype waveforms associated with low and high IBI.

###  4.5 PPG features

In this section the PPG markers and their relation to the ECG timing is
considered. Fig. 9, the various PPG features (M, F, D and ZH) are plotted as a
function of IBI. We see that there is a general trend over the various
volunteers. For the M, F and D features, this can be modeled as a low-order
polynomial (dashed lines). There occur however considerable deviations from
these models per volunteer: deviations up to 100 ms can be observed for
individual cases.

Figure 9: Scatter plots of distance of PPG markers to the R-peak as determined
from prototypes.

For ECG-blind triggering, it would be convenient if the position of the R-peak
can be predicted from the PPG signal itself. Essentially, this means that it
is desired that the offset between any of the features and the R-peak can be
accurately predicted. This is similar to the attempts of pulse arrival time
(PAT) prediction. Various combinations of PPG features were tried. Based on
the current data, the most successful prediction $\hat{d}$ was that of the
distance $d$ between the R-peak and M. This distance $d$ can be modeled as

$d=\hat{d}+e=K_{1}+K_{2}d_{ZM}+e$ (3)

where $e$ is the prediction error, $K_{1}$ and $K_{2}$ are constants and
$d_{ZM}$ is the distance between the zero-crossing of the phase of the
analytic signal and the maximum M of the PPG signal, where the distance
between the PPG markers Z and M is derived from the ECG-blind prototype
construction and the R-peak to Z distance $d$ from the prototype using the ECG
segmentation. The experimental cumulative density function of the prediction
error $e$ error for this model and also for the case that the model is just a
constant. These curves highlight that using $d_{ZM}$ reduces the outliers: the
90% confidence interval reduces from 108 to 61 ms. The error spread is
quantified by its standard deviation $\sigma=20$ ms. Taking more predictors
into account (e.g., distances between M and F or F and D, etc.) did not
improve the prediction error significantly.

Figure 10: Experimental cumulative density function of $e$.The blue curve is
the distribution of $d-\mu_{d}$, with $\mu_{d}$ the mean value of $d$. The red
curve is the error distribution of $e$ with optimized parameters $K_{1}$ and
$K_{2}$.

##  5 Discussion

For the development of a PPG-based trigger mechanism for MR, it is encouraging
that there is much qualitative agreement between the PPG signals between
different subjects. However, the data pool is still limited and deviant
behavior may have been missed. The data reveals variation between PPGs from
different participants as well as variations per person during the recording
period where these intraperson variations included arrythmias for some
participants. The recruitment only targeted healthy patients and therefore
differentiation of the PPG waveform for cardiac patients could not be
considered.

We studied intraperson variability where a subdivision of the data according
to IBI was considered. Not addressed were long-term intraperson variation,
e.g., variation of the PPG signal of a specific subject at different days. The
constructed PPG waveforms suggest that various phases in the PPG vary
differently with a change in IBI. The phase of blood volume increase appears
to be as good as constant while the phase of blood volume decrease clearly
decelerates with increased IBI. This agrees qualitatively with ECG knowledge:
there the systolic phase is less affected by a change in IBI than the
diastolic phase. The models in [8], [9] predict that at an IBI of 500 ms, the
systolic and diastolic duration are about equal (250 ms), while at an IBI of
1500 ms this changes to a 450 and 1150 ms subdivision, i.e., a much larger
increase in the diastolic phase than in the systolic one.

The PPG signal at the forehead can be described well by a low order harmonic
model. For the prototypes, a Fourier series with two harmonics provides
already an accurate description of the measurement data. The relative signal
strength of the harmonic components appears dependent on the IBI.

A model was created for blind prediction of the R-peak position from the PPG
prototype. This model was trained only for PPG data from forehead and for a
limited number of subjects. Further model validation is recommended. We also
note that the model is specific to the (used) camera setting: changing the
shutter profile may change the data and thus the parameterization of the
model. The model is extremely simple, depending on a single PPG feature only.
In contrast, in many blood pressure research, multiple PPG features have been
used and tested for blood velocity, PAT and PPT extraction. Similar to our
finding, a single feature was proposed in [10] for PTT of contact PGG at
different body sites (finger, forehead, ear). This feature (the maximum PPG
slope) was tested for our purpose as well but provided poor prediction
results.

The PPG signal was determined from a camera operating at frame rate of 40 fps.
This is much lower than the sampling rate of a typical ECG signal, yet it is
considered as a sufficiently high rate for accurate triggering for the
following reason. Accurate triggering means that the produced images are not
suffering from blurring due to trigger jitter. In cine imaging, the cardiac
cycles are binned into phases. The number of phases dictated by image quality:
a too low number would lead to blurring the constructed image. A number of
$N_{p}=40$ phases is considered being definitely at the safe side: lower
numbers are common in practice ($N_{p}=10$ to $30$). For a heart rate of 60
bpm, $N_{p}=40$ effectively corresponds to a sampling frequency of 40 Hz.
Since the cine sequences do not suffer from motion artifacts, the timing
accuracy provided when operating the camera at 40 fps is presumed sufficient.

##  6 Conclusions

Application of video-based PPG monitoring for MR triggering requires detailed
knowledge of the PPG signal. The PPG signal at the forehead has been studied
for inter- and intra-person variability by constructing prototype PPG signals
per participant. In qualitative terms, the PPG prototype is a biphasic signal
with a relatively fast decrease and a slower PPG rise. Low-order Fourier
series allow accurate modeling of the prototype. The prototype is dependent on
subject and IBI, but is not dependent on whether or not the participant is
executing a breath-hold procedure. Variation of the prototype with IBI is
mainly in the rising part of the PPG signal. Person-specific prototype signals
with equal accuracy can be constructed with or without timing information from
the ECG signal. Simple models allow to predict the R-peak position within the
prototype from prototype features with an STD of around 20 ms for the
prediction error.

## Acknowledgement

Thanks for support to E. Coezijn, A. Kok, and P. Ramamurthy.

## References

  * [1] J. Allen.  Photoplethysmography and its application in clinical physiological measurement.  Physiol. Meas., 28(3):R1–39, 2007.
  * [2] W. Verkruysse, L.O. Svaasand, and J.S. Nelson.  Remote plethysmographic imaging using ambient light.  Opt. Express., 16(26):21434–21445, 2008.
  * [3] A.C. den Brinker and W. Wang.  Model-based camera-PPG.  In W. Wang and X. Wang, editors, Contactless Vital Signs Monitoring. Academic Press, Cambridge, USA, 2021.
  * [4] T. Togawa, O. Okai, and M. Oshima.  Observation of blood flow E.M.F. in externally applied strong magnetic field by surface electrodes.  Medical and Biological Engineering and Computing, 5(2):169–170, 1967.
  * [5] R. Abächerli, C. Pasquier, F. Odille, M. Kraemer, J. Schmid, and J. Felblinger.  Suppression of MR gradient artefacts on electrophysiological signals based on an adaptive real-time filter with lms coefficient updates.  Magn. Reson. Med., 18(1):41–50, 2005.
  * [6] N. Spicher, M. Kukuk, S. Maderwald, and M. E. Ladd.  Initial evaluation of prospective cardiac triggering using photoplethysmography signals recorded with a video camera compared to pulse oximetry and electrocardiography at 7T MRI.  Biomedical Engineering Online, 15(1):126, 2016.
  * [7] W. Wang, S. Weiss, A.C. den Brinker, J.H. Wuelbern, A. Garcia i Tormo, I. Pappous, and J. Sénégas.  Fundamentals of camera-PPG based magnetic resonance imaging.  IEEE Journal of Biomedical and Health Informatics, 26(9):4378–4389, 2022.
  * [8] R. Ashman and E. Hull.  The normal human electrocardiogram and changes produced by disease. the T wave, the Q-T interval, and the U wave.  In Essentials of electrocardiography for the student and practitioner of medicine, page 162–163. Macmillan, New York, 2 edition, 1945\.
  * [9] S.D. Roes, G. Korosoglou, M. Schär, J.J. Westenberg, M.J.P. van Osch, A. de Roos, and M. Stuber.  Correction for heart rate variability during 3D whole heart MR coronary angiography.  J. of Magnetic Resonance Imaging, 27:1046–1053, 2008.
  * [10] P.S. Addison.  Slope transit time (STT): A pulse transit time proxy requiring only a single signal fiducial point.  IEEE Trans. Biomed. Eng., 63 (11):2441–2444, 2016.


[style=chinese,orcid=0000-0003-1816-1024] [1]

[style=chinese,orcid=0000-0002-8990-094X] [1]

[cor1]Corresponding authors

# Data-driven, multi-moment fluid modeling of Landau damping

Wenjie Cheng  Haiyang Fuenvelope Liang Wang  Chuanfei Dongenvelope Yaqiu Jin
Mingle Jiang  Jiayu Ma  Yilan Qin  Kexin Liu  School of Information Science
and Engineering, Fudan University, Shanghai, 200433, China  Key Laboratory for
Information Science of Electromagnetic Waves (MoE), Fudan University,
Shanghai, 200433, China  Department of Astrophysical Sciences, Princeton
University, Princeton, New Jersey, 08544, USA  Princeton Plasma Physics
Laboratory, Princeton University, Princeton, New Jersey, 08540, USA

###### Abstract

Deriving governing equations of complex physical systems based on first
principles can be quite challenging when there are certain unknown terms and
hidden physical mechanisms in the systems. In this work, we apply a deep
learning architecture to learn fluid partial differential equations (PDEs) of
a plasma system based on the data acquired from a fully kinetic model. The
learned multi-moment fluid PDEs are demonstrated to incorporate kinetic
effects such as Landau damping. Based on the learned fluid closure, the data-
driven, multi-moment fluid modeling can well reproduce all the physical
quantities derived from the fully kinetic model. The calculated damping rate
of Landau damping is consistent with both the fully kinetic simulation and the
linear theory. The data-driven fluid modeling of PDEs for complex physical
systems may be applied to improve fluid closure and reduce the computational
cost of multi-scale modeling of global systems.

###### keywords:

Data-driven modeling Multi-moment fluid closure Machine learning Kinetic model
data Landau damping

##  1 Introduction

A plasma system can be described by either microscopic kinetic theory that
involves the dynamics of individual particles or macroscopic fluid equations.
The kinetic theory can describe the distribution of particles in phase space.
Although the kinetic model is accurate, the computational cost is high for
large-scale problems. Most fluid models assume that the velocity distribution
of particles obeys the Maxwellian distribution, which has a low computational
cost with macroscopic averaged quantities such as density, velocity, and
pressure from kinetic models. However, kinetic effects are not negligible in
certain plasma problems such as ion-temperature-gradient instability [1],
collisionless magnetic reconnection [2] and ionospheric modification
experiment [3, 4]. Therefore, it is desired that macroscopic fluid models can
integrate kinetic physics to resemble the kinetic model but with a cheap
computational cost especially for global systems [5, 6, 7, 8, 9, 10, 11].

With progress in data-driven modeling of complex dynamical systems during the
past decade, it is possible to extract physical laws and partial differential
equations (PDEs) from real data. _Schmidt et al._ [12] proposed distilling
natural laws from data using evolutionary symbolic regression to discover
analytic relations automatically from motion tracking data without any
physical priors. Symbolic regression based on neural networks has also been
investigated for multiple physical data sets [13, 14, 15]. Then, _Brunton et
al._ [16] and _Rudy et al._ [17] proposed a sparse regression method PDE-FIND
to discover the governing physical PDEs via Pareto analysis from
spatiotemporal data. Recently, _Raissi et al._ [18, 19] proposed a physics-
informed neural network (PINN) that integrates the physical model to solve and
learn Navier Stokes fluid equations, which are applied to retrieve velocity
and pressure field from density data [20]. _Long et al._ [21, 22] proposed the
PDE-NET method, which approximates the spatial differentiation through the
convolutional neural network, and can realize long-time prediction even in
noisy data. For electromagnetic problems, _Xiong et al._ [23] and _Fu et al._
[24] combined the sparse regression with the PDE-FIND to build a time network
for learning Maxwell’s equations and electromagnetic wave equations,
respectively.

It has been pursued that accurate fluid models with kinetic physics can be
obtained. The earliest work could be traced back to Grad’s work [25]. Based on
Grad’s work, _Henning et al._ [26] derived the well-known 13-moment equations
by truncating the third-order distribution function in Hermite polynomials.
Recently, machine learning has been adopted to learn multi-moment fluid models
from kinetic data. _Han et al._ [27] learned fluid closure for the
hydrodynamic equations to approximate the kinetic model with high accuracy in
the absence of scale separation. _Zhang et al._ [28] adopts the PDE-FIND
algorithm to learn the governing equations of fluid dynamics from simulation
data of the Direct Simulation Monte Carlo (DSMC) method. Distinct from
learning explicit PDEs as mentioned above, _Ma et al._ [29] applied machine
learning to learn the surrogate models for the well-known Hammett-Perkins(H-P)
closure [1] in plasma physics. In addition, _Wang et al._ [30] also accurately
reproduces the Landau damping rate and is more accurate than H-P closure, and
_Maulik et al._ [31] reproduced the crucial physics inherent in known
magnetized plasma closures. Thus, machine learning has shown promising
potential to learn physical PDEs and accurate models from data.

Although the Hammett-Perkins (H-P) closure is simple in Fourier space, it is
numerically challenging to implement in fluid simulation in configuration
space. The closure by neural networks learned in the Ref. [29] [30] [31] are
the black box that can not be physically interpretable. The limitation in Ref.
[29] is that only closure has been trained, and the training data comes from
the well-known Hammett-Perkins (H-P) closure itself.

The purpose of our work is motivated by discovering the whole multi-moment
fluid model PDEs directly from kinetic data based on machine learning. The
closure form can not only be physically interpretable to capture kinetic
effects but also is expected to be easily applied in configuration space for
large-scale fluid simulations. Therefore, we attempt to construct a data-
driven neural network architecture called _mPDE-Net (Multiple Partial
Differential Equations - Network)_ to learn multi-moment fluid PDEs with
closure from kinetic data. The _mPDE-Net_ method is extended to learn
multiple-moment fluid model based on the work of _Xiong et al._ [23] and the
original method in _Long et al._ [21]. We specially develop a deep learning
neural network in conjunction with sparse regression and Pareto analysis to
deal with undefined terms. The _mPDE-Net_ method attempts to obtain a refined
equivalent form of Vlasov equations based on machine learning to mimic the
theoretical derivation of multiple moments by human beings to some extent.

The proposed _mPDE-Net_ has been tested for 1-D Landau damping as the proof-
of-concept demonstration. A multi-moment fluid with closure model has been
successfully learned from kinetic Vlasov data with a small amount of data
sampling. Based on the learned multi-moment fluid PDEs, the time evolution of
multiple physical parameters have been correctly predicted with initial and
periodic boundary condition, which agrees with the kinetic landau damping rate
well. This paper is organized as follows. Sec. 2 describes methodology of
_mPDE-Net_ and _Time prediction_. Sec. 3 describes simulation models of Landau
damping of Langmuir waves, including kinetic model, learned fluid PDEs and
multiple parameter prediction. The Landau damping of Langmuir waves learned by
_mPDE-Net_ is compared with results from theoretical predictions and kinetic
models. Finally, the conclusion and discussion will be summarized in Sec. 4.

Figure 1: The schematic diagram of data-driven modeling and prediction
architecture for Landau damping of Langmuir waves: learn the macroscopic
multi-moment fluid model from kinetic data (numerical solution of Vlasov-
Ampére equations), and then make predictions based on learned fluid PDEs with
initial and boundary conditions.

##  2 Methodology

The method proposed (_mPDE-Net_) in this paper is an extension of the work of
_Xiong et al._ [23]. As we know, the multi-moment fluid model consists of
three parts: the moment fluid equations, the field equations, and the fluid
closure equations. In this work, we will focus on the one-dimensional (1-D)
electrostatic case, and high-dimensional studies can be extended with the same
method. The physical quantities involved in the 1-D electrostatic problem
include density $n(x,t)$, average velocity $u(x,t)$, pressure $p(x,t)$,
temperature $T(x,t)$, electric field $E_{x}(x,t)$ and heat flux$q(x,t)$. We
assume that these physical quantities vary with time $t$ and the spatial
coordinate $x$, which can be expressed in the form of a set of partial
differential equations (PDEs),

$\left\\{\begin{array}[]{l}\vspace{1.5ex}\frac{{\partial n}}{{\partial
t}}={f_{n}}(n,u,{E_{x}},\frac{{\partial n}}{{\partial x}},\frac{{\partial
u}}{{\partial x}},\frac{{\partial{E_{x}}}}{{\partial x}})\\\
\vspace{1.5ex}\frac{{\partial u}}{{\partial
t}}={f_{u}}(n,u,p,{E_{x}},\frac{{\partial n}}{{\partial x}},\frac{{\partial
u}}{{\partial x}},\frac{{\partial p}}{{\partial
x}},\frac{{\partial{E_{x}}}}{{\partial x}})\\\ \vspace{1.5ex}\frac{{\partial
p}}{{\partial t}}={f_{p}}(n,u,p,q,{E_{x}},\frac{{\partial n}}{{\partial
x}},\frac{{\partial u}}{{\partial x}},\frac{{\partial p}}{{\partial
x}},\frac{{\partial q}}{{\partial x}},\frac{{\partial{E_{x}}}}{{\partial
x}})\\\ \vspace{1.5ex}\frac{{\partial T}}{{\partial
t}}={f_{T}}(n,u,T,q,{E_{x}},\frac{{\partial n}}{{\partial x}},\frac{{\partial
u}}{{\partial x}},\frac{{\partial T}}{{\partial x}},\frac{{\partial
q}}{{\partial x}},\frac{{\partial{E_{x}}}}{{\partial x}})\\\
\vspace{1.5ex}\frac{{\partial{E_{x}}}}{{\partial
t}}={f_{{E_{x}}}}(n,u,{E_{x}},\frac{{\partial n}}{{\partial
x}},\frac{{\partial u}}{{\partial x}},\frac{{\partial{E_{x}}}}{{\partial
x}})\\\ \vspace{1.5ex}q={f_{q}}(n,u,T,\frac{{\partial n}}{{\partial
x}},\frac{{\partial u}}{{\partial x}},\frac{{\partial T}}{{\partial
x}})\end{array}\right.$ (1)

where $f_{i},i\in\\{n,u,p,T,E_{x},q\\}$ are the nonlinear functions of all
candidate terms in parentheses and their cross-terms. For the equations of
$f_{i}$, where $i\in\\{n,u,p,T\\}$, the candidate terms include itself,
higher-order moment, the external force, and their spatial differential terms.
For $f_{{E_{x}}}$, the candidate terms include itself, fluid moment, and their
spatial differential terms. For the third-order moment, $q$, the candidate
terms in the closure equation ${f_{q}}$ include lower-order moments and their
spatial differential terms. It is noteworthy that we will learn the PDEs of
the Vlasov-Ampére system in this study given that the Vlasov-Ampére and
Vlasov-Poisson systems are equivalent if the initial conditions of the Vlasov-
Ampére system satisfy Poisson’s equation [32].

Fig. 1 depicts the overall framework of the proposed method, including
spatiotemporal sampling, training, and time prediction. First of all, samples
observed or calculated from spatio-temporal data in the small-scale kinetic
region are fed into _mPDE-Net_. Secondly, _mPDE-Net_ is activated to learn a
set of PDEs which are capable of describing the fluid model with Landau
damping. Finally, _Time prediction_ is adopted to predict multiple parameters
based on learned multi-moment fluid models with closure.

Figure 2: The schematic diagram of the _mPDE-Net_ architecture for data-driven
modeling of the multi-moment fluid model with Landau damping

###  2.1 mPDE-Net

The _mPDE-Net_ architecture is shown in Fig. 2. The derivative of quantities
with respect to space is calculated based on convolution [33, 34]. The central
difference with the second-order accuracy is used to calculate the first-order
and second-order spatial derivatives, respectively. The expressions of the
convolution operator with the second order accuracy are as follows:

$\begin{split}h_{x}=\frac{1}{2}\begin{pmatrix}-1&0&1\end{pmatrix}\\\
h_{xx}=\begin{pmatrix}1&-2&1\end{pmatrix}\end{split}$ (2)

Then, the circular convolution (labeled as $\otimes$ ) between filters and
function $f$ for instance is expressed as:

$\begin{split}h_{x}\otimes f\approx\delta_{x}\dfrac{\partial f}{\partial x}\\\
h_{xx}\otimes f\approx\delta_{x}^{2}\dfrac{\partial^{2}f}{\partial
x^{2}}\end{split}$ (3)

where $\delta_{z}$ is the spatial size of $z$ direction of function $u$.

Taking velocity $u(x,t)$ as $f$ for instance, its spatial derivative with
respect to $x$ can be calculated by the following formula:

$\begin{split}\dfrac{\partial u}{\partial
x}\approx\dfrac{1}{2\delta_{x}}\begin{pmatrix}-1&0&1\end{pmatrix}\otimes u\\\
\dfrac{\partial^{2}u}{\partial
x^{2}}\approx\dfrac{1}{\delta_{x}^{2}}\begin{pmatrix}1&-2&1\end{pmatrix}\otimes
u\end{split}$ (4)

Eq.(4) is also applicable for calculating spatial derivatives of other
variables $u$, $p$, $T$, $E_{x}$, $q$. In order to ensure accuracy, we do not
calculate the spatial differentiation on the edges.

The time feedforward network is constructed by the forward Euler numerical
scheme. We define $n(\cdot,t_{i})$ as all spatial value of density $n$ at
$t=t_{i}$, and $\tilde{n}(\cdot,t_{i+1})$ is the predicted value at the next
time $t_{i+1}$. According to the forward Euler method, the set of Eq (1) can
be discretized for density $n$ as an example

$\tilde{n}(\cdot,t_{i+1})=n(\cdot,t_{i})+(\Delta t)\cdot f_{n}(\cdot,t_{i})\\\
$ (5)

where $\tilde{n}(\cdot,t_{i+1})$ is the approximation value at $t_{i+1}$, and
$\Delta t$ is the time step.

The forward Euler $\Delta t$ block is similar to [21] as a layer of the neural
network to advance variables based on equations in time. In principle, the
structure of the $\Delta t$ block is an interpretation of Eq (5). The input of
this network is possible candidate terms with $f$. The prior knowledge will
help reduce the number of possible candidate functions associated with PDEs.
The output of $\Delta t$ block is the predicted value at the next time
$\tilde{u}(\cdot,t_{i+1})$. In general, the entire network consists of several
$\Delta t$ blocks by continuously connecting each block $\Delta t$, where the
number of $\Delta t$ block is determined by certain convergence criteria.

In the neural network, the weights are shared in each layer, which serves as
an output vector for the model. We define the weight vector W as

$\textbf{W}=[w_{1},w_{2},w_{3},w_{4},w_{5},...,w_{m}],$ (6)

where the non-zero weight component represents the term of the objective
equation, and $m$ represents the number of candidate terms. To reduce
unnecessary candidates, we use a sparse regression method [17], which uses
$L_{1}$ regularization. The cost function $L$ of the density $n$ equation is
defined as

$L=\sum_{i=1}^{N}\parallel\tilde{n}(\cdot,t_{i})-n(\cdot,t_{i})\parallel_{2}+\lambda_{n}\parallel\textbf{W}\parallel_{1}$
(7)

where $\lambda$ is the regularization parameter that affects the weight of the
$L_{1}$ regularization term, and $N$ is the time step. If a reasonable
inversion result is not obtained initially with manually selected $\lambda$,
we will automatically select the hyperparameter $\lambda$ using Grid Search
(GS) optimization.

Pareto analysis is a technique that can help identify top priority candidates
of the concerned problem. The application of Pareto analysis is the final
guarantee for the correct discovery of complex systems [17]. By cutting the
term with the minimum parameter out and training again, we compare the later
loss with the previous one. If loss reduces, then we repeat the process until
loss increases.

Based on the above modules, the construction procedure of the _mPDE-Net_
architecture can be described as following

(1) Determine the number of required equations $m$. If the number of
categories of physical variables input by the network is $m_{0}$, then $m$
must satisfy $m\leq{m_{0}}$. When $m\leq{m_{0}}$, the system that needs to be
learned is not closed; for $m={m_{0}}$, the system is closed.

(2) Determine the order of the equation and the sampling data, calculate the
spatial differentiation, and build the feedforward network. If equations are
discretized in first-order, the first-order forward Euler is used to build the
feedforward time network as shown in Eq. (5). If the equation is zeroth order,
the network is built directly according to the prior formula.

(3) Use sparse regression to optimize the weight W. Add $L_{1}$ regularization
and use $Adam$ optimizer to train the network.

(4) Use Pareto analysis to segment the weight W. Based on each equation,
Pareto analysis is used to optimize the weight W to achieve the most
simplified form under the premise of ensuring accuracy.

(5) Repeat iteration until convergence, and finally output the identified
equations.

###  2.2 Time prediction

Based on the learned multi-moment fluid model and initial conditions and
periodic boundary conditions, the time prediction of the physical quantities
can be written as follows,

$\displaystyle\tilde{n}(\cdot,{t_{i+1}})=\Delta
t\cdot{f_{n}}+n(\cdot,{t_{i}})$ (8)

$\tilde{u}(\cdot,{t_{i+1}})=\Delta t\cdot{f_{u}}+u(\cdot,{t_{i}})$ (9)

$\tilde{p}(\cdot,{t_{i+1}})=\Delta t\cdot{f_{p}}+p(\cdot,{t_{i}})$ (10)

$\tilde{T}(\cdot,{t_{i+1}})=\Delta t\cdot{f_{T}}+T(\cdot,{t_{i}})$ (11)

${{\tilde{E}}_{x}}(\cdot,{t_{i+1}})=\Delta
t\cdot{f_{E_{x}}}+{E_{x}}(\cdot,{t_{i}})$ (12)

$\tilde{q}(\cdot,{t_{i+1}})={f_{q}}$ (13)

Figure 3: The illustration of data sampling area for training from the kinetic
model. The transparent box (white shaded) in the top figure depict the space-
time sampling interval for $t=[2.5\omega_{pe}^{-1},3\omega_{pe}^{-1}]$ with
$500$ points, $x=[0.75\pi/k,1.25\pi/k]$ with $48$ points. The bottom figure
shows the zoomed sampling area for variables $n$, $u$, $p$, $T$, $q$, $E_{x}$
from kinetic data.

##  3 Simulation Test and Analysis

###  3.1 Kinetic model data for training

For this test, we study the Landau damping of a Langmuir wave in 1X1V by
evolving the plasma phase-space distribution function $f\left(x,v,t\right)$.
Landau damping is one of the most fundamental plasma processes, and accurate
modeling this process is an important capability of any kinetic plasma model.
Numerically, we consider an immobile ion background with density $n_{i}=n_{0}$
and a perturbed electron plasma with density

$n_{e}\left(x,t=0\right)=n_{0}\left(1+A\cos\left(kx\right)\right),$ (14)

where $k$ and $A$ are the perturbation wavenumber and amplitude, respectively,
and the electric field is perturbed to satisfy Gauss’s law. Unit normalization
$\varepsilon_{0}=1$, $q_{e}=-1$, $m_{e}=1$, $n_{0}=1$, and $T_{e}=1$ is used
so that the electron plasma frequency
$\omega_{pe}=\sqrt{n_{e}q_{e}^{2}/m_{e}\varepsilon_{0}}=1$ and Debye length
$\lambda_{De}=\sqrt{\varepsilon_{0}T_{e}/n_{e}q_{e}^{2}}=1$. The perturbation
wavenumber $k$ and amplitude $A$ are set to $0.35\lambda_{De}^{-1}$ and 0.05.
The periodic configuration domain of length $2\pi/k$ is discretized by 64
cells. The velocity space domain $\left[-6v_{{\rm th},e},\ 6v_{{\rm
th},e}\right]$ is discretized by 64 cells, where $v_{{\rm
th},e}=\sqrt{T_{e}/m_{e}}=1$ is the background electron thermal velocity.
Second-order serendipity polynomial bases are used in space and the time
integration uses the strong-stability preserving Runge-Kutta scheme with a
fixed time step $\Delta t=0.001$ and stops at $t=30\omega_{pe}^{-1}$.

We obtain fluid quantities by calculating the moments of ${f}(x,v,t)$ in the
velocity space, including electron density $n(x,t)$, average velocity
$u(x,t)$, pressure $p(x,t)$, temperature $T(x,t)$ and heat flux $q(x,t)$.

$n(x,t)=\int{f(x,v,t)dv}$ (15)

$u(x,t)=\frac{1}{{n(x,t)}}\int{vf(x,v,t)dv}$ (16)

$p(x,t)=m_{e}\int{{{(v-u(x,t))}^{2}}f(x,v,t)dv}$ (17)

$q(x,t)=m_{e}\int{{{(v-u(x,t))}^{3}}f(x,v,t)dv}$ (18)

$T(x,t)=\frac{{p(x,t)}}{{n(x,t)}}$ (19)

The Vlasov-Ampére equation is solved to obtain the kinetic data with the
numerical code Gkeyll [35], including electric field ${E_{x}}(x,t)$, electron
distribution function $f(x,v,t)$ and its moments in Eqs. (15-18). The electron
temperature $T(x,t)$ is calculated from Eq. (19).

###  3.2 Data-Driven Multi-Moment Fluid Simulation Results

Table 1: Comparison between ten-moment fluid model with gradient closure (Reference) and learned multi-moment fluid model by _mPDE-Net_ (Learned) Reference |  Learned
---|---
$\frac{{\partial n}}{{\partial t}}+u\frac{{\partial n}}{{\partial x}}+n\frac{{\partial u}}{{\partial x}}=0$ |  $\frac{{\partial n}}{{\partial t}}+1.0008n\frac{{\partial u}}{{\partial x}}+0.9972u\frac{{\partial n}}{{\partial x}}=0$
$\frac{{\partial u}}{{\partial t}}+u\frac{{\partial u}}{{\partial x}}+\frac{1}{{m_{e}n}}\frac{{\partial p}}{{\partial x}}=\frac{q_{e}}{m_{e}}{E_{x}}$ |  $\frac{{\partial u}}{{\partial t}}+1.0035u\frac{{\partial u}}{{\partial x}}+1.0038\frac{1}{{m_{e}n}}\frac{{\partial p}}{{\partial x}}=0.9978\frac{{{q_{e}}}}{{{m_{e}}}}{E_{x}}$
$\frac{{\partial p}}{{\partial t}}+u\frac{{\partial p}}{{\partial x}}+3p\frac{{\partial u}}{{\partial x}}+\frac{{\partial q}}{{\partial x}}=0$ |  $\frac{{\partial p}}{{\partial t}}+0.9942u\frac{{\partial p}}{{\partial x}}+2.9990p\frac{{\partial u}}{{\partial x}}+1.0017\frac{{\partial q}}{{\partial x}}=0$
$\frac{{\partial T}}{{\partial t}}+2T\frac{{\partial u}}{{\partial x}}+u\frac{{\partial T}}{{\partial x}}+\frac{1}{n}\frac{{\partial q}}{{\partial x}}=0$ |  $\frac{{\partial T}}{{\partial t}}+1.9991T\frac{{\partial u}}{{\partial x}}+0.9971u\frac{{\partial T}}{{\partial x}}+1.0013\frac{1}{n}\frac{{\partial q}}{{\partial x}}=0$
$\frac{{\partial{E_{x}}}}{{\partial t}}+\frac{{{q_{e}}}}{{{\varepsilon_{0}}}}nu=0$ |  $\frac{{\partial E_{x}}}{{\partial t}}+1.0005\frac{{{q_{e}}}}{{{\varepsilon_{0}}}}nu=0$
$q={\kappa}n\sqrt{T}\frac{{\partial T}}{{\partial x}}$ |  $q=-4.2655n\sqrt{T}\frac{{\partial T}}{{\partial x}}+\underbrace{\cancel{{\kappa_{1}}\frac{{{\partial^{2}}T}}{{\partial{x^{2}}}}}}_{0}+\underbrace{\cancel{{{\kappa_{2}}\frac{{\partial n}}{{\partial x}}\frac{{\partial T}}{{\partial x}}}}}_{0}$
Figure 4: The result of heat flux $q$ with learned closure form in Table 1 .
The left figure is the heat flux from kinetic data and the right is the heat
flux calculated by the learned heat flux closure in Table 1 with temperature
and density data from kinetic simulations. The time in
$[0,30{\omega_{pe}^{-1}}]$ has $30000$ points, and the space in $[0,2\pi/k]$
has $192$ points. Figure 5: Physical quantities calculated from the kinetic
data (left) and predicted with the learned multi-moment fluid model (right).
Each panel shows density $n$, velocity $u$, pressure $p$, temperature $T$,
electric field $E_{x}$ and heat-flux $q$. The time is in
$[0,30{\omega_{pe}^{-1}}]$ with $30000$ points and in space in $[0,2\pi/k]$
with $192$ points.

####  3.2.1 Constructing the data-Driven Multi-Moment Fluid Model

Based on our $1X1V$ kinetic simulation, we construct a machine learning-based
neural network to learn a multi-moment fluid model from kinetic data as shown
in Fig.2.

We train the kinetic data $n$, $u$, $p$, $T$, $q$ and $E_{x}$ from Gkeyll in
the sampling area by the _mPDE-Net_ algorithm to learn the form of the multi-
moment fluid model equations. The sampling data is shown on the top panel of
Fig. 3, where the transparent box are temporal sampling interval
$[2.5\omega_{pe}^{-1},3\omega_{pe}^{-1}]$ with $500$ points and spatial
sampling interval $[0.75\pi/k,1.25\pi/k]$ with $48$ points. For clarity, the
sampling scheme is zoomed out on the bottom panel of Fig. 3 for all variables
$n,u,p,T,q,E_{x}$ as Fig. 3 .

According to previous studies on the heat flux closure [1, 2], it is generally
believed that there exists a relationship between heat flux and temperature
derivatives. The original Hammett-Perkins closure [1] in Fourier space can be
expressed as

$\tilde{q}_{k}=-n_{\rm 0}\chi\frac{2^{1/2}v_{t}}{{|k|}}ik\tilde{T}_{k}$ (20)

where $\chi$ is dimensionless coefficients, $v_{t}$ is the electron thermal
speed, $\tilde{T}_{k}=\tilde{p}-T_{\rm 0}\tilde{n}/n_{\rm 0}$ is the Fourier
transformation of the perturbed temperature.

As mentioned in [2], the resulting heat flux for single-mode can be expressed
in configuration space as

$q=-n\chi\frac{2^{1/2}v_{t}}{|k_{s}|}\frac{\partial{T}}{\partial{x}}$ (21)

where $k_{\rm s}$ is a constant or spatially varying parameter.

It is important to note that the _mPDE-Net_ algorithm is capable of not only
obtaining proper coefficients of terms for known PDEs but also removing
trivial terms for unknown PDEs for the whole multiple-moment fluid system of
equations. Table 1 shows the reference (left-column) and learned (right-
column) multi-moment-Ampere equations, including the fluid closure. It should
be noted that our training data $(n,u,p,T,E_{x},q)$ is not constrained by the
exact reference equations. Instead, we supplied a number of additional
candidate terms to the neural network to allow it to identify the
correct/expected terms based on the nonlinear physics contained in fully
kinetic simulation data, eliminating the irrelevant candidates. Particularly,
for the heat flux $q$ closure equation, we supplied candidate terms such as
$\frac{{{\partial^{2}}T}}{{\partial{x^{2}}}}$, $\frac{{\partial n}}{{\partial
x}}\frac{{\partial T}}{{\partial x}}$, _mPDE-Net_ was able to discover the
expected term $\kappa n\sqrt{T}\frac{\partial{T}}{\partial{x}}$ for $q$ with
proper coefficients $\kappa$ and remove the redundant terms in Table 1.

Figure 6: Comparison of the evolution of electric field energy
$\int{|E_{x}|^{2}}dx$ calculated from the kinetic data (blue curve) and
predicted with the learned multi-moment fluid model (dashed yellow). The
straight lines correspond to the exponential fit to the envelope of kinetic
data (red) and learned fluid model (dashed purple) which provide the damping
rate $\gamma/\omega_{pe}$.

In order to evaluate performance, we define the relative error ${\delta}$ of
each coefficient for the known fluid PDEs as follows:

$\delta=|\dfrac{\tilde{c}-c}{c}|\times 100~{}\%$ (22)

where $\tilde{c}$ is the learned coefficient and $c$ is the theoretical
coefficient. First of all, the relative coefficient errors of the continuity
density equation are $0.08~{}\%$ and $0.28~{}\%$, respectively. Secondly, the
relative coefficient errors of the velocity equation are $0.35~{}\%$,
$0.38~{}\%$ and $0.22~{}\%$, respectively. Thirdly, the relative coefficient
errors of the pressure equation are $0.58~{}\%$, $0.03~{}\%$ and $0.17~{}\%$,
respectively. The relative coefficient errors of the temperature equation are
$0.04~{}\%$, $0.29~{}\%$ and $0.13~{}\%$, respectively. Finally, the relative
coefficient of the electric field equation is $0.05~{}\%$.

Since the theoretical coefficient of the heat flux equation is unknown, in
order to verify the heat flux closure equation in Table 1, we compare the heat
flux (right) based on learned fluid data with kinetic data (left) in Fig.4.
The comparison of the heat flux $q$ shows a similar result of heat flux
between those calculated from kinetic data and learned from the neural
network.

####  3.2.2 Time prediction

Figure 7: Sample profiles of density $n$, velocity $u$, energy $\xi$ and
electric field $E_{x}$ (from top to bottom) at $t=0,3,6$ (from left to right)
in 1-D multi-moment fluid model, obtained from the kinetic equation, the five-
moment equations, and learned multi-moment fluid model. Figure 8: The
comparison of the Landau damping rates from the linear theory, the rates
fitted to the kinetic simulation data and learned multiple-moment fluid
equations.

Now we have a multi-moment fluid model learned from kinetic data. Next, we
will activate the _Time Prediction_ block to predict all physical quantities
globally based on the learned model with initial conditions. The time
prediction can verify whether the learned multi-moment fluid model is
accurate, and it is also important for potential plasma applications.

In our numerical scheme of time prediction, the time differentiation adopts
the forward Euler and the space differentiation adopts the central difference
with periodic boundary conditions. Based on the learned multi-moment fluid
model (Table 1) and the initial conditions, the time prediction of all six
physical quantities are shown in Eq. (8) - Eq. (13) with periodic boundary
conditions.

Fig.5 shows the prediction results in comparison with kinetic data, including
density $n$, velocity $u$, pressure $p$, temperature $T$, electric field
$E_{x}$ and heat-flux $q$. Within the time interval
$[0,30{\omega_{pe}^{-1}}]$, the prediction of all six physical quantities is
in agreement with kinetic data. The errors gradually increase for a long time.
One reason is that the coefficients of the learned multi-moment fluid
equations contain small errors. The other one can be the numerical errors of
the prediction schemes.

The most critical consequence of Landau damping is the exponential damping of
the longitudinal charge waves. In Fig.6, we draw the time evolution of the
total electric field energy predicted by the kinetic model versus that
predicted by the learned fluid model. The electric field energy bounces and
decreases over time due to the damping. The local maximum values are linearly
damped as expected. The learned model’s bouncing frequency and damping rate
both agree with the kinetic data. A linear fitting of the maxima, in
particular, yields the damping rates $\gamma=-0.0384\omega_{pe}^{-1}$ and
$-0.0382\omega_{pe}^{-1}$ for the kinetic model and machine learning
prediction, respectively, with a relative error of about 0.5 percent. The
bouncing frequency differs slightly between the two models. We discovered that
the learned fluid model result contains more less-damped modes with nontrivial
contributions to the total energy after performing Fourier analysis on the
time series. These are the byproducts of the initial perturbation. These
undesirable "noise" modes are quickly damped in the fully kinetic model.
Because we assume only one mode dominates in constructing the closure in the
learned fluid model, the learned closure does not provide sufficient damping
to those noise modes that remain in the system and cause very small but
persistent phase error. Our numerical experiments demonstrated that using a
smaller time step size effectively reduces phase errors.

Fig.7 compares the 1-D profile of density $n$, velocity $u$, energy $\xi$ and
electric field $E_{x}$ profiles between the kinetic model, five-moment model
(without heat flux $q$) and learned multi-moment fluid model (with heat flux
$q$) at $t=0{\omega_{pe}^{-1}}$, $3{\omega_{pe}^{-1}}$, $6{\omega_{pe}^{-1}}$.
The data from the learned multi-moment fluid model is consistent with the
kinetic data, and shows improved accuracy than the well-known five-moment
model.

In addition, we have also test _mPDE-Net_ for different wavenumber $k$. As
shown in Fig. 8, the damping rates from the learned fluid model agree well
with the kinetic data and the results from the linear dispersion relation
solver. The absolute errors are less than $0.01$ with an optimized sampling
scheme. Therefore, the learned multi-moment fluid model by _mPDE-Net_ can
describe the kinetic Landau damping well.

##  4 Conclusion and Discussion

In this paper, we proposed a data-driven network architecture to learn the
multi-moment fluid model with Landau damping of Langmuir waves from kinetic
data. A prediction scheme of the learned multi-moment fluid model has also
been developed to achieve the prediction of Landau damping. The _mPDE-Net_
method has been verified for the 1-D electrostatic Vlasov-Ampére system to
well reproduce the Landau damping of the Langmuir wave. The Landau damping
learned by _mPDE-Net_ has shown agreement with linear theory and kinetic
simulations. The _mPDE-Net_ architecture has the capability to achieve the
physically interpretable PDEs set through machine learning.

In summary, our preliminary work adequately demonstrates the potential of
learning and predicting plasma kinetic data by multi-moment fluid equations
based on machine learning. Future work could consider a more general closure
and extend the current architecture to high-dimensional plasma systems. The
data-driven method of PDEs may pave the way to efficiently learning governing
equations of complex, multi-scale, multi-physics systems.

## Acknowledgments

The authors from Fudan University were supported by the National Key Research
and Development Program of China (2021YFA0717300) and China National Science
Foundation (42074189). C.D. was partially supported by the U.S. Department of
Energy under contract DE-AC02-09CH11466.

## References

  * [1] G. W. Hammett, F. W. Perkins, Fluid moment models for landau damping with application to the ion-temperature-gradient instability, Phys. Rev. Lett. 64 (1990) 3019–3022.
  * [2] J. Ng, A. Hakim, L. Wang, A. Bhattacharjee, An improved ten-moment closure for reconnection and instabilities, Physics of Plasmas 27 (8) (2020) 082106.
  * [3] H. Y. Fu, W. A. Scales, P. A. Bernhardt, Y. Q. Jin, S. J. Briczinski, Asymmetry in stimulated emission polarization and irregularity evolution during ionospheric electron gyroharmonic heating, Geophysical Research Letters 45 (2018).
  * [4] H. Y. Fu, M. L. Jiang, K. N. Wang, J. Wu, Y. Jin, Electron temperature inversion by stimulated brillouin scattering during electron gyro-harmonic heating at eiscat, Geophysical Research Letters (2020).
  * [5] L. Wang, A. H. Hakim, A. Bhattacharjee, K. Germaschewski, Comparison of multi-fluid moment models with particle-in-cell simulations of collisionless magnetic reconnection, Physics of Plasmas 012108 (1) (2015) 012108.
  * [6] C. Dong, L. Wang, A. Bhattacharjee, A. Hakim, Y.-M. Huang, K. Germaschewski, Magnetic reconnection in multispecies plasmas investigated by a kinetic fluid code, in: APS Division of Plasma Physics Meeting Abstracts, Vol. 2016 of APS Meeting Abstracts, 2016, p. NO7.009.
  * [7] L. Wang, K. Germaschewski, A. Hakim, C. Dong, J. Raeder, A. Bhattacharjee, Electron Physics in 3-D Two-Fluid 10-Moment Modeling of Ganymede’s Magnetosphere, Journal of Geophysical Research (Space Physics) 123 (4) (2018) 2815–2830.
  * [8] C. Dong, L. Wang, A. Hakim, A. Bhattacharjee, J. A. Slavin, G. A. DiBraccio, K. Germaschewski, Global Ten-Moment Multifluid Simulations of the Solar Wind Interaction with Mercury: From the Planetary Conducting Core to the Dynamic Magnetosphere, Geophysical Research Letters 46 (21) (2019) 11,584–11,596.
  * [9] S. Jarmak, E. Leonard, A. Akins, E. Dahl, D. R. Cremons, S. Cofield, A. Curtis, C. Dong, E. T. Dunham, B. Journaux, D. Murakami, W. Ng, M. Piquette, A. P. Girija, K. Rink, L. Schurmeier, N. Stein, N. Tallarida, M. Telus, L. Lowes, C. Budney, K. L. Mitchell, QUEST: A New Frontiers Uranus orbiter mission concept study, Acta Astronautica 170 (2020) 6–26.
  * [10] L. Wang, A. H. Hakim, J. Ng, C. Dong, K. Germaschewski, Exact and locally implicit source term solvers for multifluid-Maxwell systems, Journal of Computational Physics 415 (2020) 109510.
  * [11] E. Rulke, L. Wang, C. Dong, Three-dimensional, ten-moment multifluid simulation of the solar wind interaction with Asteroid 16 Psyche, in: AGU Fall Meeting Abstracts, Vol. 2021, 2021, pp. SM53C–08.
  * [12] M. Schmidt, H. Lipson, Distilling free-form natural laws from experimental data, Science 324 (5923) (2009) 81–85.
  * [13] S. M. Udrescu, M. Tegmark, Ai feynman: A physics-inspired method for symbolic regression, Science Advances 6 (16) (2020) eaay2631.
  * [14] S. M. Udrescu, A. Tan, J. Feng, O. Neto, M. Tegmark, Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity, arXiv e-prints (2020).
  * [15] S. Kim, P. Y. Lu, S. Mukherjee, M. Gilbert, M. Soljacic, Integration of neural network-based symbolic regression in deep learning for scientific discovery, IEEE Transactions on Neural Networks and Learning Systems PP (99) (2020) 1–12.
  * [16] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data: Sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences of the United States of America 113 (15) (2015) 3932.
  * [17] S. H. Rudy, S. L. Brunton, J. L. Proctor, J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances 3 (4) (2016).
  * [18] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations (2017).
  * [19] M. Raissi, G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, Journal of Computational Physics 357 (2018) 125–141.
  * [20] R. Maziar, Y. Alireza, K. G. Em, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science (New York, N.Y.) 367 (2020) 1026–1030.
  * [21] Z. Long, Y. Lu, X. Ma, B. Dong, Pde-net: Learning pdes from data, arXiv:1710.09668v1 (2017).
  * [22] Z. Long, Y. Lu, B. Dong, Pde-net 2.0: Learning pdes from data with a numeric-symbolic hybrid deep network, Journal of Computational Physics 399 (2019) 108925.
  * [23] B. Xiong, H. Fu, F. Xu, Y. Jin, Data-driven discovery of partial differential equations for multiple-physics electromagnetic problem, arXiv preprint arXiv:1910.13531 (2019).
  * [24] H. Fu, W. Cheng, Y. Qin, Exploration of data-driven methods for multiphysics electromagnetic partial differential equations, in: 2020 IEEE MTT-S International Conference on Numerical Electromagnetic and Multiphysics Modeling and Optimization (NEMO), 2020.
  * [25] H. Grad, On the kinetic theory of rarefied gases, Communications on Pure & Applied Mathematics 2 (4) (2010).
  * [26] H. Struchtrup, Macroscopic Transport Equations for Rarefied Gas Flows, Macroscopic Transport Equations for Rarefied Gas Flows, 2005.
  * [27] J. Han, C. Ma, Z. Ma, E. Weinan, Uniformly accurate machine learning-based hydrodynamic models for kinetic equations, Proceedings of the National Academy of Sciences 116 (44) (2019) 201909854.
  * [28] J. Zhang, W. Ma, Data-driven discovery of governing equations for fluid dynamics based on molecular simulation, Journal of Fluid Mechanics 892 (2020).
  * [29] C. Ma, B. Zhu, X. Q. Xu, W. Wang, Machine learning surrogate models for landau fluid closure, Physics of Plasmas 27 (2020).
  * [30] L. Wang, X. Q. Xu, B. Zhu, C. Ma, Y. A. Lei, Deep learning surrogate model for kinetic landau-fluid closure with collision, AIP Advances 10 (7) (2020) 075108\.
  * [31] R. Maulik, N. A. Garland, X. Z. Tang, P. Balaprakash, Neural network representability of fully ionized plasma fluid model closures, Physics of Plasmas 27 (7) (2020) 072106.
  * [32] H.-s. Xie, Constant residual electrostatic electron plasma mode in Vlasov-Ampere system, Physics of Plasmas 20 (11) (2013) 112108.
  * [33] J. F. Cai, B. Dong, S. Osher, A. Z. Shen, Image restoration: Total variation, wavelet frames and beyond, Journal of the American Mathematical Society 25 (4) (2012) 1033–1089.
  * [34] B. Dong, Q. Jiang, Z. Shen, Image restoration: Wavelet frame shrinkage, nonlinear evolution pdes, and beyond, Math.arizona.edu (2017).
  * [35] J. Juno, A. Hakim, J. Tenbarge, E. Shi, W. Dorland, Discontinuous galerkin algorithms for fully kinetic plasmas, Journal of Computational Physics (2017).