# Click-Through Rate Prediction Using Graph Neural Networks and Online
Learning

Farzaneh Rajabi  Stanford University
<EMAIL_ADDRESS>Jack Siyuan He  Google LLC
<EMAIL_ADDRESS>
##  I Introduction

Recommendation systems have been extensively studied by many literature in the
past and are ubiquitous in online advertisement, shopping industry/e-commerce,
query suggestions in search engines, and friend recommendation in social
networks. Moreover, restaurant/music/product/movie/news/app recommendations
are only a few of the applications of a recommender system. A small percent
improvement on the CTR prediction accuracy has been mentioned to add millions
of dollars of revenue to the advertisement industry. Click-Through-Rate (CTR)
prediction is a special version of recommender system in which the goal is
predicting whether or not a user is going to click on a recommended item.

A content-based recommendation approach takes into account the past history of
the user’s behavior, i.e. the recommended products and the users reaction to
them. So, a personalized model that recommends the right item to the right
user at the right time is the key to building such a model. On the other hand,
the so-called collaborative filtering approach incorporates the click history
of the users who are very similar to a particular user, thereby helping the
recommender to come up with a more confident prediction for that particular
user by leveraging the wider knowledge of users who share their taste in a
connected network of users.

Using any of these two approaches, an accurate prediction demands framing the
problem in such a way that the nonlinear high-order interaction between
features of the dataset are taken into account. Graph neural Networks are
proposed to assist with this purpose. However, even by considering the high-
order feature interactions, framing the recommendation as a static procedure
and ignoring the dynamic interactive nature between users and the recommender
systems is yet a limiting factor, resulting in a discrepancy between offline
metrics for supervised learning and the online performance of the newly
proposed models. As a solution, reinforcement learning treats the
recommendation systems as a sequential decision making process that models the
interactions between the user and the recommender engine as a dynamic
adaptation scheme. That is, the agent observes the environment and acts upon
that in order to change its state towards better states (states with higher
rewards).

In this project, we are interested in building a CTR predictor using Graph
Neural Networks complemented by an online learning algorithm that models such
dynamic interactions. By framing the problem as a binary classification task,
we have evaluated this system both on the offline models (GNN, Deep
Factorization Machines) with test-AUC of 0.7417 and on the online learning
model with test-AUC of 0.7585 using a sub-sampled version of Criteo public
dataset consisting of 10,000 data points.

##  II Literature Review

Machine learning based CTR predictor has been thoroughly discussed in the past
literature. These models started with content-based to recommend items by
considering the content similarity between items. Later, collaborative
filtering was introduced with the rationale behind that the users with similar
behaviors tend to prefer the same items, and the items consumed by similar
users tend to have the same rating. However, conventional CF-based methods
tend to suffer from the data scarcity, and the similarity calculated from
sparse data can be very unreliable. The simplest and somewhat trivial model
involves Logistic Regression (LR) [1] that can only model first-order
interactions by linearly combining extracted features. Such a model works if
features are linearly dependent. However, when features have second-order
interactions with each other, this model may fail to provide accurate results.
Naturally, we would like to create joint features using two or more raw
features, thereby taking the correlation between the features into account.
However, this leads to a quadratic explosion of feature space, rendering
learning task infeasible.

Then we found that Factorization Machines (FM) [11] came to rescue. FM is able
to represent second-order features as factorized parameters whose space grows
linearly with respect to the total number of features. Using FM, researchers
were able to improve the CTR prediction where second-order feature
interactions were present [7].

However, natural features may interact beyond second-order, making FM not
sufficient for certain use cases. The initial approach taking by the
researchers was to put those features into a deep neural network (DNN). Past
works like DeepFM [6] were able to successfully demonstrate CTR prediction
improvement over FM. In order to de-bias the system towards the low- or high-
order feature interactions, DeepFM proposes an end-to-end learning model that
emphasizes on both low- and high order feature interactions. They combine the
power of factorization machines for recommendation and deep learning for
feature learning in a new neural network architecture. The FM part models the
linear (1st order) and pairwise (2nd order) feature interactions as inner
product of respective feature latent vectors. Whereas, the deep component is a
feed-forward neural network, which is used to learn high-order feature
interactions. Compared to the Wide & Deep model from Google [3], DeepFM has a
shared input to its “wide” and “deep” parts, with no need of feature
engineering besides raw features, making them a powerful practical model.

While being better than previous models, DNN-based models use bitwise factors
and the learning results are internal to the models themselves and therefore
they lack explainability. This can be bad for business critical models that
require logical explanation from their users. Moreover, these implicit models
do not allow manual input of existing models that could speed-up the training
process or improve the metrics. As the result, researchers turned to more
explicit models such as Deep&Cross [16] that takes outer product of features
at bit level.

Very recently, a graph-based neural network model, Fi-GNN [8], was developed
to explicitly model relations among features in a graph. In this work, the
multi-field features are represented in a graph structure, where each node
corresponds to a feature filed, and different fields interact with each other
through edges. Therefore, the task of modeling feature interactions is
converted to modeling node interactions on the corresponding graph. Generally,
nodes in GNNs interact with neighbors by aggregating information from
neighborhoods and updating their hidden states. There have been many variants
of GNN with various kinds of aggregators and updaters. For instance, Gated
Graph Neural Networks use GRUs as updater. Graph Convolutional Networks
utilize the convolutional aggregators. Graph Attention Networks incorporate
the attention mechanism into the propagation step. In [8], the input sparse
m-field feature vector is first mapped into sparse one-hot embedding vectors
and then embedded to dense field embedding vectors via the embedding layer and
the multi-head self-attention layer. The field embedding vectors are then
represented as a feature graph, where each node corresponds to a feature field
and different feature fields can interact through edges. An attention scoring
layer is applied on the output of Fi-GNN to estimate the CTR. A log-loss
function as below is used as the loss function,

$\displaystyle\mathcal{L}=-\dfrac{1}{N}{\sum}_{i=1}^{N}\bigg{(}y_{i}log(\hat{y}_{i})+(1-y_{i})log(1-\hat{y}_{i})\bigg{)}$
(1)

Also, the edge-wise interactions are achieved via attentional edge weights and
edge-wise transformations, and extra residual connection along with GRU is
introduced to update states, which can help regain the low-order information.
They have reported superior accuracy boost over the sate-of-the-art models for
CTR predictions evaluated on the public benchmark dataset (in terms of their
evaluation metrics: AUC, Logloss(cross entropy) and Relative Improvement
(improvement with respect to AUC, Logloss)).

That being said, yet most CTR predictors we found so far are offline models
that require re-training when new data is available and cannot improvise to
explore different recommendation strategies. There is a big issue regarding
offline models like such. Since data is created at enormous rate in today’s
world, we may not have enough data or time to re-train our offline model. For
example, a natural disaster like COVID-19 may suddenly change users’
preference without giving enough data and time to retrain the models.
Moreover, before retraining the model, offline models will predict CTR using
previously trained parameters, resulting in an unsatisfactory period between
the change and new model update. Moreover, offline models often fail to
address the internal state of users. For example, a user who has already
bought a washing machine would most likely not buy another one. In this case,
stateful models came to rescue. We read about works that use adversarial
sampling [17] for CTR prediction. Although such works exploited real-time
feedback into their pipeline, it requires prior knowledge of the model itself
and this knowledge about the model cannot be updated in real-time.

Subsequently, Reinforcement Learning (RL) models were introduced. Past work
[18] and [9] demonstrated the ability to use RL for news article
recommendation and movie recommendation, respectively. In the former
literature, they use a combination of online-offline models to predict users’
likelihood of clicking on a news article. Updates generated by RL is
incorporated into the model through a combination of minor and major updates
to avoid frequent disturbances to the model. A similar approach can be applied
in an advertisement CTR prediction task.

Model-based RL techniques are proposed to model recommendation procedure, such
as POMDP [13] and Q-learning [15]. However, these methods are inapplicable to
complicated recommendation scenarios when the number of candidate items is
large, because a time-consuming dynamic programming step is required to update
the model. Later, model-free RL techniques were introduced, divided into two
categories: value-based and policy-based. Value-based approaches compute
Q-values of all available actions for a given state and the one with the
maximum Q-value is selected as the best action. However, due to the evaluation
on overall actions, the approaches may become very inefficient if the action
space is too large. Policy-based approaches generate a continuous parameter
vector as the representation of an action, which can be utilized in generating
the recommendation and updating the Q-value evaluator. The inefficiency
drawbacks can be overcome due to the continuous representations. Yet, these
models have one common limitation: the user state is learnt via a conventional
fully connected neural network, which does not explicitly and carefully model
the interactions between users and items.

