$\displaystyle{\textstyle\left.\left(\prod_{\sigma=\pm}\prod_{i=1}^{2}\Gamma\left(\frac{h-t_{2i-1\,2i}+\sigma
i\nu_{i}}{2}\right)\right)\,_{3}F_{2}\left({}_{\frac{2+h-t_{12}-\Delta_{5}-i\nu_{2}}{2}\,,\frac{h-t_{12}+\Delta_{5}-i\nu_{2}}{2}}^{\frac{h-t_{12}-i\nu_{1}}{2}\,,\frac{h-t_{12}+i\nu_{1}}{2}\,,\frac{h-t_{34}-i\nu_{2}}{2}};1\right)\right]}$

where we use the notation
$\delta_{ij}=\frac{\Delta_{i}+\Delta_{j}-t_{ij}}{2}$.

Let us finish this appendix by noting that a similar computation can be
performed for spinning exchanges using Schwinger parametrization (139). To do
so, at each moment, we multinomially expand the integrand decomposing it into
sums over integrands of similar form to the ones encountered for scalar
exchanges. In the end, one finds a spinning Mellin partial wave written as
several sums over scalar-type Mellin partial waves. In particular the sums are
bounded by the values of spin of the exchanged operators. This is no-good for
an analytic continuation in spin that we want to consider here. For that
reason and due to its length we do not write that result here.

##  Appendix C Explicit examples in position space

In this appendix, we single-out a conformal block contribution in position
space and compute its Regge-limit behavior.

We start with five-point conformal block lightcone limit in its integral
representation

$\displaystyle
G_{k_{1}k_{2}\ell}(u_{i})=u_{1}^{\frac{\tau_{1}}{2}}u_{3}^{\frac{\tau_{2}}{2}}(1-u_{2})^{\ell}u_{5}^{\frac{\Delta_{\phi}}{2}}\int_{0}^{1}[dt_{1}][dt_{2}]$
(147)
$\displaystyle\frac{\big{(}1-t_{1}(1-u_{2})u_{4}-u_{2}u_{4}\big{)}^{J_{2}-\ell}\big{(}1-t_{2}(1-u_{2})u_{5}-u_{2}u_{5}\big{)}^{J_{1}-\ell}}{\big{(}1-(1-u_{4})t_{2}\big{)}^{\frac{h_{2}-\tau_{1}-2\ell+\Delta_{\phi}}{2}}\big{(}1-(1-u_{5})t_{1}\big{)}^{\frac{h_{1}-\tau_{2}-2\ell+\Delta_{\phi}}{2}}\big{(}1-(1-t_{1})(1-t_{2})(1-u_{2})\big{)}^{\frac{h_{1}+h_{2}-\Delta_{\phi}}{2}}}\,,$

where $\tau_{i}=\Delta_{i}-J_{i}$ is the twist and $h_{i}=\Delta_{i}+J_{i}$
the conformal spin of the $i$-th exchanged operator. The measure is given by
$[dt]=\frac{\Gamma(\Delta_{i}+J_{i})}{\Gamma^{2}(\frac{\Delta_{i}+J_{i}}{2})}(t(1-t))^{\frac{\Delta_{i}+J_{i}}{2}-1}$.

Generically, we do not know how to evaluate these integrals in terms of known
analytic functions. However, when the exponents in the denominator of the
integrand are integers, this is no longer the case.222222The package HyperInt
[67] is particularly useful to evaluate these integrals. As a matter of
example we consider the simple case of $\Delta_{i}=\Delta_{\phi}=2$ and
$J_{1}=J_{2}=\ell=0$. Note that this is just a choice and spinning cases would
also have a similar discussion but with longer explicit expressions. In this
case, equation (147) can be integrated and yields (apart from an overall
constant)

$\displaystyle\frac{u_{1}u_{3}u_{5}}{1-u_{5}+u_{4}\left(u_{2}u_{5}-1\right)}\left[\text{Li}_{2}\left(u_{2}u_{4}\right)-\text{Li}_{2}\left(u_{4}\right)+\text{Li}_{2}\left(u_{2}u_{5}\right)-\text{Li}_{2}\left(u_{5}\right)-\text{Li}_{2}\left(u_{2}\right)\right.$
(148)
$\displaystyle\left.-\log\left(1-u_{2}\right)\log\left(u_{2}\right)-\log\left(1-u_{4}\right)\log\left(u_{4}\right)-\log\left(1-u_{5}\right)\log\left(u_{5}\right)\right.$
$\displaystyle\left.+\log\left(u_{4}\right)\log\left(u_{5}\right)+\log\left(u_{2}u_{5}\right)\log\left(1-u_{2}u_{5}\right)+\log\left(u_{2}u_{4}\right)\log\left(1-u_{2}u_{4}\right)+\zeta(2)\right]\,.$

As we perform the analytic continuation from an Euclidean to double-Reggeon
exchange kinematics that we presented in (78), we cross block branch-cuts and
it mixes with other solutions of the Casimir equations. In particular, the
discontinuities of the block contain the leading contributions in the Regge
limit. Having an explicit expression to work with we can tell the full story.

As we perform the analytic continuation and as the lightcones are crossed,
pairs of operators become timelike separated and cross-ratios $u_{2},u_{4}$
and $u_{5}$ go around 0. Note then that we are indeed crossing branch-cuts of
the expression (148). In particular, we observe that only $\log$ terms in
(148) can contribute to the discontinuity as $u_{i}$ goes around 0 with
$\log(x)\to\log(x)\pm 2\pi i$. The actual sign one picks is determined by how
one moves around branch-cuts. As we reviewed in the main text, this depends on
the ordering of operators of the Wightman function we consider. As before,
here we take an ordering compatible with the time-ordering of Regge
kinematics, i.e. $\langle\phi_{4}\phi_{1}\phi_{2}\phi_{5}\phi_{3}\rangle$.
Taking this ordering and the associated $i\epsilon$-prescription, we can
perform the path continuation to Regge kinematics in our explicit-lightcone-
block contribution and observe its discontinuities concretely. This is plotted
in figure 12.232323In this plot we only considered the terms within the
brackets in (148). Note that only this part is relevant for the
discontinuities we want to study.

Figure 12: Discontinuities of lightcone block under analytic continuation
(78). In blue, the real part of the stripped-off lightcone block. In orange,
the real part of the block with $\log(u_{2})\to\log(u_{2})+2\pi i$. In green,
the previous with $\log(u_{4})\to\log(u_{4})-2\pi i$ and in red, the latter
with $\log(u_{5})\to\log(u_{5})-2\pi i$. On the right, a zoomed-in version of
the same plot. The plots are obtained with $\delta_{2}=0.73\delta_{1}$.

As we move according to the chosen path for analytic continuation, we observe
that the lightcone block (blue) has discontinuities. The first one can be
removed if one replaces $\log(u_{2})\to\log(u_{2})+2\pi i$ as shown by the
orange line. Clearly, this shows that the discontinuity of the lightcone block
is due to a logarithmic discontinuity in $u_{2}$. Similarly, when the orange
line has a discontinuity, there is a continuation provided by the green line.
The latter is defined from the former with the replacement
$\log(u_{4})\to\log(u_{4})-2\pi i$. We conclude that a discontinuity in
$u_{4}$ has taken place. The same is true for the red line which provides the
continuation of the green line once we take $\log(u_{5})\to\log(u_{5})-2\pi i$
and once again a discontinuity, this time in $u_{5}$, has to be considered.
This simple example shows in practice what we had already guessed: the
lightcone block has discontinuities associated with $u_{2},u_{4}$ and $u_{5}$
going around 0 and all of them are important. Let us then study the
discontinuities of (148) on these variables.

It is possible to use the integral representation of the lightcone block to
argue that there are no sequential discontinuities involving $u_{2}$, i.e.

$\text{Disc}_{u_{2}}\text{Disc}_{u_{4}\,\text{or}\,u_{5}}G_{k_{1},k_{2},\ell}=\text{Disc}_{u_{4}\,\text{or}\,u_{5}}\text{Disc}_{u_{2}}G_{k_{1},k_{2},\ell}=0\,.$
(149)

In the expression (148) this is straightforward to see as there are no
products of the type $\log(u_{2})\log(u_{4})$ or $\log(u_{2})\log(u_{5})$. As
it was stated in the main text and as we will see below, it is actually the
sum
$\text{Disc}_{u_{2}}G_{k_{1},k_{2},\ell}+\text{Disc}_{u_{5}}\text{Disc}_{u_{4}}G_{k_{1},k_{2},\ell}$
that dominates the Regge behavior of the correlation function.

The discontinuity of expression (148) as $u_{2}$ goes around 0 with fixed
$u_{4},u_{5}>0$ is given by

$\pm 2\pi
i\frac{u_{1}u_{3}u_{5}}{1-u_{5}+u_{4}(u_{2}u_{5}-1)}\log\left(\frac{1-u_{2}}{(1-u_{2}u_{4})(1-u_{2}u_{5})}\right)\,,$
(150)

which in the limit $u_{4},u_{5}\to 1$ with
$\chi_{2}=\frac{1-u_{2}}{(1-u_{4})(1-u_{5})}$ fixed simplifies to

$\pm 2\pi
i\frac{\sqrt{u_{1}}\sqrt{u_{3}}}{\left(\chi_{2}-1\right)\chi_{4}\chi_{5}}\log\left(\chi_{2}\right)\,,$
(151)

where we also use $\chi_{4}=\frac{1-u_{4}}{\sqrt{u_{3}}}$ and
$\chi_{5}=\frac{1-u_{5}}{\sqrt{u_{1}}}$ which approach infinity due to the
order of limits considered. This order of limits does not correspond to the
actual Regge limit: indeed, we will call this ordered limit a boundary
condition for Regge limit. The name simply follows from the fact that we use
it below as a boundary condition for a set of recursion relations where we
compute the Regge limit of a conformal block starting from the lightcone.
Note, moreover, that the scaling in both $u_{1}$ and $u_{3}$ in the expression
above agrees with the expected $u_{i}^{(1-J_{i})/2}$ of Regge limit. As stated
above we are indeed describing a double Reggeon exchange. This clearly
contrasts with the Euclidean OPE scaling, $u_{i}^{\Delta_{i}/2}$, manifesting
the difference between Regge and Euclidean kinematics. Perhaps a more striking
example would follow from considering a spinning case from the beginning. The
story is no different in those cases but the expressions grow considerably in
size. We also note the existence of a $\log$ term in the case at hand. We
point out that some other examples where the lightcone block can be integrated
do not have these contributions in the above limit. Its existence in this case
suggests however that a generic function for the discontinuity of the
lightcone block as $u_{2}$ goes around 0 must contain $\log$ terms when the
representation labels of the external and exchanged operators conspire in a
certain way.

