# Deep Stochastic Optimization in Finance

A. Max Reppen111Questrom School of Business, Boston University, Boston, MA,
02215, USA, email<EMAIL_ADDRESS>H. Mete Soner222Department of Operations
Research and Financial Engineering, Princeton University, Princeton, NJ,
08540, USA, email<EMAIL_ADDRESS>Valentin Tissot-Daguette333Department
of Operations Research and Financial Engineering, Princeton University,
Princeton, NJ, 08540, USA, email<EMAIL_ADDRESS>
###### Abstract

This paper outlines, and through stylized examples evaluates a novel and
highly effective computational technique in quantitative finance. Empirical
Risk Minimization (ERM) and neural networks are key to this approach. Powerful
open source optimization libraries allow for efficient implementations of this
algorithm making it viable in high-dimensional structures. The free-boundary
problems related to American and Bermudan options showcase both the power and
the potential difficulties that specific applications may face. The impact of
the size of the training data is studied in a simplified Merton type problem.
The classical option hedging problem exemplifies the need of market generators
or large number of simulations.

Key words: ERM, Neural Networks, Hedging, American Options.
Mathematics Subject Classification: 91G60, 49N35, 65C05.

##  1 Introduction

Readily available and effective optimization libraries such as Tensorflow or
Pytorch now make previously intractable regression type of algorithms over
hypothesis spaces with large number of parameters computationally feasible. In
the context of stochastic optimal control and nonlinear parabolic partial
differential equations which have such representations, these exciting
advances allow for a highly efficient computational method. This algorithm,
which we call _deep empirical risk minimization_ , proposed by E & Han [21]
and E, Jentzen & Han [22], uses artificial neural networks to approximate the
feedback actions which are then trained by empirical risk minimization. As
stochastic optimal control is the unifying umbrella for almost all hedging,
portfolio or risk management problems, and many models in financial economics,
this method is also highly relevant for quantitative finance.

Although artificial neural networks as approximate controls are widely used in
optimal control and reinforcement learning [6], deep empirical risk
minimization simulates directly the system dynamics and does not necessarily
use dynamic programming. It aims to construct optimal actions and values
offline by using the assumed dynamics and the rewards structure, and often
uses market generators to simulate large training data sets. This key
difference between reinforcement learning and the proposed algorithm ushers in
essential changes to their implementations and analysis as well.

Our goal is to outline this demonstrably effective methodology, assess its
strengths and potential shortcomings, and also showcase its power through
representative examples from finance. As verified in its numerous
applications, deep empirical risk minimization is algorithmically quite
flexible and handles well a large class of high-dimensional models, even non-
Markovian ones, and adapts to complex structures with ease. To further
illustrate and evaluate these properties, we also study three classical
problems of finance with this approach. Additional examples from nonlinear
partial differential equations and stochastic optimal control are given in the
recent survey articles of Fecamp, Mikael & Warin [15] and Germain, Pham &
Warin [19]. They also provide an exhaustive literature review.

Our first class of examples is the American and Bermudan options. The analysis
of these instruments offer many-faceted complex experiments through which one
appreciates the potentials and the challenges. In a series of papers, Becker
_et al._ , [4, 5] bring forth a complete analysis with computable theoretical
upper bounds through its known convex dual. They also obtain inspiring
computational results in high dimensional problems such as Bermudan max-call
options with $500$ underlyings. Akin to deep empirical risk minimization is
the seminal regression on Monte-Carlo methods that were developed for the
American options by Longstaff & Schwartz [29] and Tsitsiklis & van Roy [38].
Many of their refinements, as delineated in the recent article of Ludkovski
[30], make them not only textbook topics but also standard industrial tools.
Still, the deep empirical risk minimization approach to optimal stopping has
some advantages over them, including its effortless ability to incorporate
market details and frictions, and to operate in high-dimensions as caused by
state enlargements needed for path-dependent claims. An example of the latter
is the American options with rough volatility models as studied by Chevalier
_at. al._ [11]. They require infinite-dimensional spaces and their numerical
analysis is given in Bayer _et al._ [3]. Other similar examples can be found
in [4, 5].

For interpretability of our results, we base the stopping decisions on a
surface separating the ‘continuation’ and ‘stopping’ regions, and approximate
directly this boundary - often called the _free boundary_ \- by an artificial
neural network. Similarly for the same reason, Ciocan & Mis̆ic [12] compute
the free boundary directly, by using tree based methods. An additional benefit
of this geometric approach to American options is to construct a tool that can
also be effectively used for financial problems with discontinuous decisions
such as regime-switching or transaction costs, as well as non-financial
applications. Indeed, the computation of the free-boundary is an interesting
problem independent of applications to finance. Recently, deep Galerkin method
[37] is used to compute the free boundary arising in the classical Stephan
problem of melting ice [39].

Our numerical results, reported in the subsections 4.5 and 4.6 below, show
that natural problem specific modifications enable the general approach to
yield excellent results comparable to the ones achieved in [4, 5]. The free
boundaries that we compute for the two-dimensional max-call options also
compare to the results by Broadie & Detemple [8] and by Detemple [14]. An
important step in our approach is to replace the stopping rule given by the
sharp interface by a relaxed stopping rule given by a _fuzzy boundary_ as
described in the subsection 4.4. Further analysis and the results of our free-
boundary methodology are given in our future manuscript [32].

Our second example of classical quadratic hedging [36] is undoubtedly one of
the most compelling benchmark for any computational technique in quantitative
finance. Thus, the evaluation of the deep empirical risk minimization
algorithm on this problem, imparts valuable insights. Readily, Bühler _et al._
, [9, 10] use this approach for multidimensional Heston type models,
delivering convincing evidence for the flexibility and the scope of the
algorithm, particularly in high-dimensions. Huré _et al._ , [25] and Bachouch
_et al._ , [2] also obtain equally remarkable results for the stochastic
optimal control using empirical minimization as well as other hybrid
algorithms partially based on dynamic programming. Extensive numerical
experimentations are also carried out by Fecamp _et. al_. [15] in an
incomplete market that models the electricity markets containing a non-
tradable volume risk [40]. Ruf & Wang [33] apply this approach to market data
of S&P 500 and Euro Stoxx 50 options. In all these applications, variants of
the quadratic hedging error is used as the loss function.

To highlight the essential features, we focus on a simple frictionless market
with Heston dynamics, and consider a vanilla Call option with quadratic loss.
In this setting, we analyze both the pure hedging problem by fixing the price
at a level lower than its known value and also the pricing and hedging problem
by training for the price as well. By the well-known results of Schweizer [35,
36] and Föllmer & Schweizer [18], we know that the minimizer of the analytical
problem in the continuous time is equal to the price obtained by Heston [23]
as the discounted expected value under the risk neutral measure with the
chosen market risk of volatility risk. Our numerical computations verify these
results as well.

As the final example, we report the results of an accompanying paper of the
first two authors [31] for a stylized Merton type problem. With simulated
data, the numerical results once again showcase the flexibility and the scope
of the algorithm, in this problem as well. We also observe that in data-poor
environments, the artificial neural networks have an amazing capability to
over-learn the data causing poor generalization. This is one of the key
results of [31] which was also observed in [28]. Despite this potential, as
demonstrated by our experiments, continual data simulation can overcome this
difficulty swiftly.

In this paper, we only discuss the properties of the algorithms that are
variants of the deep empirical risk minimization. The use of artificial neural
networks or statistical machine learning is of course not limited to this
approach. Indeed, starting from [26] and especially recently, artificial
neural networks have been extensively employed in quantitative finance. In
particular, kernel methods are applied to portfolio valuation in [7], and to
the density estimation in [16]. Gonon _et. al_. [20] use the methodology to
study an equilibrium problem in a market with frictions. For further results
and more information, we refer to the recent survey of Ruff & Wang [34] and
the references therein.

The paper is organized as follows. The next section formulates the control
problem abstractly covering many important financial applications. The
description of the algorithm follows. Section 4 is about the American and
Bermudan options. The quadratic hedging problem is the topic of Section 5.
Finally, the numerical examples related to the simple Merton problem are
discussed in Section (6).

Acknowledgements. Research of the second and the third authors was partially
supported by the National Science Foundation grant DMS 2106462.

##  2 Abstract problem

Following the formulation of [31], we start with a
${\mathcal{Z}}\subset\mathbb{R}^{d}$ valued stochastic process $Z$ on a
probability space $\Omega$. This process drives the dynamics of the problem,
and in all financial examples that we consider it is the related to the stock
returns. For that reason, in the sequel, we refer to $Z$ as the _returns
process_ , although they may be logarithmic returns in some cases. Investment
or hedging decisions are made at $N$ uniformly spaced discrete time points
labeled by $k=0,1,\ldots,N$ and let

${\mathcal{T}}:=\\{0,1,\ldots,N-1\\},\qquad\widehat{{\mathcal{T}}}:=\\{0,1,\ldots,N\\}.$

We use the notation $Z=(Z_{1},\ldots,Z_{N})$ and set $Z_{0}=0$. We further let
$\mathbb{F}=({\mathcal{F}}_{t})_{t=0,\ldots,N}$ be the filtration generated by
the process $Z$. The $\mathbb{F}$-adapted controlled state process $X$ takes
values in another Euclidean space ${\mathcal{X}}$ and it may include all or
some components of the uncontrolled returns process $Z$.

In the financial examples, the state includes the marked-to-market value of
the portfolio and maybe other relevant quantities. In a path-dependent
structure, we would be forced to include not only the current value of the
portfolio and the return, but also some past values as well (theoretically, we
need to keep all past values but in practice one stops at a finite point). In
illiquid markets, the portfolio composition is also included into the state
and even the order-book might be considered. We assume that the state is
appropriately chosen so that _the relevant decisions are feedback functions of
the state alone_ and we optimize over feedback decisions or controls. Thus,
even if the original problem is low dimensional but non-Markov, one is forced
to expand the state resulting in a high-dimensional problem.

We denote the set of possible actions or decisions by ${\mathcal{A}}$. While
the main decision variable is the portfolio composition, several other
quantities such as the speed of the change of the portfolio could be included.
Then, a feedback decision is a continuous function

$\pi:{\mathcal{T}}\times{\mathcal{X}}\mapsto{\mathcal{A}}.$

We let ${\mathcal{C}}$ be the set of all such functions. Given
$\pi\in{\mathcal{C}}$, the time evolution of the state vector is then
completely described as a function of the returns process $Z$. Hence, all
optimization problems that we consider have the following form,

$\text{minimize}\quad
v(\pi):=\mathbb{E}\left[\,\ell(\pi,Z)\,\right],\qquad\text{over all}\
\pi\in{\mathcal{C}},$

where $\ell$ is a nonlinear function. We refer the reader to [31] for a
detailed derivation of the above formulation and several examples. Although
the cost function $\ell$ could be quite complex to express analytically, it
can be easily evaluated by simply mimicking the dynamics of the financial
market. Hence, computationally they are straight-forward to compute and all
details of the markets can be easily coded into it.

The goal is to compute the optimal feedback decision, $\pi^{*}$, and the
optimal value $v^{*}$,

$\pi^{*}\in\text{argmin}_{\pi\in{\mathcal{C}}}\,v(\pi),\qquad
v^{*}:=\inf_{\pi\in{\mathcal{C}}}\,v(\pi)\ =\ v(\pi^{*}).$

When the underlying dynamics is Markovian and the cost functional has an
additive structure, the above formulation of optimization over feedback
controls is equivalent to the standard formulation which considers the larger
class of all adapted processes, sometimes called _open loop_ controls [17].
However, even without this equivalence, the minimization over the smaller
class of feedback controls is a consistent and a well-defined problem, and due
to their tractability, feedback controls are widely used. In this manuscript,
we implicitly assume that the problem is well chosen and the goal is to
construct the best feedback control.

##  3 The algorithm

In this section, we describe the _deep empirical minimization_ algorithm
proposed by _Weinan E, Jiequn Han_ , and _Arnulf Jentzen_ in [21, 22].

A _batch_ $B:=\\{Z^{1},\ldots,Z^{m}\\}$, with a size of $m$, is an i.i.d.
realization of the returns process $Z$, where
$Z^{i}=(Z^{i}_{1},\ldots,Z^{i}_{N})$ for each $i$. We set

$L(\pi,B):=\frac{1}{m}\,\sum_{i=1}^{m}\,\ell(\pi,Z^{i}),$

and consider a set of artificial neural networks parametrized by,

${\mathcal{N}}=\left\\{\,\Phi(\cdot;\theta):{\mathcal{T}}\times{\mathcal{X}}\mapsto{\mathcal{A}}\
\ :\ \theta\in\Theta\ \right\\}\,\subset\,{\mathcal{C}}.$

