11institutetext: Ultromics Ltd, Oxford, UK King’s College London, UK
22institutetext: Sherbrooke University, Canada 33institutetext: 33email:
<EMAIL_ADDRESS>
# Left Ventricle Contouring of Apical Three-Chamber Views on 2D
Echocardiography††thanks: Supported by Ultromics Ltd.

Alberto Gomez 11 2 2 Mihaela Porumb 11 Angela Mumith 11 Thierry Judge 33 Shan
Gao 11 Woo-Jin Cho Kim 11 Jorge Oliveira 11 Agis Chartsias 11

###### Abstract

We propose a new method to automatically contour the left ventricle on 2D
echocardiographic images. Unlike most existing segmentation methods, which are
based on predicting segmentation masks, we focus at predicting the endocardial
contour and the key landmark points within this contour (basal points and
apex). This provides a representation that is closer to how experts perform
manual annotations and hence produce results that are physiologically more
plausible.

Our proposed method uses a two-headed network based on the U-Net architecture.
One head predicts the 7 contour points, and the other head predicts a distance
map to the contour. This approach was compared to the U-Net and to a point
based approach, achieving performance gains of up to 30% in terms of landmark
localisation ($<4.5mm$) and distance to the ground truth contour ($<3.5mm$).

###### Keywords:

Echocardiography Segmentation Distance map.

##  1 Introduction

Left ventricular (LV) contouring in echocardiography (echo) is a crucial step
towards deriving cardiac measures including strain, LV volume, ejection
fraction, cardiac output, and others. This task is typically carried out by
manual delineation of B-mode images by experts. Such manual work is time
consuming and prone to high inter- and intra- operator variability [2]. As a
result, automation of LV contouring has been investigated thoroughly over the
past decades. Since most of the aforementioned measures can be derived from
contours obtained in two-chamber and four-chamber views, the three-chamber
view (which is crucial for strain measures and accurate LV volumes) has been
mostly neglected. In this paper, we focus our evaluation on contouring of the
three-chamber view.

Most related work has aimed at solving the _LV segmentation_ problem by
computing the closed region of the input image that includes the LV blood
pool. The best performing methods to date use convolutional neural networks
(CNNs) to predict the binary mask representing the ventricular blood pool from
a frame or for every frame in a sequence.

Leclerc et al. [2] proposed multiple encoder-decoder architectures to predict
LV segmentation masks in the CAMUS dataset, which includes two-chamber and
four-chamber views (but not three-chamber), in a frame-by-frame fashion. They
found that U-Net [9] performed best for this task. Ouyang et al. [5] proposed
a two step process towards estimating the ejection fraction (EF). First, they
trained a segmentation model on end-diastole (ED) and end-systole (ES) frames,
and applied it to all frames in the cycle. Then, they combined the frame-wise
predicted segmentation with the original sequence to regress EF using a
spatio-temporal convolution model in the clip. Their spatio-temporal model
requires resampling each cycle to a fixed length of 32 frames, losing
potentially useful information about cardiac dynamics. An alternative way of
exploiting the temporal domain is with the use of recurrent blocks. For
example, Li et al. [3] proposed an encoder-decoder architecture were the
decoder was made of convolutional LSTM blocks, passing frame information to
the next frame, but without any skip connections between the encoder and the
decoder. Before them, Poudel et al. [8] incorporated a recurrent block into
the bottleneck of a U-Net, to leverage spatial consistency across slices in a
cardiac magnetic resonance stack. In this paper, we incorporate a recurrent
block into a U-Net inspired architecture but to leverage temporal consistency
across consecutive frames. Recently, Painchaud et al. [6] proposed a post-
processing method to take frame-wise segmentations and correct temporal
inconsistencies.

All previous methods aimed at improving the quantification of area-based (or
volume based) quantities, such as EF. This makes focusing on the prediction of
closed segmentation masks appropriate. However, in many other applications,
such as regional strain, cardiac motion, and shape analysis, predicting a
endocardial contour open at the mitral valve plane would be a more sensible
choice. Moreover, retrieving the location of anatomical landmarks along the
contour, such as the basal points and the apex, are often required to split
the ventricle into segments and to track specific anatomy throughout the
cardiac cycle. Previous efforts aiming at directly predicting anatomical
points along the contour are limited to the work by Porumb et al. [7], where
points along the contour of a two-chamber or a four-chamber view were
regressed using a modified U-Net architecture using differentiable spatial to
numerical transform (DSNT) layers [4]. The major drawbacks of this approach
are that it assumes that all the landmarks can be consistently annotated by
clinicians, when in fact only the 2 basal points and the apex are anatomically
meaningful and can be identified consistently, which is insufficient to
accurately define a contour; and that it uses one frame at a time, being
unable to leverage temporal consistency of the sequence.

In this paper, we propose a contouring method for 2D cardiac ultrasound
sequences. The novel contributions are: 1) design of U-Net based multi-decoder
for predicting the endocardial LV contour (instead of a segmentation mask) and
three consistent anatomical landmarks in it: the basal points and the apex, 2)
exploit frame interdependence by using a recurrent module and temporal
regularization, and 3) propose contour evaluation metrics based on distances
(rather than overlap measures) that are intuitive and clinically relevant.

##  2 Methods

We designed a neural network that takes as input a fixed-length sequence of
three-chamber ultrasound images. This network predicts a distance-to-contour
map and the location of $N_{p}$ points on the contour for each frame. The
contour curve is derived as the shortest geodesic path passing by the
predicted points, using the predicted distance map to define geodesic
distance. An overview of the method is illustrated in Fig. 1, and each step is
described in turn in this section.

Figure 1: Overview of the method. A dual-decoder U-Net processes the input
sequence predicting, for each frame, the location of 7 contour points (as
heatmaps, on the top decoder) and a distance-to-contour map (on the bottom
decoder). Using DSNT layers [4], point coordinates are calculated for each
heatmap. The contour is calculated as the spline passing by all 7 predicted
points through the shortest geodesic distance on the distance map.

###  2.1 Data preparation

Data belongs to an institutional database and include 525 scans from 460
patients that are part of EVAREST [10], a UK multi-site, multi-vendor clinical
trial, obtained using GE, Phillips or Siemens ultrasound systems. The dataset
contains echocardiograms from patients diagnosed with coronary artery disease,
COVID patients or healthy volunteers. The EVAREST trial specifically, contains
echocardiograms obtained from patients at rest or from patients subjected to
exercise or pharmacological stress. Heart rates vary between 40 and 280 and
the number of frames per scan is between 6 and 130. For each patient there was
at least one 3 chamber sequence stored as DICOM, containing one or more
cardiac cycles. The patients were split into $T=440$ for training and 85 for
validation. The data were anonymised and pre-processed so that most text and
annotation information outside the ultrasound region is removed.

The data was annotated by expert physiologists using in-house software, as
follows. For every clip, the operator identified the end-diastole (ED) and
end-systole (ES) frames in a cycle of their choice. Both frames were in turn
annotated by selecting 7 points, of which 3 corresponded to the two basal
points (BPs) and the apex (AP). The other 4 points were picked, two on each
side of the AP, along the contour and roughly equidistant between the AP and
the corresponding BP, as exemplified in Fig. 1. Operators can add additional
points if needed. As a result, only the AP and BPs correspond to anatomical
landmarks and the other points serve to define the contour. The ground truth
contour is computed as a smooth spline that passes by the selected points.

We converted the ground truth contour information into distance maps by
defining a pixel grid, of the same size as the input images, and computing at
each pixel the shortest distance, in pixels, to the contour (e.g. Fig. 1). The
ground truth labels $L_{i}=(\\{P\\}_{j},D)_{i}$ for each image $i$ are tuples
including 7 points $\\{P\\}_{j}$ (of which the first and last are BPs and the
mid point is the apex), and the distance-to-contour map $D$. $D$ is normalised
by the half-diagonal of the image.

During training, a sequence of length 15 frames was extracted from each clip.
The sequence was extracted at random (sampling uniformly) at each epoch, while
ensuring that at least 1 frame within was annotated. The data component of the
loss (c.f. Sec 2.2) was computed only on the annotated frames.

###  2.2 Prediction of the distance-to-contour map and contour points

We propose a dual encoder U-Net based architecture to jointly predict contour
points $\\{\tilde{P}\\}_{j}$ and distance-to-contour maps $\tilde{D}$. An
overview of the neural network architecture is illustrated in Fig. 2.

Figure 2: Architecture of the contour prediction network, on the left, and
detail of the encoder (Enc), bottleneck and decoder (Dec) blocks, on the
right. The architecture is based on a U-Net but with two parallel decoder
branches, the first one used to predict the distance-to-contour map and the
second one (with $NP=7$ output channels) to predict the location of $NP=7$
contour points. On the right, detail of the model blocks.

The input to the network is an image clip of $N=15$ frames and size of
$64\times 64$ pixels. The output of the network is a tensor of size $15\times
64\times 64\times 8$, where the 8 output channels correspond to one distance-
to-contour map and 7 heat-maps predicting the location of the 7 points. These
are named, respectively $\tilde{D}$ and $\\{\tilde{H}\\}_{j=0,...,6}$. The
predicted location of the points can be computed using a Differentiable
Spatial to Numerical Transform [4] layer (DSNT). This layer computes point
coordinates from a heat-map as a spatial average weighted by normalised
coordinates, as $\tilde{P}_{j}=DSTN(\tilde{H}_{j})$. The loss function to
train the model was composed of a composite data fit term, $L_{D}$, and a
regularization term, $L_{R}$. The data term for $N_{p}=7$ points is:

$\begin{array}[]{c}L_{D}=L_{m}+L_{p}\textnormal{, where}\\\
L_{m}=\frac{1}{T}\sum_{i}MAE(D_{i},\tilde{D}_{i})\\\
L_{p}=\frac{1}{T}\sum_{i}\frac{1}{N_{p}}\sum_{j}||P_{j,i},\tilde{P}_{j,i}||^{2}.\end{array}$
(1)

The regularization term penalises the temporal gradient of the predicted
sequence as follows:

$L_{R}=\frac{1}{T}\sum_{i}\left(|\nabla_{t}\tilde{D}_{i}|+\frac{1}{N_{p}}\sum_{j}|\nabla_{t}\tilde{H}_{j,i}|\right),$
(2)