We now consider the discontinuity in $u_{4}$ with fixed and positive
$u_{2},u_{5}$. This gives

$\pm 2\pi
i\frac{u_{1}u_{3}u_{5}}{1-u_{5}+u_{4}(u_{2}u_{5}-1)}\log\left(\frac{1-u_{4}}{(1-u_{2}u_{4})u_{5}}\right)\,,$
(152)

which yields a boundary condition for Regge limit

$\pm 2\pi i\frac{u_{1}\sqrt{u_{3}}}{\chi_{4}}\,.$ (153)

From the symmetry of (148) between $u_{4}$ and $u_{5}$ we immediately see that
a similar result follows for the discontinuity in $u_{5}$. Note that these
terms are subdominant in the limit $u_{1},u_{3}\to 0$ when compared to (151).
In particular, in expression (153) $u_{1}$ scales as $u_{1}^{\Delta_{1}/2}$
whereas $u_{3}$ scales as $u_{3}^{(1-J_{2})/2}$. The converse happens in the
discontinuity in $u_{5}$ complex plane. This behavior should correspond to
single Reggeon exchanges. Notably, the sequential discontinuity in $u_{4}$ and
$u_{5}$ produces a dominant contribution for the double Reggeon kinematics. To
see this, consider (152) and take the sequential discontinuity in $u_{5}$.
This gives

$\pm
4\pi^{2}\frac{u_{1}u_{3}u_{5}}{1-u_{5}+u_{4}\left(u_{2}u_{5}-1\right)}\,,$
(154)

which fixes the boundary condition for Regge limit

$4\pi^{2}\frac{\sqrt{u_{1}}\sqrt{u_{3}}}{\left(\chi_{2}-1\right)\chi_{4}\chi_{5}}\,,$
(155)

that is as dominant as (151). We conclude in this simple example, the generic
statement we have made in the main text that
$\text{Disc}_{u_{2}}G_{k_{1},k_{2},\ell}+\text{Disc}_{u_{5}}\text{Disc}_{u_{4}}G_{k_{1},k_{2},\ell}$
provide the dominant contributions of the correlation function in Regge limit.

Even though the existence of $\log$ terms in the boundary condition for the
Regge limit is not generic, we should however show how to deal with them when
we compute the conformal blocks at Regge limit. If there are no $\log$ terms
in your case of interest, simply set those terms to zero in the procedure
below. We consider the Casimir equations in the limit of $u_{1},u_{3}\to 0$
with a block that scales as

$\mathcal{G}_{k_{1}^{\prime}k_{2}^{\prime}\ell}(x_{i})\propto
u_{1}^{\frac{1-J_{1}}{2}}u_{3}^{\frac{1-J_{2}}{2}}\mathcal{H}(\chi_{2},\chi_{4},\chi_{5})\,\,.$
(156)

In this limit, the Casimir equations for $\mathcal{H}$ simplify and read

$\displaystyle\left[\chi_{4}^{2}\left(4(2\chi_{2}-1)(\partial_{\chi_{2}}-\chi_{5}\partial_{\chi_{2}}\partial_{\chi_{5}})-(d-1)\chi_{5}^{3}\partial_{\chi_{5}}-\left(\chi_{5}^{2}-4\right)\chi_{5}^{2}\partial_{\chi_{5}}^{2}\right)\right.$
$\displaystyle\left.+4\left((\chi_{2}-1)\chi_{2}\chi_{4}^{2}+1\right)\partial_{\chi_{2}}^{2}+(\Delta_{1}-1)(\Delta_{1}-d+1)\chi_{4}^{2}\chi_{5}^{2}\right]\mathcal{H}(\chi_{2},\chi_{4},\chi_{5})=0$
(157)

with an entirely similar second equation obtained from the above by replacing
$\Delta_{1}$ by $\Delta_{2}$ and permuting the roles of $\chi_{4}$ and
$\chi_{5}$. In our particular case of study, in the limit of $u_{1},u_{3}\to
0$, large $\chi_{4},\chi_{5}$ and fixed $\chi_{2}$, the leading Regge
contribution of the block behaves as

$\frac{\sqrt{u_{1}}\sqrt{u_{3}}}{(\chi_{2}-1)\chi_{4}\chi_{5}}\left(a+b\log(\chi_{2})\right)$
(158)

where $a$ and $b$ are constants. We can thus further impose in (C)

$\mathcal{H}(\chi_{2},\chi_{4},\chi_{5})=\frac{\mathbb{H}(\chi_{2},\chi_{4},\chi_{5})}{\chi_{4}\chi_{5}}\,.$
(159)

According to (158) and considering a small-$\chi_{2}$ limit, we can look for
solutions of the Casimir equations of the form

$\mathbb{H}(\chi_{2},\chi_{4},\chi_{5})=\sum_{n_{1},n_{2},n_{3}}a_{n_{1},n_{2},n_{3}}\chi_{2}^{n_{1}}\chi_{4}^{-n_{2}}\chi_{5}^{-n_{3}}+b_{n_{1},n_{2},n_{3}}\log(\chi_{2})\chi_{2}^{n_{1}}\chi_{4}^{-n_{2}}\chi_{5}^{-n_{3}}$
(160)

where the coefficients $a_{n_{1},n_{2},n_{3}}$ and $b_{n_{1},n_{2},n_{3}}$
reduce to $a$ and $b$, respectively, when all $n_{i}$ are 0. The remaining
expansion coefficients $a_{n_{1},n_{2},n_{3}}$ and $b_{n_{1},n_{2},n_{3}}$ are
fixed by the Casimir equations. It is easy to see that this ansatz gives rise
to terms in the Casimir equations of the form

$\chi_{2}^{c_{1}+n_{1}}\chi_{4}^{c_{2}-n_{2}}\chi_{5}^{c_{3}-n_{3}}\times\begin{cases}a_{n_{1},n_{2},n_{3}}\\\
b_{n_{1},n_{2},n_{3}}\\\ b_{n_{1},n_{2},n_{3}}\log(\chi_{2})\end{cases}\,.$
(161)

Clearly, the terms that depend on $\log$ should cancel among each other in
order to satisfy the Casimir equation. This leads to two constraints per
Casimir equation, one for the $\log$-dependent terms and one for the
remaining. For the isolated $\log$ terms, we find recursion relations for the
coefficients by removing the $\chi$-dependence from the equations. To do so,
we shift each term accordingly, i.e. $n_{1}\to n_{1}-c_{1},n_{2}\to
n_{2}+c_{2}$ and $n_{3}\to n_{3}+c_{3}$. This leads to the following recursion
relations

$\displaystyle b_{n_{1},n_{2},n_{3}}=\frac{1}{n_{3}(d-n_{3}-4)}$
$\displaystyle\left[4(n_{1}+1)\left((n_{1}+n_{3})b_{n_{1}+1,n_{2},n_{3}-2}-(n_{1}+2)b_{n_{1}+2,n_{2}-2,n_{3}-2}\right)\right.$
$\displaystyle\left.-4(n_{1}+n_{3}-1)(n_{1}+n_{3})b_{n_{1},n_{2},n_{3}-2}\right]$
(162)

with a similar one where we exchange the roles of $n_{3}$ and $n_{2}$.
Clearly, the above recursion relation cannot be used whenever $n_{3}=0$. In
such case, the other recursion relation can be used instead (and vice-versa).
An entirely similar argument follows for the non-$\log$-dependent terms. We
find the recursion relations

$\displaystyle
a_{n_{1},n_{2},n_{3}}=\frac{4}{n_{2}(4-d+n_{2})}\left[(n_{1}+n_{2}-1)(n_{1}+n_{2})a_{n_{1},n_{2}-2,n_{3}}-(n_{1}+1)(n_{1}+n_{2})a_{n_{1}+1,n_{2}-2,n_{3}}\right.$
$\displaystyle\left.+(n_{1}+1)(n_{1}+2)a_{n_{1}+2,n_{2}-2,n_{3}-2}+(2n_{1}+2n_{2}-1)b_{n_{1},n_{2}-2,n_{3}}-(2n_{1}+n_{2}+1)b_{n_{1}+1,n_{2}-2,n_{3}}\right.$
$\displaystyle\left.+(2n_{1}+3)b_{n_{1}+2,n_{2}-2,n_{3}-2}\right]$ (163)

with another equivalent relation where the roles of $n_{2}$ and $n_{3}$ are
swapped.

These recursion relations are only meaningful once one prescribes a boundary
condition. We impose that

$\displaystyle
a_{n_{1},n_{2},n_{3}}=b_{n_{1},n_{2},n_{3}}=0\qquad\text{if}\qquad
n_{1}<0\,\lor\,n_{2}<0\,\lor\,n_{3}<0\,,$ $\displaystyle
a_{0,0,0}=a\qquad\text{and}\qquad b_{0,0,0}=b\,.$ (164)

It is easy to check that these recursion relations fix the behavior of all the
coefficients up to those of the form $a_{n_{1},0,0}$ and $b_{n_{1},0,0}$ but
note that these can be read from (158) by expanding it on small $\chi_{2}$
limit.

##  Appendix D Other Regge kinematics

In this short appendix, we detail other possible Regge kinematics that we did
not explore in detail in this paper but that might be worth studying in the
future.

#### Single Reggeon exchange

Within a five-point function, one can consider a single Reggeon exchange. In
terms of Mandelstam invariants $s_{25}$ or $s_{45}$ of figure 1 only one of
the two becomes large. In the context of CFTs, this translates to having only
two operators, one in the first and one in the second Poincare patches,
approaching each other in such a way that there is only one cross-ratio going
to $0$ rather than two.

One possible analytic continuation that describes single Reggeon exchange is
given by