Instead of searching for a minimizer in ${\mathcal{C}}$, we look for a
computable solution in the smaller set ${\mathcal{N}}$. That is, numerically
we approximate the following quantities:

$\displaystyle\theta^{*}$
$\displaystyle:=\theta^{*}_{\mathcal{N}}\in\text{argmin}_{\theta\in\Theta}\,v(\Phi(\cdot;\theta)),$
$\displaystyle v_{{\mathcal{N}}}$
$\displaystyle:=\inf_{\theta\in\Theta}\,v(\Phi(\cdot;\theta))=v(\Phi(\cdot;\theta^{*})).$

The classical universal approximation result for artificial neural networks
[13, 24] imply, under some natural structural assumptions on the function
$\ell$, that $v_{{\mathcal{N}}}$ approximates $v^{*}$ as the networks gets
larger as proved in [31](Theorem 5). This also implies that the performance of
the trained feedback control $\Phi(\cdot;\theta^{*})$ is almost optimal.

The pseudocode of the algorithm to compute $\theta^{*}$ and $v^{*}$ is the
following,

  * •

_Initialize_ $\theta\in\Theta$;

  * •

_Optimize by stochastic gradient descent:_ for $n=0,1,\ldots$:

    * –

Generate a batch $B:=\\{Z^{1},\ldots,Z^{m}\\}$,

    * –

Compute the derivative $d:=\nabla_{\theta}\,L(\Phi(\cdot;\theta),B)$;

    * –

_Update_ $\theta\,\leftarrow\,\theta-\kappa d$.

  * •

_Stop_ if $n$ is sufficiently large and the improvement of the value is
‘small’.

In the above $\kappa$ is the learning rate and the stochastic gradient step is
done through an optimization library.

The data generation can be done through either an assumed and calibrated
model, namely a market generator, or by random samples from a fixed financial
market data when sufficient and relevant historical data is available.
Although these two settings look similar, one may get quite different results
in these two cases, even when the fixed data set is large. One of our goals is
to better understand this dichotomy between these two data regimes and the
size of the data needed for reliable results. Theoretically, when the
simulation capability is not limited and data is continually generated, the
above algorithm should yield the desired minimizer $\theta^{*}$ and the
corresponding optimal feedback decision $\Phi(\cdot,\theta^{*})$. However,
with a fixed data set, the global minimum over ${\mathcal{N}}$ is almost
always strictly less than $v^{*}$, and the large enough networks will
eventually gravitate towards this undesirable extreme point which would be
over-learning the data as already observed and demonstrated in [31].

##  4 Exercise boundary of American type options

American and Bermudan options are particularly central to any computational
study in quantitative finance as they pose difficult and deep challenges, and
they serve as an important benchmark for any new numerical approach. Methods
successful in this setting often generalize to other problems as well. Indeed,
the seminal regression on Monte-Carlo methods that were developed for the
American options by Longstaff & Schwartz [29] and Tsitsiklis & van Roy [38]
have not only become industry standards in few years, but they have also shed
insight into other problems as well. Together with rich improvements developed
over the past decades, they can now handle many Markovian problem with ease.
However, the key feature of these algorithms is a projection onto a linear
subspace, and this space must grow exponentially with the dimension of the
ambient space, making high-dimensional problems out of reach of this otherwise
powerful technique. Examples of such high-dimensional problems are financial
instruments on many underlyings modeled with many parameters, path-dependent
options, or non-Markovian models, all requiring state enlargements and
resulting in vast state spaces.

###  4.1 Problem

As well known the problem is to decide when to stop and collect the pay-off of
a financial contract. Mathematically, for
$t\in\widehat{{\mathcal{T}}}=\\{0,\ldots,N\\}$, let
$S_{t}\in\mathbb{R}_{+}^{d}$ be the stock value at the $t$-th trading date and
$\varphi:\mathbb{R}_{+}^{d}\mapsto\mathbb{R}$ be the pay-off function. With a
given interest rate $r>0$, the problem is

$\text{maximize}\ \
v(\tau):=\mathbb{E}\left[\,e^{-r\tau}\,\varphi(S_{\tau})\,\right],$

over all $\widehat{{\mathcal{T}}}$-valued stopping times $\tau$. We use the
filtration generated by the stock price process to define the stopping times.
It is classical that the expectation is taken under the risk neutral measure.

We assume that $S$ is Markov and the pay-off is a function of the current
stock value. When it is not, then we need to enlarge the state space. In
factor models like Heston or SABR, factor process is included. In non-
Markovian models like the fractional Brownian motion, past values the stock
are added as in [3, 4, 5]. In look-back type options, the minimum or the
maximum of the stock process must be included in the state. We refer to [32]
for the details of these extensions.

We continue by defining the price at all future points. Recall that the
filtration $\mathbb{F}$ is generated by the stock price process. Let $\Xi_{t}$
be the set of all $\mathbb{F}$-stopping times with values in
$\\{t,\ldots,N\\}$. At any $t\in\widehat{{\mathcal{T}}}$,
$s\in\mathbb{R}_{+}^{d}$, let $v(t,s)$ be the maximum value or the price of
this option when $S_{t}=s$, i.e.,

$v(t,s):=\max_{\tau\in\Xi_{t}}\
\mathbb{E}[\,e^{-r(\tau-t)}\,\varphi(S_{\tau})\ \mid\ S_{t}=s\,].$

Then, $v(N,\cdot)=\varphi$ and the the _stopping region_ is given by

${\mathcal{S}}:=\left\\{\,(t,s)\ :\ v(t,s)=\varphi(s)\,\right\\}.$ (4.1)

Then the optimal stopping time is the first time to enter the region
${\mathcal{S}}$, i.e., the following stopping time in $\Xi_{t}$ is a maximizer
of the above problem:

$\tau^{*}:=\min\\{\ u\in\\{t,\ldots,N\\}\ :\ (u,S_{u})\in{\mathcal{S}}\ \\}.$

Notice that as $v(N,\cdot)=\varphi$, we always have
$(N,S_{N})\in{\mathcal{S}}$. This implies that $\tau^{*}$ is well-defined and
is bounded by $N$.

Clearly, standard call or put options are the main examples. Many other
examples that are also covered in the above abstract setting, including the
max-call option discussed below.

######  Example 4.1 (Max-Call).

Let $S_{t}=(S^{(1)}_{t},\ldots,S^{(d)}_{t})\in\mathbb{R}_{+}^{d}$ be a the
stock process of $d$ many dividend bearing stocks. We model it by a
$d$-dimensional geometric Brownian motion with constant mean-return rate and a
covariance matrix. The pay-off of the max-call is given by,

$\varphi(S_{t})=\left(\,\max_{i=1,\ldots,d}\,S^{(i)}_{t}\,-K\,\right)^{+},$

where the strike $K$ is a given constant. We study this example numerically in
subsection 4.6 below. One can also consider max-call options with factor
models with an extended state-space.

###  4.2 Relaxed stopping

Quite recently, in a series of papers, Becker _et al._ , [4, 5] use deep
empirical risk minimization in this context. As the control variable is
discrete (i.e., at any point in space, the decision is either ‘stop’ or ‘go’)
and as the training or optimization is done through a stochastic gradient
method, one has to relax the problem before applying the general procedure. We
continue by first outlining this relaxation.

In the relaxed version, we consider an adapted control process
$p=(p_{0},\ldots,p_{N})$ with values in $[0,1]$ which is the probability of
stopping at that time conditioned on the event that the process has not
stopped before $t$. Because one has to stop at maturity, we have $p_{N}=1$.
Given the process $p$, let $\xi_{t}^{p}$ be the probability of stopping
strictly before $t$. Clearly, $\xi_{0}=0$ and at other times it is defined
recursively by,

$\xi^{p}_{t+1}=\xi^{p}_{t}+p_{t}(1-\xi^{p}_{t})\,=\,p_{t}+(1-p_{t})\,\xi^{p}_{t},\qquad
t\in{\mathcal{T}}.$

It is immediate that $\xi^{p}_{t}\in[0,1]$ and is non-decreasing. Also, if
$p_{t}=1$, then $\xi^{p}_{s}=1$ for all $s>t$. The quantity $(1-\xi^{p}_{t})$
is the unused “stopping budget”, and the relaxed stopping problem is defined
by,

$\text{maximize}\ \
v_{r}(p):=\mathbb{E}\left[\sum_{t=0}^{N}\,p_{t}\,(1-\xi^{p}_{t})\,e^{rt}\,\varphi(S_{t})\,\right],$
(4.2)

over all $[0,1]$-valued, adapted processes $p$. The original problem of
stopping is included in the relaxed one, as for any given stopping time
$\tau$, $p^{\tau}_{t}:=\chi_{\\{t=\tau\\}}$ yields
$\xi^{\tau}_{t}=\chi_{\\{t>\tau\\}}$ and consequently,
$v(\tau)=v_{r}(p^{\tau})$. It is also known that this relaxation does not
change the value.

Becker _et al._ , [4, 5] study the problem through this relaxation and
implement the deep empirical risk minimization exactly as described in the
earlier section. Additionally, using the known convex dual of the stopping
problem, they are able to obtain _computable_ upper-bounds. For many financial
products of interest, they obtain remarkable results in very high-dimensions.
They also consider a fractional Brownian motion model for the stock price. As
for this example there is no Markovian structure, in their calculations the
state is all the past yielding an enormous state space. Still the algorithm is
tractable with computable guarantees.

###  4.3 The free boundary

In most examples, the optimal stopping rule is derived from a surface called
the _free boundary_. For instance, the continuation region of a one-
dimensional American Put option is the epigraph of a function of time. The
stopping region of an American max-call option on the other hand, is obtained
by comparing the maximum of the stock values to a scalar-valued function as
proved in Proposition 4.4 below. These stopping rules have the advantage of
being _interpretable_ [12] and easy to implement. Additionally, free-boundary
problems of this type appear often in financial economics as well as problems
from other disciplines. Thus numerical methods developed for the free-boundary
of an American option could have implications elsewhere as well.

To be able to apply this method, we assume that the stopping region
${\mathcal{S}}$ has a certain structure. Namely, we assume that there exists
two functions

$\alpha:\mathbb{R}_{+}^{d}\mapsto\mathbb{R},\qquad\text{and}\qquad
F:\widehat{{\mathcal{T}}}\times\mathbb{R}_{+}^{d}\mapsto\mathbb{R},$

(recall that $\widehat{{\mathcal{T}}}=\\{0,\ldots,N\\}$) so that the stopping
region of (4.1) is given by,

${\mathcal{S}}=\\{\,(t,s)\ :\ \alpha(s)\leq F(t,s)\,\\}.$

More importantly, we also assume that $\alpha$ is given by the problem and we
only need to determine $F$ which we call the _free boundary_. The following
examples clarifies this assumption which holds in a large class of problems.

######  Example 4.2.

It is known that the stopping region of an American Put option with a
Markovian stock process is given by

${\mathcal{S}}=\\{(t,s)\ :\ s\leq f(t)\ \\},$

for some function $f:[0,T]\mapsto\mathbb{R}_{+}$. In this case, $\alpha(s)=s$
and $F(t,s)=f(t)$.

In the case of the max-call option, we show in Proposition 4.4 below that for
any $s=(s_{1},\ldots,s_{d})\in R_{+}^{d}$ with
$\alpha(s)=\max\\{s_{1},\ldots,s_{d}\\}$, there exists a free boundary $F$. ∎

Given the above structure of the stopping region through the pair $(\alpha,F)$
the optimal stopping time is given by $\tau^{*}=\tau_{F}$, where for any free
boundary $F$,

$\tau_{F}\ :=\ \min\ \\{\ t\in\widehat{{\mathcal{T}}}:\ \alpha(S_{t})\leq
F(t,S_{t})\ \\}.$

In this approach, the output of the artificial neural network is a scalar
valued function $\Phi(\cdot;\theta)$ of time and the state values, and it
approximates the free boundary $F$. Then for any parameter $\theta$, the
stopping time is

$\tau_{\theta}:=\tau_{\Phi(\cdot;\theta)}=\ \min\\{\
t\in\widehat{{\mathcal{T}}}\ :\ \alpha(S_{t})\leq\Phi(t,S_{t}\,;\,\theta)\
\\}.$

###  4.4 Fuzzy boundary

A sharp free-boundary has the same problem of zero-gradients as the original
problem and its remedy is again a relaxation to allow for partial stopping.
Indeed, given a free-boundary $\Phi(\cdot;\theta)$ and a tuning-parameter
$\epsilon>0$, we define a fuzzy boundary region given by,

$F_{\Phi,\epsilon}:=\\{\ (t,s)\
:\;-\epsilon\leq\Phi(t,s;\theta)-\alpha(s)\leq\epsilon\,\\}.$