The approach presented in [9] overcomes this constraint using an “Actor-
Critic” type framework incorporated with a state representation module, which
explicitly models the complex dynamic user-item interactions, thereby treating
recommendation as a sequential decision making problem, in which the
recommender (i.e., agent) interacts with users (i.e., environment) to suggest
a list of items sequentially over the timesteps, by maximizing the cumulative
rewards of the whole recommendation procedure. The former literature is an
interesting model that could be applied on the Criteo dataset. However, given
our timeline for this class, we are inspired by the approach in [18] and will
follow a similar approach.

##  III Dataset

We used the Criteo dataset111https://www.kaggle.com/c/criteo-display-ad-
challenge/data originally released by Criteo Labs as part of the Criteo
Display Ads Challenge. The benefits of using such a dataset are:

  * •

It is freely available online for immediate download

  * •

The data format is very standardized

  * •

It’s one of the few publicly available benchmark datasets, used by many
researchers including [8], meaning that we can compare our results with theirs

The dataset consists of a portion of Criteo’s traffic over a 7 day period,
which is split into training and test dataset. Each row represents an AD shown
through Criteo’s display ADs network where the first column represents whether
it was clicked (label 1) or not (label 0) by that specific user. There are 39
anonymized feature fields; the first 13 columns are integer features, followed
by 26 categorical feature columns. The training dataset contains about 45
million rows and the test dataset contains about 6 million rows. Due to the
limiting size of the entire dataset, we used a small portion of it consisting
of 10000 user’s click records, and we divided it into 80% train, 10%
validation and 10% test portions.

The data comes in CSV format, which needs to be converted into NumPy format
first. Then the data needs to be properly scaled if log or quadratic feature
values are needed. AutoInt has provided an open source implementation of the
above prepossessing
steps222https://github.com/DeepGraphLearning/RecommenderSystems/tree/master/featureRec
[14]. In this process, the data is also segmented into 10 smaller parts so
that it fits into a normal computer’s memory.

Unfortunately, due to the time constraints, we could not investigate the data-
imbalance for the data portion that we used. However, we are aware that
investigating the data in terms of imbalance is the key to a reliable
accuracy. That is, if the number of examples with positive and negative labels
are not almost equally distributed, we have to do either majority-class down-
sampling or minority-class up-sampling as a potential solution to resolve this
issue. SMOTEBoost and RUSBoost have been proposed to alleviate the class
imbalance problem [2], [12].

That being said, we are suspect that such an issue could be the underlying
reason for not getting a higher AUC from our implemented model (and also, the
reason behind getting a lower accuracy than what [8] has reported). This came
into our mind after running the experiments and looking at the evaluation
metrics, but is definitely something of interest to inspect, and see whether
the performance improves or not.

##  IV Baseline

As for our baseline model, we have used the Fi-GNN’s open source
implementation. The Feature-interaction Graph Neural Network implementation
has modeled the low- and high- order feature interactions by modeling the
features as nodes and interactions as edges. Therefore, we used it as a pre-
training step for our online learning algorithm, and ran it with our
parameters. This GNN model has been described thoroughly in the literature
review section.

However, initially, after getting Fi-GNN’s code from GitHub, we obtained a
local test AUC of 0.5433 (meaning that model has almost no class separation
capacity) and a test LogLoss of 0.6999 as mentioned in our progress report.
This is very far away from their claimed AUC of around 0.8061 and LogLoss of
around 0.4455 in their paper [8]. We initially suspected that this discrepancy
was due to our usage of a sub-sampled smaller version of the Criteo dataset.
This is partially true. In their GitHub repository, they included a sub-
sampled version of the Criteo dataset that includes 10,000 data points,
whereas the full dataset has 45 million data points. However, this was not the
main reason for the difference.

Through our exploration, we found that Fi-GNN first splits the dataset into 10
equal partitions. For the 10,000 data point sample, each partition contains
1,000 data points. During training phase, their model loads each partition
into memory and train on its data in batches. There is one problem here: their
default batch size was set to 1024. This batch size is OK if the number of
data points in each data partition is large. For example, for the 45 million
data points full dataset, each partition would have 4.5 million data points
and training in 1024 data points per batch would still result in around 4,500
training iterations per data partition. However, for the sub-sampled dataset,
there is only 1000 data points per data partition. As the result, a batch size
of 1024 would result in only one training iteration per data partition. Since
Fi-GNN train its model on only eight data partitions, the total number of
training iteration when batch size is 1024 is only 80. The lack of sufficient
training iterations make the model not fully updated and therefore performs
very poorly.

After realizing that, we changed the batch size to 100 and the result improved
significantly. With 100 data points per batch, each 1000 data partition would
trigger 10 training iterations, summing up to 800 training iterations for all
8 training data partitions. With this change, we were able to get a test AUC
of 0.7417 and test log loss of 0.4687 using only the sub-sampled dataset. We
took this as our new baseline when evaluating our online learning models.

On a side note, we have also tried to use Deep-Factorization Machine
implementation as our second baseline. The main repository contains a lot of
models including Deep and Cross Neural Networks, however, we have only
employed whatever is relevant to the DeepFM model. This model is composed of a
deep component and a FM component that share the same input. The FM component
is a factorization machine, which is proposed in [11]. The deep component is a
feed-forward neural network that models the higher-order (higher than second-
order) feature interactions, resulting in a higher accuracy. The FM component
and deep component share the same feature embedding, which brings two
important benefits: It learns both low- and high-order feature interactions
from raw features. Also, there is no need for expertise feature engineering of
the input, as required in [3]. We got test-AUC of 0.6428 and test-Logloss of
1.5371 from this model without any hyper-parameter tuning.

##  V Main approach

Our goal in this project is not to implement a new CTR predictor, but to seek
improvements on an existing CTR predictor using some online learning
algorithm. Therefore, we took Fi-GNN as our base implementation. By default,
this model is only trained offline and cannot be updated during online
predictions. In our experiment configuration, we made several improvements to
this model to make it an online model that can be updated with more inputs.

Ideally, we want to be able to update the weight vector in the model’s neural
network using Reinforcement Learning algorithms like Q-Learning. However, this
model has a complicated seven layer graphical neural network that is buried in
a TensorFlow session. Getting the weights out and updating them via Q-Learning
will be extremely challenging and hard to achieve in the span of this quarter-
long class project. As an alternative, we opt for online training, in which we
feed test input and ground truth back to the Fi-GNN model and train it on the
fly.

In particular, as shown in our pipeline design in Figure 1, we firstly train
Fi-GNN using eight partition offline dataset of 8,000 data points. The Fi-GNN
model is then validated against another 1,000 data points. Then we throw in
the test dataset to simulate a real world ADs serving environment. In
particular we let the current Fi-GNN model predict the CTR given a test AD
input. Then after M test AD inputs, we accumulate all test ADs input in a
batch and use them to re-train the Fi-GNN model so that the model evolve into
a new generation. We then use the model to predict CTR for future ADs input.
This cycle continues until we exhaust all online ADs input. In our test
environment, we only have 1,000 online ADs input. In reality, there should be
infinite test ADs input when we try to serve user with some ADs from our
database.

Figure 1: Pipeline Design

We opt for this batched online training approach as we were inspired by a
similar work using Reinforcement Learning instead [18]. After getting the
pipeline running, we started tweaking the online learning batch size, M, and
observing how our evaluation metrics change with respect to changes in batch
size, M, especially the Area Under Curve (AUC) metric that would be explained
in the next section.

##  VI Evaluation Metric

The main evaluation metrics we used were similar to the ones in relevant
literature; AUC (Area Under the ROC curve) [5] and LogLoss under the same
definition as Fi-GNN [8]. In particular:

  * •

AUC measures the probability that a clicked item will be chosen over a
randomly selected unclicked one. Therefore, the higher the AUC, the better. It
tells how much model is capable of distinguishing between classes.

  * •

LogLoss measures the distance between the predicted CTR and the true label for
each AD, the lower the better.

  * •

Precision is the ratio (true positive)/(true positive + false positive). It’s
the ability of the classifier not to label as positive a sample that is
negative. It ranges from 0 (worst) to 1 (best), and is a measure of the
models’ relevancy.

  * •

Recall is the ratio (true positive)/(true positive + false negative). It’s the
ability of the classifier to find all the positive samples. It ranges from 0
(worst) to 1 (best), is a measure of the model’s completeness.

  * •