$\displaystyle
x_{1}=-r\left(\sinh(\delta_{1}),\cosh(\delta_{1}),\textbf{0}_{d-2}\right)\,$
$\displaystyle x_{3}=\left(0,1,\textbf{0}_{d-2}\right)\,$ $\displaystyle
x_{5}=(0,h_{1},h_{2},\textbf{0}_{d-3})$ $\displaystyle
x_{2}=r\left(\sinh(\delta_{2}),\cosh(\delta_{2}),\textbf{0}_{d-2}\right)\,$
$\displaystyle x_{4}=\left(0,-1,\textbf{0}_{d-2}\right)\,.$ (165)

with positive rapidities $\delta_{i}$ and $\textbf{0}_{d}$ denoting a
$d$-dimensional vector of zeros. In the large-rapidities limit, one can check
that $u_{1}\to 0$ and $u_{2},u_{5}\to 1$ with unfixed $u_{3}$ and $u_{4}$.
This agrees with the Euclidean OPE limit in the (12) channel. Again, we
emphasize that this limit is attained after branch-cuts are crossed and thus
in an intrinsically Lorentzian Regge sheet.

#### Six-point snowflake

The six-point conformal block of external scalars is known in the lightcone
limit in the snowflake topology [15]. Even though we did not attempt in this
paper to analyze the cut-structure of this block, we nonetheless write down an
analytic continuation prescription to achieve a Regge limit configuration that
is consistent at the level of the cross-ratios with the OPE on channels (12),
(34) and (56).

We use the set of 9 cyclic cross-ratios

$\displaystyle u_{1}=\frac{x_{12}^{2}x_{35}^{2}}{x_{13}^{2}x_{25}^{2}}\quad
u_{i+1}=u_{i}|_{x_{i}\to x_{i+1}}\quad\text{mod 6}\,$ $\displaystyle
U_{1}=\frac{x_{13}^{2}x_{46}^{2}}{x_{14}^{2}x_{36}^{2}}\quad
U_{i+1}=U_{i}|_{x_{i}\to x_{i+1}}\quad\text{mod 3}\,.$ (166)

In the snowflake OPE limit $u_{1},u_{3},u_{5}\to 0$ and the remaining all go
to 1. In the Regge limit one should reobtain the same limiting values of the
cross-ratios after some lightcones are crossed. We start with a totally
spacelike configuration and perform the analytic continuation

$\displaystyle
x_{1}=-r_{1}\left(\sinh(\delta_{1}),\cosh(\delta_{1}),\textbf{0}_{d-2}\right)\,$
$\displaystyle
x_{4}=\left(\sinh(\delta_{1}),-\cosh(\delta_{1}),\textbf{0}_{d-2}\right)\,$
$\displaystyle
x_{2}=r_{1}\left(\sinh(\delta_{2}),\cosh(\delta_{2}),\textbf{0}_{d-2}\right)\,$
$\displaystyle
x_{5}=(r_{2}\sinh(\delta_{3}),r_{3},h,r_{2}\cosh(\delta_{3}),\textbf{0}_{d-4})\,$
$\displaystyle
x_{3}=\left(-\sinh(\delta_{2}),\cosh(\delta_{2}),\textbf{0}_{d-2}\right)\,$
$\displaystyle
x_{6}=(-r_{2}\sinh(\delta_{3}),r_{4},h,-r_{2}\cosh(\delta_{3}),\textbf{0}_{d-4})\,.$
(167)

where one can see that we use 9 degrees of freedom. Note as well that for six-
point functions one can at most use the conformal symmetry to state that any
generic correlation function is related to one that lives in some half-
subspace in 4 dimensions [68]. Perhaps the most notorious difference in this
case is the need to boost a pair of points along some different plane. It is
easy to check however that this prescription indeed leads to the expected OPE
behavior for a snowflake six-point function.