where the temporal gradient operator is approximated by first order time
difference. The final loss is $L=L_{D}+\lambda L_{R}$, where the optimal value
of $\lambda=0.001$ was found by hyper-parameter search in the interval
$[0,10^{-4},10^{-3},10^{-2},10^{-1},0.5,1]$.

###  2.3 Contour curve from distance-to-contour and points

The contour is computed as a parametric smooth spline curve from the predicted
points and the predicted distance-to-contour map by finding the shortest
geodesic path from predicted basal point to basal point, passing by all
predicted points, and using the distance map as a cost for computing the
geodesic distance.

We solve this by using the Dijkstra algorithm [1] between every pair of
consecutive points in the contour. Finally, we sample the resulting parametric
curve to 21 points (the 2 BPs, the AP, and 9 points uniformly sampled at each
side of the AP).

##  3 Implementation and Experiments

###  3.1 Implementation details

The method was implemented in Python 3, using Tensorflow 2.7 for the deep
learning pipeline. Training was carried out using a NVidia GeForce RTX 2080 Ti
GPU. All models were trained using the Adam optimizer with a learning rate of
$10^{-4}$ and a batch size of 4 for 150 epochs, which was in all cases well
beyond convergence. The model with the lowest validation loss was used for
evaluation.

###  3.2 Experiments

The aim of our experiments was to assess the accuracy of the predicted
landmark points (using the Euclidean distance between predicted and ground
truth points), and the accuracy of the predicted contour as a whole (as the
distance between 18 points sampled between each basal point and the apex in
the predicted contour, to the ground truth contour). We report average and
standard deviation values using the point-based method in [7] as a baseline,
and the proposed method without (Proposed) and with (Proposed GRU) a Gated
Recurrent Unit -GRU module in the bottleneck.

Additionally, we provided experts with pairs of animations showing the
contours over the cardiac cycle, and asked them to pick the most temporally
consistent, while being blinded to the method used for each animation (which
was randomized between baseline, proposed method without recurrent module, and
whole proposed method). We then reported the frequency at which each method
was picked best of the two.

##  4 Results

The results on accuracy of predicted landmarks and contours are shown in Table
1. The proposed method provides best performance both in accuracy of predicted
landmarks (up to a 30% better than the baseline) and in terms of distance to
contour (up to 4.7% better than the baseline). Largest performance gains are
in the landmark points, and particularly at the basal points. Using a
recurrent bottleneck does match or improve the results further in all cases.

Table 1: Performance summary using two metrics: Point distance, average $\pm$ standard deviation of the Euclidean distance reported for the landmark points (basal point 1 and 2 -BP1 and BP2- and apex), and average values over those three landmarks; and distance to contour, distance from 21 uniformly sampled points in the predicted contour to the ground truth contour. The relative improvement in terms of average accuracy is presented in green. |  | PointReg | Proposed | Proposed GRU
---|---|---|---|---
Point distance | BP1 | $5.44\pm 1.50$ |  $4.20\pm 1.34$ +22.8% |  $4.22\pm 1.34$ +22.5%
Apex | $5.66\pm 1.20$ |  $4.88\pm 1.36$ +13.8% |  $4.80\pm 1.14$ +15.2%
BP2 | $6.36\pm 1.90$ | $5.29\pm 1.38\color[rgb]{0,0.6,0}{+16.8\%}$ |  $4.45\pm 1.21$ +30.0%
Average | $5.82\pm 1.53$ |  $4.79\pm 1.36$ +17.7% |  $4.49\pm 1.23$ +22.9%
Distance to contour | $3.64\pm 1.49$ |  $3.57\pm 1.45$ +1.9% |  $3.47\pm 1.28$ +4.7%

The qualitative results on temporal consistency are reported in Fig 3. The box
and whiskers plots show the compiled scores from 5 experts, indicating the
fraction of the assessed cases where each method was better or equal, in terms
of temporal consistency, than its alternative (indicated $>=$ in the figure),
or the fraction where it was better than the alternative (indicated $>$). The
proposed method, when using a recurrent bottleneck, performs better or equal
than any of the other two methods in over 50% of the cases, compared to the
baseline method which was better or equal to the alternative method in just
over 20% of the cases.

Figure 3: Relative temporal consistency of the three competing methods. From
left to right, cases where experts judged: 1) the proposed method equal or
better; 2) or just better; 3) the proposed method without the GRU bottleneck
equal or better; 4) or just better; 5) the baseline (PointReg) equal or
better; 6) or just better. In all cases, comparisons are performed against the
other two methods.

Last, we present 8 contour examples, using the our model with a GRU block, in
Fig. 4. Each of the 8 rows represents a frame from 8 different subjects. The
first column shows the input frame; Columns 2 to 8 show the predicted heatmaps
corresponding to the 7 points, where it can be seen that the three landmarks
(cols. 2 and 8: BPs and col 5: AP) are consistently predicted with much
greater accuracy than the non landmark points (p2, p3, p5, p6). The last 3
columns show, respectively, the shortest geodesic path passing by the 7 points
(pts), the smooth spline passing by the landmarks on the predicted distance to
contour map (d2cm), and the resulting contour overlaid onto the original frame
(result).

Figure 4: Example intermediate and final outputs for a frame from eight
subjects of different quality representative of the validation set. Over the
columns: the original B-mode frame, the point location heat-maps, and the
predicted contour (shortest geodesic path, fitted spline and result). Refer to
the text for further detail.

##  5 Discussion and Conclusions

We have presented a method to predict endocardial LV contours from B-mode
sequences acquired on the 3-chamber view. The proposed method computes the
contour from three predicted landmarks and a predicted distance-to-contour
map, as opposed to related work where a binary mask of the LV is predicted.
The contour based approach has the advantage of consistent retrieval of
anatomical landmarks (the basal points and the apex) that are normally needed
to compute clinically relevant metrics that involve splitting the LV into
segments, or having an accurate delimitation of the LV at the mitral valve
plane.

The proposed method leverages temporal information using a GRU module at the
bottleneck, which improves landmark detection and provides more time-
consistent contours. Specifically, our method outperforms the baseline point
regression method in [7], even without the GRU. Results shown in Fig. 4
suggest that a reason for this is the lack of consistency in annotations of
the non-landmark points, which in turn cannot be predicted accurately, and the
model seems unsure of where to locate these intermediate points along the
contour.

In this paper, we have also proposed to move from shape overlap metrics, such
as Dice and intersection-over-union, to more intuitive geometric measures,
specifically Euclidean distance to landmark points and average contour to
contour distance. The results reported by these metrics seem to correlate well
with the qualitative assessment of the resulting contours over time.

Our model achieves an average accuracy of up to 4.49mm for the landmarks and
3.47 mm for the contour itself. Leclerc et al. [2] found that inter-user
variability was up to $4.0\pm 2.0$ mm (restricting the analysis to images of
good and medium quality) which suggests that our results are of clinical use
quality. Future work will involve more experiments with clinical users to
confirm this.

## References

  * [1] Dijkstra, E.W.: A note on two problems in connexion with graphs. Numerische mathematik 1(1), 269–271 (1959)
  * [2] Leclerc, S., Smistad, E., Pedrosa, J., Østvik, A., Cervenansky, F., Espinosa, F., Espeland, T., Berg, E.A.R., Jodoin, P.M., Grenier, T., et al.: Deep learning for segmentation using an open large-scale dataset in 2D echocardiography. IEEE transactions on medical imaging 38(9), 2198–2210 (2019)
  * [3] Li, M., Zhang, W., Yang, G., Wang, C., Zhang, H., Liu, H., Zheng, W., Li, S.: Recurrent aggregation learning for multi-view echocardiographic sequences segmentation. In: International conference on medical image computing and computer-assisted intervention. pp. 678–686. Springer (2019)
  * [4] Nibali, A., He, Z., Morgan, S., Prendergast, L.: Numerical coordinate regression with convolutional neural networks. arXiv preprint arXiv:1801.07372 (2018)
  * [5] Ouyang, D., He, B., Ghorbani, A., Yuan, N., Ebinger, J., Langlotz, C.P., Heidenreich, P.A., Harrington, R.A., Liang, D.H., Ashley, E.A., et al.: Video-based ai for beat-to-beat assessment of cardiac function. Nature 580(7802), 252–256 (2020)
  * [6] Painchaud, N., Duchateau, N., Bernard, O., Jodoin, P.M.: Echocardiography segmentation with enforced temporal consistency. IEEE Transactions on Medical Imaging (2022)
  * [7] Porumb, M., Mumith, A., Gao, S., Parker, A., Beqiri, A., Sarwar, R., Upton, R., Leeson, P., Woodward, G.: Site-specific automated contouring model generalisibiliy enhancement. European Heart Journal-Cardiovascular Imaging 22(Supplement_1), jeaa356–430 (2021)
  * [8] Poudel, R.P., Lamata, P., Montana, G.: Recurrent fully convolutional neural networks for multi-slice mri cardiac segmentation. In: MICCAI RAMBO. vol. 10129, pp. 83–94. Springer Verlag (8 2016)
  * [9] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [10] Woodward, W., Dockerill, C., McCourt, A., Upton, R., O’Driscoll, J., Balkhausen, K., Chandrasekaran, B., Firoozan, S., Kardos, A., Wong, K., et al.: Real-world performance and accuracy of stress echocardiography: the EVAREST observational multi-centre study. European Heart Journal-Cardiovascular Imaging 23(5), 689–698 (2022)


# “Beam à la carte”: laser heater shaping for attosecond pulses in a
multiplexed x-ray free-electron laser