F1-score is the weighted average of precision and recall, where it has its
best value at 1 and worst value at 0. So, F1= 2*(precision*recall)/(precision
+ recall). As the output from our model ($\hat{y}$) is in fact a probability,
while the label vector $y$ has only 0 and 1s, we have plotted this measure for
different thresholds above which we take $\hat{y}$ to be 1, and otherwise 0.

We have also included the ROC curve (receiver operating characteristic curve)
which is a graph showing the performance of a classification model at all
classification thresholds. Confusion matrices for the baseline and main models
are included as well.

If we had more time, we were also going to look at another metric called
Relative Improvement (RI) of AUC and Logloss, defined as,

$\displaystyle RI-X=\dfrac{|X(model)-X(base)|}{X(base)}*100\%$ (2)

Where $X$ is AUC or Logloss. This is a useful metric because a small
improvement with respect to AUC is regarded significant for real-world CTR
tasks. But, it cold be an insightful metric if we had used a larger/more
complete version of this dataset. So, we planned to postpone it to future
work. Even though there are a few qualitative metrics for evaluating
recommender systems in general (we have mentioned some of them in the Future
Work section), all of the CTR prediction literature we studied did suffice to
AUC and Logloss, without really reporting a qualitative evaluation metric.
That’s because those qualitative metrics are more suited for a multi-label
recommendation system e.g. Netflix’s movie recommender, which has a little
different nature from what we are trying to predict here. Nevertheless, we are
reporting the recall, precision and F1-score on top of AUC and Logloss, in
order to offer a complete analysis for our classification problem.

Moreover, we also tried to evaluate the effect of
online_learning_to_predict_ratio on AUC value.

##  VII Results and Analysis

There has been several interesting discoveries from our experiments. We
initially set our batch size to 100, meaning that it learns after every 100
test input. Surprisingly, the resultant AUC actually went down to 0.7273. We
explored further by increasing the learning batch size to 200. To our delight,
the result actually went up to 0.7495 AUC. This observation challenges the
naive view that the more frequently we update our model, the better it tracks
the test input. Instead, if we update our model too closely to the test
inputs, our training becomes greedy and the model become overfitted with the
previous test input, resulting in worse model performance.

However, if the online learning batch size goes to high, it will also hinder
our model’s online performance. The larger the online learning batch size, the
fewer learning iteration would the model go through during online operations.
As result, the model wouldn’t learn very frequently and may not be able to
capture some sudden changes in user preferences. Moreover, since we have a
limited size test input of 1000 data points, a large online learning batch
size might leave no further ADs input after just a few learning batches. For
example, when the batch size goes to 500 and above, the model will only be
retrained once before exhausting all test input. This is clearly shown in
Figure 2 and Table I, which shows the changes in AUC with respect to the
baseline result mentioned in earlier section. From the results, we can see
that the AUC is actually maxed out at 400 batch size and we improved it by
about RI-AUC= $1.7\%$.

Figure 2: AUC vs learning batch size Models | AUC | Log loss
---|---|---
Baseline-DeepFM | 0.6428 | 1.5371
| Fi-GNN
---
(training batch size 1024,
no online learning)
0.4918 | 0.8142
| Fi-GNN
---
(batch size 10,
no online learning)
0.7533 | 0.4591
| Baseline Fi-GNN
---
(batch size 100,
no online learning)
0.7418 | 0.4688
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 100)
0.7273 | 0.4939
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 200)
0.7495 | 0.4722
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 300)
0.7512 | 0.4622
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 400)
0.7585 | 0.4666
| Fi-GNN +OL
---
(batch size 100,
online learning batch size 500)
0.7542 | 0.4589
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 600)
0.7472 | 0.4696
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 700)
0.7473 | 0.4916
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 800)
0.7486 | 0.4711
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 900)
0.7475 | 0.4730
| Fi-GNN + OL
---
(batch size 100,
online learning batch size 1000)
0.7427 | 0.4639
TABLE I: Results with Different Batch Sizes

The improvement is still rather small because our offline model was trained on
8,000 data points, whereas we only have 1,000 test inputs. In reality, we
could see more test inputs that will likely improve our online learning
model’s performance.

##  VIII Error Analysis

We have reported our analysis of this study both in terms of the AUC and
Logloss values for every experiments we have performed. In Table I, as well as
the ROC curve (Figure 3), F1-score curve (Figure 4) for different threshold
values (the reported $\hat{y}$ from models is a probability, and thus
converting it to a binary label to make it comparable to our true label
requires a threshold value). We have also reported precision vs recall in
Figure 5, precision vs threshold in Figure 6, and recall vs threshold in
Figure 7. All the reported results are the average of three runs for each
case; because results may vary given the stochastic nature of the algorithm or
evaluation procedure, or differences in numerical precision.

  * •

All of the plots correspond to two models, one baseline, and one baseline plus
online batch learning. The baseline Fi-GNN model has batch-size 100. The Fi-
GNN plus online learning model has batch size of 100 and learning batch size
of 400. All results are reported for the optimum threshold value of 0.3 as
represented by the F1-score vs threshold plot. Note that the two models being
plotted have very close AUC, so they are not far from each other on the plots.
If we had used the model with AUC of 0.4819 (second row in the table), we
could clearly see higher improvements.

  * •

According to the ROC curve, the experiment (Fi-GNN plus online learning) has a
slightly higher area under curve than that of the baseline with no online
learning, meaning that experiment has done a better job than the baseline in
classifying the positive class in the dataset. However, as we will see in the
confusion matrix, the experiment has done a perfect job in classifying the
negative class, but has done poorly on classifying the positive class. That’s
because ROC curve is used when there are roughly equal numbers of observations
for each class. Whereas, when there is a moderate to large class imbalance,
Precision-Recall curve gives a better insight. The reason for this
recommendation is that ROC curves present an optimistic picture of the model
on skewed datasets. Therefore, the Precision-Recall Plot is more informative
than the ROC Plot when evaluating binary classifiers on imbalanced datasets.
We did investigate our dataset after running models, and it was indeed
imbalanced, with the majority being negative class. This makes the classifier
highly biased

  * •

The Precision vs. threshold curve reveals that higher thresholds result in
higher precisions and, high precision relates to a low false positive rate.
The experiment model shows slightly higher false positive rates than the
baseline, for most thresholds.

  * •

The Recall vs. threshold curve reveals that lower thresholds result in higher
recall and, high recall relates to a low false negative rate. also, the online
learning recall curve vs the baseline shows that online learning has slightly
less false negative rate.

  * •

The sweet spot between precision abd recall threshold values is where F1-score
is maximized. We see from F1-score vs. threshold curve that optimal value of
0.3 maximizes F1-score. So, at a threshold of 0.3 for separating the
probabilities of the positive class from negative class, F1-score has its
maximum value.

  * •

According to the confusion matrices for both baseline and our online learning
experiment, while the baseline model does a reasonable job in classifying the
true positives and true negatives (Precision= 0.41, Recall=0.57,
Accuracy=0.73, negative class accuracy=0.77), the online learning model does a
perfect job on predicting negative class (96% correct), and does a poor job on
predicting the positive class (20% correct). For this experiments we got
Precision= 0.55, Recall=0.2, Accuracy=0.8, negative class accuracy=0.96. The
high false negative rate again confirms how skewed our dataset is. The recall
for this model is much lower than the baseline’s. So, any improvements on this
model (making the dataset balanced and then improving the model) could start
by making changes and watching how recall is impacted.

Figure 3: ROC Curve Figure 4: F1 Score vs Threshold Figure 5: Precision vs
Recall Figure 6: Precision Score vs Threshold Figure 7: Recall Score vs
Threshold

Figure 8: Confusion matrices for the baseline model and the best experiment,
both with threshold 0.3 and training batch size of 100. The experiment matrix
has a learning batch size of 400

##  IX Future Work

Due to the time constraint, we could not finish a lot of ideas that we had
planned for this problem. If we had more time, we were going to finish our
implementation of DRL approach proposed by [9], and compare it against the
online-learning algorithm we implemented on top of the Graph-interaction
Network, as well as our other baseline Deep Factorization Machines which is
still widely being used as a promising model in the practical/industrial
applications.

Also, we did not have enough time to perform the experiments on larger version
of the Criteo dataset. We believe that using a larger version could
extensively improve our reported accuracies. Also, we could not do much hyper-
parameter tuning on our experiments. This is a major step in any ML model, but
we just postponed it to a future attempt.