If $\Phi-\alpha\leq-\epsilon$ we stop, and if $\Phi-\alpha\geq\epsilon$ we
continue, and we do these with probability one in each case. But if the
process falls into the fuzzy region $F_{\Phi,\epsilon}$, then as in the
relaxed problem, we assign a stopping probability as a function of the
normalized distance $d_{t}^{\theta}$ to the sharp boundary
$\\{\Phi-\alpha=0\\}$, i.e.,

$p^{\theta}_{t}:=g(d_{t}^{\theta}),\qquad\text{where}\qquad
d_{t}^{\theta}=\frac{\Phi(t,S_{t};\theta)-\alpha(S_{t})}{\epsilon},$

and $g:[-1,1]\mapsto[0,1]$ is a fixed increasing, onto function. Linear or
sigmoid-like functions are the obvious choices. Once we compute the process
$p^{\theta}_{t}$, the value corresponding to the parameter $\theta$ is
$v_{r}(p^{\theta})$ with $v_{r}$ as in (4.2). Hence, the relaxed free boundary
problem is to train the network to

$\text{minimize}\ \ \theta\in\Theta\ \mapsto\ v_{r}(p^{\theta}).$

The resulting trained artificial neural network is an approximation of the
optimal free boundary.

###  4.5 American Put in one-dimension

As in [4, 5] we run the algorithm for an American put on a non-dividend paying
stock whose price process is modeled by a standard geometric Brownian motion
with parameters

$S_{0}=K=40,T=1,\sigma=0.4,r=0.06,$

where as usual $S_{0}$ is the initial stock value, $K$ is the strike, $\sigma$
is the volatility, and $r$ is the risk-free rate. In this example, the state
process is simply the stock process.

We are able to obtain accurate results for the value as well as for the free
boundary. One typical result is given in Figure 4.5 below. As the free
boundary has a large derivative near maturity, we use a non-uniform mesh near
maturity. Figure 4.5 uses 500 time points. We also employ important sampling
to ensure more crossings of the free boundary. After the training is
completed, the value corresponding to this trained free boundary is computed
by using the corresponding sharp interface. Accurate price values are obtained
rather easily. All of these calculations are implemented by python in a
personal laptop.

Figure 1: Left figure is a random initialization and the right one is the
final trained boundary. Dashed line is the optimal calculated through a
finite-difference scheme. The price is $5.311$.

###  4.6 Max-Call Options

In this subsection, we consider the max-call option studied in the seminal
paper by Broadie & Detemple [8] and also in the book by Detemple [14]. Let
$S_{t}\in\mathbb{R}^{d}$ be the price process of a dividend bearing stock. The
pay-off the max-call option at time $\tau$ is

$\varphi(S_{\tau})=\left(\,m(S_{\tau})\,-\,K\,\right)^{+},\qquad$

where the function $m:\mathbb{R}_{+}^{d}\mapsto\mathbb{R}_{+}$ is given by,

$m(s):=\max_{i=1,\ldots,d}\,\,s_{i},\qquad
s=(s_{1},\ldots,s_{d})\in\mathbb{R}_{+}^{d}.$

The main structural assumption needed is the natural sub-linear dependence of
the stock prices on their initial values.

######  Assumption 4.3 (Sublinearity).

For $t\in{\mathcal{T}}$, $s\in\mathbb{R}_{+}^{d}$, non-decreasing function
$\phi:\mathbb{R}_{+}^{d}\mapsto\mathbb{R}$, $\lambda\geq 1$ and a stopping
time $\tau\geq t$,

$\mathbb{E}\left[\phi(S_{\tau})\mid S_{t}=\lambda
s\right]\,\leq\,\mathbb{E}\left[\phi(\lambda\,S_{\tau})\mid S_{t}=s\right].$

Above assumption is satisfied in all examples. In fact, in most models the
dependency on the initial data is linear. Although in our numerical
calculations, we use a geometric Brownian motion model for the stock price
process, the method also applies more generally to all factor models.

We use this assumption to show that the stopping region has a certain
geometric structure which we exploit. The following result is already proved
in [8] and more generaly in [32]. We provide its proof for completeness. Let
${\mathcal{S}}$ be as in (4.1) and set

${\mathcal{K}}:=\left\\{\,s\in\mathbb{R}_{+}^{d}\ :\ m(s)=1\,\right\\}.$

Note that for any $s\in\mathbb{R}_{+}^{d}$, $\frac{s}{m(s)}\in{\mathcal{K}}$.

######  Proposition 4.4.

Consider the max-call option in a market satisfying the Assumption 4.3. Then,
if $(t,s)\in{\mathcal{S}}$, then $(t,\lambda s)\in{\mathcal{S}}$ for any
$\lambda\geq 1$. In particular,

${\mathcal{S}}=\\{\ (t,s)\ :\ m(s)\geq F(t,s/m(s))\ \\},$

where $F:\widehat{{\mathcal{T}}}\times{\mathcal{K}}\mapsto\mathbb{R}_{+}$ is
given by,

$F(t,k):=\inf\left\\{\,\rho>0\ :\ (t,\rho k)\in{\mathcal{S}}\ \right\\},\qquad
m\in{\mathcal{M}}.$

Above result can be equivalently stated as the $t$-section $S_{t}:=\\{\
s\in\mathbb{R}_{+}^{d}\ :\ (t,s)\in{\mathcal{S}}\ \\}$ of the continuation
region being _star-shaped_ for every $t$.

###### Proof.

Suppose that $(t,s)\in{\mathcal{S}}$ and $\lambda\geq 1$. As $\\{(N,s)\ :\
s\in\mathbb{R}^{d}_{+}\\}\subset{\mathcal{S}}$, if $t=N$, clearly $(t,\lambda
s)=(N,\lambda s)\in{\mathcal{S}}$. So we assume that $t<N$. Then, a point
$(t,s^{\prime})$ is in ${\mathcal{S}}$ if and only if $s^{\prime}>K$ and the
following inequality is satisfied for every $\tau\in\Xi_{t}$:

$\mathbb{E}\left[\,e^{-r(\tau-t)}\,(S_{\tau}-K)^{+}\,|\,S_{t}=s^{\prime}\,\right]\leq
s^{\prime}-K.$

By Assumption 4.3,

$\displaystyle\mathbb{E}\left[\,e^{-r(\tau-t)}\,(S_{\tau}-K)^{+}\,|\,S_{t}=\lambda
s\,\right]$ $\displaystyle\leq\mathbb{E}\left[\,e^{-r(\tau-t)}\,(\lambda
S_{\tau}-K)^{+}\,|\,S_{t}=s\,\right]$
$\displaystyle=\mathbb{E}\left[\,e^{-r(\tau-t)}\,(\lambda[S_{\tau}-K]+(\lambda-1)K)^{+}\,|\,S_{t}=s\,\right]$
$\displaystyle\leq\lambda\,\mathbb{E}\left[\,e^{-r(\tau-t)}\,(S_{\tau}-K)^{+}\,|\,S_{t}=s\,\right]+(\lambda-1)K$
$\displaystyle\leq\lambda(s-K)+(\lambda-1)K$ $\displaystyle=(\lambda s-K).$

Hence, we conclude that $(t,\lambda s)\in{\mathcal{S}}$. ∎

####  4.6.1 Numerical Experiments

We consider a max-Call option and in a geometric Brownian motion model under
the risk neutral measure,

$S_{t}\ =\ S_{0}\,\exp\left((r-div)t+\sigma
W_{t}-\frac{1}{2}\sigma^{2}t\right),$

with parameters

$K=100,\ S_{0}=90,100,110,\ \sigma=0.2,\ r=0.05,\ div=0.1,$

where the notation is as in the previous subsection and $div$ is the dividend
rate. We take the maturity to be $3$ years and $N=9$. Thus, each time interval
corresponds to four months. All these parameters are taken from [4, 5] to
allow for comparison. We also make qualitative comparison to the results of
[8].

Runs | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
Price | 8.0747 | 8.0757 | 8.0710 | 8.0684 | 8.0670 | 8.0731
Stdev | 0.00305 | 0.00315 | 0.00310 | 0.00310 | 0.00309 | 0.00311
Runs | 7 | 8 | 9 | 10 | Mean | Stdev
Price | 8.0686 | 8.0707 | 8.0620 | 8.0679 | 8.0699 | 0.0040
Stdev | 0.00306 | 0.00308 | 0.00307 | 0.00311 | - | -
Table 1: Ten experiments with $S_{0}=90$, batch size $2^{13}$, $7000$
iterations. Prices are calculated with $2^{23}$ simulations. Stdev in the
third and sixth rows refer to the standard deviations of the Monte-Carlo
simulations, while Stdev at the end is the standard deviation of the
calculated ten prices.

Table 1 shows the results with $d=2$, $S_{0}=90$, batch size of $B=2^{13}$ and
$7,000$ iterations. The corresponding price is computed after the training is
completed with $2^{23}$ Monte-Carlo simulations using the sharp boundary
instead of the fuzzy one. Important sampling is used with a $1.4\%$ downward
drift. We repeated the experiment ten times in a personal computer. All of the
results are within the $95\%$ confidence interval $[8.053\,,\,8.082]$ computed
in [1]. The standard deviation of each price computation is quite low. Hence,
the maximum of the values is a lower bound for the price.

We also repeated the experiments of [4, 5] in space dimensions $d=5,10,100$
with the above parameters. For each parameter set, we computed ten prices
exactly as described above. The results reported in Table 2 below are in
agreement with the results of [5] (Table 9). We should also note when $d$ is
large, the maximum of many stocks have a very strong upward drift making the
standard deviation of the rewards quite high.

Dim. | $S_{0}$ | Price | Std | Price in [5] | Max Price | Its Std
---|---|---|---|---|---|---
2 | 90 | 8.0699 | 0.0031 | 8.068 | 8.0757 | 0.0040
2 | 100 | 13.9086 | 0.0059 | 13.901 | 13.9128 | 0.0033
2 | 110 | 21.3434 | 0.0059 | 21.341 | 21.3541 | 0.0104
5 | 90 | 16.6187 | 0.0040 | 16.631 | 16.6238 | 0.0045
5 | 100 | 26.1194 | 0.0259 | 26.147 | 26.1644 | 0.0057
5 | 110 | 36.7176 | 0.0078 | 36.774 | 36.7408 | 0.0078
10 | 90 | 26.2130 | 0.0182 | 26.196 | 26.2362 | 0.0069
10 | 100 | 38.2735 | 0.0538 | 38.272 | 38.3351 | 0.0089
10 | 110 | 50.8350 | 0.0397 | 50.812 | 50.8685 | 0.0081
100 | 90 | 66.2460 | 0.4946 | 66.359 | 66.6163 | 0.0223
100 | 100 | 82.5475 | 0.6463 | 83.390 | 83.6563 | 0.0272
100 | 110 | 98.9868 | 0.0366 | 100.421 | 99.0575 | 0.0353
Table 2: Each price is the mean of ten experiments with parameters as in Table
1. Max price is the maximum of ten experiments with a standard deviation of
the price calculation with $2^{23}$ Monte-Carlo simulations.

The above table reports the average values for ten runs to be able to asses
the possible variations. However, the maximum value among these ten runs is in
fact a lower bound the actual price. As we computed these values with $2^{23}$
(roughly eight million) simulations, the standard devision of this price value
is small.

In two dimensions, the stopping region can be visualized effectively. Figures
2, 4 are stopping regions in two space dimensions obtained with initial data
$S_{0}=90$ and $S_{0}=100$. Clearly the free boundary is independent of the
initial condition and the below numerical results verify it. Also they are
similar to those obtained in [8].

Figure 2: Evolution of the Free Boundary with $S_{0}=90$ Figure 3: Evolution
of the Free Boundary with $S_{0}=100$

##  5 Valuation and Hedging

We consider a European option with stock process $S$ and pay-off
$\varphi(S_{T})$. We consider the Heston dynamics,

$\displaystyle\mathrm{d}S_{t}$
$\displaystyle=S_{t}\,(\mu\mathrm{d}t+\sqrt{v_{t}\,}\,\mathrm{d}W_{t}),$
$\displaystyle\mathrm{d}v_{t}$ $\displaystyle=\left(\kappa(\theta-
v_{t})-\lambda\,v_{t}\right)\,\mathrm{d}t+\sigma\,v_{t}\,\mathrm{d}\tilde{W}_{t},$

where $W,\tilde{W}$ are one-dimensional Brownian motions with constant
correlation of $\rho$, and the five Heston parameters
$(\mu,\kappa,\theta,\sigma,\rho)$ are chosen satisfying the Feller condition.
In particular, we choose the market price of volatility risk parameter
$\lambda$.