Siqi Li<EMAIL_ADDRESS>SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, USA  Zhen Zhang<EMAIL_ADDRESS>SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, USA  Shawn Alverson  SLAC
National Accelerator Laboratory, Menlo Park, CA 94025, USA  David Cesar  SLAC
National Accelerator Laboratory, Menlo Park, CA 94025, USA  Taran Driver  SLAC
National Accelerator Laboratory, Menlo Park, CA 94025, USA  Stanford PULSE
Institute, SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
Paris Franz  Applied Physics Department, Stanford University, Stanford, CA,
94305, USA  SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
Stanford PULSE Institute, SLAC National Accelerator Laboratory, Menlo Park, CA
94025, USA  Erik Isele  Applied Physics Department, Stanford University,
Stanford, CA, 94305, USA  SLAC National Accelerator Laboratory, Menlo Park, CA
94025, USA  Stanford PULSE Institute, SLAC National Accelerator Laboratory,
Menlo Park, CA 94025, USA  Joseph P. Duris  SLAC National Accelerator
Laboratory, Menlo Park, CA 94025, USA  Kirk Larsen  SLAC National Accelerator
Laboratory, Menlo Park, CA 94025, USA  Stanford PULSE Institute, SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, USA  Ming-Fu Lin  SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, USA  Razib Obaid  SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, USA  Jordan T O’Neal  Physics
Department, Stanford University, Stanford, CA, 94305, USA  SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, USA  Stanford PULSE Institute,
SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA  River Robles
Applied Physics Department, Stanford University, Stanford, CA, 94305, USA
SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA  Stanford
PULSE Institute, SLAC National Accelerator Laboratory, Menlo Park, CA 94025,
USA  Nick Sudar  SLAC National Accelerator Laboratory, Menlo Park, CA 94025,
USA  Zhaoheng Guo  Applied Physics Department, Stanford University, Stanford,
CA, 94305, USA  SLAC National Accelerator Laboratory, Menlo Park, CA 94025,
USA  Stanford PULSE Institute, SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, USA  Sharon Vetter  SLAC National Accelerator Laboratory,
Menlo Park, CA 94025, USA  Peter Walter  SLAC National Accelerator Laboratory,
Menlo Park, CA 94025, USA  Anna L. Wang  Applied Physics Department, Stanford
University, Stanford, CA, 94305, USA  SLAC National Accelerator Laboratory,
Menlo Park, CA 94025, USA  Stanford PULSE Institute, SLAC National Accelerator
Laboratory, Menlo Park, CA 94025, USA  Joseph Xu  Argonne National Laboratory,
Lemont, IL, USA  Sergio Carbajo  SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, USA  University of California at Los Angeles, Los Angeles, CA
90095, USA  James P. Cryan<EMAIL_ADDRESS>SLAC National Accelerator
Laboratory, Menlo Park, CA 94025, USA  Stanford PULSE Institute, SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, USA  Agostino Marinelli
<EMAIL_ADDRESS>SLAC National Accelerator Laboratory, Menlo Park,
CA 94025, USA  Stanford PULSE Institute, SLAC National Accelerator Laboratory,
Menlo Park, CA 94025, USA

###### Abstract

Electron beam shaping allows the control of the temporal properties of x-ray
free-electron laser pulses from femtosecond to attosecond timescales. Here we
demonstrate the use of a laser heater to shape electron bunches and enable the
generation of attosecond x-ray pulses. We demonstrate that this method can be
applied in a selective way, shaping a targeted subset of bunches while leaving
the remaining bunches unchanged. This experiment enables the delivery of
shaped x-ray pulses to multiple undulator beamlines, with pulse properties
tailored to specialized scientific applications.

††preprint: APS/123-QED

Free-electron lasers have emerged as the brightest sources of x-rays available
[1, 2, 3, 4, 5], enabling the study of ultrafast phenomena with atomic
resolution and femtosecond timing. Two noteworthy advances in x-ray free-
electron lasers (XFEL) technology include: the generation of high power
attosecond pulses [6, 7, 8, 9] and the demonstration of high repetition rate
operation [5]. Attosecond x-ray free-electron laser represent a groundbreaking
advancement in the field of ultrafast science, improving the peak brightness
of attosecond light sources by more than six orders of magnitude [6]. The
advent of high-repetition rate XFELs, based on superconducting accelerator
technology, represents a new significant leap forward, increasing the average
power of XFEL sources by several orders of magnitude [5]. The potential to
simultaneously support multiple users by efficiently distributing a high-
repetition rate electron beam across numerous undulator lines will play an
indispensable role in enabling sustained discovery science at XFEL facilities.
The extension of attosecond capabilities to this new generation of XFELs will
unlock the full potential of attosecond x-ray science by massively increasing
statistics and signal to noise ratio in nonlinear spectroscopy experiments.

The generation of attosecond XFEL pulses involves manipulations of the
electron beam phase space in order to generate a short high current spike,
with a duration of$\sim 1$ fs and a peak current on the order of $10$ kA [6,
7]. As we make progress in extending short-pulse capabilities to the next
generation of XFELS, it is important to be able to deliver unique pulses “on
demand” to a specific beamline while leaving the x-ray pulses delivered to
other beamlines unaffected. In this article we demonstrate such selective
shaping of electron bunches for the generation of attosecond pulses. We show
that by using a shaped laser pulse in the FEL laser heater we can selectively
shape the bunches transported to the soft x-ray beamline of the LCLS-II free-
electron laser, while leaving the operation of the hard x-ray beamline
unperturbed. This is a key step towards the next generation of XFELs, with
specialized undulator beamlines capable of delivering shaped pulses “on
demand”.

Laser heaters are commonly used in XFELs to increase the slice energy spread
of the electron beam in a controlled way. The increased energy spread
suppresses collective instabilities that degrade the electron beam quality and
decrease the performance of the XFEL [10]. This is accomplished by the
resonant interaction of the electrons with a laser pulse in a magnetic
undulator placed in the middle of a magnetic chicane. [11, 12] The energy
spread induced by the heater is proportional to the laser electric field, thus
varying the temporal shape of the laser pulse can generate a time-dependent
slice energy spread in the electron beam. Temporal shaping of the laser heater
was shown to be an effective method for selective beam spoiling for the
generation of femtosecond pulses, or for the generation of complex spectral
modulations in seeded XFELs [13, 14]. This temporal modulation of the slice
energy spread can also be used to modulate the beam current profile. The
perturbation of the current profile can be amplified by collective
instabilities, resulting in one or more high-current spikes that can be used
for the emission of attosecond pulses in an XFEL [15]. This method can also be
used to generate stable short pulses in cavity-based XFELs [16]. An important
feature of this laser-based shaping mechanism is that it can be arbitrarily
timed with respect to the electron bunch. Therefore one can selectively shape
a sub-set of the electron bunches produced by the linac, while leaving the
other bunches unchanged. In an XFEL system with multiple beamlines, this
selectivity can be used to generate a high-current spike in the bunches
directed to a desired beamline, therefore delivering short pulses to a
specific user while leaving the remaining XFEL beamlines unperturbed.

Figure 1:  Experimental layout of the laser heater shaping scheme at the LCLS-
II. Figure 2:  (a) Longitudinal phase space measured after injector by TCAV0
with different IR laser timing w.r.s. electron beam. Beam head lies to top.
(b) Measurements of FEL lasing suppression at hard x-ray and soft x-ray
undulators with laser heater shaping on all electron beams (120 Hz). (c)
Measurements of FEL lasing suppression with selective laser heater shaping for
soft x-ray line (the last 10 pulses for every 120 pulses per second). Figure
3:  Measurements of the longitudinal phase space of electron beam ($a$) w/o
and ($b$) with later heater shaping. The beam current profiles ($c$) and slice
energy spread (SES) profiles ($d$) are presented for comparisons. Figure 4:
Example pulse reconstruction using angular streaking. The left panel shows the
single shot reconstructed pulse profile, and the right panel shows the single
shot reconstructed spectrum, with the measured spectra in red. The shaded blue
curve for each shot indicates the solution to the reconstruction algorithm
with a different random seed.

##  I Experimental Results

This experimental demonstration was conducted using the LCLS-II normal-
conducting linac simultaneously feeding two undulator lines: the hard x-ray
and the soft x-ray undulators. Two indepedent Ti:Sa laser systems are used in
the laser heater, which we will term the heating laser and shaping laser. The
IR pulse from the heating laser is stretched to the duration of 20 ps to
provide a uniform laser field over the electron bunch duration. The beam-
shaping laser pulse is compressed to a shorter duration ($\sim$1.5 ps), and
overlapped with the pulse from the standard heating laser. The shaping pulse
is shorter than the electron bunch and creates a short low-density region in
the longitudinal phase-space. As the electron beam is accelerated and
compressed, this shaped phase space seeds the microbunching instability,
leading to the formation of a high-current spike [15]. It is worth emphasizing
that because we use two independent lasers it is possible to achieve
multiplexing by pulse picking. I.e. we can adjust the time-structure of the
shaping laser so it only interacts with a sub-set of the electron bunches
generated by the accelerator.

Figure 1 presents the schematic layout of the accelerator setup for this
experiment. The electron beam is generated in a high-brightness s-band
photoinjector. The laser heater is placed at the end of the injector at an
energy of approximately 135 MeV. The electron beam phase-space generated by
the injector can be diagnosed using an s-band transverse deflecting cavity
(TCAV0) placed in a dedicated diagnostic beamline at the end of the LCLS
photoinjector. After the laser heater, the electron beam is accelerated in the
main linac and compressed in a two-stage compression system (the two bunch
compressor systems are labeled as BC1 and BC2 in Fig. 1). Downstream of the
linac, the electron beam can be sent to either undulator (SXR or HXR) at the
desired beam repetition rate by using a fast kicker. The longitudinal phase-
space is diagnosed with two X-band transverse deflecting cavities (XTCAV)
placed at the end of both beamlines [17]. We measure the attosecond x-ray
pulses generated by the shaped bunches in the soft x-ray line using photo-
electron angular streaking.

The effect of the shaping laser on the longitudinal phase-space can be
observed in Fig. 2, which shows the longitudinal phase-space measured after
the photoinjector. We scan the relative timing of the shaping pulse and the
electron bunch and observe a short overheated region in the beam phase-space
moving from the head to the tail of the bunch. To demonstrate the ability to
selectively shape the bunch profile, electrons bunches are sent to both the
hard and soft x-ray undulators. The electron bunches are produced at a
repetition rate of $120$ Hz, and the first 110 bunches are sent to the hard
x-ray undulators with the rest sent to the soft x-ray undulators using a fast
kicker. Figures 2 (b) and (c) show the emitted pulse energy as a function of
shaping laser timing in both undulators. In the first data set, the shaping
laser is operating at 120 Hz and interacting with all the bunches. In the
second data set, the shaping laser interacts only with the last 10 electron
bunches per second feeding the SXR undulator, while leaving the other 110
electron bunches untouched. When the shaping pulse is set to interact with all
the electron bunches, the shaping mechanism affects the performance of both
undulators. However, when the shaping pulse is only interacting with the
bunches directed to the SXR line, the performance of the HXR line is not
affected. We note that the shaping method reduces the overall pulse energy.
This is because of the complexity of the resulting longitudinal phase space,
in which different parts of the bunch have different time-energy correlation,
requiring different undulator configurations for efficient lasing [18].