## References

  * [1] Riccardo Rattazzi, Vyacheslav S. Rychkov, Erik Tonni and Alessandro Vichi  “Bounding scalar operator dimensions in 4D CFT”  In _JHEP_ 12, 2008, pp. 031  DOI: 10.1088/1126-6708/2008/12/031
  * [2] David Poland, Slava Rychkov and Alessandro Vichi  “The Conformal Bootstrap: Theory, Numerical Techniques, and Applications”  In _Rev. Mod. Phys._ 91, 2019, pp. 015002  DOI: 10.1103/RevModPhys.91.015002
  * [3] Sheer El-Showk et al.  “Solving the 3D Ising Model with the Conformal Bootstrap”  In _Phys. Rev. D_ 86, 2012, pp. 025022  DOI: 10.1103/PhysRevD.86.025022
  * [4] Filip Kos, David Poland, David Simmons-Duffin and Alessandro Vichi  “Precision Islands in the Ising and $O(N)$ Models”  In _JHEP_ 08, 2016, pp. 036  DOI: 10.1007/JHEP08(2016)036
  * [5] Luis F. Alday and Alexander Zhiboedov  “An Algebraic Approach to the Analytic Bootstrap”  In _JHEP_ 04, 2017, pp. 157  DOI: 10.1007/JHEP04(2017)157
  * [6] David Simmons-Duffin  “The Lightcone Bootstrap and the Spectrum of the 3d Ising CFT”  In _JHEP_ 03, 2017, pp. 086  DOI: 10.1007/JHEP03(2017)086
  * [7] Zohar Komargodski and Alexander Zhiboedov  “Convexity and Liberation at Large Spin”  In _JHEP_ 11, 2013, pp. 140  DOI: 10.1007/JHEP11(2013)140
  * [8] Vladimir Rosenhaus  “Multipoint Conformal Blocks in the Comb Channel”  In _JHEP_ 02, 2019, pp. 142  DOI: 10.1007/JHEP02(2019)142
  * [9] Vasco Gonçalves, Raul Pereira and Xinan Zhou  “$20^{\prime}$ Five-Point Function from $AdS_{5}\times S^{5}$ Supergravity”  In _JHEP_ 10, 2019, pp. 247  DOI: 10.1007/JHEP10(2019)247
  * [10] Ilija Buric et al.  “From Gaudin Integrable Models to $d$-dimensional Multipoint Conformal Blocks”  In _Phys. Rev. Lett._ 126.2, 2021, pp. 021602  DOI: 10.1103/PhysRevLett.126.021602
  * [11] Ilija Buric et al.  “Gaudin models and multipoint conformal blocks III: comb channel coordinates and OPE factorisation”  In _JHEP_ 06, 2022, pp. 144  DOI: 10.1007/JHEP06(2022)144
  * [12] Ilija Buric et al.  “Gaudin models and multipoint conformal blocks. Part II. Comb channel vertices in 3D and 4D”  In _JHEP_ 11, 2021, pp. 182  DOI: 10.1007/JHEP11(2021)182
  * [13] Ilija Buric et al.  “Gaudin models and multipoint conformal blocks: general theory”  In _JHEP_ 10, 2021, pp. 139  DOI: 10.1007/JHEP10(2021)139
  * [14] David Poland and Valentina Prilepina  “Recursion relations for 5-point conformal blocks”  In _JHEP_ 10, 2021, pp. 160  DOI: 10.1007/JHEP10(2021)160
  * [15] António Antunes, Miguel S. Costa, Vasco Goncalves and Joao Vilas Boas  “Lightcone bootstrap at higher points”  In _JHEP_ 03, 2022, pp. 139  DOI: 10.1007/JHEP03(2022)139
  * [16] Carlos Bercini, Vasco Gonçalves and Pedro Vieira  “Light-Cone Bootstrap of Higher Point Functions and Wilson Loop Duality”  In _Phys. Rev. Lett._ 126.12, 2021, pp. 121603  DOI: 10.1103/PhysRevLett.126.121603
  * [17] Luis F. Alday, Vasco Gonçalves and Xinan Zhou  “Supersymmetric Five-Point Gluon Amplitudes in AdS Space”  In _Phys. Rev. Lett._ 128.16, 2022, pp. 161601  DOI: 10.1103/PhysRevLett.128.161601
  * [18] Vasco Gonçalves et al.  “Kaluza-Klein Five-Point Functions from $\textrm{AdS}_{5}\times S_{5}$ Supergravity”, 2023  arXiv:2302.01896 [hep-th]
  * [19] David Poland, Valentina Prilepina and Petar Tadić  “The five-point bootstrap”, 2023  arXiv:2305.08914 [hep-th]
  * [20] Lorenzo Cornalba, Miguel S. Costa, Joao Penedones and Ricardo Schiappa  “Eikonal Approximation in AdS/CFT: Conformal Partial Waves and Finite N Four-Point Functions”  In _Nucl. Phys. B_ 767, 2007, pp. 327–351  DOI: 10.1016/j.nuclphysb.2007.01.007
  * [21] Lorenzo Cornalba, Miguel S. Costa and Joao Penedones  “Eikonal approximation in AdS/CFT: Resumming the gravitational loop expansion”  In _JHEP_ 09, 2007, pp. 037  DOI: 10.1088/1126-6708/2007/09/037
  * [22] Miguel S. Costa, Vasco Goncalves and Joao Penedones  “Conformal Regge theory”  In _Jhep_ 12, 2012, pp. 091  DOI: 10.1007/jhep12(2012)091
  * [23] Simon Caron-Huot  “Analyticity in Spin in Conformal Theories”  In _Jhep_ 09, 2017, pp. 078  DOI: 10.1007/jhep09(2017)078
  * [24] António Antunes et al.  “The perturbative CFT optical theorem and high-energy string scattering in AdS at one loop”  In _JHEP_ 04, 2021, pp. 088  DOI: 10.1007/JHEP04(2021)088
  * [25] Simon Caron-Huot et al.  “Detectors in weakly-coupled field theories”  In _JHEP_ 04, 2023, pp. 014  DOI: 10.1007/JHEP04(2023)014
  * [26] Nima Afkhami-Jeddi, Thomas Hartman, Sandipan Kundu and Amirhossein Tajdini  “Einstein gravity 3-point functions from conformal field theory”  In _JHEP_ 12, 2017, pp. 049  DOI: 10.1007/JHEP12(2017)049
  * [27] Miguel S. Costa, Tobias Hansen and João Penedones  “Bounds for OPE coefficients on the Regge trajectory”  In _JHEP_ 10, 2017, pp. 197  DOI: 10.1007/JHEP10(2017)197
  * [28] Manuela Kulaxizi, Andrei Parnachev and Alexander Zhiboedov  “Bulk Phase Shift, CFT Regge Limit and Einstein Gravity”  In _JHEP_ 06, 2018, pp. 121  DOI: 10.1007/JHEP06(2018)121
  * [29] Daliang Li, David Meltzer and David Poland  “Conformal Bootstrap in the Regge Limit”  In _JHEP_ 12, 2017, pp. 013  DOI: 10.1007/JHEP12(2017)013
  * [30] Simon Caron-Huot, Dalimil Mazac, Leonardo Rastelli and David Simmons-Duffin  “AdS bulk locality from sharp CFT bounds”  In _JHEP_ 11, 2021, pp. 164  DOI: 10.1007/JHEP11(2021)164
  * [31] Pulkit Agarwal, Richard C. Brower, Timothy G. Raben and Chung-I Tan  “Embedding Space Approach to Lorentzian CFT Amplitudes and Causal Spherical Functions”, 2023  arXiv:2302.06469 [hep-th]
  * [32] T. Regge  “Introduction to complex orbital momenta”  In _Nuovo Cim._ 14, 1959, pp. 951  DOI: 10.1007/BF02728177
  * [33] P. Goddard and A. R. White  “Complex helicity and the sommerfeld-watson transformation of group-theoretic expansions”  In _Nuovo Cim. A_ 1, 1971, pp. 645–679  DOI: 10.1007/bf02734390
  * [34] A. R. White  “The signatured froissart-gribov continuation of multiparticle amplitudes to complex helicity and angular momentum”  In _Nucl. Phys. B_ 39, 1972, pp. 432–460  DOI: 10.1016/0550-3213(72)90381-1
  * [35] A. R. White  “The analytic continuation of multiparticle unitarity in the complex angular momentum plane”  In _Nucl. Phys. B_ 39, 1972, pp. 461–478  DOI: 10.1016/0550-3213(72)90382-3
  * [36] A. R. White  “Signature, factorization and unitarity in multi-regge theory \- the five-point function”  In _Nucl. Phys. B_ 67, 1973, pp. 189–231  DOI: 10.1016/0550-3213(73)90325-8
  * [37] Alan R. White  “The Analytic Foundations of Regge Theory”  In _Institute on Structural Analysis of Multiparticle Collision Amplitudes in Relativistic Quantum Theory_ , 1976
  * [38] R. C. Brower, Carleton E. DeTar and J. H. Weis  “Regge Theory for Multiparticle Amplitudes”  In _Phys. Rept._ 14, 1974, pp. 257  DOI: 10.1016/0370-1573(74)90012-x
  * [39] J. H. Weis  “Singularities in complex angular momentum and helicity”  In _Phys. Rev. D_ 6, 1972, pp. 2823–2841  DOI: 10.1103/PhysRevD.6.2823
  * [40] H. D. I. Abarbanel and A. Schwimmer  “Analytic structure of multiparticle amplitudes in complex helicity”  In _Phys. Rev. D_ 6, 1972, pp. 3018–3031  DOI: 10.1103/PhysRevD.6.3018
  * [41] Henry P. Stapp and Alan R. White  “An Asymptotic Dispersion Relation for the Six Particle Amplitude”  In _Phys. Rev. D_ 26, 1982, pp. 2145  DOI: 10.1103/PhysRevD.26.2145
  * [42] I. T. Drummond  “Multi-Reggeon Behavior of Production Amplitudes”  In _Phys. Rev._ 176 American Physical Society, 1968, pp. 2003–2013  DOI: 10.1103/PhysRev.176.2003
  * [43] Alan White  “Analytic Multi-regge Theory And The Pomeron In Qcd - Part 1” [Online; accessed 15-November-2022], https://www.osti.gov/servlets/purl/7002400
  * [44] C. E. DeTar et al.  “Helicity Poles, Triple-Regge Behavior, and Single-Particle Spectra in High-Energy Collisions”  In _Phys. Rev. Lett._ 26 American Physical Society, 1971, pp. 675–676  DOI: 10.1103/PhysRevLett.26.675
  * [45] Steinmann  “Ueber den Zusammenhang zwischen den Wightmanfunktionen und den retardierten Kommutatoren” [Online; accessed 15-November-2022], https://www.research-collection.ethz.ch/handle/20.500.11850/135473
  * [46] John H. Schwarz  “Dual resonance theory”  In _Phys. Rept._ 8, 1973, pp. 269–335  DOI: 10.1016/0370-1573(73)90003-3
  * [47] K. Bardakci and H. Ruegg  “Meson resonance couplings in a five-point veneziano model”  In _Phys. Lett. B_ 28, 1969, pp. 671–675  DOI: 10.1016/0370-2693(69)90219-6
  * [48] Carlos Bercini, Vasco Gonçalves, Alexandre Homrich and Pedro Vieira  “The Wilson loop — large spin OPE dictionary”  In _Jhep_ 07, 2022, pp. 079  DOI: 10.1007/jhep07(2022)079
  * [49] Juan Maldacena, David Simmons-Duffin and Alexander Zhiboedov  “Looking for a bulk point”  In _JHEP_ 01, 2017, pp. 013  DOI: 10.1007/JHEP01(2017)013
  * [50] Matthijs Hogervorst and Slava Rychkov  “Radial Coordinates for Conformal Blocks”  In _Phys. Rev. D_ 87, 2013, pp. 106004  DOI: 10.1103/PhysRevD.87.106004
  * [51] Simon Caron-Huot and Joshua Sandor  “Conformal Regge Theory at Finite Boost”  In _Jhep_ 05, 2021, pp. 059  DOI: 10.1007/jhep05(2021)059
  * [52] Petr Kravchuk, Jiaxin Qiao and Slava Rychkov  “Distributions in CFT. Part I. Cross-ratio space”  In _Jhep_ 05, 2020, pp. 137  DOI: 10.1007/jhep05(2020)137
  * [53] Petr Kravchuk, Jiaxin Qiao and Slava Rychkov  “Distributions in CFT. Part II. Minkowski space”  In _Jhep_ 08, 2021, pp. 094  DOI: 10.1007/jhep08(2021)094
  * [54] M. Luscher and G. Mack  “Global Conformal Invariance in Quantum Field Theory”  In _Commun. Math. Phys._ 41, 1975, pp. 203–234  DOI: 10.1007/bf01608988
  * [55] G. Mack  “Convergence of Operator Product Expansions on the Vacuum in Conformal Invariant Quantum Field Theory”  In _Commun. Math. Phys._ 53, 1977, pp. 155  DOI: 10.1007/bf01609130
  * [56] R. F. Streater and A. S. Wightman  “PCT, spin and statistics, and all that”, 1989
  * [57] H. Casini  “Wedge reflection positivity”  In _J. Phys. A_ 44, 2011, pp. 435202  DOI: 10.1088/1751-8113/44/43/435202
  * [58] Juan Maldacena, Stephen H. Shenker and Douglas Stanford  “A bound on chaos”  In _Jhep_ 08, 2016, pp. 106  DOI: 10.1007/jhep08(2016)106
  * [59] Murat Kologlu, Petr Kravchuk, David Simmons-Duffin and Alexander Zhiboedov  “Shocks, Superconvergence, and a Stringy Equivalence Principle”  In _Jhep_ 11, 2020, pp. 096  DOI: 10.1007/jhep11(2020)096
  * [60] David Simmons-Duffin, Douglas Stanford and Edward Witten  “A spacetime derivation of the Lorentzian OPE inversion formula”  In _Jhep_ 07, 2018, pp. 085  DOI: 10.1007/jhep07(2018)085
  * [61] Gerhard Mack  “D-independent representation of Conformal Field Theories in D dimensions via transformation to auxiliary Dual Resonance Models. Scalar amplitudes”, 2009  arXiv:0907.2407 [hep-th]
  * [62] Vasco Gonçalves, João Penedones and Emilio Trevisani  “Factorization of Mellin amplitudes”  In _Jhep_ 10, 2015, pp. 040  DOI: 10.1007/jhep10(2015)040
  * [63] Miguel S. Costa, Vasco Gonçalves and João Penedones  “Spinning AdS Propagators”  In _Jhep_ 09, 2014, pp. 064  DOI: 10.1007/jhep09(2014)064
  * [64] Petr Kravchuk and David Simmons-Duffin  “Light-ray operators in conformal field theory”  In _Jhep_ 11, 2018, pp. 102  DOI: 10.1007/jhep11(2018)102
  * [65] Joao Penedones  “Writing CFT correlation functions as AdS scattering amplitudes”  In _Jhep_ 03, 2011, pp. 025  DOI: 10.1007/jhep03(2011)025
  * [66] K. Symanzik  “On Calculations in conformal invariant field theories”  In _Lett. Nuovo Cim._ 3, 1972, pp. 734–738  DOI: 10.1007/BF02824349
  * [67] Erik Panzer  “Algorithms for the symbolic integration of hyperlogarithms with applications to Feynman integrals”  In _Comput. Phys. Commun._ 188, 2015, pp. 148–166  DOI: 10.1016/j.cpc.2014.10.019
  * [68] Petr Kravchuk and David Simmons-Duffin  “Counting Conformal Correlators”  In _Jhep_ 02, 2018, pp. 096  DOI: 10.1007/jhep02(2018)096
  * [69] Thomas Hartman, Sachin Jain and Sandipan Kundu  “Causality Constraints in Conformal Field Theory”  In _JHEP_ 05, 2016, pp. 099  DOI: 10.1007/JHEP05(2016)099
  * [70] Jean-François Fortin, Wenjie Ma and Witold Skiba  “Higher-Point Conformal Blocks in the Comb Channel”  In _JHEP_ 07, 2020, pp. 213  DOI: 10.1007/JHEP07(2020)213
  * [71] Jean-François Fortin, Wen-Jie Ma and Witold Skiba  “Six-point conformal blocks in the snowflake channel”  In _JHEP_ 11, 2020, pp. 147  DOI: 10.1007/JHEP11(2020)147
  * [72] Jean-François Fortin, Wen-Jie Ma and Witold Skiba  “Seven-point conformal blocks in the extended snowflake channel and beyond”  In _Phys. Rev. D_ 102.12, 2020, pp. 125007  DOI: 10.1103/PhysRevD.102.125007
  * [73] Sarthak Parikh  “A multipoint conformal block chain in $d$ dimensions”  In _JHEP_ 05, 2020, pp. 120  DOI: 10.1007/JHEP05(2020)120
  * [74] Till Bargheer, Vsevolod Chestnov and Volker Schomerus  “The Multi-Regge Limit from the Wilson Loop OPE”  In _JHEP_ 05, 2020, pp. 002  DOI: 10.1007/JHEP05(2020)002
  * [75] Richard C. Brower, Horatiu Nastase, Howard J. Schnitzer and Chung-I. Tan  “Implications of multi-Regge limits for the Bern-Dixon-Smirnov conjecture”  In _Nucl. Phys. B_ 814, 2009, pp. 293–326  DOI: 10.1016/j.nuclphysb.2009.02.009
  * [76] V. Del Duca et al.  “All-order amplitudes at any multiplicity in the multi-Regge limit”  In _Phys. Rev. Lett._ 124.16, 2020, pp. 161602  DOI: 10.1103/PhysRevLett.124.161602
  * [77] Vittorio Del Duca, Claude Duhr, Falko Dulat and Brenda Penante  “All two-loop MHV remainder functions in multi-Regge kinematics”  In _JHEP_ 01, 2019, pp. 162  DOI: 10.1007/JHEP01(2019)162
  * [78] Vittorio Del Duca et al.  “Amplitudes in the Multi-Regge Limit of $\mathcal{N}$=4 SYM”  In _Acta Phys. Polon. Supp._ 12.4, 2019, pp. 961–966  DOI: 10.5506/APhysPolBSupp.12.961
  * [79] Vittorio Del Duca et al.  “The seven-gluon amplitude in multi-Regge kinematics beyond leading logarithmic accuracy”  In _JHEP_ 06, 2018, pp. 116  DOI: 10.1007/JHEP06(2018)116
  * [80] Vittorio Del Duca et al.  “Multi-Regge kinematics and the moduli space of Riemann spheres with marked points”  In _JHEP_ 08, 2016, pp. 152  DOI: 10.1007/JHEP08(2016)152
  * [81] Simon Caron-Huot et al.  “Multi-Regge Limit of the Two-Loop Five-Point Amplitudes in $\mathcal{N}=4$ Super Yang-Mills and $\mathcal{N}=8$ Supergravity”  In _JHEP_ 10, 2020, pp. 188  DOI: 10.1007/JHEP10(2020)188
  * [82] S. Donnachie, Hans Gunter Dosch, O. Nachtmann and P. Landshoff  “Pomeron physics and QCD”  Cambridge University Press, 2004
  * [83] V. N. Gribov  “Possible Asymptotic Behavior of Elastic Scattering”  In _JETP Lett._ 41, 1961, pp. 667–669
  * [84] J. Bartels, V. Schomerus and M. Sprenger  “Multi-Regge Limit of the n-Gluon Bubble Ansatz”  In _Jhep_ 11, 2012, pp. 145  DOI: 10.1007/jhep11(2012)145
  * [85] N. E. J. Bjerrum-Bohr, Poul H. Damgaard, Thomas Sondergaard and Pierre Vanhove  “The Momentum Kernel of Gauge and Gravity Theories”  In _Jhep_ 01, 2011, pp. 001  DOI: 10.1007/jhep01(2011)001
  * [86] Richard C. Brower, Joseph Polchinski, Matthew J. Strassler and Chung-I Tan  “The Pomeron and gauge/string duality”  In _Jhep_ 12, 2007, pp. 005  DOI: 10.1088/1126-6708/2007/12/005
  * [87] Deeksha Chandorkar, Subham Dutta Chowdhury, Suman Kundu and Shiraz Minwalla  “Bounds on Regge growth of flat space scattering from bounds on chaos”  In _Jhep_ 05, 2021, pp. 143  DOI: 10.1007/jhep05(2021)143
  * [88] Kevin E. Cahill and Henry P. Stapp  “Optical Theorems And Steinmann Relations”  In _Annals Phys._ 90, 1975, pp. 438  DOI: 10.1016/0003-4916(75)90006-8
  * [89] Denis Karateev, Petr Kravchuk and David Simmons-Duffin  “Harmonic Analysis and Mean Field Theory”  In _Jhep_ 10, 2019, pp. 217  DOI: 10.1007/jhep10(2019)217
  * [90] Simon Caron-Huot, Lance J. Dixon, Andrew McLeod and Matt Hippel  “Bootstrapping a Five-Loop Amplitude Using Steinmann Relations”  In _Phys. Rev. Lett._ 117.24, 2016, pp. 241601  DOI: 10.1103/PhysRevLett.117.241601
  * [91] David Simmons-Duffin  “TASI lectures on Conformal Field Thoery” [Online; accessed 15-November-2022], https://gitlab.com/davidsd/lorentzian-cft-notes/-/blob/master/TASI-Lorentzian-CFT.pdf
  * [92] Miguel F. Paulos, Marcus Spradlin and Anastasia Volovich  “Mellin Amplitudes for Dual Conformal Integrals”  In _Jhep_ 08, 2012, pp. 072  DOI: 10.1007/jhep08(2012)072
  * [93] I. T. Drummond, P. V. Landshoff and W. J. Zakrzewski  “Signature in production amplitudes”  In _Phys. Lett. B_ 28, 1969, pp. 676–678  DOI: 10.1016/0370-2693(69)90220-2
  * [94] I. O. Moen  “Factorization of the Six Particle Multi-Regge Amplitude”  In _Phys. Rev. D_ 12, 1975, pp. 2486  DOI: 10.1103/PhysRevD.12.2486
  * [95] Dean Carmi and Simon Caron-Huot  “A Conformal Dispersion Relation: Correlations from Absorption”  In _JHEP_ 09, 2020, pp. 009  DOI: 10.1007/JHEP09(2020)009
  * [96] Joao Penedones, Joao A. Silva and Alexander Zhiboedov  “Nonperturbative Mellin Amplitudes: Existence, Properties, Applications”  In _JHEP_ 08, 2020, pp. 031  DOI: 10.1007/JHEP08(2020)031