Let $p^{*}$ be the price of this claim, and $Z$ be the return process, i.e.,

$Z_{t+1}:=\frac{S_{t+1}-S_{t}}{S_{t}},\qquad t\in{\mathcal{T}}.$ (5.1)

Further, let the feedback actions be the continuous functions

$\pi:{\mathcal{T}}\times\mathbb{R}_{+}\times\mathbb{R}\ \mapsto\ \mathbb{R},$

representing the dollar amount invested in the stock. The corresponding wealth
process is given by,

$X^{\pi,x}_{t+1}=(1+r)\,X^{\pi,x}_{t}+\pi(t,X^{\pi,x}_{t},Z_{t})\,\left(Z_{t+1}-r\,\right),\quad
t\in{\mathcal{T}},$ (5.2)

with initial data $X^{\pi,x}_{0}=x$.

We first fix an initial wealth of $x<p^{*}$ and consider the following _pure-
hedging problem_ of minimizing the square hedging error, i.e.,

$v^{*}(x):=\min_{\pi\in{\mathcal{C}}}\,v(x,\pi),\qquad\text{where}\qquad
v(x,\pi):=\mathbb{E}\left[(\varphi(S_{T})-X^{\pi,x}_{T})^{2}\ \right].$ (5.3)

In the second problem, we minimize over $x$ as well, i.e,

$v^{*}:=\min_{x\in\mathbb{R}}\,v^{*}(x)=\min_{(x,\pi)\in\mathbb{R}\times{\mathcal{C}}}\,v(x,\pi).$
(5.4)

As proved by Föllmer & Schweizer [18], it is well-known that in continuous
time the solution to the second problem, $v^{*}$, is equal to the Heston
price. Thus, for sufficiently fine discretization $v^{*}$ is close to zero,
$x^{*}$ is close to the known continuous-time Heston price. Also the numerical
hedge $\pi^{*}$ must be equal to the continuous time hedge.

If $r=0$, then, $X^{\pi,x}_{t}=x+X^{\pi,0}_{t}$ and the initial wealth $x$
only influences the mean of the hedging error. Therefore, we expect that after
an initial adjustment to minimize the mean, the networks would minimize the
variance which is independent of the initial wealth. This approximate
reasoning indicates that after an initial transient region, both minimization
problems may behave similarly when there is large data.

###  5.1 Numerical results

We implemented the above hedging problem in Julia’s Flux [27] by
parameterizing the portfolio at each time point, including the initial wealth
level. In particular, we hedge a call option with strike $K$, i.e.,
$\varphi(x)=(x-K)_{+}=\max\\{x-K,0\\}$. Our implementation follows the scheme
in Section 3, which we here describe in greater detail for this particular
problem.

We see in (5.4) that the two quantities we optimize over are $x$ and $\pi$. As
$x$ is a scalar, we directly parameterize it with a 1-element tensor, which
after optimization is the option price. The policy $\pi$, however, can be
approximated in various ways. We here opt for a very direct method in which we
represent it by a single neural network with time and stock data as inputs.
This contrasts [9], where the authors discretize time and design one neural
network per time point. As we shall see, our implementation of a single neural
network also performs well, with the additional benefit of allowing changes to
the time discretization during training. There are also other training
differences between the two parameterizations, as, for instance, the one used
here accomplishes a large degree of parameter sharing. Nevertheless, a
thorough account of these differences is outside the scope of the present
paper.

Another detail of our implementation is that we write $\pi$ as a function of
$t$ and $S_{t}$ instead of the formulation in (5.2). It is clear that the two
are mathematically equivalent, although they could differ in training
performance. Ours is a naïve choice and we make it because we find it more
natural, not because it necessarily leads to better performance. The neural
network is designed with two hidden layers of width 20 and with ReLU
activation. In-between layers, batch normalization is employed.444Although we
believe that the following parameters are not crucial for replicating our
results (because they were not tuned), we list them here for completeness:
batch size: 512; optimizer: Adam with the Flux default parameters
$(\eta,\beta_{1},\beta_{2})=(0.001,0.9,0.999)$; and the number of epochs was a
fixed value for which the training error of a typical run had plateaued.

The results of our computations are presented in Table 3. We compare our
numerical solution to the Heston prices from https://www.quantlib.org/. No
significant tuning has gone into producing our values, and it is nevertheless
clear that accurate prices are consistently attained. We see, for instance,
that the absolute error is approximately the same for all three strikes, which
we argue is a consequence of (i) not tuning the training parameters to each
individual problem and (ii) our hedging is in discrete time, which introduces
a time discretization error. Although this is only a one-dimensional problem,
it gives credence to the method’s effectiveness, effectiveness that does
translate into higher-dimensional performance, as we illustrated for the
American options problems.

$K$ | QuantLib price | Price | Avg. abs. error | Error std. dev.
---|---|---|---|---
90 | 10.076508 | 10.078163 | 0.001869 | 0.001174
100 | 2.295405 | 2.295211 | 0.002018 | 0.001065
110 | 0.128136 | 0.127069 | 0.001971 | 0.001793
Table 3:  Hedging performance of a call option with strike $K$ in a Heston
model with parameters $S_{0}=100$, $v_{0}=0.04$, $\kappa=0.9$, $\theta=0.04$,
$r=\lambda=0$ and $\sigma=0.2$. The maturity $T=1/12$ and the time interval is
discretized in 22 steps. Each row lists the deep hedging price average over
100 runs along with the standard deviation over the same 100 runs.

Theoretically, in continuous-time, the optimal hedge is independent of the
initial wealth. We also studied this by fixing the initial portfolio value to
the price and also to its half value. One simulation of the trained hedge is
given in the Figure 4 shows that the dependence is minimal.

Figure 4: Optimal Hedges for the Heston model. The orange curve is the trained
feedback hedge with an initial wealth half the option price, while the red
curve is trained with initial wealth equal to the option price. The stock
price is rescaled to start at $S_{0}=1$.

##  6 Merton Problem and Overlearning

In this section, we summarize the results of [31] by the first two authors. As
in that paper, to emphasize the essential features of the algorithm, a simple
financial market without any frictions and constant coefficients is
considered. Additionally, consumption is not taken into account. All these
details can incorporated into the model and problems with complex market
structures have already been studied extensively by Buehler _et.al._ [9, 10].

Consider a stock price process $S_{t}\in\mathbb{R}_{+}^{d}$ in discrete time
and assume a constant interest rate of $r$. Let the return process $Z$ be as
in (5.1) and $X^{\pi}=X^{\pi,x}$ be as in (5.2). We suppress the dependence of
the initial wealth $x$ for simplicity. Then, the classical investment problem
is to maximize $v(\pi):=\mathbb{E}[U(X^{\pi}_{T})]$ with a given utility
function $U$.

In [31] it is proved that the deep empirical risk minimization algorithm
converges as the size of the training data gets larger. On the other hand, it
is also shown that for fixed training data sets, larger and deeper neural
networks have the capability of overlearning the data, however large it might
be. In such situations, the trained neural networks while substantially over-
perform the theoretical optimum on the training set, they do not generalize
and perform poorly on other data sets.

These theoretical results are demonstrated in the following stylized example
with an explicit solution in [31](Section 8). In that example, the utility
function is taken to be the exponential with parameter one, and as the
decisions are independent of the initial value for these class of utilities,
the initial value is fixed as one dollar. To simplify even further, for one
period this amount is invested uniformly on all stocks. Then, with
${\bf{1}}:=(1,\ldots,1)$, $\pi_{0}={\bf{1}}/d$ and
$X_{1}=(Z_{1}\cdot{\bf{1}})/d-r$ are uncontrolled, and the investment problem
is to choose the feedback portfolio $\pi_{1}(Z_{1})\in\mathbb{R}^{d}$ so as to
maximize

$v(\pi)=\mathbb{E}\big{[}1-\exp(-X^{\pi}_{2})\big{]},$

where $X^{\pi}_{2}=(1+r)X_{1}+a(Z_{1})\cdot(Z_{2}-r{\bf{1}})$. The _certainty
equivalent_ of a utility value $v<1$ given by

$\mathrm{ce}(v):=\ln(1-v)\quad\Longleftrightarrow\quad v=U(\mathrm{ce}(v))$

is a more standard way of comparing different utility values. Indeed, agents
with expected utility preferences would be indifferent between an action $\pi$
and a cash amount of $\mathrm{ce}(v(\pi))$ because the utilities of both
positions are equal to each other. Thus, for these agents the cash equivalent
of the action $\pi$ is $\mathrm{ce}(v(\pi))$.

The following table [31](Table 1) clearly demonstrates overlearning. In this
experiments the training data of size $N=100,000$ and an artificial neural
network with three hidden layers of width 10 is trained on this set for four
or five epochs. For each dimension the algorithm is run thirty times and Table
4 below reports the mean and the standard. deviation. Although conservative
stopping rules are employed in [31], there is substantial overperformance
increasing with dimension.

dims |  $p_{in}$ (%) |  $p_{in}-p_{out}$ (%)
---|---|---
$\mu\ \quad$ | $\sigma\ \quad$ | $\mu\ \quad$ | $\sigma\ \quad$
100 | 10.12820 | 1.09290 | 23.67080 | 2.01177
85 | 8.38061 | 1.35575 | 20.16440 | 2.30489
70 | 7.32720 | 0.86458 | 15.62060 | 1.94043
55 | 5.05783 | 0.81518 | 10.93950 | 1.54431
40 | 3.74648 | 0.62588 | 7.91105 | 1.32581
25 | 2.11501 | 0.43845 | 4.58954 | 0.88461
10 | 0.53982 | 0.34432 | 1.46138 | 0.39078
Table 4: Average relative in-sample performance, and its comparison to the
out-of-sample performance with the above described conservative stopping rule.
Everything is in % with training size of $N=100,000$ and three hidden layers
of width 10. The $\mu$ value is the average of 30 runs and $\sigma$ is the
standard deviation.

##  7 Conclusion

The deep empirical risk minimization proposed by E, Han & Jentzen [21, 22]
provides a flexible and a highly effective tool for stochastic optimization
problems arising in computational finance. Recent development of optimization
libraries make this algorithm tractable in very high dimensions allowing to
include important market details such as factors and frictions, as well as
models with long memory. Once a large training set is given, the algorithm
mimics the market dynamics with all its details. This simple description
together with powerful new computational tools are keys to the power of the
algorithm. We have demonstrated above properties in three different classes of
problems. As it is always the case, each requires problem specific but natural
modifications. Moreover, the output can be designed to be exactly the decision
rule that is under investigation.

The method on the other hand needs large data sets for reliable results. In
the financial setting this essentially limits its scope to model driven
markets with an unlimited simulation capability. However, due to its seamless
transition to more complex structures, more interesting parametric models are
now feasible. Thus, on-going research on market generators will be an
important factor on further developments.