Inspecting data-imbalance was another point mentioned before. We did not check
this here, but we believe it would considerably boost our accuracies if we
make sure the dataset in not unreasonably imbalanced in terms of positive and
negative classes.

One future direction could be building a recommendation engine using a
powerful and efficient dynamic adaptation algorithm (like an Actor-Critic RL
scheme), which is able to recommend the top-K items that users are most likely
to click based on each user’s history of clicking on the past recommended
items. This is similar to how the Netflix’s Movie recommendation algorithm
does an offline testing, followed by an online A/B testing, resulting in the
top-K movies to recommend 333https://netflixtechblog.com/netflix-
recommendations-beyond-the-5-stars-part-1-55838468f429,
444https://netflixtechblog.com/netflix-recommendations-beyond-the-5-stars-
part-2-d9b96aa399f5. Having such a top-K recommender, we could analyze a wider
range of RS-specific quantitative and qualitative metrics like: top-N hit-
rate, average reciprocal rank, cumulative hit-rate, rating hit-rate, coverage,
diversity, novelty, churn, responsiveness, and so on. Yet, besides all of
these, the online A/B testing has been mentioned as the most reliable metric
for making better recommendations.

To this end, by simplifying the problem into a CTR predictor, we did not have
to deal with the K nearest neighbor pruning that can be technically
challenging and was beyond our time limit. Nevertheless, by recasting the CTR
problem into a multi-label binary classification approach (and applying an
algorithm like Binary Relevance for the multi-label problem if there’s no
correlation between the different recommendations, or using a label powerset
approach if we believe there is correlation between the recommended
campaigns), and outputting the probability of the user clicking on a specific
recommendation, followed by sorting all the outputted probabilities, and then
picking the top-K (e.g. top-10), we can unambiguously extend our models into a
top-K recommender system.

Another idea we could try to improve AUC was performing an ablative study
similar to the one in [8]. Their idea was adding an edge-wise node
interactions via attentional edge weights and edge-wise transformation, as
well as extra residual connections to update state along with GRU in a GGNN.
We could remove any of these additions (e.g. edge-wise interactions) and rerun
the online learning to see whether we get any improvements.

Another very interesting approach that we would like to investigate in the
future, is implementing a Meta-learning approach [4], [10]. The motivation
behind such model is that the recommender sytem’s dataset are usually very
heterogeneous and highly personalized, and no single model could give the best
recommendation for every user. Therefore, a user-level adaptive model
selection will train a collection of different recommendation models
(including DFM, Fi-GNN, DRL, and so on) with data from all users, on top of
which a model-selector is trained via meta-learning to select the best single
model for each user with the user-specific historical data. Such an approach
has been reported to achieve improvements over the single-model baseline.

## References

  * [1] Olivier Chapelle, Eren Manavoglu and Romer Rosales  “Simple and Scalable Response Prediction for Display Advertising”  In _ACM Trans. Intell. Syst. Technol._ 5.4 New York, NY, USA: Association for Computing Machinery, 2015  DOI: 10.1145/2532128
  * [2] Nitesh V Chawla, Aleksandar Lazarevic, Lawrence O Hall and Kevin W Bowyer  “SMOTEBoost: Improving prediction of the minority class in boosting”  In _European conference on principles of data mining and knowledge discovery_ , 2003, pp. 107–119  Springer
  * [3] Heng-Tze Cheng et al.  “Wide & deep learning for recommender systems”  In _Proceedings of the 1st workshop on deep learning for recommender systems_ , 2016, pp. 7–10
  * [4] Tiago Cunha, Carlos Soares and André CPLF Carvalho  “Metalearning and Recommender Systems: A literature review and empirical study on the algorithm selection problem for Collaborative Filtering”  In _Information Sciences_ 423 Elsevier, 2018, pp. 128–144
  * [5] Tom Fawcett  “An introduction to ROC analysis” ROC Analysis in Pattern Recognition  In _Pattern Recognition Letters_ 27.8, 2006, pp. 861–874  DOI: https://doi.org/10.1016/j.patrec.2005.10.010
  * [6] Huifeng Guo et al.  “DeepFM: A Factorization-Machine based Neural Network for CTR Prediction”, 2017  arXiv:1703.04247 [cs.IR]
  * [7] Yuchin Juan, Yong Zhuang, Wei-Sheng Chin and Chih-Jen Lin  “Field-Aware Factorization Machines for CTR Prediction”  In _Proceedings of the 10th ACM Conference on Recommender Systems_ , RecSys ’16  Boston, Massachusetts, USA: Association for Computing Machinery, 2016, pp. 43–50  DOI: 10.1145/2959100.2959134
  * [8] Zekun Li et al.  “Fi-gnn: Modeling feature interactions via graph neural networks for ctr prediction”  In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_ ACM, 2019  DOI: 10.1145/3357384.3357951
  * [9] Feng Liu et al.  “Deep reinforcement learning based recommendation with explicit user-item interactions modeling”  In _arXiv preprint arXiv:1810.12027_ , 2018
  * [10] Mi Luo et al.  “MetaSelector: Meta-Learning for Recommendation with User-Level Adaptive Model Selection”  In _Proceedings of The Web Conference 2020_ , 2020, pp. 2507–2513
  * [11] Steffen Rendle  “Factorization machines”  In _2010 IEEE International Conference on Data Mining_ , 2010, pp. 995–1000  IEEE
  * [12] Chris Seiffert, Taghi M Khoshgoftaar, Jason Van Hulse and Amri Napolitano  “RUSBoost: A hybrid approach to alleviating class imbalance”  In _IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans_ 40.1 IEEE, 2009, pp. 185–197
  * [13] Guy Shani, David Heckerman and Ronen I Brafman  “An MDP-based recommender system”  In _Journal of Machine Learning Research_ 6.Sep, 2005, pp. 1265–1295
  * [14] Weiping Song et al.  “Autoint: Automatic feature interaction learning via self-attentive neural networks”  In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_ , 2019, pp. 1161–1170  ACM
  * [15] Nima Taghipour and Ahmad Kardan  “A hybrid web recommender system based on q-learning”  In _Proceedings of the 2008 ACM symposium on Applied computing_ , 2008, pp. 1164–1168
  * [16] Ruoxi Wang, Bin Fu, Gang Fu and Mingliang Wang  “Deep & Cross Network for Ad Click Predictions”  In _CoRR_ abs/1708.05123, 2017  arXiv: http://arxiv.org/abs/1708.05123
  * [17] Yikai Wang et al.  “Regularized Adversarial Sampling and Deep Time-Aware Attention for Click-Through Rate Prediction”  In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_ , CIKM ’19  Beijing, China: Association for Computing Machinery, 2019, pp. 349–358  DOI: 10.1145/3357384.3357936
  * [18] Guanjie Zheng et al.  “DRN: A Deep Reinforcement Learning Framework for News Recommendation”  In _Proceedings of the 2018 World Wide Web Conference_ , WWW ’18  Lyon, France: International World Wide Web Conferences Steering Committee, 2018, pp. 167–176  DOI: 10.1145/3178876.3185994


# High energy density plasma mediated by collisionless resonance absorption
inside dielectrics

Kazem Ardaneh<EMAIL_ADDRESS>Remi Meyer  Mostafa Hassan  Remo Giust
FEMTO-ST Institute, Univ. Bourgogne Franche-Comté, CNRS,15B avenue des
Montboucons,
25030, Besançon Cedex, France  Chen Xie  Key Laboratory of Opto-electronic
Information Technology of Ministry of Education,
School of Precision Instruments and Opto-electronics Engineering,
Tianjin University, 300072 Tianjin, China
Benoit Morel  Ismail Ouadghiri-Idrissi  Luca Furfaro  Luc Froehly  FEMTO-ST
Institute, Univ. Bourgogne Franche-Comté, CNRS,15B avenue des Montboucons,
25030, Besançon Cedex, France  Arnaud Couairon  CPHT, CNRS, Ecole
Polytechnique, Institut Polytechnique de Paris,
Route de Saclay, F-91128 Palaiseau, France
Guy Bonnaud  CEA, Centre de Paris-Saclay, DRF, Université Paris-Saclay,
91191 Gif-sur-Yvette, France
Francois Courvoisier<EMAIL_ADDRESS>FEMTO-ST Institute,
Univ. Bourgogne Franche-Comté, CNRS,15B avenue des Montboucons,
25030, Besançon Cedex, France

###### Abstract