Ginwidth=,height=,keepaspectratio mm[#1]

# Caustics: A Python Package for Accelerated Strong Gravitational Lensing
Simulations

Connor Stone  Ciela Institute - Montréal Institute for Astrophysical Data
Analysis and Machine Learning, Montréal, Québec, Canada  Department of
Physics, Université de Montréal, Montréal, Québec, Canada  0000-0002-9086-6398
Alexandre Adam  Ciela Institute - Montréal Institute for Astrophysical Data
Analysis and Machine Learning, Montréal, Québec, Canada  Department of
Physics, Université de Montréal, Montréal, Québec, Canada  0000-0001-8806-7936
Adam Coogan Work done while at UdeM, Ciela, and Mila Ciela Institute -
Montréal Institute for Astrophysical Data Analysis and Machine Learning,
Montréal, Québec, Canada  Department of Physics, Université de Montréal,
Montréal, Québec, Canada  0000-0002-0055-1780 M. J. Yantovski-Barth  Ciela
Institute - Montréal Institute for Astrophysical Data Analysis and Machine
Learning, Montréal, Québec, Canada  Department of Physics, Université de
Montréal, Montréal, Québec, Canada  Mila - Québec Artificial Intelligence
Institute, Montréal, Québec, Canada  0000-0001-5200-4095 Andreas Filipp  Ciela
Institute - Montréal Institute for Astrophysical Data Analysis and Machine
Learning, Montréal, Québec, Canada  Department of Physics, Université de
Montréal, Montréal, Québec, Canada  Mila - Québec Artificial Intelligence
Institute, Montréal, Québec, Canada  0000-0003-4701-3469 Landung Setiawan
eScience Institute Scientific Software Engineering Center, 1410 NE Campus
Pkwy, Seattle, WA 98195, USA  0000-0002-1624-2667 Cordero Core  eScience
Institute Scientific Software Engineering Center, 1410 NE Campus Pkwy,
Seattle, WA 98195, USA  0000-0002-3531-3221 Ronan Legin  Ciela Institute -
Montréal Institute for Astrophysical Data Analysis and Machine Learning,
Montréal, Québec, Canada  Department of Physics, Université de Montréal,
Montréal, Québec, Canada  Mila - Québec Artificial Intelligence Institute,
Montréal, Québec, Canada  0000-0001-9459-6316 Charles Wilson  Ciela Institute
- Montréal Institute for Astrophysical Data Analysis and Machine Learning,
Montréal, Québec, Canada  Department of Physics, Université de Montréal,
Montréal, Québec, Canada  Mila - Québec Artificial Intelligence Institute,
Montréal, Québec, Canada  0000-0001-7071-5528 Gabriel Missael Barco  Ciela
Institute - Montréal Institute for Astrophysical Data Analysis and Machine
Learning, Montréal, Québec, Canada  Department of Physics, Université de
Montréal, Montréal, Québec, Canada  Mila - Québec Artificial Intelligence
Institute, Montréal, Québec, Canada  0009-0008-5839-5937 Yashar Hezaveh  Ciela
Institute - Montréal Institute for Astrophysical Data Analysis and Machine
Learning, Montréal, Québec, Canada  Department of Physics, Université de
Montréal, Montréal, Québec, Canada  Mila - Québec Artificial Intelligence
Institute, Montréal, Québec, Canada  Center for Computational Astrophysics,
Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA
0000-0002-8669-5733 Laurence Perreault-Levasseur  Ciela Institute - Montréal
Institute for Astrophysical Data Analysis and Machine Learning, Montréal,
Québec, Canada  Department of Physics, Université de Montréal, Montréal,
Québec, Canada  Mila - Québec Artificial Intelligence Institute, Montréal,
Québec, Canada  Center for Computational Astrophysics, Flatiron Institute, 162
5th Avenue, 10010, New York, NY, USA  0000-0003-3544-3939

(19 March 2024)

## Summary

Gravitational lensing is the deflection of light rays due to the gravity of
intervening masses. This phenomenon is observed in a variety of scales and
configurations, involving any non-uniform mass such as planets, stars,
galaxies, clusters of galaxies, and even the large scale structure of the
universe. Strong lensing occurs when the distortions are significant and
multiple images of the background source are observed. The lens objects must
align on the sky of order $\sim 1$ arcsecond for galaxy-galaxy lensing, or
10’s of arcseonds for cluster-galaxy lensing. As the discovery of lens systems
has grown to the low thousands, these systems have become pivotal for
precision measurements and addressing critical questions in astrophysics.
Notably, they facilitate the measurement of the Universe’s expansion rate
(e.g. ref-holycowK. C. Wong et al., 2020), dark matter (e.g. ref-
Hezaveh2016Hezaveh et al., 2016; ref-Vegetti2014Vegetti & Vogelsberger, 2014),
supernovae (e.g. ref-Rodney2021Rodney et al., 2021), quasars (e.g. ref-
Peng2006Peng et al., 2006), and the first stars (e.g. ref-Welch2022Welch et
al., 2022) among other topics. With future surveys expected to discover
hundreds of thousands of lensing systems, the modelling and simulation of such
systems must occur at orders of magnitude larger scale then ever before. Here
we present caustics, a Python package designed to handle the extensive
computational demands of modeling such a vast number of lensing systems.

## Statement of need

The next generation of astronomical surveys, such as the Large Synoptic Sky
Telescope survey and Euclid wide survey, are expected to uncover hundreds of
thousands of gravitational lenses (ref-Collett2015Collett, 2015), dramatically
increasing the scientific potential of gravitational lensing studies.
Currently, analyzing a single lensing system can take several days or weeks,
which will be infeasible as the number of known lenses increases by orders of
magnitude. Thus, advancements such as computational acceleration via GPUs
and/or algorithmic advances such as automatic differentiation are needed to
reduce the analysis timescales. Machine learning will be critical to achieve
the necessary speed to process these lenses, it will also be needed to meet
the complexity of strong lens modelling. caustics is built with the future of
lensing in mind, using PyTorch (ref-pytorchPaszke et al., 2019) to accelerate
the low level computation and enable deep learning algorithms which rely on
automatic differentiation. Automatic differentiation also benefits classical
algorithms such as Hamiltonian Monte Carlo (ref-hmcBetancourt, 2017). With
these tools available, caustics will provide greater than two orders of
magnitude acceleration to most standard operations, enabling previously
impractical analyses at scale.

Several other simulation packages for strong gravitational lensing are already
publicly available. The well established lenstronomy package has been in use
since 2018 (ref-lenstronomyBirrer et al., 2021); GLAMER is a C++ based code
for modelling complex and large dynamic range fields (ref-GLAMERMetcalf &
Petkova, 2014); PyAutoLens is also widely used (ref-PyAutoLensNightingale et
al., 2021); and GIGA-Lens is a specialized JAX (ref-JAXBradbury et al., 2018)
based gravitational lensing package (ref-GIGALensGu et al., 2022); among
others (ref-Galan2021Galan et al., 2021; ref-Keeton2011Keeton, 2011; ref-
Kneib2011Kneib et al., 2011; ref-Wagner2024Wagner-Carena et al., 2024). There
are also several “in house” codes developed for specialized analysis which are
then not publicly released (e.g. ref-Suyu2010Suyu & Halkola, 2010). Caustics
development has been primarily focused on three aspects: processing speed,
user experience, and flexibility. The code is optimized to fully exploit
PyTorch’s capabilities, significantly enhancing processing speed. The user
experience is streamlined by providing three interfaces to the code:
configuration file, object-oriented, and functional, where each interface
level requires more expertise but allows more capabilities. In this way, users
with all levels of gravitational lensing simulation experience may effectively
engage with the software. Flexibility is achieved by a determined focus on
minimalism in the core functionality of caustics. All of these elements
combine to make Caustics a capable lensing simulator to support machine
learning applications, and classical analysis.

Caustics fulfills a timely need for a differentiable lensing simulator.
Several other fields have already benefitted from such simulators, for
example: gravitational wave analysis (ref-Coogan2022Coogan et al., 2022; ref-
Edwards2023Edwards et al., 2023; ref-Wong2023K. W. K. Wong et al., 2023);
astronomical image photometry (ref-Stone2023Stone et al., 2023); point spread
function modelling (ref-Desdoigts2023Desdoigts et al., 2023); light curves
(ref-Millon2024Millon et al., 2024); and even generic optimization for
scientific problems (ref-Nikolic2018Nikolic, 2018). Caustics is built on
lessons from other differentiable codes, with the goal of enabling machine
learning techniques in the field of strong lensing. With caustics it will now
be possible to analyze over 100,000 lenses in a timely manner (ref-
Hezaveh2017Hezaveh et al., 2017; ref-Perreault2017Perreault Levasseur et al.,
2017).

## Scope

Caustics is a gravitational lensing simulator. The purpose of the project is
to streamline the simulation of strong gravitational lensing effects on the
light of a background source. This includes a variety of parametric lensing
profiles such as: Singular Isothermal Ellipsoid (SIE), Elliptical Power Law
(EPL), Pseudo-Jaffe, Navarro-Frenk-White (NFW), and External Shear.
Additionally, it offers non-parametric representations such as a gridded
convergence or potential field. For the background source caustics provides a
Sérsic light profile, as well as a pixelized light image. Users may easily
extend these lens and source lists using templates provided in the
documentation.

Once a lensing system has been defined, caustics may then perform various
computational operations on the system such as raytracing through the lensing
system, either forwards and backwards. Users may compute the lensing
potential, convergence, deflection field, time delay field, and magnification.
All of these operations may readily be performed in a multi-plane setting to
account for interlopers or multiple sources. Since the code is differentiable
(via PyTorch), one may easily also compute the derivatives of these
quantities, such as when finding critical curves or computing the Jacobian of
the lens equation. For example, one may project interlopers from multiple
lensing planes to a single plane by computing the effective convergence,
obtained with the trace of the Jacobian.

With these building blocks in place, one may construct fast and accurate
simulators used to produce training sets for machine learning models or for
inference on real-world systems. Neural networks have become a widespread tool
for amortized inference of gravitational lensing parameter (ref-
Hezaveh2017Hezaveh et al., 2017) or in the detection of gravitational lenses
(ref-Huang2021Huang et al., 2021; ref-Petrillo2017Petrillo et al., 2017), but
they require large and accurate training sets which can be created quickly
with caustics. A demonstration of such a simulator is given in Figure 1 which
also demonstrates the importance of sub-pixel sampling. This involves
raytracing through the lensing mass and extracting the brightness of the
background source. Further, the image then must be convolved with a PSF for
extra realism. All of these operations are collected into a single simulator
which users may access and use simply as a function of the relevant lensing
and light source parameters. Since caustics is written in PyTorch, its
simulators are differentiable, thus one can compute gradients through the
forward model. This enables machine learning algorithms such as recurrent
inference machines (ref-Adam2023Adam et al., 2023) and diffusion models (ref-
Adam2022Adam et al., 2022).

Figure 1: Example simulated gravitational lens system defined by a Sérsic
source, SIE lens mass, and Sérsic lens light. Left, the pixel map is sampled
only at the midpoint of each pixel. Middle, the pixel map is supersampled and
then integrated using gaussian quadrature integration for greater accuracy.
Right, the fractional difference between the two is shown. We can see that in
this case the midpoint sampling is inaccurate by up to 30% of the pixel value
in areas of high contrast. The exact inaccuracy depends greatly on the exact
configuration.

The current scope of caustics does not include weak lensing, microlensing, or
cluster scale lensing simulations. While the underlying mathematical
frameworks are similar, the specific techniques commonly used in these areas
are not yet implemented, though they represent an avenue for future
development.

Caustics’ defined scope ends at lensing simulation, thus it does not include
functionality to optimize or sample the resulting functions. Users are
encouraged to use already existing optimization and sampling codes like
scipy.optimize (ref-scipyVirtanen et al., 2020), emcee (ref-emceeForeman-
Mackey et al., 2013), dynesty (ref-dynestySpeagle, 2020), Pyro (ref-
pyroBingham et al., 2019), and torch.optim (ref-pytorchPaszke et al., 2019).

## Performance

Here we discuss the performance enhancements enabled by caustics. Via PyTorch,
the code allows operations to be batched, multi-threaded on CPUs, or offloaded
to GPUs to optimize computational efficiency. In Figure 2 we demonstrate this
by sampling images of a Sérsic with an SIE model lensing the image (much like
Figure 1). In the two subfigures we show performance for simply sampling a
128x128 image using the pixel midpoint (left), and sampling a “realistic”
image (right) which is upsampled by a factor of 4 and convolved with a PSF.
All parameters are randomly resampled for each simulation to avoid caching
effects. This demonstrates a number of interesting facts about numerical
performance in such scenarios.

We compare the performance with that of Lenstronomy as our baseline. The most
direct comparison between the two codes can be observed by comparing the
Lenstronomy line with the “caustics unbatched cpu” line. Lenstronomy is
written using the numba (ref-numbaLam et al., 2015) package which compiles
python code into lower level C code. The left plot shows that caustics suffers
from a significant overhead compared with Lenstronomy, which is nearly twice
as fast as the “caustics unbatched cpu” line. This occurs because the pure
Python (intepreted language) elements of caustics are much slower than the
C/Cuda PyTorch backends (compiled language). This is most pronounced when
fewer computations are needed to perform a simulation. Despite this overhead,
caustics showcases a strong performance when using the batched GPU setting,
especially in the more realistic scenario with extra computations in the
simulator including 4x oversampling of the raytracing and the PSF convolution.

Figure 2: Runtime comparisons for a simple lensing setup. We compare the
amount of time taken (y-axis) to generate a certain number of lensing
realizations (x-axis) where a Sérsic model is lensed by an SIE mass
distribution. For CPU calculations we use Intel Gold 6148 Skylake and for the
GPU we use a NVIDIA V100, all tests were done at 64 bit precision. On the
left, the lensing system is sampled 128 pixel resolution only at pixel
midpoints. On the right, a more realistic simulation includes upsampled pixels
and PSF convolution. From the two tests we see varying performance
enhancements from compiled, unbatched, batched, multi-threaded, and GPU
processing setups.

Comparing the “caustics unbatched cpu” and “caustics batched cpu” lines we see
that batching can provide more efficient use of the same, single CPU,
computational resources. However, in the realistic scenario the batching has
minimal performance enhancement, likely because the python overhead of a for-
loop is no longer significant compared to the large number of numerical
operations being performed.

Comparing “caustics batched cpu” and “caustics batched 4cpu” we see that
PyTorch’s automatic multi-threading capabilities can indeed provide
performance enhancements. However, the enhancement is not a direct multiple of
the number of CPUs due to overhead. For tasks that are “embarrassingly
parallel,” such as running multiple MCMC chains, it is more effective to
parallelize at the job level rather than at the thread level to avoid these
overheads.

The most dramatic improvements are observed when comparing any CPU operations
with “caustics batched gpu”. Although communication between the CPU and GPU
can be slow, consolidating calculations into fewer, larger batches allows
caustics to fully exploit GPU capabilities. In the midpoint sampling, the GPU
never “saturates” meaning that it runs equally fast for any number of samples.
In the realistic scenario we reach the saturation limit of the GPU memory at
100 samples and could no longer simultaneously model all the systems, we thus
entered a linear regime in runtime just like the CPU sampling does for any
number of simulations. Nonetheless, it is possible to easily achieve over 100X
speedup over CPU performance, making GPUs by far the most efficient method to
perform large lensing computations such as running many MCMC chains or
sampling many lensing realizations (e.g. for training machine learning
models).

## User experience

Caustics offers a tiered interface system designed to cater to users with
varying levels of expertise in gravitational lensing simulation. This section
outlines the three levels of interfaces that enhance user experience by
providing different degrees of complexity and flexibility.

Configuration file interface: The most accessible level of interaction is
through configuration files. Users can define simulators in .yaml format,
specifying parameters such as lens models, light source characteristics, and
image processing details like PSF convolution and sub-pixel integration. The
user may then load such a simulator in a single line of Python and carry on
using that simulator as a pure function f(x) which takes in parameters such as
the Sérsic index, position, SIE Einstein radius, etc. and returns an image.
This interface is straightforward for new users and for simplifying the
sharing of simulation configurations between users.

Object oriented interface: This intermediate level allows users to manipulate
lenses and light sources as objects. The user may build simulators just like
the configuration file interface, or they may interact with the objects in a
number of other ways accessing further details about each lens. Each lensing
object has (where meaningful) a convergence, potential, time delay, and
deflection field. We provide examples to visualize all of these. Users may
apply the full flexibility of Python with these lensing objects and construct
customized analysis code, though there are many default routines which enable
one to quickly perform typical analysis tasks.

For both the object oriented and .yaml interfaces, the final simulator object
can be analyzed in a number of ways, Figure 3 demonstrates how one can
investigate the structure of a simulator in the form of a directed acyclic
graph of calculations. Note that one may also fix a subset of parameter
values, making them “static” instead of the default which is “dynamic”. Users
can produce such a graph representation in a single line of Python for any
caustics simulator.

Figure 3: Example directed acyclic graph representation of the simulator from
Figure 2. Ellipses are caustics objects and squares are parameters; open
squares are dynamic parameters and greyed squares are static parameters.
Parameters are passed at the top level node (Lens_Source) and flow down the
graph automatically to all other objects which require parameter values to
complete a lensing simulation.

Functional interface: The functional interface eschews the object oriented
caustics code, instead giving the user access to individual mathematical
operations related to lensing, most of which are drawn directly from
gravitational lensing literature. All such functions include references in
their documentation to the relevant papers and equation numbers from which
they are derived. These equations have been tested and implemented in a
reasonably efficient manner. This interface is ideal for researchers and
developers looking to experiment with novel lensing techniques or to modify
existing algorithms while leveraging robust, pre-tested components.

Each layer is in fact built on the one below it, making the transition from
one to the other a matter of following documentation and code references. This
makes the transition easy since one may very clearly observe how their current
analysis can be reproduced in the lower level. Caustics thus provides a
straightforward pipeline for users to move from beginner to expert. Users at
all levels are encouraged to investigate the documentation as the code
includes extensive docstrings for all functions, including units for most
functions. This transparency not only aids in understanding and utilizing the
functions correctly but also enhances the reliability and educational value of
the software.

## Flexibility

The flexibility of caustics is fundamentally linked to its design philosophy,
which is focused on providing a robust yet adaptable framework for
gravitational lensing simulations. A focus on minimalism in the core
functionality means that research-ready analysis routines must be built by the
user. To facilitate this, our Jupyter notebook tutorials include examples of
many typical analysis tasks, with the details laid out for the user so they
may simply copy and modify to suit their particular analysis task. Thus, we
achieve flexibility both by allowing many analysis paradigms, and by
supporting the easy development of production code.

Research is an inherently dynamic process and gravitational lensing is an
evolving field. Designing flexible codes for such environments ensures long
lasting relevance. Caustics is well poised to grow and evolve with the needs
of the community.

## Machine Learning

One of the core purposes of caustics is to advance the application of machine
learning to strong gravitational lensing analysis. This is accomplished
through two avenues. First, as demonstrated in Figure 2, caustics efficiently
generates large samples of simulated mock lensing images by leveraging GPUs.
Since many machine learning algorithms are “data hungry”, this translates to
better performance with more examples to learn from. Literature on machine
learning applications in strong gravitational lensing underscores the benefits
of this generation capacity (ref-Brehmer2019Brehmer et al., 2019; ref-
Chianese2020Chianese et al., 2020; ref-Coogan2020Coogan et al., 2020; ref-
Karchev2022Karchev, Coogan, et al., 2022; ref-Karchev2022bKarchev, Anau
Montel, et al., 2022; ref-Mishra2022Mishra-Sharma & Yang, 2022). Second, the
differentiable nature of caustics allows it to be integrated directly into
machine learning workflows. This could mean using caustics as part of a loss
function. Alternatively, it could be through a statistical paradigm like
diffusion modelling, in which caustics would be directly integrated in the
sampling procedure. It has already been shown that differentiable lensing
simulators, coupled with machine learning and diffusion modelling, can
massively improve source reconstruction in strong gravitational lenses (ref-
Adam2022Adam et al., 2022) and in weak lensing (ref-Remy2023Remy et al.,
2023).

## Conclusions

Here we have presented caustics, a gravitational lensing simulator framework
which allows for greater than 100 times speedup over traditional CPU
implementations by leveraging GPU resources. Caustics is fully-featured,
meaning one can straightforwardly model any strong lensing system with state-
of-the-art techniques. The code and documentation facilitate users transition
from beginner to expert by providing three interfaces which allow increasingly
more flexibility in how one wishes to model a lensing system. Caustics is
designed to be the gravitational lensing simulator of the future and to meet
the hundreds of thousands of lenses soon to be discovered with modern
computational resources.

## Acknowledgements

This research was enabled by a generous donation by Eric and Wendy Schmidt
with the recommendation of the Schmidt Futures Foundation. We acknowledge the
generous software support of the University of Washington Scientific Software
Engineering Center (SSEC) at the eScience Institute, via matching through the
Schmidt Futures Virtual Institute for Scientific Software (VISS). CS
acknowledges the support of a NSERC Postdoctoral Fellowship and a CITA
National Fellowship. This research was enabled in part by support provided by
Calcul Québec and the Digital Research Alliance of Canada. The work of A.A.
and R.L. were partially funded by NSERC CGS D scholarships. Y.H. and L.P.
acknowledge support from the National Sciences and Engineering Council of
Canada grants RGPIN-2020-05073 and 05102, the Fonds de recherche du Québec
grants 2022-NC-301305 and 300397, and the Canada Research Chairs Program.
Thanks to Simon Birrer for communications regarding benchmarking of
lenstronomy.

## References

## References

  * Adam, A., Coogan, A., Malkin, N., Legin, R., Perreault-Levasseur, L., Hezaveh, Y., & Bengio, Y. (2022). Posterior samples of source galaxies in strong gravitational lenses with score-based priors. _Machine Learning and the Physical Sciences Workshop_ , E1. https://doi.org/10.48550/arXiv.2211.03812
  * Adam, A., Perreault-Levasseur, L., Hezaveh, Y., & Welling, M. (2023). Pixelated Reconstruction of Foreground Density and Background Surface Brightness in Gravitational Lensing Systems Using Recurrent Inference Machines. _951_(1), 6. https://doi.org/10.3847/1538-4357/accf84
  * Betancourt, M. (2017). A Conceptual Introduction to Hamiltonian Monte Carlo. _arXiv e-Prints_ , arXiv:1701.02434. https://doi.org/10.48550/arXiv.1701.02434
  * Bingham, E., Chen, J. P., Jankowiak, M., Obermeyer, F., Pradhan, N., Karaletsos, T., Singh, R., Szerlip, P., Horsfall, P., & Goodman, N. D. (2019). Pyro: Deep universal probabilistic programming. _The Journal of Machine Learning Research_ , _20_(1), 973–978.
  * Birrer, S., Shajib, A. J., Gilman, D., Galan, A., Aalbers, J., Million, M., Morgan, R., Pagano, G., Park, J. W., Teodori, L., Tessore, N., Ueland, M., Vyvere, L. V. de, Wagner-Carena, S., Wempe, E., Yang, L., Ding, X., Schmidt, T., Sluse, D., … Amara, A. (2021). Lenstronomy II: A gravitational lensing software ecosystem. _Journal of Open Source Software_ , _6_(62), 3283. https://doi.org/10.21105/joss.03283
  * Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., & Zhang, Q. (2018). _JAX: Composable transformations of Python+NumPy programs_ (Version 0.3.13). http://github.com/google/jax
  * Brehmer, J., Mishra-Sharma, S., Hermans, J., Louppe, G., & Cranmer, K. (2019). Mining for Dark Matter Substructure: Inferring Subhalo Population Properties from Strong Lenses with Machine Learning. _886_(1), 49. https://doi.org/10.3847/1538-4357/ab4c41
  * Chianese, M., Coogan, A., Hofma, P., Otten, S., & Weniger, C. (2020). Differentiable strong lensing: uniting gravity and neural nets through differentiable probabilistic programming. _496_(1), 381–393. https://doi.org/10.1093/mnras/staa1477
  * Collett, T. E. (2015). THE POPULATION OF GALAXY–GALAXY STRONG LENSES IN FORTHCOMING OPTICAL IMAGING SURVEYS. _The Astrophysical Journal_ , _811_(1), 20. https://doi.org/10.1088/0004-637X/811/1/20
  * Coogan, A., Edwards, T. D. P., Chia, H. S., George, R. N., Freese, K., Messick, C., Setzer, C. N., Weniger, C., & Zimmerman, A. (2022). Efficient gravitational wave template bank generation with differentiable waveforms. _106_(12), 122001. https://doi.org/10.1103/PhysRevD.106.122001
  * Coogan, A., Karchev, K., & Weniger, C. (2020). Targeted Likelihood-Free Inference of Dark Matter Substructure in Strongly-Lensed Galaxies. _arXiv e-Prints_ , arXiv:2010.07032. https://doi.org/10.48550/arXiv.2010.07032
  * Desdoigts, L., Pope, B. J. S., Dennis, J., & Tuthill, P. G. (2023). Differentiable optics with $\partial$Lux: I—deep calibration of flat field and phase retrieval with automatic differentiation. _Journal of Astronomical Telescopes, Instruments, and Systems_ , _9_(2), 028007. https://doi.org/10.1117/1.JATIS.9.2.028007
  * Edwards, T. D. P., Wong, K. W. K., Lam, K. K. H., Coogan, A., Foreman-Mackey, D., Isi, M., & Zimmerman, A. (2023). ripple: Differentiable and Hardware-Accelerated Waveforms for Gravitational Wave Data Analysis. _arXiv e-Prints_ , arXiv:2302.05329. https://doi.org/10.48550/arXiv.2302.05329
  * Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. (2013). emcee: The MCMC Hammer. _125_(925), 306. https://doi.org/10.1086/670067
  * Galan, A., Peel, A., Joseph, R., Courbin, F., & Starck, J.-L. (2021). SLITRONOMY: Towards a fully wavelet-based strong lensing inversion technique. _647_ , A176. https://doi.org/10.1051/0004-6361/202039363
  * Gu, A., Huang, X., Sheu, W., Aldering, G., Bolton, A. S., Boone, K., Dey, A., Filipp, A., Jullo, E., Perlmutter, S., Rubin, D., Schlafly, E. F., Schlegel, D. J., Shu, Y., & Suyu, S. H. (2022). GIGA-Lens: Fast Bayesian Inference for Strong Gravitational Lens Modeling. _935_(1), 49. https://doi.org/10.3847/1538-4357/ac6de4
  * Hezaveh, Y., Dalal, N., Holder, G., Kisner, T., Kuhlen, M., & Perreault Levasseur, L. (2016). Measuring the power spectrum of dark matter substructure using strong gravitational lensing. _2016_(11), 048. https://doi.org/10.1088/1475-7516/2016/11/048
  * Hezaveh, Y., Perreault Levasseur, L., & Marshall, P. J. (2017). Fast automated analysis of strong gravitational lenses with convolutional neural networks. _548_(7669), 555–557. https://doi.org/10.1038/nature23463
  * Huang, X., Storfer, C., Gu, A., Ravi, V., Pilon, A., Sheu, W., Venguswamy, R., Banka, S., Dey, A., Landriau, M., Lang, D., Meisner, A., Moustakas, J., Myers, A. D., Sajith, R., Schlafly, E. F., & Schlegel, D. J. (2021). Discovering New Strong Gravitational Lenses in the DESI Legacy Imaging Surveys. _909_(1), 27. https://doi.org/10.3847/1538-4357/abd62b
  * Karchev, K., Anau Montel, N., Coogan, A., & Weniger, C. (2022). Strong-Lensing Source Reconstruction with Denoising Diffusion Restoration Models. _arXiv e-Prints_ , arXiv:2211.04365. https://doi.org/10.48550/arXiv.2211.04365
  * Karchev, K., Coogan, A., & Weniger, C. (2022). Strong-lensing source reconstruction with variationally optimized Gaussian processes. _512_(1), 661–685. https://doi.org/10.1093/mnras/stac311
  * Keeton, C. R. (2011). _GRAVLENS: Computational Methods for Gravitational Lensing_. Astrophysics Source Code Library, record ascl:1102.003.
  * Kneib, J.-P., Bonnet, H., Golse, G., Sand, D., Jullo, E., & Marshall, P. (2011). _LENSTOOL: A Gravitational Lensing Software for Modeling Mass Distribution of Galaxies and Clusters (strong and weak regime)_. Astrophysics Source Code Library, record ascl:1102.004.
  * Lam, S. K., Pitrou, A., & Seibert, S. (2015). Numba: A llvm-based python jit compiler. _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_ , 1–6.
  * Metcalf, R. B., & Petkova, M. (2014). GLAMER - I. A code for gravitational lensing simulations with adaptive mesh refinement. _445_(2), 1942–1953. https://doi.org/10.1093/mnras/stu1859
  * Millon, M., Michalewicz, K., Dux, F., Courbin, F., & Marshall, P. J. (2024). Image deconvolution and PSF reconstruction with STARRED: a wavelet-based two-channel method optimized for light curve extraction. _arXiv e-Prints_ , arXiv:2402.08725. https://doi.org/10.48550/arXiv.2402.08725
  * Mishra-Sharma, S., & Yang, G. (2022). Strong Lensing Source Reconstruction Using Continuous Neural Fields. _Machine Learning for Astrophysics_ , 34. https://doi.org/10.48550/arXiv.2206.14820
  * Nightingale, James. W., Hayes, R. G., Kelly, A., Amvrosiadis, A., Etherington, A., He, Q., Li, N., Cao, X., Frawley, J., Cole, S., Enia, A., Frenk, C. S., Harvey, D. R., Li, R., Massey, R. J., Negrello, M., & Robertson, A. (2021). ‘PyAutoLens‘: Open-source strong gravitational lensing. _Journal of Open Source Software_ , _6_(58), 2825. https://doi.org/10.21105/joss.02825
  * Nikolic, B. (2018). Acceleration of Non-Linear Minimisation with PyTorch. _arXiv e-Prints_ , arXiv:1805.07439. https://doi.org/10.48550/arXiv.1805.07439
  * Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., … Chintala, S. (2019). PyTorch: An imperative style, high-performance deep learning library. In _Advances in neural information processing systems 32_ (pp. 8024–8035). Curran Associates, Inc. http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * Peng, C. Y., Impey, C. D., Rix, H.-W., Kochanek, C. S., Keeton, C. R., Falco, E. E., Lehár, J., & McLeod, B. A. (2006). Probing the Coevolution of Supermassive Black Holes and Galaxies Using Gravitationally Lensed Quasar Hosts. _649_(2), 616–634. https://doi.org/10.1086/506266
  * Perreault Levasseur, L., Hezaveh, Y. D., & Wechsler, R. H. (2017). Uncertainties in Parameters Estimated with Neural Networks: Application to Strong Gravitational Lensing. _850_(1), L7. https://doi.org/10.3847/2041-8213/aa9704
  * Petrillo, C. E., Tortora, C., Chatterjee, S., Vernardos, G., Koopmans, L. V. E., Verdoes Kleijn, G., Napolitano, N. R., Covone, G., Schneider, P., Grado, A., & McFarland, J. (2017). Finding strong gravitational lenses in the Kilo Degree Survey with Convolutional Neural Networks. _472_(1), 1129–1150. https://doi.org/10.1093/mnras/stx2052
  * Remy, B., Lanusse, F., Jeffrey, N., Liu, J., Starck, J.-L., Osato, K., & Schrabback, T. (2023). Probabilistic mass-mapping with neural score estimation. _672_ , A51. https://doi.org/10.1051/0004-6361/202243054
  * Rodney, S. A., Brammer, G. B., Pierel, J. D. R., Richard, J., Toft, S., O’Connor, K. F., Akhshik, M., & Whitaker, K. E. (2021). A gravitationally lensed supernova with an observable two-decade time delay. _Nature Astronomy_ , _5_ , 1118–1125. https://doi.org/10.1038/s41550-021-01450-9
  * Speagle, J. S. (2020). DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences. _493_(3), 3132–3158. https://doi.org/10.1093/mnras/staa278
  * Stone, C. J., Courteau, S., Cuillandre, J.-C., Hezaveh, Y., Perreault-Levasseur, L., & Arora, N. (2023). ASTROPHOT: fitting everything everywhere all at once in astronomical images. _525_(4), 6377–6393. https://doi.org/10.1093/mnras/stad2477
  * Suyu, S. H., & Halkola, A. (2010). The halos of satellite galaxies: the companion of the massive elliptical lens SL2S J08544-0121. _524_ , A94. https://doi.org/10.1051/0004-6361/201015481
  * Vegetti, S., & Vogelsberger, M. (2014). On the density profile of dark matter substructure in gravitational lens galaxies. _442_(4), 3598–3603. https://doi.org/10.1093/mnras/stu1284
  * Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., … SciPy 1.0 Contributors. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. _Nature Methods_ , _17_ , 261–272. https://doi.org/10.1038/s41592-019-0686-2
  * Wagner-Carena, S., Lee, J., Pennington, J., Aalbers, J., Birrer, S., & Wechsler, R. H. (2024). A Strong Gravitational Lens Is Worth a Thousand Dark Matter Halos: Inference on Small-Scale Structure Using Sequential Methods. _arXiv e-Prints_ , arXiv:2404.14487. https://doi.org/10.48550/arXiv.2404.14487
  * Welch, B., Coe, D., Zackrisson, E., de Mink, S. E., Ravindranath, S., Anderson, J., Brammer, G., Bradley, L., Yoon, J., Kelly, P., Diego, J. M., Windhorst, R., Zitrin, A., Dimauro, P., Jiménez-Teja, Y., Abdurro’uf, Nonino, M., Acebron, A., Andrade-Santos, F., … Vikaeus, A. (2022). JWST Imaging of Earendel, the Extremely Magnified Star at Redshift z = 6.2. _940_(1), L1. https://doi.org/10.3847/2041-8213/ac9d39
  * Wong, K. C., Suyu, S. H., Chen, G. C.-F., Rusu, C. E., Million, M., Sluse, D., Bonvin, V., Fassnacht, C. D., Taubenberger, S., Auger, M. W., Birrer, S., Chan, J. H. H., Courbin, F., Hilbert, S., Tihhonova, O., Treu, T., Agnello, A., Ding, X., Jee, I., … Meylan, G. (2020). H0LiCOW - XIII. A 2.4 per cent measurement of H0 from lensed quasars: 5.3$\sigma$ tension between early- and late-Universe probes. _498_(1), 1420–1439. https://doi.org/10.1093/mnras/stz3094
  * Wong, K. W. K., Isi, M., & Edwards, T. D. P. (2023). Fast Gravitational-wave Parameter Estimation without Compromises. _958_(2), 129. https://doi.org/10.3847/1538-4357/acf5cd