Figure 3 shows the longitudinal phase-space measured at the end of the HXR
undulator with the XTCAV. Panel (a) shows the phase-space with no shaping
laser, while panel (b) shows the phase-space with shaping laser on and
interacting with the HXR bunches. The current and slice energy spread profiles
with and without laser heater shaping are presented in Figs. 3 (c) and (d)
respectively. A high current spike appears at the core part of the electron
beam when the shaping laser is present. Note that the temporal resolution of
the measurement ($\sim 1fs$, root mean square) is insufficient to fully
resolve the temporal shape of the spike, leading to an underestimate of the
peak current. A distinct feature of the shaped electron beam phase space is
the presence of a large time-energy correlation within the high-current spike
induced by the longitudinal space charge force. This correlation is lost in
the measurement due to the limited time resolution, and appears as a large
uncorrelated energy spread in the center of the bunch.

##  II Attosecond x-ray FEL pulses

Table 1:  Beam and machine parameters. parameters | values (Fig 3) | (Fig. 4)
---|---|---
charge (pC) | 165 | 180
BC1 current (kA) | 0.195 | 0.195
BC2 current (kA) | 2 | 3.1
electron beam energy (GeV) | 5 | 6.5
central photon energy (eV) | N/A | 715 - 725

The high-current spike generated with heater shaping is employed for the
generation of attosecond pulses in the SXR undulator, using the methods
described in [15]. The details of the undulator and accelerator configuration
are shown in the supplementary information. We operated the SXR undulator with
shaped laser heater while delivering to the HXR line during user operation
with the regular laser heater. Critically, the FEL performance of the HXR line
was unaffected by the shaping experiment in the SXR line (See the
Supplementary Information). We observe single-spike x-ray pulses at a photon
energy of 720 eV. The entire data set shows a median pulse energy of 24
$\mu$J. The spectrum of the x-ray pulses is measured by a variable line
spectrometer (VLS) [19, 20] and demonstrates a median bandwidth of 7.5 eV full
width half maximum (FWHM).

We use an angular streaking diagnostic to measure the pulse profile and
duration of the pulses on a single-shot basis [6, 21, 22, 23]. This
measurement was conducted in the TMO hutch of the LCLS facility. The
measurement and pulse reconstruction were performed using the methods
described in [6, 21]. We reconstruct the temporal shape of the x-ray pulse for
a subset of the data based on the measured x-ray pulse energy. The shots used
for reconstruction have a median pulse energy of 57 $\mu$J and the
corresponding x-ray spectra have a median bandwidth of 9.3 eV FWHM. Figure 4
shows the result of 10 individual pulse reconstructions, together with the
measured spectrum. The reconstructed full width half max pulse duration has a
median of 368 as with a standard deviation of +/- 150 as. In this data set,
the heating laser pulse energy was kept at 10 $\mu$J and the shaping laser
pulse energy was kept at 46 $\mu$J. In a different data set where we
continuously scanned the pulse energy of the shaping laser from 18 $\mu$J to
712 $\mu$J, the reconstruction consistently showed a full width half max pulse
duration of 320 as with a standard deviation of 186 as.

In conclusion we have demonstrated the generation of attosecond x-ray pulses
using laser shaping of electron bunches. The pulses were characterized in the
time domain with an angular streaking diagnostic. The median measured pulse
duration was 368 as, with a median FWHM bandwdith of 9.3 eV. This method can
be used to selectively shape a subset of electron bunches, with application to
multiplexed XFEL facilities. Specifically, we applied laser heater shaping to
the soft x-ray line of the LCLS-II XFEL while leaving the bunches used in the
hard x-ray line unperturbed. This method paves the way to shaped pulse
generation on demand in high-repetition rate XFEL facilities serving several
users at a time.

###### Acknowledgements.

Use of the Linac Coherent Light Source (LCLS), SLAC National Accelerator
Laboratory, is supported by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences under Contract No. DE-AC02-76SF00515. A.M.,
D.C., P.F., Z.G. and Z.Z. acknowledge support from the Accelerator and
Detector Research Program of the Department of Energy, Basic Energy Sciences
division. Z.G., R.R. and P.F. also acknowledge support from Robert Siemann
Fellowship of Stanford University. The effort from T.D.D., J.W., M.F.K, T.W.,
and J.P.C. is supported by DOE, BES, Chemical Sciences, Geosciences, and
Biosciences Division (CSGB). J. X. acknowledges support from the Department of
Energy, Office of Science, under contract number DE-AC02-06CH11357.

## References

  * Emma  _et al._ [2010] P. Emma, R. Akre, J. Arthur, R. Bionta, C. Bostedt, J. Bozek, A. Brachmann, P. Bucksbaum, R. Coffee, F.-J. Decker, _et al._ , First lasing and operation of an ångstrom-wavelength free-electron laser, nature photonics 4, 641 (2010).
  * Ishikawa  _et al._ [2012] T. Ishikawa, H. Aoyagi, T. Asaka, Y. Asano, N. Azumi, T. Bizen, H. Ego, K. Fukami, T. Fukui, Y. Furukawa, _et al._ , A compact x-ray free-electron laser emitting in the sub-ångström region, nature photonics 6, 540 (2012).
  * Prat  _et al._ [2020] E. Prat, R. Abela, M. Aiba, A. Alarcon, J. Alex, Y. Arbelo, C. Arrell, V. Arsov, C. Bacellar, C. Beard, _et al._ , A compact and cost-effective hard x-ray free-electron laser driven by a high-brightness and low-energy electron beam, Nature Photonics 14, 748 (2020).
  * Kang  _et al._ [2017] H.-S. Kang, C.-K. Min, H. Heo, C. Kim, H. Yang, G. Kim, I. Nam, S. Y. Baek, H.-J. Choi, G. Mun, _et al._ , Hard x-ray free-electron laser with femtosecond-scale timing jitter, Nature Photonics 11, 708 (2017).
  * Decking  _et al._ [2020] W. Decking, S. Abeghyan, P. Abramian, A. Abramsky, A. Aguirre, C. Albrecht, P. Alou, M. Altarelli, P. Altmann, K. Amyan, _et al._ , A mhz-repetition-rate hard x-ray free-electron laser driven by a superconducting linear accelerator, Nature photonics 14, 391 (2020).
  * Duris  _et al._ [2020] J. Duris, S. Li, T. Driver, E. G. Champenois, J. P. MacArthur, A. A. Lutman, Z. Zhang, P. Rosenberger, J. W. Aldrich, R. Coffee, _et al._ , Tunable isolated attosecond x-ray pulses with gigawatt peak power from a free-electron laser, Nature Photonics 14, 30 (2020).
  * Zhang  _et al._ [2020] Z. Zhang, J. Duris, J. P. MacArthur, A. Zholents, Z. Huang, and A. Marinelli, Experimental demonstration of enhanced self-amplified spontaneous emission by photocathode temporal shaping and self-compression in a magnetic wiggler, New Journal of Physics 22, 083030 (2020).
  * Marinelli  _et al._ [2017] A. Marinelli, J. MacArthur, P. Emma, M. Guetg, C. Field, D. Kharakh, A. Lutman, Y. Ding, and Z. Huang, Experimental demonstration of a single-spike hard-x-ray free-electron laser starting from noise, Applied Physics Letters 111 (2017).
  * Prat  _et al._ [2023] E. Prat, A. Malyzhenkov, C. Arrell, P. Craievich, S. Reiche, T. Schietinger, and G. Wang, Coherent sub-femtosecond soft x-ray free-electron laser pulses with nonlinear compression, APL Photonics 8 (2023).
  * Saldin  _et al._ [2004] E. L. Saldin, E. A. Schneidmiller, and M. Yurkov, Longitudinal space charge-driven microbunching instability in the tesla test facility linac, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 528, 355 (2004).
  * Huang  _et al._ [2010] Z. Huang, A. Brachmann, F.-J. Decker, Y. Ding, D. Dowell, P. Emma, J. Frisch, S. Gilevich, G. Hays, P. Hering, _et al._ , Measurements of the linac coherent light source laser heater and its impact on the x-ray free-electron laser performance, Physical Review Special Topics-Accelerators and Beams 13, 020703 (2010).
  * Ratner  _et al._ [2015] D. Ratner, C. Behrens, Y. Ding, Z. Huang, A. Marinelli, T. Maxwell, and F. Zhou, Time-resolved imaging of the microbunching instability and energy spread at the linac coherent light source, Physical Review Special Topics-Accelerators and Beams 18, 030704 (2015).
  * Marinelli  _et al._ [2016] A. Marinelli, R. Coffee, S. Vetter, P. Hering, G. West, S. Gilevich, A. Lutman, S. Li, T. Maxwell, J. Galayda, _et al._ , Optical shaping of x-ray free-electron lasers, Physical review letters 116, 254801 (2016).
  * Roussel  _et al._ [2015] E. Roussel, E. Ferrari, E. Allaria, G. Penco, S. Di Mitri, M. Veronese, M. Danailov, D. Gauthier, and L. Giannessi, Multicolor high-gain free-electron laser driven by seeded microbunching instability, Physical review letters 115, 214801 (2015).
  * Cesar  _et al._ [2021] D. Cesar, A. Anakru, S. Carbajo, J. Duris, P. Franz, S. Li, N. Sudar, Z. Zhang, and A. Marinelli, Electron beam shaping via laser heater temporal shaping, Physical Review Accelerators and Beams 24, 110703 (2021).
  * Tang  _et al._ [2022] J. Tang, J. P. Duris, G. Marcus, and A. Marinelli, Electron beam shaping for actively outcoupling radiation from an x-ray regenerative amplifier free-electron laser, Physical Review Accelerators and Beams 25, 080701 (2022).
  * Behrens  _et al._ [2014] C. Behrens, F.-J. Decker, Y. Ding, V. A. Dolgashev, J. Frisch, Z. Huang, P. Krejcik, H. Loos, A. Lutman, T. Maxwell, _et al._ , Few-femtosecond time-resolved measurements of x-ray free-electron lasers, Nature communications 5, 3762 (2014).
  * Saldin  _et al._ [2006] E. L. Saldin, E. A. Schneidmiller, and M. V. Yurkov, Self-amplified spontaneous emission fel with energy-chirped electron beam and its application for generation of attosecond x-ray pulses, Phys. Rev. ST Accel. Beams 9, 050702 (2006).
  * Hettrick  _et al._ [1988] M. C. Hettrick, J. H. Underwood, P. J. Batson, and M. J. Eckart, Resolving power of 35,000 (5 ma) in the extreme ultraviolet employing a grazing incidence spectrometer, Applied optics 27, 200 (1988).
  * Obaid  _et al._ [2018] R. Obaid, C. Buth, G. L. Dakovski, R. Beerwerth, M. Holmes, J. Aldrich, M.-F. Lin, M. Minitti, T. Osipov, W. Schlotter, _et al._ , Lcls in—photon out: Fluorescence measurement of neon using soft x-rays, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 034003 (2018).
  * Li  _et al._ [2018] S. Li, Z. Guo, R. N. Coffee, K. Hegazy, Z. Huang, A. Natan, T. Osipov, D. Ray, A. Marinelli, and J. P. Cryan, Characterizing isolated attosecond pulses with angular streaking, Optics express 26, 4531 (2018).
  * Hartmann  _et al._ [2018] N. Hartmann, G. Hartmann, R. Heider, M. Wagner, M. Ilchen, J. Buck, A. Lindahl, C. Benko, J. Grünert, J. Krzywinski, _et al._ , Attosecond time–energy structure of x-ray free-electron laser pulses, Nature Photonics 12, 215 (2018).
  * Zhao  _et al._ [2022] X. Zhao, S. Li, T. Driver, V.-H. Hoang, A.-T. Le, J. P. Cryan, A. Marinelli, and C. Lin, Characterization of single-shot attosecond pulses with angular streaking photoelectron spectra, Physical Review A 105, 013111 (2022).