## References

  * Andersen and Broadie [2004] L. Andersen and M. Broadie.  Primal-dual simulation algorithm for pricing multidimensional American options.  _Management Science_ , 50(9):1222–1234, 2004\.
  * Bachouch et al. [2018] A. Bachouch, C. Huré, N. Langrené, and H. Pham.  Deep neural networks algorithms for stochastic control problems on finite horizon, part 2: Numerical applications.  _arXiv:1812.05916_ , 2018.
  * Bayer et al. [2020] C. Bayer, R. Tempone, and S. Wolfers.  Pricing American options by exercise rate optimization.  _Quantitative Finance_ , 20(11):1749–1760, 2020\.
  * Becker et al. [2019] S. Becker, P. Cheridito, and A. Jentzen.  Deep optimal stopping.  _Journal of Machine Learning Research_ , 20(4):1–25, 2019.
  * Becker et al. [2021] S. Becker, P. Cheridito, A. Jentzen, and T. Welti.  Solving high-dimensional optimal stopping problems using deep learning.  _European Journal of Applied Mathematics_ , 32(3):470–514, 2021.
  * Bertsekas and Tsitsiklis [1996] D. Bertsekas and J. Tsitsiklis.  _Neuro-dynamic programming_.  Athena Scientific, 1996.
  * Boudabsa and Filipović [2021] L. Boudabsa and D. Filipović.  Machine learning with kernels for portfolio valuation and risk management.  _Finance and Stochastics_ , pages 1–42, 2021.
  * Broadie and Detemple [1997] M. Broadie and J. Detemple.  The valuation of American options on multiple assets.  _Mathematical Finance_ , 7(3):241–286, 1997.
  * Buehler et al. [2019a] H. Buehler, L. Gonon, J. Teichmann, and B. Wood.  Deep hedging.  _Quantitative Finance_ , 19(8):1271–1291, 2019a.
  * Buehler et al. [2019b] H. Buehler, L. Gonon, J. Teichmann, B. Wood, and B. Mohan.  Deep hedging: hedging derivatives under generic market frictions using reinforcement learning.  Technical report, Swiss Finance Institute, 2019b.
  * Chevalier et al. [2021] E. Chevalier, S. Pulido, and E. Zúñiga.  American options in the Volterra Heston model.  _arXiv:2103.11734_ , 2021.
  * Ciocan and Mišić [2022] D. Ciocan and V. Mišić.  Interpretable optimal stopping.  _Management Science_ , 68(3):1616–1638, 2022\.
  * Cybenko [1989] G. Cybenko.  Approximations by superpositions of a sigmoidal function.  _Mathematics of Control, Signals and Systems_ , 2:183–192, 1989.
  * Detemple [2005] J. Detemple.  _American-style derivatives: Valuation and computation_.  CRC Press, 2005.
  * Fecamp et al. [2020] S. Fecamp, J. Mikael, and X. Warin.  Deep learning for discrete-time hedging in incomplete markets.  _Journal of computational Finance_ , 2020.
  * Filipovic et al. [2021] D. Filipovic, M. Multerer, and P. Schneider.  Adaptive joint distribution learning.  _arXiv:2110.04829_ , 2021.
  * Fleming and Soner [2006] W. H. Fleming and H. M. Soner.  _Controlled Markov processes and viscosity solutions_ , volume 25.  Springer Science & Business Media, 2006.
  * Föllmer and Schweizer [1991] H. Föllmer and M. Schweizer.  Hedging of contingent claims under incomplete information.  _Applied stochastic analysis_ , 5(389-414):19–31, 1991.
  * Germain et al. [2020] M. Germain, H. Pham, and X. Warin.  Deep backward multistep schemes for nonlinear pdes and approximation error analysis.  _arXiv:2006.01496_ , 2020.
  * Gonon et al. [2021] L. Gonon, J. Muhle-Karbe, and X. Shi.  Asset pricing with general transaction costs: Theory and numerics.  _Mathematical Finance_ , 31(2):595–648, 2021\.
  * Han and E [2016] J. Han and W. E.  Deep learning approximation for stochastic control problems.  In _Deep Reinforcement Learning Workshop, NIPS_ , 2016.
  * Han et al. [2018] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, 2018.
  * Heston [1993] S. L. Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  _The Review of Financial Studies_ , 6(2):327–343, 1993.
  * Hornik [1991] K. Hornik.  Approximation capabilities of multilayer feedforward networks.  _Neural networks_ , 4(2):251–257, 1991.
  * Huré et al. [2018] C. Huré, H. Pham, A. Bachouch, and N. Langrené.  Deep neural networks algorithms for stochastic control problems on finite horizon, part I: convergence analysis.  _arXiv:1812.04300_ , 2018.
  * Hutchinson et al. [1994] J. M. Hutchinson, A. W. Lo, and T. Poggio.  A nonparametric approach to pricing and hedging derivative securities via learning networks.  _The Journal of Finance_ , 49(3):851–889, 1994\.
  * Innes [2018] M. Innes.  Flux: Elegant machine learning with julia.  _Journal of Open Source Software_ , 3(25):602, 2018.
  * Laurière et al. [2021] M. Laurière, G. Pagès, and O. Pironneau.  Performance of a Markovian neural network versus dynamic programming on a fishing control problem.  _arXiv:2109.06856_ , 2021.
  * Longstaff and Schwartz [2001] F. A. Longstaff and E. S. Schwartz.  Valuing American options by simulation: a simple least-squares approach.  _The Review of Financial Studies_ , 14(1):113–147, 2001.
  * Ludkovski [2020] M. Ludkovski.  mlosp: Towards a unified implementation of regression Monte Carlo algorithms.  _arXiv:2012.00729_ , 2020.
  * Reppen and Soner [2020] A. M. Reppen and H. M. Soner.  Deep empirical risk minimization in finance: an assessment.  _in preparation_ , 2020.
  * Reppen et al. [2022] A. M. Reppen, H. M. Soner, and V. Tissot-Daguette.  Neural optimal stopping boundary.  _in preparation_ , 2022.
  * Ruf and Wang [2020a] J. Ruf and W. Wang.  Hedging with neural networks.  _arXiv:2004.08891_ , 2020a.
  * Ruf and Wang [2020b] J. Ruf and W. Wang.  Neural networks for option pricing and hedging: a literature review.  _Journal of Computational Finance_ , 24(1), 2020b.
  * Schweizer [1991] M. Schweizer.  Option hedging for semimartingales.  _Stochastic Processes and Their Applications_ , 37(2):339–363, 1991.
  * Schweizer [1999] M. Schweizer.  A guided tour through quadratic hedging approaches.  Technical report, SFB 373 Discussion Paper, 1999.
  * Sirignano and Spiliopoulos [2018] J. Sirignano and K. Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  _Journal of Computational Physics_ , 375:1339–1364, 2018\.
  * Tsitsiklis and Van Roy [2001] J. N. Tsitsiklis and B. Van Roy.  Regression methods for pricing complex American-style options.  _IEEE Transactions on Neural Networks_ , 12(4):694–703, 2001.
  * Wang and Perdikaris [2021] S. Wang and P. Perdikaris.  Deep learning of free boundary and Stefan problems.  _Journal of Computational Physics_ , 428:109914, 2021.
  * Warin [2019] X. Warin.  Variance optimal hedging with application to electricity markets.  _Journal of Computational Finance_ , 23(3):33–59, 2019.


$\sum_{\lambda}\epsilon_{\lambda}^{\mu}(p_{c},\,n)\epsilon_{\lambda}^{\nu*}(p_{c},\,n)=-g^{\mu\nu}+\frac{p_{c}^{\mu}n^{\nu}+p_{c}^{\nu}n^{\mu}}{p_{c}n}-n^{2}\frac{p_{c}^{\mu}p_{c}^{\nu}}{(p_{c}n)^{2}}\,,$
(218)

and obtain the gauge-invariant result121212 While it is nice that this result
is manifestly gauge-invariant, without Ward identities the polarization sum
rules of Eq. (218) would lead to a final $\left|\mathcal{M}\right|^{2}$ that
depends on the arbitrary gauge vector $n^{\mu}$. Instead, when computing the
matrix element, only the physical polarizations (corresponding to the gauge
choice $n^{\mu}=(1,0,0,-1)$) should be used, as prescribed long ago by
Altarelli and Parisi [88], and as we do in App. A.

$\left|\mathcal{M}_{a\to
bc}^{(0)}\right|^{2}=4E_{a}^{2}|g|^{2}\left(\frac{2p_{a,s}p_{b,h}}{p_{a,s}p_{c}~{}p_{b,h}p_{c}}-\frac{p_{a,s}^{2}}{(p_{a,s}p_{c})^{2}}-\frac{p_{b,h}^{2}}{(p_{b,h}p_{c})^{2}}\right).$
(219)

The perturbative splitting probability associated to Eq. (219) then reads

$dP^{{}^{\prime}}_{E}=\frac{d^{3}p_{c}}{(2\pi)^{3}2E_{c}}\frac{\left|\mathcal{M}_{a\to
bc}^{(0)}\right|^{2}}{4E_{a}}\simeq\zeta_{a}\frac{dk_{\perp}^{2}}{k_{\perp}^{2}}dx\,x\,\left(\frac{k_{\perp}^{2}+xm_{b,h}^{2}}{k_{\perp}^{2}+x^{2}m_{b,h}^{2}}\right)^{2}.$
(220)

The splitting probability in Eq. (220) does not vanish in the limit where the
order parameter of the phase transition goes to zero, either $m_{b,h}\ll
k_{\perp}$ or $m_{b,h}\to m_{b,s}$. Therefore the friction pressure is left
intact in the limit in which the wall disappears $\left<\phi\right>\to 0$
[58], thus signaling a possible problem with Eq. (220). We suggest that the
root of the problem may lie in the effective assumptions discussed earlier,
that lead to unexpectedly satisfying the Ward identities.

In contrast, our splitting probability does vanish when $m_{c,h}\ll k_{\perp}$
or $m_{b,h}\ll k_{\perp}$ (but also when $m_{i,h}\to m_{i,s}$), for vector
bosons $c$ which are respectively either massive, see Eq. (41), or massless,
see Eq. (189).

###  F.3 Splitting at all orders

##### The average exchange momentum.

In [1], the average exchanged momentum is computed according to

$\left<\Delta
p\right>=\int_{k_{\perp}>\mu}dP^{{}^{\prime}}_{E}\,\frac{k_{\perp}^{2}}{2xE_{a}}\,\exp\left[-\int_{\tilde{k}_{\perp}>k_{\perp}}dP^{{}^{\prime}}_{E}(\tilde{k}_{\perp})\right],$
(221)

where the perturbative splitting probability $dP^{{}^{\prime}}_{E}$ is given
by Eq. (220). The integrand in Eq. (221) is dominated by the region where
$x\sim 1$ and $k_{\perp}\sim E_{a}$ which leads the authors of [1] to conclude
that

$\left<\Delta p\right>\sim\zeta_{a}\,E_{a}.$ (222)

##### Our comments.

The Sudakov resummation operated in Eq. (221) can be obtained from the master
formula in Eq. (71) if we assume that the exchanged momentum in Eq. (70) is
dominated by the largest $k_{\perp}$

$\Delta
p=\underset{k_{\perp}}{\textrm{Max}}\left[\frac{k_{\perp,1}^{2}}{2x_{1}E_{a}},\,\frac{k_{\perp,2}^{2}}{2x_{2}E_{a}},\cdots,\frac{k_{\perp,n}^{2}}{2x_{n}E_{a}}\right].$
(223)

In that case, Eq. (71) becomes

$\displaystyle\left<\Delta p\right>$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\left[\prod_{j=1}^{n}\int
dP^{{}^{\prime}}_{E,\,j}\right]\Delta p~{}\exp\left[-\int
dP^{{}^{\prime}}_{E}\right]$ $\displaystyle=\left[\int
dP^{{}^{\prime}}_{E}(k_{\perp})\,\frac{k_{\perp}^{2}}{2xE_{a}}\right]\sum_{n=1}^{\infty}\frac{1}{(n-1)!}\left[\int
dP^{{}^{\prime}}_{E}(\tilde{k}_{\perp})\,\theta\big{(}k_{\perp}-\tilde{k}_{\perp}\big{)}\right]^{n-1}\,\exp\left[-\int
dP^{{}^{\prime}}_{E}\right]$ $\displaystyle=\int
dP^{{}^{\prime}}_{E}(k_{\perp})\,\frac{k_{\perp}^{2}}{2xE_{a}}\,\exp\left[\int_{\tilde{k}_{\perp}<k_{\perp}}dP^{{}^{\prime}}_{E}(\tilde{k}_{\perp})-\int_{\tilde{k}_{\perp}>\mu}dP^{{}^{\prime}}_{E}(\tilde{k}_{\perp})\right]$
$\displaystyle=\int_{k_{\perp}>\mu}dP^{{}^{\prime}}_{E}(k_{\perp})\,\frac{k_{\perp}^{2}}{2xE_{a}}\,\exp\left[-\int_{\tilde{k}_{\perp}>k_{\perp}}dP^{{}^{\prime}}_{E}(\tilde{k}_{\perp})\right],$
(224)

which coincides with Eq. (221). We have checked that to take the max in Eq.
(223) instead of the sum of the $\Delta p_{i}$ in Eq. (70) can lead to the
underestimation of the final $\left<\Delta p\right>$ by a factor
$\sim\mathcal{O}(30)$ in the limit of large supercooling $T_{\rm nuc}/T_{\rm
start}\sim 10^{-4}$ and large $\zeta_{a}\sim 0.1$.

##### Conclusion.

We conclude that the assumption $p_{a}=p_{b}+p_{c}$ at the level of the vertex
leads the authors of [1] to Eq. (220) instead of Eq. (41) for massive vector
boson, or instead of Eq. (189) for massless vector boson,131313The authors of
[1] consider the case of a massless vector boson $m_{c}=0$. which
overestimates the final $\left<\Delta p\right>$ by an amount $E_{a}/m_{c}$ or
$(E_{a}/m_{b})^{2}$ respectively, and does not vanish in the limit
$\left<\phi\right>\to 0$. Additionally, the approximation in Eq. (223) instead
of Eq. (70) leads the authors of [1] to underestimate $\left<\Delta p\right>$
by an amount $\sim\mathcal{O}(30)$.

###  F.4 Finite wall thickness