We demonstrate for the first time to our knowledge the generation of
overcritical plasma densities inside transparent solids over long distances
using femtosecond laser pulses. This opens new avenues for high energy density
physics in confined geometry such as warm dense matter study or the synthesis
of new material phases. We show both with experiments and first-principles
simulations, that femtosecond conical interference via a Bessel beam creates a
dense plasma rod with typically 100 nm diameter in sapphire. The interaction
is in ideal conditions to trigger collisionless resonance absorption. This
mechanism plays a primary role in the energy deposition process, yielding a
plasma with an energy density on the order of MJ/cm3 and a length that can
reach several cm using only tabletop femtosecond lasers.

Extreme focusing of a femtosecond laser pulse with numerical apertures in
excess of unity inside transparent materials generate a nanosphere of dense
plasma, but the volume of the excited material is too small for X-ray
investigation of warm dense matter [1, 2] or for the synthesis of new material
phases [3, 4]. Reducing the focusing strength enlarges the focal volume, but –
in contrast with laser interaction at the surface of solids – the plasma
density build-up inside the solid defocuses a Gaussian beam so that high
plasma densities cannot be reached, even with higher pulse energy. The problem
is illustrated in Fig. 1(a): after plasma generation on the optical axis by
the initial part of a laser pulse, the trailing part of the pulse impinges on
the plasma with a high angle of incidence $i$ with respect to the plasma
surface normal. Instead of heating the plasma core, the optical pulse is
deviated outwards to the regions with higher refractive index, even for very
low plasma densities. This is the wave-turning phenomenon [5, 6] which
presents a fundamental challenge. Importantly, the most intense components of
a Gaussian beam, at the center of its far-field distribution, have largest
incidence angles. They turn away at the largest distance from the plasma
center [see circled inset in Fig. 1(a)]. Due to the wave-turning, plasma
generation continues outside the initial volume and consequently the plasma
remains sub-critical, i.e., below the critical density at which the plasma
permittivity turns to zero, and the deposited energy density remains low [7,
8].

Figure 1: Laser-plasma interaction within dielectrics. (a) Gaussian beams are
conventionally limited to create sub-critical plasmas in dielectrics because
of the wave-turning phenomenon. Most of the waves composing the Gaussian beam
impinge on the plasma with a very large incidence angle: these are deviated
away well before they reach the plasma center. The plasma buildup expands the
plasma volume but the density remains sub-critical. The radius of the turning
point $r_{\rm t}$ is defined for an incidence angle $i$ by the density reached
$n(r_{\rm t})=n_{\rm c}\cos^{2}i$, where $n_{\rm c}$ is the plasma critical
density. (b) In a Bessel beam, the far-field k-space spectrum is narrow, the
whole pulse energy impinges on the plasma with a relatively low incidence
angle. The steep plasma density profile generated by the pulse onset makes
resonance absorption possible, which transfers energy from the laser beam into
the plasma. The energy transfer and ionization are therefore extremely
localized on a sub-wavelength rod. The geometry of the interaction is
longitudinally invariant.

In this article, we report the creation of over-critical plasmas along quasi-
arbitrary distances with femtosecond Bessel beams in sapphire, using only a
moderately high numerical aperture (0.4). The plasma absorbs high amount of
the pulse energy as a result of the resonance absorption mechanism, observed
here for the first time in the bulk of solids.

A zeroth-order Bessel beam is a circularly-symmetric interference field, so-
called diffraction free [9], whose length is independent from the central lobe
diameter [Fig. 1(b)]. Its conical structure uniformly distributes energy by a
transverse flow along the plasma created in the central core. Resonance
absorption has been widely studied within the context of inertial confinement
fusion, by illuminating the surface of solids [10, 11, 12]. It occurs when a
laser pulse obliquely impinges on a plasma ramp and the laser electric field
has a component parallel to the plasma density gradient. In this condition,
mode conversion resonantly drives electron plasma waves at the critical
surface. This can transfer up to 70% of the laser energy into the plasma [13,
5]. Mode conversion is efficient for the moderate incidence angle of Bessel
beams on the plasma density ramp and generally inefficient in Gaussian beams
because most of the energy is located on components at grazing incidence.

We determine the typical plasma profile by comparing complementary
experimental diagnostics to first-principles Particle-In-Cell (PIC)
simulations using the 3D massively parallel code EPOCH [14] because direct
measurement of the plasma density inside the bulk cannot be performed via
optical techniques when the plasma reaches the critical density at sub-
wavelength scales. We demonstrate that in sapphire, over-critical densities
are reached over a typical radius of 75-200 nm with typically 50% absorption
of the laser pulse energy. This is additionally confirmed by the observation
of second harmonic from the bulk of a centrosymmetric dielectric for the first
time to our knowledge. While conventional models of laser energy deposition
hitherto have relied on collisional absorption, we show here that most of the
absorption is due to collisionless resonance absorption. Our findings open a
new route to develop high energy density physics in a confined geometry inside
solids, but also open new perspectives for nonlinear optics with epsilon-near-
zero materials [15], synthesis of new material phases and laser material
processing.

Figure 2: Near-field fluence distribution of the Bessel beam from PIC
simulations. Fluence distribution in $xy-$plane at $z=7$ µm (peak of the axial
fluence profile), (a) in the absence of plasma and (b) in the presence of a
circular plasma rod with a transverse Gaussian density profile. (c) Zoom-in
view of the central part of (b).

The interaction of an ultrashort pulse with a dielectric medium can be viewed
as a two-stage process. In the first stage, field ionization and electron-
impact ionization create a plasma, with which the trailing part of the laser
pulse interacts during the second stage. Our experiments show that the
femtosecond Bessel pulse generates a plasma over a length of 18 µm with a
typical diameter of a few hundreds of nanometers and absorbs approximately
$50$% of the 1 µJ pulse energy (similar plasma was produced over a length
close to 1 cm with mJ pulse energy [16]). We point out that the ionization
energy absorbed to produce an over-critical plasma of this size is $n_{\rm
c}VU_{\rm g}$ where $n_{\rm c}$ is the critical density at 800 nm, $V\approx
1.5$ µm3 the plasma volume and $U_{\rm g}\approx 10$ eV the bandgap of
sapphire. The ionization energy is less than 0.04 µJ which is negligible in
comparison with the total absorbed energy. We therefore hypothesize that the
first stage is very short compared to the pulse duration and we focus our
modeling efforts on the second stage of the interaction. Our model does not
take into account ionization in order to isolate the laser-plasma interaction.

We first qualitatively describe the physics of Bessel beam interaction with
plasma using numerical simulations. We have performed 3D PIC simulations to
model the interaction of a 100 fs Bessel beam with a plasma rod. The pulse,
with a central wavelength $\lambda=800$ nm, propagates along $z-$ direction
and is polarized along the $x-$ direction [17]. In the absence of plasma, the
intensity reaches $6\times 10^{14}\,{\rm W/cm^{2}}$, as in the experiments.
The beam cone angle is $\theta=25^{\circ}$, which corresponds to an incidence
angle $i=65^{\circ}$ on the cylindrical surface of the plasma [Fig. 1(b)]. In
the central lobe the amplitude of the $z-$component of electric field is 25%
of the $x-$component, while the $y-$component remains negligible (2%). We have
presented in detail our numerical scheme in the Supplemental Material [18]. In
Fig. 2(a), we show the fluence distribution, i.e., the time-integration of the
Poynting vector of the Bessel beam, in $xy-$plane at a propagation distance
corresponding to the peak of axial fluence profile, without plasma. In Fig.
2(b), the beam interacts with an over-critical circular plasma rod with a
turning point radius of 274 nm. All lobes are shifted outwards by a distance
of about 300 nm due to the wave-turning phenomenon (the slight discrepancy
stems from the fact that the estimate of the turning point radius is made for
a monochromatic plane wave).

In the zoom-in view presented in Fig. 2(c), two high-intensity structures
appear inside the region of the initial central lobe. One is composed of two
electromagnetic lobes outside the plasma ($y=\pm 0.5$ µm), oriented
perpendicular to the pulse polarization. They originate from the difference in
permittivity between the plasma rod and the external medium [19, 20]. This
effect occurs even for sub-critical plasmas. It reduces the amplitude of the
field in the polarization direction and amplifies it in the perpendicular
direction. The second structure is much more localized in the vicinity of the
critical surface and shows two lobes along the polarization direction. It
corresponds to resonantly driven electrostatic plasma waves at the critical
surface. Therefore, this structure cannot be observed experimentally by beam
imaging. This is why, in the following figures, we exclude the electrostatic
field components from the simulated fluence distributions.