# The effect of the width of the incident pulse to the dielectric transition
layer in the scattering of an electromagnetic pulse - a quantum lattice
algorithm simulation

George Vahala  Department of Physics, William & Mary, Williamsburg, VA23185
Linda Vahala  Department of Electrical & Computer Engineering, Old Dominion
University, Norfolk, VA 12319  Abhay K. Ram  Plasma Science and Fusion Center,
MIT, Cambridge, MA 02139  Min Soe  Department of Mathematics and Physical
Sciences, Rogers State University, Claremore, OK 74017

$\bf{Abstract}$: The effect of the thickness of the dielectric boundary layer
that connects a material of refractive index $n_{1}$ to another of index
$n_{2}$ is considered for the propagation of an electromagnetic pulse. For
very thin boundary layer the scattering properties of the pulse mimics that
found from the Fresnel jump conditions for a plane wave - except that the
transmission to incident amplitudes are augmented by a factor of
$\sqrt{n_{2}/n_{1}}$. As the boundary layer becomes thicker one finds
deviations away from the Fresnel conditions and eventually one approaches WKB
propagation. However there is found a small but unusual dip in part of the
transmitted pulse that persists in time. The quantum lattice algorithm (QLA)
used recovers the Maxwell equations to second order in a small parameter – but
QLA still recovers Maxwell equations when this parameter is unity. The
expansion parameter is the speed of the pulse in medium $n_{1}$.

##  1 Introduction

###  1.1 Maxwell equations in a vacuum

There has been considerable interest in connecting the vacuum equations of
electromagnetism with a photonic wave function [1-6]. For example, Majorana
[5] postulated that the vacuum Maxwell equations could be rewritten in the
form of a Dirac equation for the ”wave function” $\mathbf{F}$ that is defined
explicitly in terms of the electric and magnetic fields $\mathbf{E}$ and
$\mathbf{B}$ by (under appropriate normalization)

$\mathbf{F}=\mathbf{E}-i\mathbf{B}.$ (1)

This function $\mathbf{F}$ had been introduced earlier [7] and is known as the
Riemann-Silberstein-Weber (RSW) vector. The vacuum Maxwell equations can then
be written in terms of this RSW function

$\frac{\partial\mathbf{F}}{\partial
t}=i\nabla\crossproduct\mathbf{F}\qquad,\qquad\nabla\cdot\mathbf{F}=0$ (2)

Using the correspondance principle

$\mathcal{E}\rightarrow i\frac{\partial}{\partial
t}\qquad,\qquad\mathbf{p}\rightarrow-i\nabla$ (3)

Eq. (2) takes a quasi-Dirac form

$(\mathcal{E}-\bm{\alpha}\cdot\mathbf{p})\mathbf{F}=0\qquad,\qquad\mathbf{p}\cdot\mathbf{F}=0$
(4)

where the $\bm{\alpha}$ are $3\crossproduct 3$ Hermitian matrices, whose
specific forms are not needed here. The first part of Eq. (4) corresponds to
the curl-parts of Maxwell equations while the second part of Eq. (4) to the
divergence equations.

Khan [8] introduced a 4-spinor $\Psi$ using certain combination of the RSW
vector

$\displaystyle{\Psi}=\begin{pmatrix}-F_{x}+iF_{y}\\\ F_{z}\\\ F_{z}\\\
F_{x}+iF_{y}\end{pmatrix},$ (5)

and showed that the two curl and the two divergence equations of Maxwell could
be written in a single 4-dimensional (4D) Dirac-like equation

$\frac{\partial\Psi}{\partial t}=-\mathbf{M}\cdot\nabla\Psi$ (6)

where the vector $4\crossproduct 4$ matrices $\mathbf{M}$ are given by the
tensor products of the Pauli spin matrices
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z}).$ with the $2\times 2$
identity matrix $\mathbf{I_{2}}$:

$\mathbf{M}=\bm{\sigma}\otimes\mathbf{I_{2}}.$ (7)

Indeed, the vacuum 4-spinor Maxwell Eq. (6), written out explicitly, takes the
form

$\frac{\partial}{\partial t}\begin{bmatrix}\psi_{0}\\\ \psi_{1}\\\ \psi_{2}\\\
\psi_{3}\\\ \end{bmatrix}=-\frac{\partial}{\partial
x}\begin{bmatrix}\psi_{2}\\\ \psi_{3}\\\ \psi_{0}\\\ \psi_{1}\\\
\end{bmatrix}+i\frac{\partial}{\partial y}\begin{bmatrix}\psi_{2}\\\
\psi_{3}\\\ -\psi_{0}\\\ -\psi_{1}\\\ \end{bmatrix}-\frac{\partial}{\partial
z}\begin{bmatrix}\psi_{0}\\\ \psi_{1}\\\ -\psi_{2}\\\ -\psi_{3}\\\
\end{bmatrix}$ (8)

while the 4-spinor Dirac equation for a massless free particle is

$\frac{\partial}{\partial t}\begin{bmatrix}\psi_{0}\\\ \psi_{1}\\\ \psi_{2}\\\
\psi_{3}\\\ \end{bmatrix}=\frac{\partial}{\partial
x}\begin{bmatrix}\psi_{3}\\\ \psi_{2}\\\ \psi_{1}\\\ \psi_{0}\\\
\end{bmatrix}+i\frac{\partial}{\partial y}\begin{bmatrix}-\psi_{3}\\\
\psi_{2}\\\ -\psi_{1}\\\ \psi_{0}\\\ \end{bmatrix}+\frac{\partial}{\partial
z}\begin{bmatrix}\psi_{2}\\\ -\psi_{3}\\\ \psi_{0}\\\ -\psi_{1}\\\
\end{bmatrix}.$ (9)

Now Yepez [9, 10] has introduced a sequence of interleaved non-commuting
unitary collision and streaming operators that perturbatively recover the
Dirac equation to second order accuracy. We [11-16] then generalized these
non-commuting unitary collision and streaming operators in order to generate a
qubit lattice algorithm (QLA) for the solution of the vacuum Maxwell Eqs. (8).
Interestingly, these unitary QLA’s can, in principle, be readily encoded onto
error-correcting qubit quantum computers to yield a quantum algorithm solving
the initial value problem for Maxwell equations.

###  1.2 Maxwell equations in scalar dielectric media

Khan[8] has also extended his matrix representation of Maxwell equations to
include propagation in a scalar dielectric medium $\varepsilon(\mathbf{x})$,
with refractive index $n(\mathbf{x})=\sqrt{\varepsilon(\mathbf{x})}$. The
medium inhomogeniety will now couple the two possible electromagnetic pulse
polarizations so one must introduce two 4-spinors

$\displaystyle{\Psi^{\pm}}=\begin{pmatrix}-F^{\pm}_{x}\pm iF^{\pm}_{y}\\\
F^{\pm}_{z}\\\ F^{\pm}_{z}\\\ F^{\pm}_{x}\pm iF^{\pm}_{y}\end{pmatrix},$ (10)

with generalized RSW vectors

$\mathbf{F^{\pm}}=\sqrt{\epsilon}\mathbf{E}\pm
i\frac{\mathbf{B}}{\sqrt{\mu_{0}}}.$ (11)

($\mu_{0}$ is the vacuum magnetic permeability). The Maxwell equations, with
no free sources, are

$\displaystyle\nabla\bm{.}\mathbf{D}=0,\qquad\nabla\bm{.}\mathbf{B}=0$
$\nabla\times\mathbf{E}=-\frac{\partial\mathbf{B}}{\partial
t},\qquad\nabla\times\mathbf{H}=\frac{\partial\mathbf{D}}{\partial t}$ (12)

with $\mathbf{D}=\varepsilon(\mathbf{x})\mathbf{E}$ and
$\mathbf{B}=\mu_{0}\mathbf{H}$. The evolution equations for the coupled
4-spinor RSW vectors [8] are