Finally, let us add that if the result $\left<\Delta
p\right>\sim\left<k_{\perp}^{2}/2xE_{a}\right>\sim\zeta_{a}\,E_{a}$ claimed by
[1] was right, the infinitely-thin-wall assumption $L_{\rm w}\Delta p\lesssim
1$ would violently break down. Indeed, the perturbative splitting probability
used by [1] and given in Eq. (220) should receive the additional factor
$W_{\rm wall}(L_{\rm w},\,k_{\perp},\,x)$ defined in Eq. (134), which
evaluated at $k_{\perp}^{2}/xE_{a}\sim\zeta_{a}\,E_{a}$, is exponentially
suppressed

$W_{\rm wall}\propto e^{-\pi L_{\rm w}E_{a}}\sim
e^{-\pi\frac{E_{a}}{m_{c,h}}}\ll 1.$ (225)

Figure 16:  We display here the average momentum transferred to the wall
$\Delta p$ as a function of the IR cut-off $\mu$, both in some units of
$m_{c}$. In contrast to the analytical result in dotted, cf. Eq (80), which
overestimates the exchanged momentum in the limit $\mu\to 0$, the MC
simulation in solid accounts for the depletion of the incoming energy-momentum
due to the successive boson emissions, see Sec. 5.2. This implies a maximal
number of boson emission and a maximal value for the exchanged momentum
compatible with energy-momentum conservation. The impact of energy depletion
is higher and higher as we decrease the incoming energy $E_{a}/m_{c}$. For
simplicity, $m_{c}=m_{c,h}$ in this plot.

Figure 17:  Comparison between analytical treatment in dashed lines and MC
treatment in solid lines. The blue and green lines show the exchanged momentum
induced by radiated bosons resummed at all leading-log orders, either
transmitted (green) or reflected (blue) by the wall boundary. The red dashed
line shows the LO contribution assuming no splitting radiation. In the run-
away regime, we have $E_{a}/m_{c,h}\sim(M_{\rm pl}/\left<\phi\right>)(T_{\rm
nuc}/T_{\rm start})^{2}$, see Eq. (60) and Eq. (105). Therefore, at very
strong supercooling the initial particle energy $E_{a}$ in the wall frame
becomes smaller than the particle masses $m_{a,h}$, $m_{b,h}$ and/or $m_{c,h}$
in the broken phase, such that $a$, $b$ and/or $c$ are reflected.

##  Appendix G Energy-momentum conservation

###  G.1 Backreaction

In the Sudakov resummation in Sec. 4.2, we have neglected the depletion of the
initial energy and momentum $(E_{b},\,p_{b})$ due to the multiple boson
emissions. When considering those effects, which we call ‘backreaction’, Eq.
(74) should instead become

$\displaystyle\left<\Delta
p\right>=\sum_{n=0}^{\infty}\frac{1}{n!}\left[\prod_{j=1}^{n}\int
dP_{E,\,j}\Theta\big{(}(1-X)^{2}E_{a}^{2}-K_{\perp}^{2}\big{)}\Theta\big{(}1-X\big{)}\right]\sum_{i=1}^{n}\Delta
p_{i}~{}\exp\left[-\int dP_{E}\right].$ (226)

where $X$ and $K_{\perp}$ are defined in Eq. (67). The $\Theta$ functions
prevent the emission of more momentum and energy than what is available, i.e
the momentum $p_{b}^{z}$ and energy $E_{b}$ in Eq. (15) must remain positive.
We account for the effect of backreaction with our Monte-Carlo simulation, cf.
Sec. 5, in which the Theta functions in Eq. (226) enter through the boundaries
of the integral in Eq. (92) and through the condition in Eq. (93).

In Fig. 16, we can see that in contrast to the analytical result (dotted
line), in the numerical study in which the effect of energy depletion due to
successive emissions is included (solid lines), the number of emissions and
the resulting $\left<\Delta p\right>$ saturates to a plateau. This could by
itself regulate the logarithmic divergence if the thermal mass as well as
other possible IR cut-offs were zero, see Sec. 3.4. However, in presence of
the thermally induced cut-off $\mu=\alpha^{1/2}T_{\rm nuc}$, we find that
backreaction can be safely neglected.

###  G.2 Snowplow region: $E_{a}\lesssim m_{i,h}$

The expressions for the friction pressure, $\mathcal{P}_{{{\text{LO}}}}$ and
$\mathcal{P}_{{{\text{LL}}}}$, induced by $1\to 1$ and $1\to n$ processes, cf.
Eq. (96) and Eq. (98), are valid in the regime $m_{a,h},\,m_{b,h}\ll E_{a}$
and $\mu\ll m_{c,h}\ll E_{a}$, respectively. We now derive their expression in
the regime where $E_{a}\lesssim m_{a,h},\,m_{b,h}$ and $E_{a}\lesssim
m_{c,h}$, in which the particles $a$, $b$ and $c$ are reflected by the wall
boundary. Since the particles then accumulate in front of the wall, we call
this the ‘snowplow region’.141414 We thank Gilad Perez for suggesting the name
‘bulldozzer’, we end up choosing ‘snowplow’ due to proximity with the
Christmas season at the time of arXiv first submission.  Assuming for the sake
of simplicity that $a$ and $b$ are identical, the exchanged momentum between
$a$ and the wall is given by $\left<\Delta p\right>=\left<2p_{a}\right>$, and
the friction pressure at LO in the relativistic limit is, cf. Eq. (1)

$E_{a}\lesssim
m_{a,h}\qquad\implies\qquad\mathcal{P}_{{{\text{LO}}}}\simeq\frac{\pi^{2}}{15}\sum_{a}g_{a}e_{a}~{}\gamma^{2}\,T_{\rm
nuc}^{4},\qquad e_{a}=1~{}(7/8)~{}\textrm{for bosons (fermions)}$ (227)

The exchanged momentum due to splitting radiation resummed at all leading-log
orders is, cf. Eq. (22), Eq. (72) and Eq. (74)

$\displaystyle\left<\Delta
p\right>\simeq\int_{0}^{\frac{E_{a}^{2}}{4}}\frac{dk_{\perp}^{2}}{k_{\perp}^{2}}\,\int_{1}^{\frac{\sqrt{k_{\perp}^{2}+m_{c,h}^{2}}}{E_{a}}}\frac{dx}{x}\,\frac{k_{\perp}^{4}}{(k_{\perp}^{2}+\mu^{2})^{2}}\,2xE_{a}\simeq
4E_{a}\log\left(\frac{E_{a}}{\mu}\right),$ (228)

which after injection into Eq. (6) with $E_{a}=\left<p_{a}\right>\simeq
2.7\,\gamma\,T_{\rm nuc}$, leads to the friction pressure151515Note that we
recover the scaling $\mathcal{P}_{{\text{LL}}}\propto\alpha\gamma^{2}\,T_{\rm
nuc}^{4}$ found in [1], but only in the limit $E_{a}\lesssim m_{\rm c,h}$.

$\displaystyle E_{a}\lesssim
m_{c,h}\qquad\implies\qquad\mathcal{P}_{{\text{LL}}}\simeq\frac{11\,\zeta(3)}{\pi^{3}}\left[\sum_{a,b,c}\nu_{a}g_{a}C_{abc}\right]\alpha\,\text{ln}\left(\dfrac{E_{a}}{\mu}\right)\,\gamma^{2}\,T_{\rm
nuc}^{4}.$ (229)

The exchanged momentum due to the Eq. (228) is shown by the dashed blue line
in Fig. 17. Note that the region $E_{a}\lesssim m_{c,h}\sim m_{\rm a,h}$ falls
in the region where interactions between reflected particles and newly
incoming particles can be neglected, cf. dotted purple and blue lines in Fig.
8, such that the possible corrections to the friction pressure discussed in
Sec. 4.4 are not expected. It is interesting that a region $E_{a}\lesssim
m_{i,h}$, where particles are too soft to penetrate inside the bubble wall,
exists in the large supercooling limit (as already pointed out in [12] for
confining PTs). We leave its study for further works.

##  Appendix H Lorentz factor in the run-away regime

###  H.1 Generic potential

##### Energy conservation.

The vacuum energy gained upon formation of a bubble of radius $R$ is

$E_{\text{bubble}}\simeq\frac{4}{3}\,\pi\,R^{3}~{}\mathcal{P},$ (230)

where $\mathcal{P}$ is the expanding pressure. As the bubble grows, the energy
stored in the bubble wall is given by

$E_{\text{wall}}\simeq 4\,\pi\,R^{2}\,\gamma~{}\sigma,$ (231)

where $\sigma$ is the surface energy of the wall in the wall frame

$\sigma\equiv\int_{0}^{\infty}dr\,\left[\frac{1}{2}(\phi^{{}^{\prime}}(r))^{2}+V(\phi(r))-V(\phi(0))\right].$
(232)

Insuring energy conservation by equating Eq. (230) and Eq. (231), we obtain
that in the run-away regime, the wall Lorentz factor grows linearly with the
bubble radius, e.g. [80, 116]

$\gamma=\frac{\mathcal{P}}{3\,\sigma}R.$ (233)

In the supercooled and run-away regime, the expanding pressure is given by the
vacuum energy difference $\mathcal{P}=c_{\rm vac}\left<\phi\right>^{4}$ while
the surface tension can be estimated as

$\sigma\approx L_{\rm w}^{\rm tot.}c_{\rm vac}\left<\phi\right>^{4},$ (234)

where $L_{\rm w}^{\rm tot.}$ is the total thickness of the wall, accounting
for the rising part plus the multiple oscillations around the true vacuum
value, until the oscillations become negligible. This should not be mistaken
with $L_{\rm w}$ introduced earlier, which is the thickness of the rising part
of the wall only, see footnote 11. The Lorentz factor of the running-away
bubble wall in Eq. (233) becomes

$\gamma=\frac{R}{3\,L_{\rm w}^{\rm tot.}}.$ (235)

##### Total wall thickness.

The wall profile can be found by sticking the space-like profile with the
time-like profile, see e.g. [117]. The former is solution of the euclidean
equation of motion

$\phi_{\rm E}^{\prime\prime}(s_{\rm E})+\frac{d-1}{s_{\rm E}}\phi_{\rm
E}^{\prime}(s_{\rm E})=\frac{dV}{d\phi_{\rm
E}},\qquad\textrm{with}\quad\phi_{\rm
E}^{{}^{\prime}}(0)=0,\quad\textrm{and}\quad\lim_{r\to\infty}\phi_{\rm
E}(r)=0.$ (236)

where $s_{\rm E}=\sqrt{\vec{r}^{2}+t_{E}^{2}}=\sqrt{\vec{r}^{2}-t^{2}}$ is the
space-like light-cone coordinate, $t_{E}=i\,t$ is the Euclidean time and $d=3$
or $4$. The latter is solution of the real-time equation of motion

$\frac{\partial^{2}\phi}{\partial
s^{2}}+\frac{3}{s}\frac{\partial\phi}{\partial s}+\frac{\partial
V}{\partial\phi}=0,\qquad\textrm{with}\quad\phi(s)=\phi_{\rm
E}(0),\quad\textrm{and}\quad\phi^{{}^{\prime}}(0)=\phi_{\rm
E}^{{}^{\prime}}(0).$ (237)

where $s=\sqrt{t^{2}-\vec{r}^{2}}$ is now the time-like light-cone coordinate.
The total wall thickness $L_{\rm w}^{\rm tot.}$ is set by the sum of the
characteristic length scales of the two profiles

$L_{\rm w}^{\rm tot.}=L_{\rm space-like}+L_{\rm time-like},$ (238)

which are themselves set by the damping terms of the two respective Eqs. (236)
and (237). Due to the matching condition $\phi^{{}^{\prime}}(0)=\phi_{\rm
E}^{{}^{\prime}}(0)$, the two damping terms are of the same order and we
conclude that

$L_{\rm time-like}\simeq L_{\rm space-like}\equiv R_{\rm nuc},$ (239)

where the last equality defines the bubble radius at nucleation. The Lorentz
factor of the running-away bubble wall in Eq. (235) becomes161616 In [80], the
wall tension is expressed in term of the bubble radius at nucleation $R_{\rm
nuc}$, obtained from minimizing the total energy
$E_{\text{bubble}}+E_{\text{wall}}$ $\sigma=R_{\rm nuc}\mathcal{P}/2,$ (240)
such that the Lorentz factor of the running-away bubble wall in Eq. (233)
becomes  $\gamma=\frac{2R}{3R_{\rm nuc}}.$ (241) However, as the bubble
expands, the scalar field undergoes damped oscillations toward the true
vacuum, such that Eq. (240) under-estimates the surface tension $\sigma$. With
Eq. (242), we claim that Eq. (240) only under-estimates $\sigma$ by an
$\mathcal{O}(1)$ factor.

$\gamma=\frac{R}{3R_{\rm nuc}}.$ (242)

###  H.2 Shallow potential