Figure 3: Near-field fluence distributions in sapphire: experiments (left) and
PIC simulation (right). (a) and (h) Slices in $xy-$plane at $z=20$ µm in (a)
and $z=7$ µm in (h). (b) and (i) Zoom-in views of (a), (h) respectively. (c)
and (f) Slices of the fluence in the $zx-$plane ($y=0$). (d) and (g) Evolution
of pulse energy as a function of propagation distance $z$ for collisional
(solid blue line) and collisionless simulation (dashed red line). (e)Scanning
electron microscopy images of void channels produced in single shot
illumination for different polarizations. The leftmost figure is for circular
polarization and the channel cross-section is circular. In all other cases,
the polarization is linear, oriented along the white arrow.

In Figs. 3(a)-3(d), we present experimental results of near field fluence
imaging in 400 µm thick C-cut sapphire Al_2O_3. They have been obtained by
collecting the 100 fs laser pulse after its interaction with the sample, using
a single shot beam sectioning technique [21]. We have detailed our
experimental techniques in the Supplemental Material [18]. Figures 3(a) and
3(b) show the $xy-$cuts of the distribution at the peak of the axial profile
($z=20$ µm). Figure 3(c) shows the fluence distribution in $zx-$plane while
Fig. 3(d) displays the evolution of the pulse energy as a function of
propagation distance, measured with the imaging technique.

We observe a significant energy drop of nearly 50$\%$ over an 18 µm segment of
propagation, which corresponds to the region where plasma is generated. This
is approximately the length of the void channel that is created in sapphire in
single shot under similar laser illumination conditions [22]. We independently
checked that the energy drop is effectively equal to the total loss within the
sample, confirming that energy deposition is localized in the high-intensity
region, i.e., not in an earlier stage of the propagation as it can be the case
for Gaussian beams.

Importantly, we observe in Fig. 3(c) that in the plasma region ($z=14-32$ µm),
the surrounding lobes negligibly shift with respect to those in the regions
where plasma is absent ($z<14$ µm, $z>32$ µm). This sets the upper limit to
the turning point radius $r_{\rm{t}}$ roughly to the experimental sensitivity
associated with the width of the Bessel lobes, i.e., typically 500 nm. Given
the incidence angle, the density at a distance $r_{\rm{t}}$ is less than
$n_{\rm c}\cos^{2}i=0.086n_{\rm{c}}$.

In Figs. 3(a) and 3(c), the shift due to wave turning is invisible within the
central lobe. Because of the limited numerical aperture, the fine details of
the fields inside the central lobe are not transferred to the imaging plane.
To summarize, our experiments show that the onset of the Bessel pulse
generates a plasma with a sub-wavelength diameter and that the trailing part
of the laser pulse is significantly absorbed.

We determine the typical plasma density profile with a series of numerical
simulations with different plasma profiles. We limited the turning radius to
$r_{\rm t}=500\,{\rm nm}$ and chose the plasma transverse profile as a
Gaussian distribution for which we adjust peak and width. Collisional PIC
simulations for sub-critical plasmas were showing an absorption of only a few
percent, even for collision times as short as $0.1-2.5\,{\rm fs}$. As an
example, with a peak density of $0.25~{}n_{\rm{c}}$, the turning radius has to
exceed $900\,{\rm nm}$ to reach an absorption of 20%. However, this would
cause a remarkable radial shift of the Bessel beam lobes which is not found in
our experiments. In contrast, only simulations with over-critical densities
were showing strong absorption in the range of $40-60$% as in our experimental
results.

Figure 4: Angular distribution of far-field fluence for the first and second
harmonics.The simulations share the same parameters as in Fig. 3.

The experimental near field profile of Figs. 3(a) and 3(b) is not
cylindrically-symmetric. With circular plasma rods, we were not able to
reproduce the orientation of the asymmetries, neither the far-field
distribution [see Fig.4(a)]. In contrast, we obtain a very good quantitative
agreement with all our experimental results using elliptical plasma rods,
having the major axis oriented perpendicular to the polarization (critical
radius $r_{\rm c}=75-190$ nm along the polarization direction and $r_{\rm
c}=380$ nm perpendicular to the polarization).

The ellipticity can be explained by the field amplification in the first stage
of field ionization, due to the difference in permittivity between the medium
and the initially sub-critical plasma [19, 20]. We experimentally confirmed
that the plasma is elliptical and elongated perpendicular to the polarization,
by high-resolution imaging of voids formed in sapphire after the laser-induced
micro-explosion, as shown for different input polarizations directions in Fig.
3(e).

In Figs. 3(f)-3(i), our simulations of the near field show negligible lobe
shift of the Bessel beam and the asymmetries of the fluence distribution are
oriented as in the experiment (we remind that the small dip observed in the
central lobe of Figs. 3(h-i) is filtered by the imaging operation). In the
collisional plasma model [Fig. 3(g) solid line], the pulse energy decreases
linearly with a slope of $-0.03$ µJ/µm and the total absorption is 54%, both
in good agreement with the experimental values. Overall, a steep plasma
profile tends to increase the absorption factor and a narrow plasma density
profile reduces the shift due to wave-turning. Importantly, the dashed curve
in Fig. 3(g) shows the result for a collisionless simulation. Since the
difference with the collisional case is only a few percent (5 points), the
main absorption mechanism is collisionless resonance absorption.

Our simulations show that the plasma diameter and its cross-section geometry
strongly affect the far-field pattern. In the linear, absorption-free, regime,
the far-field fluence of a Bessel beam is a homogeneous ring [Fig. 1(b)].
Supplementary Fig. S.1(a) shows the experimental far-field fluence
distribution collected after the interaction. Its angular distribution is
illustrated in Fig. 4(a) (solid blue). It shows two lobes oriented in the
direction of the input pulse polarization. The far-field from the simulation
[Fig. 4(a), red dashed line] is in quantitative agreement with our
experimental data for the elliptical plasma rod.

Overall, our simulations allowed us, up to here, to identify that the plasma
is over-critical and with an elliptic cross-section.

We experimentally observed second harmonic emission from the bulk of sapphire
despite its zero second-order nonlinear susceptibility [23]. For this, we
imaged the far-field fluence distribution in the spectral range 400$\pm$20 nm,
as shown in Supplementary Fig. S.1(c). Our experimental technique is described
in the Supplemental Material [18]. The distribution is composed of a low-
intensity background originating from black-body radiation, and of a thin
second harmonic emission at an angle nearly identical to the reflection of the
pump. This reflection is superimposed to the transmission of the pump because
of the conical interaction. We show in Fig. 4(b) that our PIC simulations are
in very good agreement with the experimental results.

In laser-plasma interaction, second-harmonic emission originates either from
resonance absorption or from parametric decay instability [24]. In the latter
case, the second harmonic polarization is perpendicular to the incident laser
one [25] while in the former case, it is parallel [26, 27], as in our
experimental results. Therefore, the second-harmonic diagnostic confirms that
resonance absorption occurs.

In the simulations, the ratio of second harmonic to the fundamental intensity
is $\sim 10^{-5}$. Using this ratio for Fig. S.1(c) and Stefan-Boltzmann law,
one can estimate the plasma temperature on the order of 10 eV. This is
consistent with our numerical framework and with the fact that Bessel beam-
induced plasma relaxation forms warm dense matter in silicon dioxide [28].

We have demonstrated, for the first time to our knowledge, i) that ultrafast
Bessel beams can generate over-critical plasma within the bulk of sapphire,
ii) that collisionless resonance absorption is a key process for energy
deposition and iii) second harmonic emission in single shot from the bulk of
sapphire. Our findings therefore open a new route to generate high energy
density matter inside solids. They answer a central question in the energy
deposition by ultrafast pulses within the bulk of solids, which is
particularly crucial in the field of laser micro-/nano- structuring of
materials. Since the conical shape is propagation-invariant, our results can
be easily extended to much longer propagation distances: a pulse energy of 1
mJ is enough to create a 1 cm long plasma in identical conditions as discussed
above [16].

Nanoplasma generation with over-critical density inside solids opens
interesting opportunities for nonlinear plasmonics [29], strong THz fields
generation [30] and extreme ultraviolet sources [31]. Controlling resonance
absorption inside dielectrics provides useful new tools to drive solids into
extreme thermodynamics pathways in a geometry where the material cannot expand
into vacuum. The energy density deposited in our experiments is on the order
of MJ/cm3, which is a typical range for driving a transformation to warm dense
matter [32, 33]. Consequently, millijoule energies are enough to generate
relatively large volumes of warm dense matter on the order of $1000$ µ$\rm
m^{3}$.