$\displaystyle\frac{\partial}{\partial t}\begin{pmatrix}\Psi^{+}\\\
\Psi^{-}\end{pmatrix}=$ $\displaystyle-
v_{ph}\begin{pmatrix}\mathbf{M}\bm{.}\nabla-\bm{\Sigma}\bm{.}\displaystyle{\frac{\nabla\epsilon}{4\epsilon}}&+iM_{z}\bm{\Sigma}\bm{.}\displaystyle{\frac{\nabla\epsilon}{2\epsilon}}\alpha_{y}\\\
+iM_{z}\bm{\Sigma}^{*}\bm{.}\displaystyle{\frac{\nabla\epsilon}{2\epsilon}}\alpha_{y}&\mathbf{M}^{*}\bm{.}\nabla-\bm{\Sigma}^{*}\bm{.}\displaystyle{\frac{\nabla\epsilon}{4\epsilon}}\end{pmatrix}\begin{pmatrix}\Psi^{+}\\\
\Psi^{-}\end{pmatrix}.$ (13)

where $v_{ph}=(\epsilon\mu_{0})^{-1/2}$ is the pulse phase velocity, with the
$4\times 4$ matrices

$\bm{\alpha}=\begin{pmatrix}0&\bm{\sigma}\\\ \bm{\sigma}&0\\\ \end{pmatrix},\
\ \ \bm{\Sigma}=\begin{pmatrix}\bm{\sigma}&0\\\ 0&\bm{\sigma}\\\
\end{pmatrix}.$ (14)

Now the time evolution of the 8-spinor $(\Psi^{+},\Psi^{-})$ is no longer
fully Hermitian [8] - but several of the matrices (which depend on
$\nabla\epsilon$) are anti-Hermitian. In particular, for 1D pulse propagation
in the z-direction, the 8-spinor system, Eqs. (13), reduces to

$\displaystyle\frac{\partial}{\partial t}\begin{bmatrix}\psi_{0}\\\
\psi_{1}\\\ \psi_{2}\\\ \psi_{3}\\\
\end{bmatrix}=-\frac{1}{n(z)}\frac{\partial}{\partial
z}\begin{bmatrix}\psi_{0}\\\ \psi_{1}\\\ -\psi_{2}\\\ -\psi_{3}\\\
\end{bmatrix}+\frac{n^{\prime}(z)}{2n^{2}(z)}\begin{bmatrix}\psi_{0}-\psi_{7}\\\
-\psi_{1}-\psi_{6}\\\ \psi_{2}+\psi_{5}\\\ -\psi_{3}+\psi_{4}\\\
\end{bmatrix}$ $\frac{\partial}{\partial t}\begin{bmatrix}\psi_{4}\\\
\psi_{5}\\\ \psi_{6}\\\ \psi_{7}\\\
\end{bmatrix}=-\frac{1}{n(z)}\frac{\partial}{\partial
z}\begin{bmatrix}\psi_{4}\\\ \psi_{5}\\\ -\psi_{6}\\\ -\psi_{7}\\\
\end{bmatrix}+\frac{n^{\prime}(z)}{2n^{2}(z)}\begin{bmatrix}\psi_{4}-\psi_{3}\\\
-\psi_{5}-\psi_{2}\\\ \psi_{6}+\psi_{1}\\\ -\psi_{7}+\psi_{0}\\\
\end{bmatrix}.$ (15)

To develop a QLA for the 1D z-propagation of a pulse in a dielectric medium we
need to determine the specific collision and streaming operators required so
that Eqs. (15) are recovered perturbatively to 2nd order accuracy. The
collision operator must couple two qubits at a given lattice site so that we
have quantum entanglement, and then that entanglement is spread throughout the
lattice by the streaming. If we look at the time evolution of any of the
8-spinor $\psi_{i}$ in Eqs. (15), one sees that it is coupled to the
z-derivative $\partial\psi_{i}/\partial z$ of the same spinor component
$\psi_{i}$. Hence for a QLA for z-propagation we are forced into a $16$-qubit
representation. This coupling of the same spinor component in
$\partial/\partial t$ with $\partial/\partial z$ can be traced to the diagonal
property of the Pauli $\sigma_{z}$ matrix. This does not occur in the QLA for
x- or y-propagation of the 8-spinor since different spinor components of
$\partial\psi_{i}/\partial t$ couple with $\partial\psi_{j}/\partial x$ and
$\partial\psi_{j}/\partial y$ with $j\neq i$ [11-14]. Thus for x- and
y-propagation we can directly work with an 8-qubit set that is just the
8-spinor set of Khan’s Maxwell representation. To develop QLA for 2D and 3D
simulations, the interleaved sequence of collide-streams in each of the
Cartesian directions will act consecutively on the chosen qubit basis. It is
thus more convenient to work with a 16-qubit basis by simply generalizing the
8-spinor set to the 16-qubit space.

###  1.3 QLA for 1D Propagation in the z-direction

We connect the 16-qubit field $Q=(\overline{q}_{0}....\overline{q}_{15})^{T}$
to the 8-spinor field $(\psi_{0}...\psi_{7})^{T}$ by

$\displaystyle\psi_{0}=\overline{q}_{0}+\overline{q}_{2},\psi_{1}=\overline{q}_{1}+\overline{q}_{3},\psi_{2}=\overline{q}_{4}+\overline{q}_{6},\psi_{3}=\overline{q}_{5}+\overline{q}_{7},$
$\psi_{4}=\overline{q}_{8}+\overline{q}_{10},\psi_{5}=\overline{q}_{9}+\overline{q}_{11},\psi_{6}=\overline{q}_{12}+\overline{q}_{14},\psi_{7}=\overline{q}_{13}+\overline{q}_{15}.$
(16)

The QLA on the 16-qubit basis consists of an interleaved sequence of collision
($C$) and streaming ($S$) operators

$\displaystyle U$
$\displaystyle=\overline{S}_{-}C\overline{S}_{+}C^{\dagger}\bm{.}\>S_{+}CS_{-}C^{\dagger},$
(17) $\displaystyle\overline{U}$
$\displaystyle=S_{+}C^{\dagger}S_{-}C\bm{.}\>\overline{S}_{-}C^{\dagger}\overline{S}_{+}C$

where $C^{\dagger}$ is the adjoint of $C$. The operator $S_{+}$ streams the
8-qubit subset
$[\overline{q}_{0},\overline{q}_{1},\overline{q}_{4},\overline{q}_{5},\overline{q}_{8},\overline{q}_{9},\overline{q}_{12},\overline{q}_{13}]$
one lattice unit in the positive z-direction. Notice that these qubits are
each the first qubits on the right hand sides of Eqs. (16). $S_{-}$ streams
these 8 qubits one lattice unit in the negative z-direction. Similarly for the
streaming operator $\overline{S}$ for the other 8 qubit which are the second
qubits on the right hand sides of Eqs. (16).

The spatial and time derivatives on the 8-spinor field are determined from
$\overline{U}.UQ(t)$. In particular, an appropriate $16\times 16$ unitary
collision operator $C$ couples the qubits
[$\overline{q}_{0}-\overline{q}_{2}$, $\overline{q}_{1}-\overline{q}_{3}$],
[$\overline{q}_{4}-\overline{q}_{6}$, $\overline{q}_{5}-\overline{q}_{7}$],
[$\overline{q}_{8}-\overline{q}_{10}$, $\overline{q}_{9}-\overline{q}_{11}$],
[$\overline{q}_{12}-\overline{q}_{14}$, $\overline{q}_{13}-\overline{q}_{15}$]
and has $4\crossproduct 4$ block diagonal structure

$C(\theta)=\begin{bmatrix}C_{4}(\theta)&0&0&0\\\ 0&C_{4}(\theta)^{T}&0&0\\\
0&0&C_{4}(\theta)&0\\\ 0&0&0&C_{4}(\theta)^{T}\\\ \end{bmatrix},$ (18)

with the $4\crossproduct 4$ submatrix

$C_{4}(\theta)=\begin{bmatrix}\cos\theta&0&\sin\theta&0\\\
0&\cos\theta&0&\sin\theta\\\ -\sin\theta&0&\cos\theta&0\\\
0&-\sin\theta&0&\cos\theta\\\ \end{bmatrix},$ (19)

and collision angle $\theta$

$\theta=\frac{\epsilon}{4n(z)}$ (20)

for some small perturbation parameter $\epsilon$. This sequence of collide-
stream operators will generate, in the continuum limit, the required first
term $\partial/\partial z$ in Eqs. (15) to $O(\epsilon^{2})$.

To recover the two $n^{\prime}(z)$ terms in the 8-spinor Eqs. (15) one
introduces a particular potential collision operator for each term. Since this
first $n^{\prime}(z)$ term has exactly the same couplings to the spinor-
components as the $\partial/\partial t$ and $\partial/\partial z$ terms, the
appropriate QLA potential collision operator will again have this
$4\crossproduct 4$ diagonal block structure

$\overline{P}_{1Z}(\gamma)=\begin{bmatrix}P_{4}(\gamma)&0&0&0\\\
0&P_{4}(\gamma)&0&0\\\ 0&0&P_{4}(\gamma)&0\\\ 0&0&0&P_{4}(\gamma)\\\
\end{bmatrix},$ (21)

with

$P_{4}(\gamma)=\begin{bmatrix}\cos\gamma&0&-\sin\gamma&0\\\
0&\cos\gamma&0&-\sin\gamma\\\ -\sin\gamma&0&\cos\gamma&0\\\
0&-\sin\gamma&0&\cos\gamma\\\ \end{bmatrix},$ (22)

for some collisional angle $\gamma$. The second potential collision operator
has diagonal structure in its two $8\crossproduct 8$ matrices

$\overline{P}_{2Z}(\gamma)=\begin{bmatrix}P_{81}(\gamma)&P_{82}(\gamma)\\\
P_{82}(\gamma)&P_{81}(\gamma)\\\ \end{bmatrix},$ (23)

where

$P_{81}(\gamma)=\begin{bmatrix}\cos\gamma&0&0&0&0&0&0&0\\\
0&\cos\gamma&0&0&0&0&0&0\\\ 0&0&\cos\gamma&0&0&0&0&0\\\
0&0&0&\cos\gamma&0&0&0&0\\\ 0&0&0&0&\cos\gamma&0&0&0\\\
0&0&0&0&0&\cos\gamma&0&0\\\ 0&0&0&0&0&0&\cos\gamma&0\\\
0&0&0&0&0&0&0&\cos\gamma\\\ \end{bmatrix}.$ (24)

and