The friction pressure in Eq. (95) is suppressed by powers of $T_{\rm nuc}/f\ll
1$, where $f\equiv\left<\phi\right>$. Therefore, a run-away regime at bubble
collision time requires supercooled phase transitions. Typically, supercooled
phase transitions are generated by shallow zero-temperature potentials, namely
potentials where the curvature close to the false vacuum

$\sqrt{V^{\prime\prime}(\phi)}\Big{|}_{\phi\ll f}\simeq
f\,\exp(-c/\epsilon),\qquad\epsilon\ll 1,\quad c=O(1),$ (243)

is much smaller than the curvature $f$ close to the true vacuum. In that case,
the tunneling exit point is very close to the false minimum and the bounce
action is only sensitive to the scale $f\,\exp(-c/\epsilon)$. We expect both
temperature and bubble radius at nucleation to be related to that scale, and
therefore171717 Note that the authors of [39] have set $R_{\rm nuc}\sim
f^{-1}$ instead of $R_{\rm nuc}\sim T_{\rm nuc}^{-1}$. For typical potentials
leading to large supercooling, e.g. Coleman-Weinberg [76] or light-dilation
[74], the latter choice is the correct one, see App. A of [12].

$R_{\rm nuc}\simeq c_{\rm w}\,T_{\rm nuc}^{-1},$ (244)