Our results therefore provide new avenues for generating relatively large
volumes of new material phases and our configuration can be considered as a
convenient new platform for the exploration of warm dense matter physics.

###### Acknowledgements.

K.A and R.M equally contributed, with the largest fraction of, respectively,
simulation work and experimental work. Technical assistance by C. Billet and
E. Dordor as well as fruitful discussions with J.M. Dudley and D. Brunner are
gratefully acknowledged. We thank the EPOCH support team for help
(https://cfsa-pmw.warwick.ac.uk), and French RENATECH network. The authors
acknowledge the financial supports of: European Research Council (ERC)
682032-PULSAR, Region Bourgogne-Franche-Comte and Agence Nationale de la
Recherche (EQUIPEX+ SMARTLIGHT platform ANR-21-ESRE-0040), Labex ACTION
ANR-11-LABX-0001-01, I-SITE BFC project (contract ANR-15-IDEX-0003), and the
EIPHI Graduate School ANR-17-EURE-0002. This work was granted access to the
PRACE HPC resources at CINECA, Italy, under the Project ”PULSARPIC”
(PRA19_4980), PRACE HPC resource at TGCC, France under the Project ”PULSARPIC”
(RA5614), HPC resource at TGCC, France under the projects A0070511001 and
A0090511001, and Mésocentre de Calcul de Franche-Comté.

## References

  * Guillot [1999] T. Guillot, Interiors of giant planets inside and outside the solar system, Science 286, 72 (1999).
  * Saiz  _et al._ [2008] E. G. Saiz, G. Gregori, D. O. Gericke, J. Vorberger, B. Barbrel, R. J. Clarke, R. R. Freeman, S. H. Glenzer, F. Y. Khattak, M. Koenig, O. L. Landen, D. Neely, P. Neumayer, M. M. Notley, A. Pelka, D. Price, M. Roth, M. Schollmeier, C. Spindloe, R. L. Weber, L. van Woerkom, K. Wünsch, and D. Riley, Probing warm dense lithium by inelastic X-ray scattering, Nat. Phys. 4, 940 (2008).
  * Vailionis  _et al._ [2011] A. Vailionis, E. G. Gamaly, V. Mizeikis, W. Yang, A. V. Rode, and S. Juodkazis, Evidence of superdense aluminium synthesized by ultrafast microexplosion, Nat. Commun 2, 445 (2011).
  * Smillie  _et al._ [2020] L. A. Smillie, M. Niihori, L. Rapp, B. Haberl, J. S. Williams, J. E. Bradby, C. J. Pickard, and A. V. Rode, Exotic silicon phases synthesized through ultrashort laser-induced microexplosion: Characterization with Raman microspectroscopy, Phys. Rev. Materials 4, 093803 (2020).
  * Kruer [1988] W. Kruer,  _The physics of laser plasma interactions_ , Frontiers in physics (Westview Press, Boulder, Colorado 80301, 1988).
  * Eliezer [2002] S. Eliezer,  _The Interaction of High-Power Lasers with Plasmas_ , Series in Plasma Physics (Institute of Physics Publishing, Bristol BS1 6BE, 2002).
  * Bulgakova  _et al._ [2015] N. M. Bulgakova, V. P. Zhukov, S. V. Sonina, and Y. P. Meshcheryakov, Modification of transparent materials with ultrashort laser pulses: What is energetically and mechanically meaningful?, J. Appl. Phys. 118, 233108 (2015).
  * Naseri  _et al._ [2020] N. Naseri, G. Dupras, and L. Ramunno, Mechanism of laser induced filamentation in dielectrics, Opt. Express 28, 26977 (2020).
  * Durnin  _et al._ [1987] J. Durnin, J. J. Miceli, and J. H. Eberly, Diffraction-free beams, Phys. Rev. Lett. 58, 1499 (1987).
  * Gibbon and Bell [1992] P. Gibbon and A. R. Bell, Collisionless absorption in sharp-edged plasmas, Phys. Rev. Lett. 68, 1535 (1992).
  * Teubner  _et al._ [1993] U. Teubner, J. Bergmann, B. van Wonterghem, F. P. Schäfer, and R. Sauerbrey, Angle-dependent X-ray emission and resonance absorption in a laser-produced plasma generated by a high intensity ultrashort pulse, Phys. Rev. Lett. 70, 794 (1993).
  * Palastro  _et al._ [2018] J. P. Palastro, J. G. Shaw, R. K. Follett, A. Colaïtis, D. Turnbull, A. V. Maximov, V. N. Goncharov, and D. H. Froula, Resonance absorption of a broadband laser pulse, Phys. Plasmas 25, 123104 (2018).
  * Denisov [1957] N. Denisov, On a singularity of the field on an electromagnetic wave propagated in an inhomogeneous plasma, J. Exp. Theor. Phys. 4, 544 (1957).
  * Arber  _et al._ [2015] T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, Contemporary particle-in-cell approach to laser-plasma modelling, Plasma Phys. Control. Fusion 57, 113001 (2015).
  * Alam  _et al._ [2016] M. Z. Alam, I. De Leon, and R. W. Boyd, Large optical nonlinearity of indium tin oxide in its epsilon-near-zero region, Science 352, 795 (2016).
  * Meyer  _et al._ [2019] R. Meyer, L. Froehly, R. Giust, J. D. Hoyo, L. Furfaro, C. Billet, and F. Courvoisier, Extremely high-aspect-ratio ultrafast Bessel beam generation and stealth dicing of multi-millimeter thick glass, Appl. Phys. Lett. 114, 201105 (2019).
  * Ardaneh  _et al._ [2020] K. Ardaneh, R. Giust, B. Morel, and F. Courvoisier, Generation of a Bessel beam in FDTD using a cylindrical antenna, Opt. Express 28, 2895 (2020).
  * [18] See supplemental material at URL for a document with the information about setups experiments and PIC simulations.
  * Bhardwaj  _et al._ [2006] V. R. Bhardwaj, E. Simova, P. P. Rajeev, C. Hnatovsky, R. S. Taylor, D. M. Rayner, and P. B. Corkum, Optically produced arrays of planar nanostructures inside fused silica, Phys. Rev. Lett. 96, 057404 (2006).
  * Liao  _et al._ [2015] Y. Liao, J. Ni, L. Qiao, M. Huang, Y. Bellouard, K. Sugioka, and Y. Cheng, High-fidelity visualization of formation of volume nanogratings in porous glass by femtosecond laser irradiation, Optica 2, 329 (2015).
  * Xie  _et al._ [2015] C. Xie, V. Jukna, C. Milián, R. Giust, I. Ouadghiri-Idrissi, T. Itina, J. M. Dudley, A. Couairon, and F. Courvoisier, Tubular filamentation for laser material processing, Sci. Rep. 5, 8914 (2015).
  * Rapp  _et al._ [2016] L. Rapp, R. Meyer, R. Giust, L. Furfaro, M. Jacquot, P. A. Lacourt, J. M. Dudley, and F. Courvoisier, High aspect ratio micro-explosions in the bulk of sapphire generated by femtosecond Bessel beams, Sci. Rep. 6, 34286 (2016).
  * Boyd and Prato [2008] R. Boyd and D. Prato,  _Nonlinear Optics_ (Academic Press, Cambridge, Massachusetts, 2008).
  * Erokhin  _et al._ [1969] N. Erokhin, V. Zakharov, and S. Moiseev, Second harmonic generation by an electromagnetic wave incident on inhomogeneous plasma, J. Exp. Theor. Phys. 29, 101 (1969).
  * Vlnogradovand and Pustovalov [1973] A. V. Vlnogradovand and V. V. Pustovalov, Second-harmonic generation in an inhomogeneous laser plasma, J. Exp. Theor. Phys. 36, 492 (1973).
  * Auer  _et al._ [1979] G. Auer, K. Sauer, and K. Baumgärtel, Second-harmonic emission at resonance absorption, Phys. Rev. Lett. 42, 1744 (1979).
  * Dragila [1982] R. Dragila, Second-harmonic generation in a laser-produced plasma, Phys. Rev. A 25, 1127 (1982).
  * del Hoyo  _et al._ [2020] J. del Hoyo, R. Meyer, L. Furfaro, and F. Courvoisier, Nanoscale confinement of energy deposition in glass by double ultrafast Bessel pulses, Nanophotonics 10, 1089 (2020).
  * Kauranen and Zayats [2012] M. Kauranen and A. V. Zayats, Nonlinear plasmonics, Nat. Photonics 6, 737 (2012).
  * Kampfrath  _et al._ [2013] T. Kampfrath, K. Tanaka, and K. A. Nelson, Resonant and nonresonant control over matter and light by intense terahertz transients, Nat. Photonics 7, 680 (2013).
  * Han  _et al._ [2016] S. Han, H. Kim, Y. W. Kim, Y.-J. Kim, S. Kim, I.-Y. Park, and S.-W. Kim, High-harmonic generation by field enhanced femtosecond pulses in metal-sapphire nanostructure, Nat. Commun 7, 13105 (2016).
  * Denoeud  _et al._ [2014] A. Denoeud, A. Benuzzi-Mounaix, A. Ravasio, F. Dorchies, P. M. Leguay, J. Gaudin, F. Guyot, E. Brambrink, M. Koenig, S. Le Pape, and S. Mazevet, Metallization of warm dense SiO_2 studied by XANES spectroscopy, Phys. Rev. Lett. 113, 116404 (2014).
  * Ernstorfer  _et al._ [2009] R. Ernstorfer, M. Harb, C. T. Hebeisen, G. Sciaini, T. Dartigalongue, and R. D. Miller, The formation of warm dense matter: Experimental evidence for electronic bond hardening in gold, Science 323, 1033 (2009).

Supplemental material

##  I Particle-In-Cell simulations

###  I.1 Pulse injection

Our numerical simulations model a plasma of electrons and ions in vacuum. We
used ion-to-electron mass ratio $102\times 1836$ as the product of the
sapphire molar mass and the proton mass. We injected from the $z=0$ boundary a
Gaussian pulse which propagates along the positive $z-$direction and is
polarized along the $x-$direction. We applied a phase
$\Phi(r)=-(2\pi/\lambda)r\sin\theta$ to the Gaussian beam, as it is the case
experimentally using the spatial light modulator. The wavelength and cone
angle are $\lambda=0.8$ µm, $25^{\circ}$, respectively. The temporal profile
of the beam is a Gaussian function with a width of $100\,{\rm fs}$ FWHM. The
Gaussian beam waist is $w=10$ µm. All simulations using a reduced cone angle,
to simulate refraction, gave very similar results in terms of absorption,
characteristics of near and far-field patterns. The Bessel beam length has
been reduced to $\approx$ 18 µm, by cropping the Gaussian beam window, to
reduce the simulation run time. The peak intensity in the Bessel zone is
$6\times 10^{14}\,{\rm W/cm^{2}}$ in absence of plasma, i.e., on the same
order of magnitude as in the experiments.

###  I.2 Plasma parameters

We have used circular and elliptical plasma rods. The axial density profile of
the plasma rod follows $\tanh(z/z_{0})$ with $z_{0}=1$ µm. It shows a
$z-$invariant density profile extending from $\approx$ 2 µm away from the
onset of the Bessel beam. The velocity distribution of both plasma species
have been initialized with a Maxwellian distribution at $T_{\rm e}=T_{\rm i}$.
We have run simulations for different initial plasma temperatures of $1\times
10^{-3}$, 1 and 10 eV and we have obtained very similar results in terms of
pulse absorption, structure of the fields, and hot electron population. All
simulations results have included electron-ion collisions, although the
collisionless case is almost indistinguishable.

###  I.3 Resolution

We have used a computational box of $15\times 15\times 30$ µ$\rm m^{3}$ volume
and run the simulations up to $t_{\rm run}=320\,{\rm fs}$. The spatial
resolution of the grid for second-order FDTD is set so that $dxk^{0}_{\rm
r}=dyk^{0}_{\rm r}=0.04$, and $dzk^{0}_{\rm z}=0.1$ where $k^{0}_{\rm
r}=k_{0}\sin\theta$, $k^{0}_{\rm z}=k_{0}\cos\theta$. Using a fourth-order
FDTD solver has not shown any significant change of the results. The temporal
resolution is constrained by the Courant stability condition. We have used the
convolutional perfectly matched layers boundary condition for the fields and
open boundary condition for the particles. After testing the simulations with
different numbers of particles per cell from 32 to 128 to ensure energy
conservation, we have chosen $32$ particles per cell per species (a total
number of 550 million particles). We have used a first-order triangle shape
function with 3 points to represent particles giving third-order particle
weighting into the grid. We have obtained similar results using a third-order
b-spline particle shape with 5 points. In the simulation of laser-plasma
interaction using PIC simulation, the Debye length with typical value in the
range $0.2-2$ nm is not resolved by the grid cell, as it is the case in most
simulations of laser-solid interaction [14]. However, the energy is conserved
using the energy conservative scheme for the integration of particle
trajectory. In the EPOCH PIC code, the numerical heating is reported to grow
as ${\rm d}T_{\rm eV}/{\rm d}t_{\rm ps}=\alpha_{\rm H}n_{23}^{3/2}\Delta
x_{\rm nm}^{2}/n_{\rm ppc}$ where $T_{\rm eV}$ is the temperature in eV,
$t_{\rm ps}$ time in ps, $n_{23}$ density in $10^{23}\,{\rm 1/cm^{3}}$,
$\Delta x_{\rm nm}$ grid cell in nm, $n_{\rm ppc}$ number of particles per
cell, and $\alpha_{\rm H}$ depends on the particle shape and current
smoothing. For the triangular particle shape with current smoothing and the
above resolutions $\alpha_{\rm H}=60$. Therefore, for the chosen resolution
${\rm d}T_{\rm eV}/{\rm d}t\sim 1\,{\rm eV/ps}$.

###  I.4 Far-field distributions for first and second harmonics

We have obtained the intensity spectrum $I(\omega,k_{\rm x},k_{\rm y})$ by
performing a discrete Fourier transform on $B_{\rm x:y:z}(t,x,y)$ at a fixed
propagation distance of $z=20$ µm. The Fourier transformation has been
performed on the recorded data for $|x|\leqslant 5$ µm and $|y|\leqslant 5$
µm. This is the reason for the width of the simulation far-fields in Figs.
(S.1). The first and second harmonic radiation spectra are then obtained by
filtering the power spectrum at the central frequency of the pulse
$\omega_{0}$, and $2\omega_{0}$, respectively.

##  II Experimental setup

###  II.1 Fluence distributions

We spatially shaped 100 fs laser pulses at a central wavelength of 800 nm.
Using a spatial light modulator combined to a 2f-2f telecentric arrangement to
produce a Bessel beam with cone angle 25∘ in air. The pulse duration of 100 fs
has been characterized at the sample position. After the interaction, the
pulse is collected with a $\times$100 microscope objective (Olympus MPLFLN) in
a 2f-2f arrangement to image the pulse on a camera (Stingray F146B, 14 bits).
The imaging system, with a numerical aperture of 0.9 is placed on an
independent translation stage. The near-field fluence distributions of Fig.
(3) have been obtained by interrupting the propagation of the Bessel beam in
the sample at the exit surface, as described in detail in reference [21]. The
evolution of energy along the propagation direction has been measured by
spatially integrating on the camera the far-field fluence at each propagation
distance. The overall absorption in the sample has been confirmed by measuring
the pulse energy on large-area photodiodes before and after the sample, as in
reference [28].

###  II.2 Imaging of second harmonic generation

Far-field images were recorded on a camera using an accumulation over 500
shots when the sample is continuously translated to separate the laser impacts
by 5 µm. The field was collected using a $\times$50 microscope objective
(Olympus MPLFLN). The far-field distributions of Figs. 4(a) and 4(b) have been
symmetrized with respect to the $y-$axis. The measurement has been performed
with the Bessel beam fully inside the bulk of sapphire. We remark that the
second harmonic signal disappears when the input pulse is stretched to ps
pulse duration. The signal is reduced when part of the beam crosses one of the
air-dielectric interfaces, since the plasma length is reduced. The
normalization of the axes of experimental measurements of Fig. S.1 (a) and (c)
has been performed using the pupil radius as the maximal numerical aperture
(NA=0.8).

Figure S.1: Far-field fluence distributions for the first and second harmonic.
Far-field distribution for first harmonic from (a) experiment and (b)
numerical simulation using the same parameters as in Fig. 3. (c) Experimental
far-field distribution of the emission in the spectral range $400\pm 20$ nm
and (d) simulated far-field distribution for the second harmonic. All axes are
normalized by the radial component of the Bessel pulse, $k_{\rm
r}^{0}=2\pi\sin\theta/\lambda$.