$P_{82}(\gamma)=\begin{bmatrix}0&0&0&0&0&0&0&-\sin\gamma\\\
0&0&0&0&0&0&-\sin\gamma&0\\\ 0&0&0&0&0&-\sin\gamma&0&0\\\
0&0&0&0&-\sin\gamma&0&0&0\\\ 0&0&0&\sin\gamma&0&0&0&0\\\
0&0&\sin\gamma&0&0&0&0&0\\\ 0&\sin\gamma&0&0&0&0&0&0\\\
\sin\gamma&0&0&0&0&0&0&0\\\ \end{bmatrix}.$ (25)

To recover the 8-spinor Eqs. (15) to $O(\epsilon^{2})$, one chooses

$\gamma=\varepsilon^{2}\frac{n^{\prime}(z)}{2n(z)^{2}}.$ (26)

The final QLA that is used for the simulations is

$Q(t+\delta
t)=\overline{P}_{2Z}(\gamma).\overline{P}_{1Z}(\gamma).\overline{U}.U.Q(t)$
(27)

###  1.4 Initial conditions for the 16 qubit QLA

We are given the initial electromagnetic field components of the pulse. From
the RSW transformation, Eq. (10), we then have the initial 8-spinor
$(\psi_{0},....\psi_{7})$. The initial condition for the 16-qubit $Q(t=0)$ is
then chosen:

$\displaystyle\overline{q}_{0}=\overline{q}_{2}=\psi_{0}/2,\overline{q}_{1}=\overline{q}_{3}=\psi_{1}/2,\overline{q}_{4}=\overline{q}_{6}=\psi_{2}/2,\overline{q}_{5}=\overline{q}_{7}=\psi_{3}/2,$
$\overline{q}_{8}=\overline{q}_{10}=\psi_{4}/2,\overline{q}_{9}=\overline{q}_{11}=\psi_{5}/2,\overline{q}_{12}=\overline{q}_{14}=\psi_{6}/2,\overline{q}_{13}=\overline{q}_{15}=\psi_{7}/2.$
(28)

##  2 Effect of Dielectric Boundary Layer Thickness on Reflection and
Transmission of a 1D Pulse

Our earlier QLA simulations [11, 12] of an electromagnetic pulse propagating
from one dielectric (with refractive index $n_{1}$) to another (with
refractive index $n_{2}$) at normal incidence assumed the boundary layer
thickness connecting the two media, $\Delta_{BL}$, was very much less than the
width of the incident pulse $\Delta_{pulse}$: $\Delta_{BL}<<\Delta_{pulse}$.
In this case we found that the ratio of the reflected to initial field
amplitudes were just those derived from the Fresnel equations for a boundary
value problem of a plane wave normally incident onto a dielectric
discontinuity [17]. However, for pulse propagation as an initial value
problem, the transmitted to incident field amplitude is augmented from the
plane wave Fresnel results by a factor $\sqrt{n_{2}/n_{1}}$ :

$\frac{E_{refl}}{E_{inc}}=\frac{n_{1}-n_{2}}{n_{1}+n_{2}}\qquad,\quad\frac{E_{trans}}{E_{inc}}=\frac{2n_{1}}{n_{1}+n_{2}}\sqrt{\frac{n_{2}}{n_{1}}}.$
(29)

As the QLA simulations repeated stamped out Eq. (29) for many different
simulations with various choices of $n_{1}$ and $n_{2}$ and different pulse
geometries, we [12] developed a theory for Gaussian pulses that validates Eq.
(29). A back-of-the-envelope argument also shows this: now the reflected pulse
preserves all the reflected Fresnel plane wave characteristics. Thus they must
have the energy. Since the total energy of the scattering is conserved, then
we must have the same conservation of transmitted energy. However, the speed
of propagation of the transmitted pulse and its width are reduced by a factor
of $n_{2}/n_{1}$. Hence the amplitude of the transmitted pulse must increase
by a factor of $\sqrt{n_{2}/n_{1}}$.

###  2.1 Case 1 : $\Delta_{BL}<<\Delta_{pulse}$ with $\epsilon=0.3$

We consider the case of $n_{1}=1$ (for $z<4990$) and $n_{2}=2$ (for $z>5010$),
with the boundary layer to pulse width (see Fig. 1) of
$\Delta_{BL}/\Delta_{pulse}=0.1$. The vertical dashed line in Fig. 1 is the
mid point of the boundary layer.

(a) dielectric boundary layer (b) initial $E_{x},B_{y}$ profile

Figure 1: (a) A dielectric boundary layer connecting medium $n_{1}=1$, for
$z<4990$ to medium $n_{2}=2$ for $z>5010$. The vertical dashed line indicates
the midpoint of the boundary layer, (b) The initial vacuum fields (in our
units) $E_{x}=B_{y}$ with $\Delta_{pulse}=200$. Initially $B_{y}=n_{1}E_{x}$
and the two field profiles overlay each other. $E_{x}$ \- blue, $B_{y}$ \-
red.

The perturbation parameter $\epsilon$, introduced into the QLA collision
angles Eq. (20) and (26), is chosen to be $\epsilon=0.3$. In the QLA
simulation units, the time increment $\delta t=1$, Eq. (27).

As the pulse propagates into the boundary layer, part of it will be
transmitted and part of it will be reflected. Away from the boundary layer,
asynptotically, the transmitted amplitudes have $B_{y}=(n_{2}/n_{1})E_{x}$,
while the reflected amplitudes have $B_{y}=-E_{x}>0$. (Fig. 2). From QLA
simulations, we find

$\frac{E_{refl}}{E_{inc}}=-0.32\qquad,\quad\frac{E_{trans}}{E_{inc}}=0.94\qquad.$
(30)

This is in excellent agreement with Eq. (29), while the Fresnel plane wave
boundary conditions would yield $E_{trans}/E_{inc}=0.67$.

(a) reflected and transmitted fields (b) blowup of the profiles

Figure 2: (a) The transmitted and reflected fields with
$\Delta_{BL}<<\Delta_{pulse}$ after 5000 time iterations. (i.e., at t = 5k).
(b) A blowup showing the Gaussian nature of the reflected and transmitted
profiles. Note there is no disturbance between the reflected and transmitted
profiles. $E_{x}$ \- blue, $B_{y}$ \- red.

Note that the width of the reflected pulse is equal to that of the incident
pulse, but twice that of the transmitted pulse
$\Delta_{refl}=2\Delta_{trans}$, and the transmitted pulse travels at half the
speed of the reflected pulse.

###  2.2 Case 2 : $\Delta_{BL}\approx\Delta_{pulse}$ with $\epsilon=0.3$

When the boundary layer thickness is on the order of the pulse width, Fig. 3a,
both the transmitted and reflected field amplitudes are significantly
affected. In Fig. 3b, at time $t=4k$ much of the pulse is within this boundary
layer $4900<z<5100$. There is developing characteristics of the $n_{2}=2$
region as the magnetic field amplitude is basically twice that of the electric
field amplitude for $z>5000$ while part of the pulse for $z<5000$ is showing
the characteristics of the $n_{1}=1$ region with a significant range having
$E_{x}<0$.

(a) thicker dielectric boundary layer (b) pulse overlapping the boundary layer

Figure 3: (a) The dielectric boundary layer now extends $4900<z<5100$, with
$\Delta_{BL}\approx\Delta_{pulse}$ (b) At $t=4k$, the pulse is overlapping the
boundary layer. The basic 2:1 ratio of the transmitted magnetic to transmitted
electric field is being established for $z>5000$. $E_{x}$ \- blue, $B_{y}$ \-
red. In the reflected pulse, a significant portion already has $E_{x}<0$ since
$n_{1}<n_{2}$.

By $t=5k$, the reflected and transmitted pulses are approaching their quasi-
asymptotic state. The reflected pulse has a lower amplitude and greater width
than that predicted by the Fresnel plane wave conditions, Fig. 4a. Of some
interest is the trailing edge of the transmitted pulse which exhibits a slight
dip with both $B_{y}<0$ and $E_{x}<0$, but with $B_{y}/E_{x}\approx 2$, Fig.
4b.

(a) pulse at $t=5k$ (b) blowup of pulse at $t=5k$

Figure 4: (a) The transmitted and reflected fields when
$\Delta_{BL}\approx\Delta_{pulse}$. The reflected amplitude is almost a factor
of 3 lower than the reflected Fresnel plane wave solution and the reflected
pulse is significantly broader. (b) A blowup showing asymmetry in both the
reflected and transmitted profiles. There is a small region in $z>5000$ in
which the transmitted pulse exhibits the unusual behavior of $B_{z}<0$ and
$E_{z}<0$. $E_{x}$ \- blue, $B_{y}$ \- red.

In Fig. 5 we see the asymptotic profiles for the reflected and transmitted
pulses. $E_{refl}/E_{inc}=-0.12$ while the Fresnel plane wave solution yields
$-0.33$.

(a) pulse at time $t=6k$ (b) Blowup of the pulse at $t=6k$.

Figure 5: (a) The transmitted and reflected fields with
$\Delta_{BL}<<\Delta_{pulse}$. $E_{x}$ \- blue, $B_{y}$ \- red. (b) A blowup
showing the asymmetric nature of the reflected and transmitted profiles. At
the back end of the transmitted pulse one finds a small region in which both
$B_{y}<0$ and $E_{x}<0$ and with $B_{y}\approx 2E_{x}$.

Around $z\approx 5210$ one finds a section of the transmitted pulse that has
$B_{y}<0,E_{x}<0$ with the magnetic to electric field peaks being in the 2:1
ratio, as might be expected in the $n_{2}=2$ dielectric region. This feature
is asymptotically stable and will continue traveling as a part of the
transmitted pulse. It is not present when $\Delta_{BL}<<\Delta_{pulse}$, Fig.
2b.

###  2.3 Case 3: $\Delta_{BL}>>\Delta_{pulse}$ with $\epsilon=0.3$

We now consider the case when the boundary layer thickness $\Delta_{BL}\approx
2000$, and $\Delta_{pulse}\approx 200$, Fig 6a.

(a) refractive index with $\Delta_{BL}>>\Delta_{pulse}$. (b) Pulse at $t=2k$

Figure 6: (a) The refractive index for the case when $\Delta_{BL}\approx 2000$
and $\Delta_{pulse}\approx 200$. (b) By $2k$ iterations the pulse has moved
into a slightly non-vacuum dielectric state, as seen by the amplitude peak in
$B_{y}$ ebing greater than that in $E_{x}$. $E_{x}$ \- blue, $B_{y}$ \- red.