where $c_{\rm w}$ is a model-dependent numerical factor. At the time of
collision, the bubble size is given by181818 The factor $8\pi$ comes from the
expected number of bubbles nucleating in a volume $V$. To derive this, one has
to take into account the fact that the nucleating bubble must be in the false
vacuum. Taking the wall velocity to be unity for simplicity, and taking the
nucleation rate per unit time and volume to be $\Gamma(t)=\Gamma_{*}e^{\beta
t}$, we sum up the differential nucleation probability for time interval
$[t_{n},t_{n}+dt_{n}]$ to get (expected $\\#$ of bubbles) $=V\times\sum$
(prob. for the nucleation point to be in the false vacuum) $\times$ (prob. for
nucleation between $[t_{n},t_{n}+dt_{n}]$)
$=V\times\int_{-\infty}^{\infty}\,e^{-\frac{4\pi}{3}\int_{-\infty}^{t_{n}}dt\,(t_{n}-t)^{3}\Gamma_{*}e^{\beta
t}}\times\Gamma_{*}e^{\beta
t_{n}}dt_{n}=V\times\int_{-\infty}^{\infty}\,e^{-8\pi\Gamma_{*}e^{\beta
t_{n}}/\beta^{4}}\times\Gamma_{*}e^{\beta t_{n}}dt_{n}=\beta^{3}V/8\pi$.

$R_{\rm coll}\simeq(8\pi)^{1/3}\,v_{w}\beta^{-1},$ (245)

where $v_{w}$ is the bubble wall velocity, $\beta$ is the expansion
coefficient of the exponent of the nucleation rate around the typical
transition time $t_{*}$

$\displaystyle\Gamma$ $\displaystyle\propto e^{\beta(t-t_{*})+\cdots},$ (246)

which gives the inverse of the bubble propagation time. Therefore, in the run-
away regime ($v_{w}=1$ and $\mathcal{P}>0$), the Lorentz factor at the time of
collision in Eq. (242) becomes

$\gamma_{\rm run}\simeq~{}\frac{T_{\rm nuc}}{c_{\rm w}\beta},$ (247)

which we have use in the main text, cf. Eq. (105).

## References

  * [1] S. Höche, J. Kozaczuk, A. J. Long, J. Turner and Y. Wang, _Towards an all-orders calculation of the electroweak bubble wall velocity_ , _JCAP_ 03 (2021) 009, [2007.10343].
  * [2] E. Witten, _Cosmic Separation of Phases_ , _Phys. Rev._ D30 (1984) 272–285.
  * [3] C. J. Hogan, _Gravitational radiation from cosmological phase transitions_ , _Mon. Not. Roy. Astron. Soc._ 218 (1986) 629–636.
  * [4] M. Kamionkowski, A. Kosowsky and M. S. Turner, _Gravitational radiation from first order phase transitions_ , _Phys. Rev._ D49 (1994) 2837–2851, [astro-ph/9310044].
  * [5] T. Konstandin and G. Servant, _Cosmological Consequences of Nearly Conformal Dynamics at the TeV scale_ , _JCAP_ 1112 (2011) 009, [1104.4791].
  * [6] A. Falkowski and J. M. No, _Non-thermal Dark Matter Production from the Electroweak Phase Transition: Multi-TeV WIMPs and ’Baby-Zillas’_ , _JHEP_ 02 (2013) 034, [1211.5615].
  * [7] T. Hambye and A. Strumia, _Dynamical generation of the weak and Dark Matter scale_ , _Phys. Rev. D_ 88 (2013) 055022, [1306.2329].
  * [8] T. Hambye, A. Strumia and D. Teresi, _Super-cool Dark Matter_ , _JHEP_ 08 (2018) 188, [1805.01473].
  * [9] Y. Bai, A. J. Long and S. Lu, _Dark Quark Nuggets_ , _Phys. Rev. D_ 99 (2019) 055047, [1810.04360].
  * [10] M. J. Baker, J. Kopp and A. J. Long, _Filtered Dark Matter at a First Order Phase Transition_ , _Phys. Rev. Lett._ 125 (2020) 151102, [1912.02830].
  * [11] A. Azatov, M. Vanvlasselaer and W. Yin, _Dark Matter production from relativistic bubble walls_ , _JHEP_ 03 (2021) 288, [2101.05721].
  * [12] I. Baldes, Y. Gouttenoire and F. Sala, _String Fragmentation in Supercooled Confinement and Implications for Dark Matter_ , _JHEP_ 04 (2021) 278, [2007.08440].
  * [13] J.-P. Hong, S. Jung and K.-P. Xie, _Fermi-ball dark matter from a first-order phase transition_ , _Phys. Rev. D_ 102 (2020) 075028, [2008.04430].
  * [14] P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer and J. Smirnov, _Thermal squeezeout of dark matter_ , _Phys. Rev. D_ 104 (2021) 095013, [2103.09827].
  * [15] I. Baldes, Y. Gouttenoire, F. Sala and G. Servant, _Supercool Composite Dark Matter beyond 100 TeV_ , 2110.13926.
  * [16] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, _On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe_ , _Phys. Lett. B_ 155 (1985) 36.
  * [17] A. G. Cohen, D. B. Kaplan and A. E. Nelson, _Progress in electroweak baryogenesis_ , _Ann. Rev. Nucl. Part. Sci._ 43 (1993) 27–70, [hep-ph/9302210].
  * [18] T. Konstandin and G. Servant, _Natural Cold Baryogenesis from Strongly Interacting Electroweak Symmetry Breaking_ , _JCAP_ 1107 (2011) 024, [1104.4793].
  * [19] G. Servant, _Baryogenesis from Strong $CP$ Violation and the QCD Axion_, _Phys. Rev. Lett._ 113 (2014) 171803, [1407.0030].
  * [20] A. Katz and A. Riotto, _Baryogenesis and Gravitational Waves from Runaway Bubble Collisions_ , _JCAP_ 1611 (2016) 011, [1608.00583].
  * [21] J. M. Cline and K. Kainulainen, _Electroweak baryogenesis at high bubble wall velocities_ , _Phys. Rev. D_ 101 (2020) 063525, [2001.00568].
  * [22] G. C. Dorsch, S. J. Huber and T. Konstandin, _On the wall velocity dependence of electroweak baryogenesis_ , _JCAP_ 08 (2021) 020, [2106.06547].
  * [23] A. Azatov, M. Vanvlasselaer and W. Yin, _Baryogenesis via relativistic bubble walls_ , _JHEP_ 10 (2021) 043, [2106.14913].
  * [24] I. Baldes, S. Blasi, A. Mariotti, A. Sevrin and K. Turbang, _Baryogenesis via relativistic bubble expansion_ , 2106.15602.
  * [25] M. J. Baker, M. Breitbach, J. Kopp, L. Mittnacht and Y. Soreq, _Filtered Baryogenesis_ , 2112.08987.
  * [26] H. Kodama, M. Sasaki and K. Sato, _Abundance of Primordial Holes Produced by Cosmological First Order Phase Transition_ , _Prog. Theor. Phys._ 68 (1982) 1979.
  * [27] S. W. Hawking, I. G. Moss and J. M. Stewart, _Bubble Collisions in the Very Early Universe_ , _Phys. Rev._ D26 (1982) 2681.
  * [28] C. Gross, G. Landini, A. Strumia and D. Teresi, _Dark Matter as dark dwarfs and other macroscopic objects: multiverse relics?_ , _JHEP_ 09 (2021) 033, [2105.02840].
  * [29] M. J. Baker, M. Breitbach, J. Kopp and L. Mittnacht, _Primordial Black Holes from First-Order Cosmological Phase Transitions_ , 2105.07481.
  * [30] M. J. Baker, M. Breitbach, J. Kopp and L. Mittnacht, _Detailed Calculation of Primordial Black Hole Formation During First-Order Cosmological Phase Transitions_ , 2110.00005.
  * [31] K. Kawana and K.-P. Xie, _Primordial black holes from a cosmic phase transition: The collapse of Fermi-balls_ , _Phys. Lett. B_ 824 (2022) 136791, [2106.00111].
  * [32] T. Vachaspati, _Magnetic fields from cosmological phase transitions_ , _Phys. Lett. B_ 265 (1991) 258–261.
  * [33] T. W. B. Kibble, _Topology of Cosmic Domains and Strings_ , _J. Phys. A_ 9 (1976) 1387–1398.
  * [34] T. W. B. Kibble and A. Vilenkin, _Phase equilibration in bubble collisions_ , _Phys. Rev. D_ 52 (1995) 679–688, [hep-ph/9501266].
  * [35] J. Borrill, T. W. B. Kibble, T. Vachaspati and A. Vilenkin, _Defect production in slow first order phase transitions_ , _Phys. Rev. D_ 52 (1995) 1934–1943, [hep-ph/9503223].
  * [36] N. Craig, N. Levi, A. Mariotti and D. Redigolo, _Ripples in Spacetime from Broken Supersymmetry_ , _JHEP_ 21 (2020) 184, [2011.13949].
  * [37] C. Caprini et al., _Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions_ , _JCAP_ 1604 (2016) 001, [1512.06239].
  * [38] R.-G. Cai, Z. Cao, Z.-K. Guo, S.-J. Wang and T. Yang, _The Gravitational-Wave Physics_ , _Natl. Sci. Rev._ 4 (2017) 687–706, [1703.00187].
  * [39] C. Caprini et al., _Detecting gravitational waves from cosmological phase transitions with LISA: an update_ , _JCAP_ 03 (2020) 024, [1910.13125].
  * [40] C. Caprini and D. G. Figueroa, _Cosmological Backgrounds of Gravitational Waves_ , _Class. Quant. Grav._ 35 (2018) 163001, [1801.04268].
  * [41] A. Mazumdar and G. White, _Review of cosmic phase transitions: their significance and experimental signatures_ , _Rept. Prog. Phys._ 82 (2019) 076901, [1811.01948].
  * [42] M. B. Hindmarsh, M. Lüben, J. Lumma and M. Pauly, _Phase transitions in the early universe_ , 2008.09136.
  * [43] J. Ignatius, K. Kajantie, H. Kurki-Suonio and M. Laine, _The growth of bubbles in cosmological phase transitions_ , _Phys. Rev. D_ 49 (1994) 3854–3868, [astro-ph/9309059].
  * [44] G. D. Moore and T. Prokopec, _How fast can the wall move? A Study of the electroweak phase transition dynamics_ , _Phys. Rev. D_ 52 (1995) 7182–7204, [hep-ph/9506475].
  * [45] G. D. Moore, _Electroweak bubble wall friction: Analytic results_ , _JHEP_ 03 (2000) 006, [hep-ph/0001274].
  * [46] P. John and M. G. Schmidt, _Do stops slow down electroweak bubble walls?_ , _Nucl. Phys. B_ 598 (2001) 291–305, [hep-ph/0002050].
  * [47] A. Megevand and A. D. Sanchez, _Detonations and deflagrations in cosmological phase transitions_ , _Nucl. Phys. B_ 820 (2009) 47–74, [0904.1753].
  * [48] A. Megevand and A. D. Sanchez, _Velocity of electroweak bubble walls_ , _Nucl. Phys. B_ 825 (2010) 151–176, [0908.3663].
  * [49] J. R. Espinosa, T. Konstandin, J. M. No and G. Servant, _Energy Budget of Cosmological First-order Phase Transitions_ , _JCAP_ 1006 (2010) 028, [1004.4187].
  * [50] S. J. Huber and M. Sopena, _An efficient approach to electroweak bubble velocities_ , 1302.1044.
  * [51] G. C. Dorsch, S. J. Huber and T. Konstandin, _Bubble wall velocities in the Standard Model and beyond_ , _JCAP_ 12 (2018) 034, [1809.04907].
  * [52] A. Friedlander, I. Banta, J. M. Cline and D. Tucker-Smith, _Wall speed and shape in singlet-assisted strong electroweak phase transitions_ , _Phys. Rev. D_ 103 (2021) 055020, [2009.14295].
  * [53] S. Balaji, M. Spannowsky and C. Tamarit, _Cosmological bubble friction in local equilibrium_ , _JCAP_ 03 (2021) 051, [2010.08013].
  * [54] M. Barroso Mancha, T. Prokopec and B. Swiezewska, _Field-theoretic derivation of bubble-wall force_ , _JHEP_ 01 (2021) 070, [2005.10875].
  * [55] W.-Y. Ai, B. Garbrecht and C. Tamarit, _Bubble wall velocities in local equilibrium_ , 2109.13710.
  * [56] D. Bodeker and G. D. Moore, _Can electroweak bubble walls run away?_ , _JCAP_ 0905 (2009) 009, [0903.4099].
  * [57] D. Bodeker and G. D. Moore, _Electroweak Bubble Wall Speed Limit_ , _JCAP_ 1705 (2017) 025, [1703.08215].
  * [58] A. Azatov and M. Vanvlasselaer, _Bubble wall velocity: heavy physics effects_ , _JCAP_ 01 (2021) 058, [2010.02590].
  * [59] E. W. Kolb and S. Wolfram, _Spontaneous Symmetry Breaking and the Expansion Rate of the Early Universe_ , _Astrophys. J._ 239 (1980) 428.
  * [60] P. Creminelli, A. Nicolis and R. Rattazzi, _Holography and the electroweak phase transition_ , _JHEP_ 03 (2002) 051, [hep-th/0107141].
  * [61] L. Randall and G. Servant, _Gravitational waves from warped spacetime_ , _JHEP_ 05 (2007) 054, [hep-ph/0607158].
  * [62] G. Nardini, M. Quiros and A. Wulzer, _A Confining Strong First-Order Electroweak Phase Transition_ , _JHEP_ 09 (2007) 077, [0706.3388].
  * [63] J. R. Espinosa, T. Konstandin, J. M. No and M. Quiros, _Some Cosmological Implications of Hidden Sectors_ , _Phys. Rev. D_ 78 (2008) 123528, [0809.3215].
  * [64] T. Konstandin, G. Nardini and M. Quiros, _Gravitational Backreaction Effects on the Holographic Phase Transition_ , _Phys. Rev._ D82 (2010) 083513, [1007.1468].
  * [65] R. Jinno and M. Takimoto, _Probing a classically conformal B-L model with gravitational waves_ , _Phys. Rev. D_ 95 (2017) 015020, [1604.05035].
  * [66] L. Marzola, A. Racioppi and V. Vaskonen, _Phase transition and gravitational wave phenomenology of scalar conformal extensions of the Standard Model_ , _Eur. Phys. J. C_ 77 (2017) 484, [1704.01034].
  * [67] S. Iso, P. D. Serpico and K. Shimada, _QCD-Electroweak First-Order Phase Transition in a Supercooled Universe_ , _Phys. Rev. Lett._ 119 (2017) 141301, [1704.04955].
  * [68] C.-W. Chiang and E. Senaha, _On gauge dependence of gravitational waves from a first-order phase transition in classical scale-invariant $U(1)^{\prime}$ models_, _Phys. Lett. B_ 774 (2017) 489–493, [1707.06765].
  * [69] B. von Harling and G. Servant, _QCD-induced Electroweak Phase Transition_ , _JHEP_ 01 (2018) 159, [1711.11554].
  * [70] K. Fujikura, Y. Nakai and M. Yamada, _A more attractive scheme for radion stabilization and supercooled phase transition_ , _JHEP_ 02 (2020) 111, [1910.07546].
  * [71] E. Megías, G. Nardini and M. Quirós, _Cosmological Phase Transitions in Warped Space: Gravitational Waves and Collider Signatures_ , _JHEP_ 09 (2018) 095, [1806.04877].
  * [72] S. Arunasalam, A. Kobakhidze, C. Lagger, S. Liang and A. Zhou, _Low temperature electroweak phase transition in the Standard Model with hidden scale invariance_ , _Phys. Lett. B_ 776 (2018) 48–53, [1709.10322].
  * [73] S. Bruggisser, B. Von Harling, O. Matsedonskyi and G. Servant, _Baryon Asymmetry from a Composite Higgs Boson_ , _Phys. Rev. Lett._ 121 (2018) 131801, [1803.08546].
  * [74] S. Bruggisser, B. Von Harling, O. Matsedonskyi and G. Servant, _Electroweak Phase Transition and Baryogenesis in Composite Higgs Models_ , _JHEP_ 12 (2018) 099, [1804.07314].
  * [75] P. Baratella, A. Pomarol and F. Rompineve, _The Supercooled Universe_ , _JHEP_ 03 (2019) 100, [1812.06996].
  * [76] I. Baldes and C. Garcia-Cely, _Strong gravitational radiation from a simple dark matter model_ , _JHEP_ 05 (2019) 190, [1809.01198].
  * [77] T. Prokopec, J. Rezacek and B. Świeżewska, _Gravitational waves from conformal symmetry breaking_ , _JCAP_ 02 (2019) 009, [1809.11129].
  * [78] V. Brdar, A. J. Helmboldt and J. Kubo, _Gravitational Waves from First-Order Phase Transitions: LIGO as a Window to Unexplored Seesaw Scales_ , _JCAP_ 02 (2019) 021, [1810.12306].
  * [79] C. Marzo, L. Marzola and V. Vaskonen, _Phase transition and vacuum stability in the classically conformal B–L model_ , _Eur. Phys. J. C_ 79 (2019) 601, [1811.11169].
  * [80] J. Ellis, M. Lewicki, J. M. No and V. Vaskonen, _Gravitational wave energy budget in strongly supercooled phase transitions_ , _JCAP_ 1906 (2019) 024, [1903.09642].
  * [81] K. Agashe, P. Du, M. Ekhterachian, S. Kumar and R. Sundrum, _Cosmological Phase Transition of Spontaneous Confinement_ , 1910.06238.
  * [82] K. Agashe, P. Du, M. Ekhterachian, S. Kumar and R. Sundrum, _Phase Transitions from the Fifth Dimension_ , _JHEP_ 02 (2021) 051, [2010.04083].
  * [83] L. Delle Rose, G. Panico, M. Redi and A. Tesi, _Gravitational Waves from Supercool Axions_ , _JHEP_ 04 (2020) 025, [1912.06139].
  * [84] B. Von Harling, A. Pomarol, O. Pujolas and F. Rompineve, _Peccei-Quinn Phase Transition at LIGO_ , 1912.07587.
  * [85] I. M. Bloch, C. Csáki, M. Geller and T. Volansky, _Crunching away the cosmological constant problem: dynamical selection of a small $\Lambda$_, _JHEP_ 12 (2020) 191, [1912.08840].
  * [86] Y. Gouttenoire, _Beyond the Standard Model Cocktail_ , 2207.01633.
  * [87] A. Caputo, M. Regis, M. Taoso and S. J. Witte, _Detecting the Stimulated Decay of Axions at RadioFrequencies_ , _JCAP_ 03 (2019) 027, [1811.08436].
  * [88] G. Altarelli and G. Parisi, _Asymptotic Freedom in Parton Language_ , _Nucl. Phys. B_ 126 (1977) 298–318.
  * [89] P. B. Arnold, G. D. Moore and L. G. Yaffe, _Effective kinetic theory for high temperature gauge theories_ , _JHEP_ 01 (2003) 030, [hep-ph/0209353].
  * [90] P. B. Arnold, G. D. Moore and L. G. Yaffe, _Transport coefficients in high temperature gauge theories. 2. Beyond leading log_ , _JHEP_ 05 (2003) 051, [hep-ph/0302165].
  * [91] R. Ellis, W. Stirling and B. Webber, _QCD and collider physics_ , vol. 8.  Cambridge University Press, 2, 2011.
  * [92] M. E. Peskin and D. V. Schroeder, _An Introduction to quantum field theory_.  Addison-Wesley, Reading, USA, 1995.
  * [93] M. D. Schwartz, _Quantum Field Theory and the Standard Model_.  Cambridge University Press, 3, 2014.
  * [94] A. J. Larkoski, _An Unorthodox Introduction to QCD_ , 1709.06195.
  * [95] E. W. Kolb and M. S. Turner, _The Early Universe_ , _Front. Phys._ 69 (1990) 1–547.
  * [96] A. Banfi, G. P. Salam and G. Zanderighi, _Principles of general final-state resummation and automated implementation_ , _JHEP_ 03 (2005) 073, [hep-ph/0407286].
  * [97] G. Marchesini and B. R. Webber, _Monte Carlo Simulation of General Hard Processes with Coherent QCD Radiation_ , _Nucl. Phys. B_ 310 (1988) 461–526.
  * [98] R. Jinno, H. Seong, M. Takimoto and C. M. Um, _Gravitational waves from first-order phase transitions: Ultra-supercooled transitions and the fate of relativistic shocks_ , _JCAP_ 1910 (2019) 033, [1905.00899].
  * [99] J. Ellis, M. Lewicki and J. M. No, _Gravitational waves from first-order cosmological phase transitions: lifetime of the sound wave source_ , _JCAP_ 07 (2020) 050, [2003.07360].
  * [100] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, _Gravitational waves from the sound of a first order phase transition_ , _Phys. Rev. Lett._ 112 (2014) 041301, [1304.2433].
  * [101] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, _Numerical simulations of acoustically generated gravitational waves at a first order phase transition_ , _Phys. Rev._ D92 (2015) 123009, [1504.03291].
  * [102] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, _Shape of the acoustic gravitational wave power spectrum from a first order phase transition_ , _Phys. Rev._ D96 (2017) 103520, [1704.05871].
  * [103] M. Hindmarsh, _Sound shell model for acoustic gravitational wave production at a first-order phase transition in the early Universe_ , _Phys. Rev. Lett._ 120 (2018) 071301, [1608.04735].
  * [104] R.-G. Cai and S.-J. Wang, _Energy budget of cosmological first-order phase transition in FLRW background_ , _Sci. China Phys. Mech. Astron._ 61 (2018) 080411, [1803.03002].
  * [105] M. Hindmarsh and M. Hijazi, _Gravitational waves from first order cosmological phase transitions in the Sound Shell Model_ , _JCAP_ 12 (2019) 062, [1909.10040].
  * [106] H.-K. Guo, K. Sinha, D. Vagie and G. White, _Phase Transitions in an Expanding Universe: Stochastic Gravitational Waves in Standard and Non-Standard Histories_ , _JCAP_ 01 (2021) 001, [2007.08537].
  * [107] X. Wang, F. P. Huang and Y. Li, _Sound velocity effects on the phase transition gravitational wave spectrum in the Sound Shell Model_ , 2112.14650.
  * [108] J. D. Jackson, _Classical Electrodynamics_.  Wiley, 1998.
  * [109] P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia and A. Urbano, _Weak Corrections are Relevant for Dark Matter Indirect Detection_ , _JCAP_ 03 (2011) 019, [1009.0224].
  * [110] J. Chen, T. Han and B. Tweedie, _Electroweak Splitting Functions and High Energy Showering_ , _JHEP_ 11 (2017) 093, [1611.00788].
  * [111] J. Chen, _On the Feynman Rules of Massive Gauge Theory in Physical Gauges_ , 1902.06738.
  * [112] I. Feige and M. D. Schwartz, _An on-shell approach to factorization_ , _Phys. Rev. D_ 88 (2013) 065021, [1306.6341].
  * [113] M. Abramowitz, I. A. Stegun and R. H. Romer, _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , 1988.
  * [114] B. C. Hall, _Quantum theory for mathematicians_ , vol. 267.  Springer, 2013.
  * [115] F. A. Berends, R. Kleiss, P. De Causmaecker, R. Gastmans and T. T. Wu, _Single Bremsstrahlung Processes in Gauge Theories_ , _Phys. Lett. B_ 103 (1981) 124–128.
  * [116] R.-G. Cai and S.-J. Wang, _Effective picture of bubble expansion_ , _JCAP_ 03 (2021) 096, [2011.11451].
  * [117] R. Jinno, T. Konstandin and M. Takimoto, _Relativistic bubble collisions—a closer look_ , _JCAP_ 09 (2019) 035, [1906.02588].