As the pulse propagates, Fig. 6b-8b, the peak amplitude ratio of $B_{y}/E_{x}$
basically scales as the value of the local refractive index at that point. In
essence, the time evolution of the pulse is similar to WKB . There is a very
low order reflected pulse - but its peaks are 2 orders of magnitude less than
the transmitted pulse (Fig. 8b).

(a) slowly varying refractive index (b) Pulse at $t=3k$

Figure 7: (a) The refractive index for the case when
$\Delta_{BL}>>\Delta_{pulse}$. (b) By $2k$ iterations the pulse has moved into
a slightly non-vacuum state, as seen by the amplitude peak in $B_{y}$ ebbing
greater than that in $E_{x}$. $E_{x}$ \- blue, $B_{y}$ \- red.

(a) slowly varying refractive index (b) Pulse at $t=3k$

Figure 8: (a) The refractive index for the case when
$\Delta_{BL}>>\Delta_{pulse}$. (b) By $2k$ iterations the pulse has moved into
a slightly non-vacuum state, as seen by the amplitude peak in $B_{y}$ ebing
greater than that in $E_{x}$. $E_{x}$ \- blue, $B_{y}$ \- red.

###  2.4 Case 4 : $\Delta_{BL}\approx\Delta_{pulse}$ with $\epsilon=1.0$

QLA is a perturbation theory, based on the parameter $\epsilon$ introduced
into the various collision operators, Eqs. (20) and (26). Moreover, QLA
recovers the Maxwell equations only to errors $O(\epsilon^{2})$. In our
earlier QLA for solitons [18] the perturbation parameter was related to the
amplitude of the wave function $\psi_{NLS}$.. Since the evolution equations
for solitons/quantum vortices involve a $|\psi_{NLS}|^{2}$ nonlinearity, we
[18] found that QLA simulations strongly deviated away from the exact
nonlinear soliton dynamics if $\epsilon$ was chosen too high - typically
$\epsilon\approx 0.45$. These deviations took the form of background turbulent
noise along with distorted to disintegrating soliton shapes.

Here we are developing a QLA for the linear Maxwell equations and find that
the basic Maxwell equations are still being modeled even up to $\epsilon=1.0$.
The essential physics is retained (Fig. 9) but there is a slight deviation in
the pulse propagation speed. Indeed $\epsilon$ is nothing but the pulse speed
in medium 1. For $\epsilon=0.15$, in 20k time steps the pulse peak location is
delayed by 7 lattice units. For $\epsilon=0.3$ , by $t=10k$, the pulse peak is
delayed by 16 lattice units. This delay increases for increasing $\epsilon$:
for $\epsilon=0.6$, the peak is off by 45 lattice units, and by $125$ lattice
units for $\epsilon=1.0$.

(a) $\epsilon=0.3,t=20k$ (b) $\epsilon=1.0,t=6k$

Figure 9: (a) A blowup of the reflection-transmission pulses for the case
$\Delta_{BL}\approx\Delta_{pulse}$ at t = 20k and with $\epsilon=0.3$ (b) The
corresponding pulses but now for $\epsilon=1.0$ and at $t=6K$. For this higher
$\epsilon$ there is a slight decrease in the QLA pulse propagation speed.
$E_{x}$ \- blue, $B_{y}$ \- red.

##  3 Summary and Conclusions

In this paper we have examined the effect of the boundary layer thickness on
the 1D normal propagation of a pulse in a dielectric medium. For very sharp
boundary layers connecting the two dielectric regions, our initial value QLA
simulations give reflected and transmitted electromagnetic fields in agreement
with the boundary value Fresnel plane wave conditions - except that the ratio
of the transmitted to incident field amplitudes is augmented by a factor
$\sqrt{(}n_{2}/n_{1})$. As the thickness of the boundary layer, $\Delta_{BL}$,
increases the reflected pulse becomes more and more modified. It becomes
smaller in amplitude as well as increasing in width till we approach the WKB-
like solution when $\Delta_{BL}>>\Delta_{pause}$. A somewhat unsuspected
feature was seen in the trailing edge of the transmitted pulse. We noticed a
region in the $n_{2}$ dielectric in which both $B_{y}<0$ and $E_{x}<0$. This
”dip” is stable and propagates away undistorted from the boundary layer along
with the rest of the transmitted pulse. We have also performed a QLA run in
which $n_{1}=2>n_{2}=1$. In region $n_{1}$ we initially have $B_{y}=2E_{x}$.
The transmitted pulse now has $B_{y}=E_{x}$, but again with a ”dip in the
transmitted fields with $B_{y}=E_{x}<0.$. The Poynting flux has been
determined for each of the QLA runs, and we find energy is conserved with
normalized variations on the order of $1.35\times 10^{-3}$.

Our interest in the effects of the boundary layer thickness on scattering
stems from some of our 2D QLA simulations [15] where we have considered a
plane 1D electromagnetic pulse scattering from a small dielectric cylinder.
For $\Delta_{BL}<<|Delta_{pulse}$ we have found quite complicated structures
being emitted from the dielectric cylinder due to the internal bouncing from
the cylinder walls. For example, in Fig. 10, we have a 1D plane pulse
propagating to along the x- axis. In the middle of the plane is a 2D
dielectric cylinder of internal refractive index $n_{2}=3$ whose diameter
$>>\Delta_{pulse}$. The pulse speed within the dielectric cylinder is reduced
by a factor of $3$ over the rest of the plane propagation as it continues
along the axis. However, within the cylinder the pulse will bounce and then
undergo reflection and transmission at the dielectric cylinder boundaries.
This accounts for the multiple circular wave fronts for the electric field
being emitted at various times in the QLA simulation. Fig. 10 is a snapshot of
$E_{y}(x,z)$ at time t = $20k$. The 1D pulse can be seen as a dark thin strip
to the right of the dielectric cylinder.

Figure 10: The scattering from a dielectric cylinder centered at $x=L/2,z=L/2$
of a plane 1D pulse propagating along the x-axis. L is the number of lattice
sites in each direction. Plotted is the electric field $E_{y}(x,z,t)$ at time
$t=20k$. The multiple circular ring disturbances arise from the internal
scattering and then transmission out of the dielectric cylinder regime.

##  4 Acknowledgments

This research was partially supported by Department of Energy grants DE-
SC0021647, DE-FG02-91ER-54109, DE-SC0021651, DE-SC0021857, and DE-SC0021653.
The 2D QLA simulation was performed on $CORI$ of the National Energy Research
Scientific Computing Center (NERSC), A U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory,
operated under Contract No. DE-AC02-05CH11231. The 1D QLA simulations were
performed on a laptop.

##  5 References

[1] O. Laporte $\&$ E. G. Uhlenbeck, “Application of spinor analysis to the
Maxwell and Dirac equations”, Phys. Rev. 37, 1380-1397 (1931)

[2] J. R. Oppenheimer, “Note on light quanta and the electromagnetic field”,
Phys. Rev. 38, 725-746 (1931).

[3] E. Moses, “Solutions of Maxwell’s equations in terms of a spinor notation:
the direct and inverse problems”, Phys. Rev. 113, 1670-1679 (1959)

[4] I. Bialynicki-Birula, “Photon Wave Function” , in Progress in Optics, Vol.
34, pp. 248-294, ed. E. Wolf (North-Holland, 1996).

[5] E. Majorana (unpublished notes) : quoted after R. Mignani, E. Recami $\&$
M. Baldo, “About a Diraclike Equation for the Photon, according to Ettore
Majorana.”, Lett. Nuovo Cimento, 11, 568-572 (1974).

[6] S. Esposito, “Covariant Majorana Formulation of Electrodynamics”,
Foundations of Physics 28, 231-244 (1998)

[7] L. Silberstein, Annalen der Physik 22, 579 (1907); 24, 783 (1907)

[8] S. A. Khan, “Maxwell Optics: I. An exact matrix representation of the
Maxwell equations in a medium”, Physica Scripta 71, 440-442 (2005); also
arXiv: 0205083v1 (2002)

[9] Y. Yepez, ”An efficient and accurate quantum algorithm for the Dirac
equation”, arXiv: 0210093 (2002).

[10] J. Yepez, ”Relativistic Path Integral as a Lattice-Based Quantum
Algorithm”, Quant. Info. Proc. 5, 471-509 (2005)

[11] G. Vahala, L. Vahala, M. Soe $\&$ A. K. Ram, ”Unitary Quantum Lattice
Simulations for Maxwell Equations in Vacuum and in Dielectric Media” , J.
Plasma Phys. 86m 905860518 (2020)

[12] A. K. Ram, G. Vahala, L. Vahala, $\&$ M. Soe. ”Reflection and
transmission of electromagnetic pulses at a planar dielectric interface -
theory and quantum lattice simulations”, AIP Advances 11, 105116 (2021)

[13] G. Vahala, L. Vahala, M. Soe $\&$ A. K. Ram, “One- and two-dimensional
quantum lattice algorithms for Maxwell equations in inhomogeneous scalar
dielectric media I: theory” Rad. Effects Defects in Solids, 176, 49-63 (2021)

[14] G. Vahala, M. Soe, L. Vahala $\&$ A. K. Ram, “One- and two-dimensional
quantum lattice algorithms for Maxwell equations in inhomogeneous scalar
dielectric media II: Simulations” Rad. Effects Defects in Solids, 176, 64-72
(2021)

[15] G. Vahala, M. Soe, L. Vahala $\&$ A. K. Ram, “Two Dimensional
Electromagnetic Scattering from Dielectric Objects using Quantum Lattice
Algorithm”, arXiv:2110.05480

[16] G. Vahala, J. Hawthorne, L. Vahala, A. K. Ram $\&$ M. Soe, “Quantum
lattice representation for the curl equations of Maxwell equations”,
arXiv:2111.09745

[17] J. D. Jackson, “Classical Electrodynamics”, 3rd Ed., (Wiley, New York,
1998)

[18] A. Oganesov, G. Vahala, L. Vahala, J. Yepez $\&$ M. Soe, “Benchmarking
the Dirac-generated unitary lattice qubit collision-stream algorithm for 1D
vector Manakov soliton collisions”, Computers Math. with Applic. 72, 386
(2016)