#  binary_c-python: A Python-based stellar population synthesis tool and
interface to binary_c

D. D. Hendriks E-mail<EMAIL_ADDRESS>(DDH) Department of Physics,
University of Surrey, Guildford, GU2 7XH, Surrey, UK  0000-0002-8717-6046 R.
G. Izzard  Department of Physics, University of Surrey, Guildford, GU2 7XH,
Surrey, UK  0000-0003-0378-4843

(11 April 2023)

## Summary

We present the software package binary_c-python which provides a convenient
and easy-to-use interface to the binary_c (Izzard et al., 2004, 2006, 2009,
2018; Izzard & Jermyn, 2022) framework, allowing the user to rapidly evolve
individual systems and populations of stars. binary_c-python is available on
Pip and on GitLab.

binary_c-python contains many useful features to control and process the
output of binary_c, like by providing binary_c-python with logging statements
that are dynamically compiled and loaded into binary_c. Moreover, we have
recently added standardised output of events like Roche-lobe overflow or
double compact-object formation to binary_c, and automatic parsing and
managing of that output in binary_c-python. binary_c-python uses
multiprocessing to utilise all the cores on a particular machine, and can run
populations with HPC cluster workload managers like HTCondor and Slurm,
allowing the user to run simulations on large computing clusters.

Recent developments in binary_c include standardised output datasets called
_ensembles_ , which contain nested histograms of binned data like supernovae
rates or chemical yields. binary_c-python easily processes these datasets and
provides a suite of utility functions to handle them. Furthermore, binary_c
now includes the _ensemble-manager_ class, which makes use of the core
functions and classes of binary_c-python to evolve a grid of stellar
populations with varying input physics, allowing for large, automated
parameter studies through a single interface.

binary_c-python is easily integrated with other Python-based codes and
libraries, e.g. sampling codes like Emcee (Foreman-Mackey et al., 2013) or
Dynesty (Speagle, 2020), or the astrophysics oriented package Astropy (Astropy
Collaboration et al., 2013; Astropy Collaboration et al., 2018). Moreover, it
is possible to provide custom system-generating functions through function
hooks, allowing third-party packages to manage the properties of the stars in
the populations and evolve them through binary_c-python.

We provide documentation that is automatically generated based on _docstrings_
and a suite of Jupyter notebooks. These notebooks consist of technical
tutorials on how to use binary_c-python and use-case scenarios aimed at doing
science. Much of binary_c-python is covered by unit tests to ensure
reliability and correctness, and the test coverage is continually increased as
the package is improved.

## Statement of need

In the current scientific climate Python is ubiquitous. While lower level code
written in, e.g., Fortran or C is still widely used, much new software is
written in Python, either entirely or as a wrapper around other code and
libraries. Education in programming also often includes Python courses because
of its ease of use and its flexibility. Moreover, Python has a large community
with many resources and tutorials. We have created binary_c-python to allow
students and scientists alike to explore current scientific issues while
enjoying the familiar syntax, and at the same time make use of the plentiful
scientific and astrophysical packages like Numpy, Scipy, Pandas, Astropy and
platforms like Jupyter.

Over time, many (binary) stellar population synthesis codes have been created.
Each one usually has a slightly different focus, like gravitational waves or
spectral synthesis. Recent ones are BPASS/HOKI (Stevance et al., 2020), COMPAS
(Compas et al., 2022), COSMIC (Breivik et al., 2019), MOBSE (Giacobbo &
Mapelli, 2019) and SeBa (Toonen & Nelemans, 2013). Most of these have Python
interfaces, which further highlights the need to develop and release a modern
Python interface to binary_c, i.e. binary_c-python.

The previous interface to binary_c, binary_grid was written in Perl, where
much of the logic and structure was developed and debugged. While much of the
code base of binary_c-python has changed significantly from its Perl
predecessor, the initial porting to Python and the development of binary_c-
python greatly benefitted from the existence of this earlier interface.

## Projects that use binary_c-python

binary_c-python has already been used in a variety of situations, ranging from
pure research to educational purposes, as well as in outreach events. In the
summer of 2021 we used binary_c-python as the basis for the interactive
classes on stellar ecosystems during the International Max-Planck Research
School summer school 2021 in Heidelberg. Students were introduced to the topic
of population synthesis and were able to use our notebooks to perform their
own calculations. binary_c-python has been used in Mirouh et al. (submitted),
implementing improvements to tidal interactions between stars and varying
initial birth parameter distributions to match to observed binary systems in
star clusters. The binary_c-python and Emcee packages were used in a Master’s
thesis project to find the birth system parameters of the V106 stellar system,
compare observations to results of binary_c, and perform uncertainty with
Bayesian uncertainty inference through Markov chain Monte Carlo sampling.

Currently, binary_c-python is used in several ongoing projects. These include
a study on the effect of birth distributions on the occurrence of carbon-
enhanced metal-poor (CEMP) stars, the occurence and properties of accretion
disks in main-sequence stars (Hendriks & Izzard, in prep), and the predicted
observable black hole distribution by combining star formation and metallicity
distributions with the output of binary_c (Hendriks et al., in prep). We also
use the _ensemble_ output structure to generate datasets for galactochemical
evolution over cosmological timescales, where we rely heavily on the utilities
of binary_c-python (Yates et al., in prep).

## Acknowledgements

We are grateful for the helpful discussions and testing efforts from A.
Aryeaipour, M. Delorme, S. Dykes, G. Mirouh, M. Matteuzzi, N. Rees, M. Renzo,
L. van Son, K. Temmink, D. Tracey and R. Yates, and the early work of J.
Andrews which inspired our Python-C interface code. DDH thanks the UKRI/UoS
for the funding grant H120341A. RGI thanks STFC for funding grants
ST/R000603/1 and ST/L003910/2.

## References

reAstropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M.,
Lim, P. L., Crawford, S. M., Conseil, S., Shupe, D. L., Craig, M. W.,
Dencheva, N., Ginsburg, A., VanderPlas, J. T., Bradley, L. D., Pérez-Suárez,
D., de Val-Borro, M., Aldcroft, T. L., Cruz, K. L., Robitaille, T. P.,
Tollerud, E. J., … Astropy Contributors. (2018). The Astropy Project: Building
an Open-science Project and Status of the v2.0 Core Package. _The Astronomical
Journal_ , _156_ , 123. https://doi.org/10.3847/1538-3881/aabc4f

preAstropy Collaboration, Robitaille, T. P., Tollerud, E. J., Greenfield, P.,
Droettboom, M., Bray, E., Aldcroft, T., Davis, M., Ginsburg, A., Price-Whelan,
A. M., Kerzendorf, W. E., Conley, A., Crighton, N., Barbary, K., Muna, D.,
Ferguson, H., Grollier, F., Parikh, M. M., Nair, P. H., … Streicher, O.
(2013). Astropy: A community Python package for astronomy. _Astronomy and
Astrophysics_ , _558_ , A33. https://doi.org/10.1051/0004-6361/201322068

preBreivik, K., Coughlin, S., Zevin, M., Rodriguez, C., Andrews, J., Kimball,
C., mcdigman, & 1nhtran. (2019). _COSMIC-PopSynth/COSMIC: COSMIC now
integrated with Globular Cluster Simulation Software ClusterMonteCarlo (CMC)_
(Version v3.1.0). Zenodo. https://doi.org/10.5281/zenodo.3482915

preCompas, T., Riley, J., Agrawal, P., Barrett, J. W., Boyett, K. N. k.,
Broekgaarden, F. S., Chattopadhyay, D., Gaebel, S. M., Gittins, F., Hirai, R.,
Howitt, G., Justham, S., Khandelwal, L., Kummer, F., Lau, M. Y. m., Mandel,
I., Mink, S. E. de, Neijssel, C., Riley, T., … Willcox, R. (2022). COMPAS: A
rapid binary population synthesis suite. _Journal of Open Source Software_ ,
_7_(69), 3838. https://doi.org/10.21105/joss.03838

preForeman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. (2013). Emcee: The
MCMC Hammer. _Publications of the Astronomical Society of the Pacific_ ,
_125_(925), 306–312. https://doi.org/10.1086/670067

preGiacobbo, N., & Mapelli, M. (2019). The impact of electron-capture
supernovae on merging double neutron stars. _482_(2), 2234–2243.
https://doi.org/10.1093/mnras/sty2848

preHendriks, D. D., & Izzard, R. G. (in prep). Disk mass loss during mass
transfer onto main sequence accretors.

preHendriks, D. D., van Son, L. A. C., Renzo, M., & Izzard, R. G. (in prep).
Pulsational pair-instability supernovae in gravitational-wave and
electromagnetic transients.

preIzzard, R. G., Dray, L. M., Karakas, A. I., Lugaro, M., & Tout, C. A.
(2006). Population nucleosynthesis in single and binary stars. I. Model.
_460_(2), 565–572. https://doi.org/10.1051/0004-6361:20066129

preIzzard, R. G., Glebbeek, E., Stancliffe, R. J., & Pols, O. R. (2009).
Population synthesis of binary carbon-enhanced metal-poor stars. _508_(3),
1359–1374. https://doi.org/10.1051/0004-6361/200912827

preIzzard, R. G., & Jermyn, A. S. (2022). Circumbinary discs for stellar
population models. _Monthly Notices of the Royal Astronomical Society_.
https://doi.org/10.1093/mnras/stac2899

preIzzard, R. G., Preece, H., Jofre, P., Halabi, G. M., Masseron, T., & Tout,
C. A. (2018). Binary stars in the Galactic thick disc. _473_(3), 2984–2999.
https://doi.org/10.1093/mnras/stx2355

preIzzard, R. G., Tout, C. A., Karakas, A. I., & Pols, O. R. (2004). A new
synthetic model for asymptotic giant branch stars. _350_(2), 407–426.
https://doi.org/10.1111/j.1365-2966.2004.07446.x

preMirouh, G. M., Hendriks, D. D., Dykes, S., Moe, M., & Izzard, R. G.
(submitted). Detailed equilibrium and dynamical tides: impact on
circularization and synchronization in open clusters.

preSpeagle, J. S. (2020). Dynesty: A dynamic nested sampling package for
estimating Bayesian posteriors and evidences. _Monthly Notices of the Royal
Astronomical Society_ , _493_(3), 3132–3158.
https://doi.org/10.1093/mnras/staa278

preStevance, H., Eldridge, J., & Stanway, E. (2020). Hoki: Making BPASS
accessible through Python. _The Journal of Open Source Software_ , _5_(45),
1987. https://doi.org/10.21105/joss.01987

preToonen, S., & Nelemans, G. (2013). The effect of common-envelope evolution
on the visible population of post-common-envelope binaries. _557_ , A87.
https://doi.org/10.1051/0004-6361/201321753

preYates, R. M., Hendriks, D. D., Vijayan, A. P., Izzard, R. G., Thomas, P.
A., & Das, P. (in prep). Modelling binary stellar populations & dust in
cosmological-scale galaxy evolution simulations.

p


# What can a detected photon with a given gravitational redshift tell us about
the maximum density of a compact star?

Shahar Hod  The Ruppin Academic Center, Emeq Hefer 40250, Israel  The Hadassah
Institute, Jerusalem 91010, Israel

###### Abstract

Far away observers can in principle bound from below the dimensionless
maximum-density parameter $\Lambda\equiv 4\pi R^{2}\rho_{\text{max}}$ of a
compact star by measuring the gravitational redshift factor
$z\equiv\nu_{\text{e}}/\nu_{\infty}-1$ of photons that were emitted from the
surface of the star: $\Lambda\geq{3\over 2}[1-(1+z)^{-2}]$ [here $R$ is the
radius of the star and $\\{\nu_{\text{e}},\nu_{\infty}\\}$ are respectively
the frequency of the emitted light as measured at the location of the emission
and by asymptotic observers]. However, if photons that were created somewhere
inside the star can make their way out and reach the asymptotic observers,
then the measured redshift parameter $z$ may not determine uniquely the
surface properties of the star, thus making the above bound unreliable. In the
present compact paper we prove that in these cases, in which the creation
depth of a detected photon is not known to the far away observers, the
empirically measured redshift parameter can still be used to set a (weaker)
lower bound on the dimensionless density parameter of the observed star:
$\Lambda\geq{3\over 2}[1-(1+z)^{-2/3}]$.

##  I Introduction

The gravitational redshift effect in general relativity Chan ; ShTe  implies
that electromagnetic waves traveling out of a gravitational well (which, for
example, may be created by the presence of a compact star) seem to lose
energy. This physically important phenomenon is reflected in the fact that the
frequencies of the detected photons as measured by flat-space asymptotic
observers are lower than the original frequencies of the corresponding photons
as determined at the source of the emission (the compact star).

In spherically symmetric spacetimes the ratio between these two frequencies
can be expressed in terms of the $tt$-component of the curved line element at
the emission site Chan ; ShTe :

$1+z\equiv{{\nu_{\text{e}}}\over{\nu_{\infty}}}={{1}\over{\sqrt{g_{00}}}}\ ,$
(1)

where $\\{\nu_{\text{e}},\nu_{\infty}\\}$ are respectively the photon
frequency as measured at the source of emission and the detected frequency as
measured by far away asymptotic observers. In particular, since the spacetime
region just outside the surface of a compact star of mass $M$ and radius $R$
is characterized by the simple relation $g_{00}=1-2M/R$ Noteunits , the
frequencies of photons that were created near the surface of the star are
gravitationally redshifted according to the simple dimensionless ratio Chan ;
ShTe

$1+z\equiv{{\nu_{\text{e}}}\over{\nu_{\infty}}}=\Big{(}1-{{2M}\over{R}}\Big{)}^{-{1\over
2}}\ .$ (2)

An immediate consequence of the simple relation (2) is that the dimensionless
maximum-density-area parameter

$\Lambda\equiv 4\pi R^{2}\rho_{\text{max}}\ $ (3)

of the emitting star can in principle be bounded from below by far away
observers who measure the gravitational redshift parameter $z$ of an emitted
photon that was created near the surface of the star:

$\Lambda\geq{3\over 2}[1-(1+z)^{-2}]\ .$ (4)

(We have used here the simple inequality $M\leq{4\over 3}\pi
R^{3}\rho_{\text{max}}$ for the mass of the observed star NoteK1 ).

However, caution should be taken in interpreting observational data of stellar
emission spectra: one has to take into account the possibility that photons
can in principle also be created in the non-vacuum region inside the star,
where the radial functional behavior of the metric component $g_{00}(r)$ [see
Eq. (1)] may not be known to the far away observers.

In particular, if a photon that was created in some unknown depth inside the
non-vacuum region of the star can make its way out, then the empirically
measured redshift parameter of the detected photon cannot be used by the
asymptotic observers to determine uniquely the dimensionless compactness $M/R$
[see Eq. (2)] of the emitting star, thus making the simple lower bound (4)
unreliable.

The main goal of the present compact paper is to present a compact theorem
that reveals the interesting fact that in these situations, in which the
detected photon was originally created in the non-vacuum region inside the
star (in some depth which may not be known to the far away asymptotic
observers), the observers can still use the empirically measured redshift
factor of the detected photon in order to derive a lower bound, which is
albeit weaker than (4), on the dimensionless maximum density parameter
$\Lambda$ of the observed star.

##  II Description of the system

We consider a compact star of mass $M$ and radius $R$ whose spherically
symmetric asymptotically flat spacetime is described, using the Schwarzschild
spacetime coordinates, by the curved line element Chan ; ShTe

$ds^{2}=-e^{-2\delta}\mu
dt^{2}+\mu^{-1}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta d\phi^{2})\ $ (5)

with the radially-dependent metric functions $\mu=\mu(r)$ and
$\delta=\delta(r)$.

The Einstein-matter field equations, $G^{\mu}_{\nu}=8\pi T^{\mu}_{\nu}$, yield
the non-linear differential relations May ; Hodt1

${{d\mu}\over{dr}}=-8\pi r\rho+{{1-\mu}\over{r}}\ $ (6)

and

${{d\delta}\over{dr}}=-{{4\pi r(\rho+p)}\over{\mu}}\ $ (7)

for the metric functions. Here

$\rho\equiv-T^{t}_{t}\ \ \ \ \text{and}\ \ \ \ p\equiv T^{r}_{r}$ (8)

are respectively the energy density and radial pressure of the matter fields,
which are assumed to be non-negative in the interior region of the star and
respect the dominant energy condition HawEl ; Bond1 :

$0\leq p\leq\rho\ \ \ \ \ \text{for}\ \ \ \ \ r\leq R\ .$ (9)

In addition, the energy density and pressure are assumed to vanish outside the
compact surface of the star:

$\rho=p=0\ \ \ \ \ \text{for}\ \ \ \ \ r>R\ .$ (10)

We assume that the spacetime of the compact star is spatially regular, which
implies the near-origin relations May ; Hodt1

$\mu(r\to 0)=1+O(r^{2})\ \ \ \ {\text{and}}\ \ \ \ \delta(0)<\infty\ .$ (11)

In addition, the metric functions of the asymptotically flat spacetime of the
star are characterized by the functional behaviors May ; Hodt1

$\mu(r\to\infty)\to 1\ \ \ \ {\text{and}}\ \ \ \ \delta(r\to\infty)\to 0\ $
(12)

at spatial infinity.

Taking cognizance of the Einstein equation (6), one finds the simple
expression

$\mu(r)=1-{{2m(r)}\over{r}}\ $ (13)

for the metric function, where the gravitational mass contained within a
sphere of radius $r$ is given by the integral relation May ; Hodt1

$m(r)=4\pi\int_{0}^{r}x^{2}\rho(x)dx\ $ (14)

and is characterized by the boundary condition [see Eq. (10)]

$m(r=R)=M\ .$ (15)

##  III Generic lower bound on the dimensionless maximum density parameter of
an optically observed star

In the present section we shall derive, using analytical techniques, a generic
lower bound on the dimensionless maximum-density-area parameter $\Lambda$ of
optically observed stars. In particular, our goal is to derive a model-
independent bound which would be valid even in situations in which the
creation depth inside the star of the asymptotically detected photon is not
known to the far away observers.

We first point out that Eqs. (1), (5), and (13) yield the functional relation

$1+z(r_{\text{cr}})={{e^{\delta(r_{\text{cr}})}}\over{\sqrt{1-{{2m(r_{\text{cr}})}\over{r_{\text{cr}}}}}}}\
$ (16)

for the asymptotically measured gravitational redshift factor,
$z=z(r_{\text{cr}})$, that characterizes a photon that originally was created
inside the star at a radial distance $r=r_{\text{cr}}$ from its center and
eventually was detected by the far away asymptotic observers NoteK2 .

In order to emphasize the fact that far away observers can in principle
determine the redshift factor $z$ of an asymptotically detected photon but, in
general, they may not have precise knowledge about the exact depth (or,
equivalently, the exact radius $r=r_{\text{cr}}$) at which the photon was
originally created inside the star, we shall henceforth use the notation $z$
[instead of $z(r_{\text{cr}})$] for the empirically determined redshift factor
of an asymptotically detected photon.

Our main goal is to derive a lower bound on the dimensionless density
parameter $\Lambda$ of the observed star. To this end, we shall first prove
that the metric function $e^{\delta(r)}/\sqrt{\mu(r)}$, which determines the
redshift parameter (16), is monotonically decreasing inside the star. Using
the Einstein equations (6) and (7) with the radial relation (13), one obtains
the gradient relation [see Eq. (9)]

${{d\Big{[}{{e^{\delta(r)}}\over{\sqrt{\mu(r)}}}\Big{]}}\over{dr}}=-{{e^{\delta(r)}}\over{\sqrt{\mu(r)}}}\cdot{{{{m(r)}\over{r}}+4\pi
r^{2}p}\over{\mu r}}<0\ ,$ (17)

which implies [see Eq. (11)]

${\text{max}}_{r}\Big{\\{}{{e^{\delta(r)}}\over{\sqrt{\mu(r)}}}\Big{\\}}=e^{\delta(0)}\
\ \ \ \ \text{for}\ \ \ \ \ r\in[0,R]\ .$ (18)

We shall next derive a useful generic upper bound on the radially-dependent
dimensionless metric function $\delta(r)$. Using the simple relation [see Eq.
(14)]

$m(r)\leq{{4\pi}\over{3}}r^{3}\cdot\rho_{\text{max}}\ $ (19)

for the gravitational mass contained within a sphere of radius $r$, one finds
from Eq. (13) the characteristic inequality

$\mu(r)\geq 1-{{8\pi}\over{3}}r^{2}\rho_{\text{max}}\ .$ (20)

We shall henceforth assume the relation

$\Lambda<{3\over 2}$ (21)

for the dimensionless density parameter of the star [otherwise we immediately
obtain the lower bound $\Lambda\geq{3\over 2}$, which is actually stronger
than the bound (29) that we shall prove in the present section], which implies
the characteristic inequalities [see Eq. (20)]

$0<\mu(r)\leq 1\ .$ (22)

Taking cognizance of Eqs. (7), (9), and (20), one can write the following
series of inequalities

$-{{d\delta}\over{dr}}\leq{{8\pi r\rho(r)}\over{\mu(r)}}\leq{{8\pi
r\rho_{\text{max}}}\over{1-{{8\pi}\over{3}}r^{2}\rho_{\text{max}}}}\ $ (23)

for the gradient of the metric function in the interior region ($r\leq R$) of
the compact star. From (23) one obtains the integral relation

$-\int^{\delta(R)}_{\delta(r)}d\delta\leq\int^{R}_{r}{{8\pi
x\rho_{\text{max}}}\over{1-{{8\pi}\over{3}}x^{2}\rho_{\text{max}}}}dx\ .$ (24)

Performing the integration in (24) with the help of the boundary relation [see
Eqs. (7), (10), and (12)]

$\delta(r=R)=0\ ,$ (25)

one finds the characteristic inequality

$\delta(r)\leq{3\over 2}\cdot\ln\Big{(}{{3-8\pi
r^{2}\rho_{\text{max}}}\over{3-8\pi R^{2}\rho_{\text{max}}}}\Big{)}\ $ (26)

for the radially-dependent metric function. In particular, one can write the
upper bound

$\delta(0)\leq{3\over 2}\cdot\ln\Big{(}{{3}\over{3-8\pi
R^{2}\rho_{\text{max}}}}\Big{)}\ $ (27)

on the value of the metric function at the center of the star.

Taking cognizance of Eqs. (3), (16), (18), and (27), one obtains the
dimensionless inequalities

$1+z\leq e^{\delta(0)}\leq\Big{(}{{3}\over{3-2\Lambda}}\Big{)}^{3\over 2}\ ,$
(28)

which imply the lower bound

$\Lambda\geq{3\over 2}[1-(1+z)^{-2/3}]\ $ (29)

on the characteristic maximum-density-area parameter of the optically observed
compact star.

We note that the bound (29) is weaker than the bound (4). It is important to
emphasize, however, that while the familiar bound (4) is only valid if one
assumes that the asymptotically detected photon was created near the surface
of the star, the analytically derived lower bound (29) is valid even in
situations in which the creation depth of the asymptotically detected photon
inside the non-vacuum region of the star is not known to us.

##  IV Summary and physical implications

Determining the value of the composed maximum-density-area parameter
$\Lambda\equiv 4\pi R^{2}\rho_{\text{max}}$ of highly compact stars is a
challenging task that may help physicists to determine the correct equation of
state of highly dense nuclear matter configurations. Interestingly, using the
bound (4), the value of this dimensionless physical parameter can in principle
be bounded from below by far away observers who measure the gravitational
redshift parameter $z$ of an asymptotically detected photon that was created
near the surface of the compact star.

In the present paper we have emphasized the fact that, in principle, photons
can also be created in the non-vacuum region inside the star, in which case
the asymptotically measured redshift parameter $z$ may not determine uniquely
the surface properties of the star, thus making the simple bound (4)
unreliable.

Motivated by this observation, we have presented a compact theorem that
reveals the physically important fact that in these situations, in which the
detected photon may have been created in the non-vacuum region inside the
compact star (in some depth which is not necessarily known to the far away
observers), the asymptotic observers can still use the empirically measured
redshift factor of the detected photon in order to set an upper bound on the
value of the dimensionless density parameter $\Lambda$ that characterizes the
observed star.

In particular, using the non-linearly coupled Einstein-matter field equations,
we have derived the relation [see Eqs. (3) and (29)]

$\Lambda\equiv 4\pi R^{2}\rho_{\text{max}}\geq{3\over 2}[1-(1+z)^{-2/3}]\ $
(30)

between the dimensionless maximum-density-area parameter of an emitting star
and the gravitational redshift factor that characterizes the detected photon
as measured by far away asymptotic observers.

Inspection of the analytically derived lower bound (30) reveals the fact that
the strongest bound on the dimensionless density parameter of the star may be
deduced by substituting in (30) the value $z_{\text{max}}$ that characterizes
the most redshifted photon which is detected in the emission spectrum of the
star by the asymptotic observers.

In order to illustrate the important physical implications of our results, let
us assume that asymptotic observers detect an emitted photon with say $z=1/2$.
According to the familiar vacuum relation (4), this gravitational redshift
factor yields the lower bound $\Lambda\geq 5/6\simeq 0.833$ on the
dimensionless density parameter of the emitting star. However, the use of the
vacuum relation (4) may not be justified if the detected photon was created in
the non-vacuum region inside the star. In these cases one should instead rely
on the analytically derived density-redshift relation (30), which implies that
the asymptotically detected photon with the property $z=1/2$ corresponds to
the corrected bound $\Lambda\gtrsim 0.355$ on the dimensionless maximum-
density-area parameter of the observed star.

Finally, we would like to emphasize again that the analytically derived lower
bound (30) on the dimensionless density parameter (3) of an optically observed
star is valid even in situations in which the creation depth of the
asymptotically detected photon inside the compact star is not known to us.

ACKNOWLEDGMENTS

This research is supported by the Carmel Science Foundation. I would like to
thank Yael Oren, Arbel M. Ongo, Ayelet B. Lata, and Alona B. Tea for
stimulating discussions.

## References

  * (1) S. Chandrasekhar, The Mathematical Theory of Black Holes, (Oxford University Press, New York, 1983).
  * (2) S. L. Shapiro and S. A. Teukolsky, Black Holes, White Dwarfs and Neutron Stars: The Physics of Compact Objects, 1st ed. (Wiley-Interscience, 1983).
  * (3) We shall use gravitational units in which $G=c=1$.
  * (4) Note that the dimensionless density-area parameter $\Lambda=4\pi R^{2}\rho_{\text{max}}$ of the compact star is related to its dimensionless compactness parameter $C=M/R$ by the simple inequality $C\leq\Lambda/3$.
  * (5) A. E. Mayo and J. D. Bekenstein, Phys. Rev. D 54, 5059 (1996).
  * (6) S. Hod, Phys. Rev. D 84, 124030 (2011) [arXiv:1112.3286]; S. Hod, Phys. Rev. D 84, 104024 (2011) [arXiv:1201.0068].
  * (7) H. Bondi, Mon. Not. Roy. Astr. Soc. 259, 365 (1992).
  * (8) S. W. Hawking and G. F. Ellis, The Large Scale Structure of Spacetime (Cambridge University Press, Cambridge 1973).
  * (9) It is interesting to note that Eqs. (3), (14), and (16) yield the series of inequalities $\Lambda/3={4\over 3}\pi R^{2}\rho_{\text{max}}\geq{4\over 3}\pi r^{2}_{\text{cr}}\cdot\rho_{\text{max}}\geq m(r_{\text{cr}})/r_{\text{cr}}={1\over 2}\big{\\{}1-[1+z(r_{\text{cr}})]^{-2}e^{2\delta(r_{\text{cr}})}\big{\\}}$ from which one deduces the lower bound $\Lambda\geq{3\over 2}\big{\\{}1-[1+z(r_{\text{cr}})]^{-2}e^{2\delta(r_{\text{cr}})}\big{\\}}$, where $r_{\text{cr}}$ is the radius at which the photon was originally created inside the star. As a consistency check we point out that, as should be the case, the relation (4) is recovered from this inequality in the $r_{\text{cr}}\to R$ limit where $\delta(R)=0$ [see Eqs. (7), (10), and (12) which yield the relation $\delta(r=R)=0$ on the surface of the compact star]. It is worth emphasizing that far away observers can in principle measure the redshift factor of a detected photon but, in general, they may not have precise knowledge about the exact radius $r=r_{\text{cr}}$ (the exact depth) at which the photon was originally created inside the star. Thus, far away observers may not know the exact value of $\delta(r_{\text{cr}})$ at the creation point of the photon inside the star [note that Eqs. (7), (9), (10), and (12) yield the generic inequality $\delta(r_{\text{cr}})\geq 0$]. This fact motivates us to derive in the present paper a robust lower bound on the density-area parameter $\Lambda$ of the observed star which does not depend on the possibly unknown value of the metric function $\delta(r_{\text{cr}})$ at the creation point of the photon inside the star.


# Optimization under rare events: scaling laws for
linear chance-constrained programs

Jose Blanchet1 1Stanford University, Department of Management Science and
Engineering  ,  Joost Jorritsma2 2Centrum voor Wiskunde & Informatica,
Amsterdam  and  Bert Zwart2,3 3Eindhoven University of Technology, Department
of Mathematics and Computer Science<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We consider a class of chance-constrained programs in which profit needs to be
maximized while enforcing that a given adverse event remains rare. Using
techniques from large deviations and extreme-value theory, we show how the
optimal value scales as the prescribed bound on the violation probability
becomes small and how convex programs emerge in the limit. We also use our
results to show that the popular CVaR approximation is asymptotically optimal
under light-tail assumptions, but it is sub-optimal in a heavy-tailed setting.
Then, we apply point process techniques and random set theory to study the
suboptimality gap in the Monte Carlo based scenario approach.

Keywords: chance-constrained optimization, rare events, extreme values,
regular variation.

MSC Class: 90C15 (primary), 60F10, 60G70 (secondary).

##  1\. Introduction

We analyze a class of optimization problems where the goal is to maximize
profit, subject to a reliability constraint that requires the probability of
an undesirable event below a small level $\delta$. These formulations, known
as chance-constrained (CC) optimization problems, arise in a wide range of
applications, including:

  * –

Financial planning. Banks and insurance companies need to optimize their
resources such that the probability of insolvency is kept small. For example,
in fixed income portfolio optimization, an investment grade portfolio has a
default rate of about $\delta=10^{-4}$, see Frank (2008).

  * –

Network design. In ultra-reliable communication system design for 5G wireless
networks, the failure probability should be no bigger than $\delta=10^{-5}$,
as is explained in Alsenwi _et al._ (2019).

  * –

Power grid operations. The European Union has set a goal that the amount of
reserves available to mitigate possible frequency disturbance should be
sufficient, up to once every 20 years. As these reserves are typically
procured on an hourly basis, this amounts to a safety constraint of $\delta$
smaller than $10^{-5}$, see Mosca (2024).

Additional examples originate in the design of bridges and control of
aircraft, cf. Shi and Holzapfel (2021), where the probability of a rare event
is sometimes required to be below $\delta=10^{-9}$. For more examples, see
Subramanyam (2020), Tong _et al._ (2022b).

Because of their direct interpretation and versatility to be applied to
various important areas, chance-constrained optimization has been a subject of
interest in Operations Research for several decades. It has become one of the
standard formulations in risk-constrained optimization and is covered as such
in several textbooks cf. Prekopa (1970), Prékopa (2003), Klein Haneveld _et
al._ (2020), Shapiro _et al._ (2021). Unfortunately, chance-constrained
problems are typically not very tractable: the set of feasible solutions may
be non-convex, and evaluating the chance constraint can be a computational
challenge. In general, CC problems are provably NP-hard, see Luedtke _et al._
(2010). As a result, researchers often resort to various approaches that
provide bounds or relaxations to approximate the solution to CC optimization
problems.

Our goal in this paper is to leverage the fact that in many applications, such
as those described earlier, $\delta$ is very close to zero. We use this
feature, together with large-deviations principles and extreme-value theory to
gain insights into the asymptotic behavior (as $\delta\downarrow 0$) of a
stylized family of CC problems and their relaxations. In order to provide a
more detailed description of our results and their implications, let us
provide a brief overview of various methods in CC optimization.

One method, introduced in Nemirovski and Shapiro (2006a), constructs a convex
restriction of the feasible region by upper bounding the reliability
constraint by the so-called conditional value at risk (abbreviated with CVaR;
henceforth we refer to this as the CVaR approach). Another option, popular
because of its ease of implementation, is the so-called scenario approach,
introduced by Calafiore and Campi (2006). The scenario approach consists of
replacing the probability constraint by several Monte Carlo-based scenarios,
_all of which_ need to satisfy the reliability constraint. If the undesirable
event needs to be kept rare, the number of required samples is large, which
can be prohibitive, though it is sometimes possible to reduce the number of
samples using importance sampling techniques, see Nemirovski and Shapiro
(2006b), Barrera _et al._ (2016), Blanchet _et al._ (2024).

Recently, several authors (see e.g. Subramanyam (2020), Nesti _et al._ (2019),
Tong _et al._ (2022a)) have suggested a different method, namely to
approximate the probability constraint using ideas from large-deviations
theory. The approach can be described as follows. Suppose that ${\bf
P}(G(x))\leq\delta$ is a safety constraint, $x$ is the optimization decision
variable, and $G(x)$ is an undesirable event. A (heuristic) large-deviations
approximation of ${\bf P}(G(x))$ can be obtained by replacing the integral of
the density of an underlying probabilistic model (say $f$) over the set $G(x)$
by a supremum, i.e.:

${\bf P}(G(x))=\int_{G(x)}f(u)\mathrm{d}u\approx\sup_{u\in
G(x)}f(u)=\exp\\{-\inf_{u\in G(x)}\log f(u)\\}.$ (1)

The previous approximation can often be made rigorous by introducing a
suitable asymptotic regime. Using (1), the chance-constrained program can be
approximated by a bi-level (or nested) optimization problem (involving $x$ at
the top level and $u$ at the bottom level).

The papers Subramanyam (2020), Nesti _et al._ (2019), and Tong _et al._
(2022a) focus on various computational aspects of this bi-level procedure. The
analysis in those papers requires sufficient tractability for the underlying
density $f$, and this is carried out in situations where (for example) the
randomness is a mixture of Gaussian distributions or where the randomness
arises from another specific light-tailed distribution. Further, this analysis
focuses on developing approximations as in (1) that result in computationally
tractable procedures.

Moving on to a more precise description of our contributions, this paper takes
a different perspective: we derive rigorous approximations like (1) allowing
both light and heavy-tailed distributions; in the latter case, the analogues
of (1) can be harder to formulate. In turn, we consider a sequence of CC
problems indexed by the safety parameter $\delta$, and complete the following
technical tasks:

  1. (a)

we show how large-deviations approximations for the solution of CC problems
emerge in the limit and analyze the scaling of the chance-constrained
program’s optimal value $v_{\delta}$ as $\delta\downarrow 0$;

  2. (b)

we compare $v_{\delta}$ both to the performance of the CVaR relaxation and the
scenario approach;

  3. (c)

we identify conditions for which the latter approaches are asymptotically
optimal or not.

An attractive feature of our analysis is that it is rooted in multivariate
extreme-value theory, allowing us to consider classes of distributions that
are either light-tailed or heavy-tailed with a rich dependence structure. We
restrict our analysis to a class of linear programs to keep it crisp.
Specifically, the insights that we obtained based on our technical development
can be summarized as follows:

  1. (1)

If the underlying random variables are light-tailed, then we show (Theorem 1
below) that the basic formulation of the chance-constrained program, given in
($\rm{CCP}_{\delta}$) below, converges to a convex optimization problem as
$\delta\downarrow 0$, and that an appropriate version of (1) emerges. While
this leads to a characterization of the growth rate of $v_{\delta}$, it is
nontrivial to construct a sequence of asymptotically optimal solutions from
the limiting optimization problem: essentially, the limit of the logarithmic
estimate (1) is rather rough to guarantee feasibility for fixed $\delta>0$.
However, the asymptotics are sharp enough to derive solutions with any
prescribed (non-negative) asymptotic optimality gap.

  2. (2)

Also in the context of light-tailed distributions, we show that both the CC
problem and its CVaR relaxation converge to the same convex program in the
regime $\delta\downarrow 0$. Thus, CVaR provides a sequence of convex
approximations that are both conservative and asymptotically optimal: their
optimal values scale in the same way as the actual profit. For a precise
formulation, see Theorem 2. In short, for light-tailed settings, the CVaR
relaxation is asymptotically exact.

  3. (3)

A similar, but more nuanced conclusion, once again in the light-tailed
environment, holds for the scenario approach: in this case, we show that the
sequence of random linear programs, indexed by the number of scenarios,
converges, after appropriate normalization, to a deterministic optimization
problem as the number of scenarios is large. However, it must be emphasized
that to match the behavior of ($\rm{CCP}_{\delta}$), the scenario approach
requires a careful choice in the number of scenarios $k$ as function of
$\delta$, which differs from the standard prescription provided by Calafiore
and Campi (2006); see also (6) below. Again, feasibility for the original
problem cannot be guaranteed in general. If we take the standard
recommendation for $k$ as given in Calafiore and Campi (2006) to guarantee
feasibility with high probability, the gap between the actual optimal profit
and the optimal profit resulting from the scenario approach diverges. We refer
to Theorem 3 and its discussion for details.

  4. (4)

If the underlying random variables are heavy-tailed (more precisely, if their
tails are of regular variation), we show (Theorem 4) that our basic chance-
constrained problem converges to a convex optimization problem after
appropriate renormalization. It turns out that (1) is not adequate here:
essentially logarithmic asymptotics as in (1) are not enough; a more precise
asymptotic estimate is needed. Fortunately, such an estimate is available in
the heavy-tailed case, using multivariate regular variation.

  5. (5)

A normalized sequence of CVaR constraints converges to a similar limit as the
chance constraints, except with a different value in a certain constant. This
makes CVaR asymptotically suboptimal in the heavy-tailed case. For example, if
the underlying random variables have a power law with index $\alpha$, CVaR
yields a value of the objective function which is a factor $1-1/\alpha$ of the
actual optimum, as $\delta\downarrow 0$. A detailed statement is Theorem 5
below.

  6. (6)

The scenario approach behaves even worse in the heavy-tailed case: the
appropriately normalized sequence of random optimization problems converges to
a limit, but this limit is still random, cf. Theorem 6 below. This is in
contrast to methods such as sample-average approximation (also known) as
empirical optimization in which, as the Monte Carlo sample size increases, the
sample-based approximation converges to the desired solution. Our results
suggest that in the heavy-tailed case the scenario approach must be modified
to obtain asymptotically optimal solutions. We utilize concepts from point-
process theory by Resnick (2007) and random set theory, see Molchanov (2017),
to establish these results.

The rest of this paper is organized as follows. We define our model and state
our assumptions in Section 2. In Section 3, we derive our results for light-
tailed distributions. Our results for heavy-tailed distributions can be found
in Section 4. In Section 5 we consider an example showing how it is possible
to model dependence between random variables in our framework for light-tailed
distributions, and how using approximations like (1) impacts feasibility.
Proofs of our results can be found in Section 7, which also develops a
general-purpose result that circumvents the use of epi-convergence.

##  2\. Problem description

Basic formulation. Consider the problem

$\displaystyle\mbox{maximize }c^{T}x\hskip 28.45274pt\mbox{ subject to }x\in
X,\mbox{ and }{\bf P}(\phi(x,L)>1)\leq\delta,$ ($\rm{CCP}_{\delta}$)

where $x\in\mathbb{R}^{m}$ is an $m$-dimensional decision vector, $c$ is a
vector with non-negative entries, $X$ is a compact convex subset of
$[0,\infty)^{m}$ which contains the set $[0,h]^{m}$ for some $h>0$, and $L\geq
0$ is an $n$-dimensional random vector. The elements of $L$ are often referred
to as risk factors and may be dependent. We assume there exist $n\times m$
matrices $A^{(i)},i=1,...,d$ with non-negative entries and set

$\phi(x,L):=\sup_{i=1,...,d}x^{T}A^{(i)}L.$ (2)

This formulation is general enough to model linear programs with random
coefficients in the constraints.

As mentioned earlier, the parameter $\delta>0$ denotes the risk tolerance. We
focus on the case where the reliability constraint is a rare event, i.e.
$\delta\downarrow 0$. We aim to understand how the optimal profit scales as
$\delta\downarrow 0$.

######  Definition 1 (Asymptotic optimality).

Let $v_{\delta}$ be the optimal value of $c^{T}x$ in ($\rm{CCP}_{\delta}$).
Let $(x_{\delta})$ be a sequence of $m$-dimensional vectors such that
$x_{\delta}\in X$ that is also asymptotically feasible, i.e.,

$\limsup_{\delta\downarrow 0}{\bf P}(\phi(x_{\delta},L)>1)/\delta\leq 1.$ (3)

We call $(x_{\delta})_{\delta>0}$ asymptotically optimal if

$c^{T}x_{\delta}\,\big{/}\,v_{\delta}\rightarrow 1,\qquad\text{as
}\delta\downarrow 0.$ (4)

We now introduce two formulations that are popular alternatives to
($\rm{CCP}_{\delta}$); we are interested in whether their solutions are
asymptotically optimal. The determination of asymptotic optimality is
facilitated by Theorems 1 and 4 below, as they establish the asymptotic
scaling of ($\rm{CCP}_{\delta}$) for light-tail and heavy-tail distributions,
respectively.

CVaR formulation. The problem ($\rm{CCP}_{\delta}$) is generally not convex.
For this reason, Nemirovski and Shapiro (2006a) introduce a convex relaxation
which typically (in particular in the case of non-atomic distributions)
coincides with the conditional value-at-risk

${\bf
CVaR}(\phi(x,L))=\inf_{\tau\in{\mathbb{R}}}\Big{[}\tau+\frac{1}{\delta}{\bf
E}[(\phi(x,L)-1-\tau)^{+}]\Big{]}.$ (5)

In particular, they propose to solve the optimization problem

$\displaystyle\mbox{maximize }c^{T}x\hskip 28.45274pt\mbox{ subject to }x\in
X,\mbox{ and }{\bf CVaR}(\phi(x,L))\leq 0.$ ($\rm{CVaR}_{\delta}$)

Let $v_{{\bf CVaR},\delta}$ be the optimal solution of the CVaR version of
($\rm{CCP}_{\delta}$). The CVaR constraint in ($\rm{CVaR}_{\delta}$) is more
conservative than the chance constraint in ($\rm{CCP}_{\delta}$), so $v_{{\bf
CVaR},\delta}\leq v_{\delta}$. However, as this optimization problem is
convex, it may be more tractable than ($\rm{CCP}_{\delta}$). Theorems 2 and 5
below show how $v_{{\bf CVaR},\delta}$ compares to $v_{\delta}$ as
$\delta\downarrow 0$.

Scenario formulation. Another approach towards solving chance-constrained
programs is to resort to Monte-Carlo methods, where the chance constraint is
replaced by constraints of the form $\phi(x,L^{({i})})\leq 1$, with
$L^{(1)},...,L^{(k)}$ a number of samples called scenarios. This results in
the optimization problem

$\displaystyle\mbox{maximize }c^{T}x\hskip 28.45274pt\mbox{subject to }x\in
X\mbox{ and }\phi(x,L^{(i)})\leq 1,\quad\forall i\in[k],$ ($\mathrm{SCP}_{k}$)

where $[k]:=\\{1,\ldots,k\\}$. To ensure feasibility of ($\rm{CCP}_{\delta}$)
with high enough probability, $k$ needs to be chosen appropriately as
$\delta\downarrow 0$. In Calafiore and Campi (2006) it is shown that the
choice

$k(\delta)=\frac{2}{\delta}\log\frac{1}{\beta}+2d+\frac{2d}{\delta}\log\frac{2}{\delta}$
(6)

gives a solution that is feasible with probability at least $1-\beta$. This
bound is valid for any function $\phi$. We show in Theorems 3 and 6 below that
it leads to conservative solutions in our setting: letting $X^{k(\delta)}$ be
a solution of ($\mathrm{SCP}_{k}$) with $k=k(\delta)$, then the random
variable $c^{T}X^{k(\delta)}/v_{\delta}$ always converges to 0 in probability
under our assumptions.

Plan. We analyze the above three families of optimization problems in the next
two sections by making appropriate assumptions on $L$. We will apply concepts
and techniques from multivariate extreme-value theory and large-deviations
theory. We begin with the case in which the elements of $L$ are light-tailed.

##  3\. Light tails

In this section, we analyze ($\rm{CCP}_{\delta}$), ($\rm{CVaR}_{\delta}$), and
($\mathrm{SCP}_{k}$) for a large class of distributions guaranteeing that all
moments of $L$ are finite. In particular, we assume that $Q(x):=-\log{\bf
P}(L>x)$ is multivariate regularly varying: there exists a constant $\beta>0$,
a slowly varying function $h$ and a non-negative continuous function
$\lambda:[0,\infty)^{n}\mapsto\mathbb{R}$ such that for $x\in[0,\infty)^{n}$,
with $q(r)=h(r)r^{\beta}$,

$Q(rx)/q(r)\rightarrow\lambda(x),\quad\text{as }r\rightarrow\infty.$ (7)

The convergence as stated above is pointwise. However, it turns out to hold
uniformly on compact sets because of regular variation, see Resnick (2007).
Note that if we divide by $|x|$ and use that $h(r|x|)/h(r)\rightarrow 1$ as
$r\rightarrow\infty$ for $|x|>0$, we obtain

$\lambda(x)=\lambda(x/|x|)\times|x|^{\beta}.$ (8)

This family of distributions includes multi-variate normal distributions (for
which $\beta=2$), exponential and gamma distributions ($\beta=1$) and Weibull
distributions with a sub-exponential tail ($\beta<1$). The dependence
structure can be rather arbitrary: the main form of dependence which is
relevant is to which extent large values in the sample are correlated, which
is encoded in $\lambda$. We give an example of $\lambda$ in Section 5.

We now apply large-deviations theory to approximate ($\rm{CCP}_{\delta}$) for
small $\delta$. Set for $b\in[0,\infty)^{n}$

$I(b):=\inf_{x}\big{\\{}\lambda(x):bx\geq 1\big{\\}}=|b|^{-\beta}I(b/|b|).$
(9)

In Section 7.1 we show that the set $\\{b:I(b)\geq 1\\}$ is convex, compact,
possesses a non-empty interior, and

$-\log{\bf P}(b^{T}L>r)\sim I(b)q(r).$ (10)

Here, $\sim$ means that the ratio of both functions converges to $1$ as
$r\rightarrow\infty$. Next, we use the principle of the largest term (see e.g.
Ganesh _et al._ (2004), Lemma 2.1) to approximate our chance constraint by

$-\log{\bf
P}(\max_{i}y^{T}A^{(i)}L>r)\sim\min_{i}I(y^{T}A^{(i)})q(r)=:J(y)q(r),$ (11)

making the change of variable $x=y/r$. Now, consider the problem

$\mbox{maximize }c^{T}y\hskip 28.45274pt{\mbox{ subject to }}y\geq 0\mbox{ and
}J(y)\geq 1.$ ($\rm{LT}$)

Formulation ($\rm{LT}$) is motivated by making the change of variable $x=y/r$
in ($\rm{CCP}_{\delta}$), replacing the probability in ($\rm{CCP}_{\delta}$)
by the approximation ${\bf
P}(\max_{i}y^{T}A^{(i)}L>r)=\exp(-J(y)q(r)(1+o(1))$, and choosing
$r=q^{-1}(\log 1/\delta)$. This ensures that the chance constraint, after
taking logarithms, has the desired behavior. To make this work, we assume that
$q$ is continuous for large $r$.

The problem ($\rm{LT}$) is tractable: using Sion’s minimax theorem, cf. Sion
(1958), strong duality holds, and since $\\{J(y)\geq 1\\}$ has non-empty
interior, ($\rm{LT}$) can be solved by standard Lagrangian methods. In Section
5, we give examples for $\lambda$ for which $I$ itself has an explicit
expression.

Our first result is that ($\rm{LT}$) determines the correct growth rate of
$v_{\delta}$:

######  Theorem 1.

Let $y^{*}$ solve ($\rm{LT}$). Then $v_{\delta}=c^{T}y^{*}/q^{-1}(\log
1/\delta)(1+o(1))$ as $\delta\downarrow 0$.

To prove this result, a verification of the pointwise convergence property
(11) is not sufficient. The standard way in the optimization literature to
prove a result of the type of Theorem 1 is through epi-convergence, see
Rockafellar and Wets (1998). In our proof, we circumvent this and follow a
more direct approach; we refer to Section 7 for details.

Based on Theorem 1, it is tempting to conjecture that $y^{*}/q^{-1}(\log
1/\delta)$ is asymptotically optimal, but this may not be the case, as the
asymptotic estimate (11) is not sharp enough to guarantee the asymptotic
feasibility condition (3), as is exemplified in Section 5. One may instead
consider $(1-\eta)y^{*}/q^{-1}(\log 1/\delta)$, which is asymptotically
feasible for any $\eta>0$ (this follows from the proof of the next theorem in
Section 7.2), and accept an optimality gap $\eta$.

In contrast, ($\rm{CVaR}_{\delta}$) produces an asymptotically optimal
sequence of solutions:

######  Theorem 2.

$v_{{\bf CVaR},\delta}/v_{\delta}\rightarrow 1$ as $\delta\downarrow 0$.

The reason why this is true is that the overshoot $\phi(y^{*}/q^{-1}(\log
1/\delta),L)-1$, conditional on $\phi(y^{*}/q^{-1}(\log 1/\delta)/,L))>1$,
converges to $0$. In addition, the optimizing value of $\tau$ in (5) converges
to $0$ as $\delta\downarrow 0$. Together, these two properties imply that the
constraint in the CVaR formulation ($\rm{CVaR}_{\delta}$) behaves
asymptotically the same as the probability in ($\rm{CCP}_{\delta}$). This
makes the CVaR formulation appealing when optimizing under rare event
constraints in a light-tail setting.

To assess ($\mathrm{SCP}_{k}$) in the region $k\rightarrow\infty$, let
$L^{(j)},j=1,...,k$ be $k$ i.i.d. samples of $L$. Using the insights from
Theorem 1, this motivates us to consider a scaled problem. We make a similar
change of variables $x=y/q^{-1}(\log k)$. We exploit that
$\phi(y/r,L)=\phi(y,L)/r$ and consider

$\displaystyle\mbox{maximize }c^{T}y\hskip 14.22636pt\mbox{subject to
}\phi(y,L^{(i)})\leq q^{-1}(\log k),\quad\forall i\in[k].$
(${\rm{SCP}\textnormal{-}\rm{LT}}_{k}$)

######  Theorem 3.

Let $Y^{(k)}$ solve (${\rm{SCP}\textnormal{-}\rm{LT}}_{k}$). Then, as
$k\to\infty$, $c^{T}Y^{(k)}/q^{-1}(\log k)\rightarrow c^{T}y^{*}$, almost
surely, with $y^{*}$ a solution of ($\rm{LT}$).

To recover the optimal growth rate for the optimal value $v_{\delta}$ of
($\rm{CCP}_{\delta}$), we see from Theorem 1 that it is required that $k$ is
proportional to $1/\delta$. As the choice $k(\delta)$ to guarantee feasibility
diverges faster than $1/\delta$, the choice $k=1/\delta$ may lead to an
asymptotically infeasible solution.

Summarizing our findings in the light-tail case, we see that the behavior of
the optimal value $v_{\delta}$ can be characterized using the rare-event
approximation ($\rm{LT}$). Because this approximation is based on logarithmic
asymptotics (11), the rescaled optimum of $y^{\ast}$ may not be feasible, even
for small $\delta$.

The sampling approach leads to a random solution $Y^{(k)}$, which, after
appropriate normalization, converges to the solution of ($\rm{LT}$). Apart
from the computational challenges that play a role if $k$ is large, this
procedure leads to a tension between feasibility and optimality: to get a
feasible solution, the choice $k(\delta)$ leads to an optimal profit which is
of smaller order than $v_{\delta}$ as $\delta\downarrow 0$.

Finally, the CVaR approach, while not as tractable as ($\rm{LT}$), and more
conservative than ($\rm{CCP}_{\delta}$), is asymptotically optimal in a light-
tail regime, reinforcing the fact that CVaR is one of the main approaches for
solving chance constraint programs. In the next section, we will examine if
these findings still hold if the underlying random variables are heavy-tailed.

##  4\. Heavy tails

Throughout this section, we suppose that $L$ has a distribution that is of
multivariate regular variation, which covers the case in which the elements of
$L$ have a power law. To describe our assumptions in a convenient form, we use
a representation which (apart from not taking logarithms) is different from
the one in Section 3, building on an equivalent characterization of regular
variation developed in Basrak _et al._ (2002)). Specifically, let $|\cdot|$ be
the ${\mathcal{L}}_{1}$ norm. Then $\bar{F}(r)={\bf
P}(|L|>r)=h(r)r^{-\alpha},\alpha>1$, with $h$ slowly varying (i.e.
$h(ar)/h(r)\rightarrow 1$ as $r\rightarrow\infty$ for $a>0$) and

${\bf P}\big{(}L/|L|\in\cdot\,\big{|}\,|L|>r\big{)}\rightarrow{\bf
P}(\Theta\in\cdot),\hskip 28.45274ptr\rightarrow\infty.$ (12)

The idea behind this representation is that, if we transform $L$ in polar
coordinates, the radius and angle become independent, conditional upon the
radius getting large. This formulation allows dependence between extreme
values of risk factors. For background on multivariate regular variation, see
Basrak _et al._ (2002), Resnick (2007).

We now construct an appropriately normalized version of ($\rm{CCP}_{\delta}$)
in the regime $\delta\downarrow 0$. Set $x=y/r_{\delta}$, with
$r_{\delta}:=\bar{F}^{-1}(\delta)=\inf\\{r:{\bf P}(|L|>r)\leq\delta\\}$. Now,
write the constraint in ($\rm{CCP}_{\delta}$) as ${\bf
P}(\phi(y,L)>r_{\delta})\leq\delta.$ We assume that $\bar{F}(t)$ is eventually
strictly monotone and continuous which makes the constraint equivalent to

${\bf P}(\phi(y,L)>r_{\delta})/{\bf P}(|L|>r_{\delta})\leq 1.$ (13)

Letting $\delta\downarrow 0$, using (12), we show in Section 7.4 that ${\bf
P}(\phi(y,L)>r_{\delta})/{\bf P}(|L|>r_{\delta})\rightarrow{\bf
E}[\phi(y,\Theta)^{\alpha}]$, and that ($\rm{CCP}_{\delta}$) can be
approximated by

$\mbox{maximize }y^{T}c\hskip 28.45274pt{\mbox{ subject to }}y\geq 0\mbox{ and
}{\bf E}[\phi(y,\Theta)^{\alpha}]\leq 1.$ ($\rm{HT}$)

An appealing feature of ($\rm{HT}$) is that the feasible region of this
problem is convex, which generally is not the case for ($\rm{CCP}_{\delta}$).
We prove the following theorem in Section 7.4:

######  Theorem 4.

If $y^{*}$ solves ($\rm{HT}$), then $y^{*}/r_{\delta}$ is asymptotically
optimal and $v_{\delta}=c^{T}y^{*}/r_{\delta}(1+o(1))$ as $\delta\downarrow
0$.

Unlike in the light-tail case, the value $y^{*}$ leads to an asymptotically
optimal solution as asymptotic feasibility is guaranteed. In contrast, in the
heavy-tailed case ($\rm{CVaR}_{\delta}$) is no longer asymptotically optimal:

######  Theorem 5.

$v_{{\bf CVaR},\delta}/v_{\delta}\rightarrow 1-1/\alpha$ as $\delta\downarrow
0$.

The reason behind the subobtimality is that the overshoot
$\phi(y,L)-r_{\delta}$, conditioned on $\phi(y,L)>r_{\delta}$ is not of
smaller order than $r_{\delta}$ when $\phi(y,L)$ has a heavy tail.

In the previous section, we observed that the sample approximation
($\mathrm{SCP}_{k}$) of ($\rm{CCP}_{\delta}$) has some optimality properties
in the light-tail case, since the space of feasible solutions converges, after
appropriate normalization, to a deterministic set which coincides with the
feasible region of the limiting optimization problem ($\rm{LT}$). In the
heavy-tail case, this is no longer true. In contrast, a random solution
emerges, making any optimality property elusive.

To see this, consider the simplest possible case, where the goal is to
maximize a scalar $x$ subject to the constraint that ${\bf}P(xL>1)\leq\delta$,
with $L$ also scalar. ($\mathrm{SCP}_{k}$) now can be formulated as maximizing
$x$ subject to $x\max_{i=1}^{k}L^{(i)}\leq 1$. If $L$ is light-tailed, we can
normalize $\max_{i=1}^{k}L^{(i)}$ such that it converges to a deterministic
constant (we have done this in a more general setting to obtain Theorem 3).
However, in the heavy-tail case, such a normalization is not possible. In
fact, $\max_{i=1}^{k}L^{(i)}/\bar{F}^{-1}(1/k)$ converges in distribution to a
random variable with a Fréchet distribution, see Resnick (2007), Nair _et al._
(2022) for background.

To describe our general technical result we use point-process machinery
commonly used in extreme-value theory, cf. Resnick (2007) as well as random
set theory as developed in (Molchanov and Strokorb (2016)). Let $N$ be a
Poisson random measure on the space $(0,\infty]\times S_{n-1}^{+}$, where
$S_{n-1}^{+}$ is the intersection of the unit sphere and the positive orthant.

Let $N$ be a Poisson random measure on $(0,\infty]\times S_{n-1}^{+}$ with
intensity measure $\zeta(\cdot)=\nu_{\alpha}(\cdot)\times{\bf
P}(\Theta\in\cdot)$, where $\nu_{\alpha}$ has density $\alpha u^{-\alpha-1}$.
This should be interpreted as follows: if $F$ is a subset of $(0,\infty]\times
S_{n-1}^{+}$, then $N(F)$ has a Poisson distribution with parameter
$\zeta(F)$.

Under our assumptions, the point process (with $\epsilon$ being the Dirac
measure)

$N_{k}(\cdot)=\sum_{i=1}^{k}\epsilon_{|L^{(i)}|/\bar{F}^{-1}(1/k),L^{(i)}/|L^{(i)}|}(\cdot)$
(14)

is converging in distribution to $N$; this is a direct application of Theorem
6.2 of Resnick (2007). Similar to (${\rm{SCP}\textnormal{-}\rm{LT}}_{k}$) we
perform a change of variables $x=y/\bar{F}^{-1}(1/k)$, exploit the homogeneity
of $\phi$ and define the random variable $V_{SCP}(k)$ to be the optimal value
of the optimization problem

$\displaystyle\mbox{maximize }c^{T}y\hskip 28.45274pt\mbox{subject to }y\geq
0\mbox{ and }\phi(y,L^{(i)})\leq\bar{F}^{-1}(1/k),\quad\forall i\in[k].$
(${\rm{SCP}\textnormal{-}\rm{HT}}_{k}$)

In Section 7.6 we apply concepts from random set theory (cf. Molchanov (2017))
to show

######  Theorem 6.

$V_{SCP}(k)$ converges weakly to the optimal value $V^{*}$ of the optimization
problem

$\displaystyle\mbox{maximize }c^{T}y\hskip 28.45274pt\mbox{subject to }y\geq
0\mbox{ and }\phi(y,z)\leq 1,\quad\forall z\in N.$ (15)

In conclusion, comparing the three main methods to (approximately) solve
($\rm{CCP}_{\delta}$) in the heavy-tailed case, the large-deviations
approximation ($\rm{HT}$) yields an asymptotically optimal solution, and is
characterized by a convex optimization problem. In contrast, the CVaR
formulation yields a solution that stays a factor $1-1/\alpha$ away from the
optimum, while the scenario approach does not converge to a deterministic
solution or approximation of ($\rm{CCP}_{\delta}$) at all.

##  5\. Example for light-tailed distributions

While the heavy-tail case directly leads to the convex program ($\rm{HT}$),
the light-tail case leads to a bi-level optimization problem involving the
function $I$ in (9), making it less straightforward to parse. In this section,
we provide a computationally tractable example for the light-tail case that
also provides insights into how the tail behavior of $L$ and its dependence
structure may impact the optimal solution of ($\rm{CCP}_{\delta}$) as
$\delta\downarrow 0$. Throughout this section, we assume $d=1,m=n$ and
consider the problem

$\displaystyle\mbox{maximize }c^{T}x\hskip 28.45274pt\mbox{ subject to }x\geq
0,\mbox{ and }{\bf P}\Big{(}\sum_{i=1}^{n}x_{i}a_{i}L_{i}\geq
1\Big{)}\leq\delta,$ (16)

for some constants $a_{i}>0,i=1,...,n$. Suppose that $L$ satisfies (7), with
$\lambda$ given by

$\lambda(x):=\ell(x^{\beta}):=\ell(x_{1}^{\beta},...,x_{n}^{\beta}):=\bigg{(}\sum_{i=1}^{n}x_{i}^{\beta\theta}\bigg{)}^{1/\theta},\hskip
28.45274pt\theta\geq 1.$ (17)

Here, $\ell(x)$ corresponds to the so-called Gumbel-Hougaard copula (see e.g.
Segers (2012)). For $\theta=1$, large coordinates in $L$ occur independently;
this is no longer true if $\theta>1$. For $\theta=\infty$, we take
$\ell(x)=\max_{i}x_{i}$, so that large values are perfectly correlated.

If $\gamma=\beta\theta\leq 1$, then
$I(b)=\min\\{(\sum_{i=1}^{n}x_{i}^{\gamma})^{1/\theta}:b^{T}x\geq 1\\}$ has as
optimal solution $x_{i^{*}}=1/b_{i^{*}}$ with $i^{*}=\arg\max b_{i}$ (with
arbitrary tie-breaking rule) and $x_{i}=0$ for $i\neq i^{*}$. Consequently,

$I(b)=|b|_{\infty}^{-\beta}.$ (18)

Thus, if $\gamma\leq 1$, the limiting approximation ($\rm{LT}$) of
($\rm{CCP}_{\delta}$) simplifies to

$\displaystyle\mbox{maximize }c^{T}y\hskip 28.45274pt\mbox{ subject to }y\geq
0,\mbox{ and }\max_{i}a_{i}y_{i}\leq 1,$ (19)

which has solution $y_{i}^{\ast}=1/a_{i}$.

If $\gamma=\beta\theta>1$, then we find using the Karush-Kuhn-Tucker (KKT)
conditions that the optimal $x_{i}$ in the definition of $I(b)$ equals
$b_{i}^{1/(\gamma-1)}/\sum_{j}b_{j}^{\gamma/(\gamma-1)}$, and after some
simplifications we obtain

$I(b)=\left(\sum_{i=1}^{n}b_{i}^{\gamma/(\gamma-1)}\right)^{-(\gamma-1)/\theta}.$
(20)

Thus, if $\gamma>1$, the limiting approximation ($\rm{LT}$) of
($\rm{CCP}_{\delta}$) simplifies to

$\displaystyle\mbox{maximize }c^{T}y\hskip 28.45274pt\mbox{ subject to }y\geq
0,\mbox{ and }\sum_{i=1}^{n}(a_{i}y_{i})^{\gamma/(\gamma-1)}\leq 1.$ (21)

Using the KKT conditions we obtain the solution

$y_{i}^{\ast}=(c_{i}/a_{i})^{\gamma-1}/\Big{(}\sum_{j=1}^{n}(c_{j}/a_{j})^{\gamma}\Big{)}.$
(22)

If the coordinates of $L$ are independent in the tail ($\theta=1)$, then the
nature of $I(b)$ changes at the point $\beta=1$ which is precisely the point
where the moment-generating function of $L$ becomes finite. For $\beta<1$, the
tails of $L_{i}$ are heavy-tailed Weibull, for which it is known that large
values of sums are caused by a single big jump, leading to the simple form of
the optimal solution.

Interestingly, the impact of the dependence structure, summarized by $\theta$,
is nontrivial, in the sense that the single-big jump intuition is not relevant
in our decision-making context if the tail of $L_{i}$ is not very heavy, and
the dependence is sufficiently strong, i.e. if $\beta>1/\theta$.

We mentioned before that asymptotic optimality in the light-tailed case is not
guaranteed as the solution of ($\rm{LT}$) after normalization may not be
feasible for the problem ($\rm{CCP}_{\delta}$). To illustrate this in more
detail in the context of our example, consider the case where $\theta=1$ and
$\beta<1$. Take also $q(u)=u^{\beta}$. In this case, the chance constraint in
(16), taking the scaled solution $x_{i}=1/(a_{i}q^{-1}(\log 1/\delta))$, leads
to the probability ${\bf P}\big{(}\sum_{i=1}^{n}L_{i}>q^{-1}(\log
1/\delta)\big{)}$. The distribution of $L$ is sub-exponential and we have the
estimate, as $\delta\downarrow 0$,

${\bf P}\Big{(}\sum_{i=1}^{n}L_{i}>q^{-1}(\log 1/\delta)\Big{)}\sim n{\bf
P}\big{(}L_{1}>q^{-1}(\log 1/\delta)\big{)}=n\delta,$ (23)

cf. Chapter 3 in Nair _et al._ (2022). Thus, our choice leads to a violation
probability of a factor $n$ larger than it should be to satisfy (16).

For $\beta=1$ we obtain a sum of i.i.d. exponential random variables, for
which

$\displaystyle{\bf P}\Big{(}\sum_{i\in[n]}L_{i}>q^{-1}(\log 1/\delta)\Big{)}$
$\displaystyle\sim{\bf P}\big{(}\mathrm{Gamma}(n,1)>(\log 1/\delta)\big{)}$
(24) $\displaystyle=\int_{u=(\log
1/\delta)}^{\infty}\frac{1}{\Gamma(n)}x^{n-1}{\mathrm{e}}^{-x}\mathrm{d}u$
(25) $\displaystyle=\frac{1}{n!}\mathrm{e}^{-(\log 1/\delta)}(\log
1/\delta)^{n}=\delta\frac{(\log 1/\delta)^{n}}{n!}.$ (26)

This behavior is even worse compared to $\beta<1$ as $\delta\downarrow 0$.

##  6\. Concluding comments

This paper provides the first systematic large-deviations analysis of chance-
constrained optimization problems. To facilitate the development of our
algorithmic and structural insights, we focus on a linear objective and a
safety constraint controlling (i.e. upper bounding by a fixed constant) the
$1-\delta$ quantile of the maxima of bi-linear functions (bi-linear both in
the multivariate risk factor and the decision parameter). We distinguish two
regimes according to the tails of the risk factor, namely, light and heavy
tailed regimes.

Under the light-tail assumption, the main takeaway from Section 3 is that the
CVaR formulation is a good approach to chance-constrained optimization in a
rare-event regime. The scenario approach and the large-deviations
approximation (both suitably modified as explained in Section 3) can be a good
approach to CC optimization, though for those methods, there is a nontrivial
trade-off between feasibility and optimality.

In the heavy-tailed regime, we explicitly show that the CVaR relaxation is no
longer asymptotically optimal: its solution remains a factor $1-1/\alpha$ away
from the optimum. The scenario approach behaves even worse, as explained in
Section 4. In contrast, the large-deviations approximation is asymptotically
optimal in this case, as it is based on a sharper version of the large-
deviations approximation (1).

There are several natural research questions that may be considered: (i) we
expect that large-deviations approximations, for both light-tails and heavy-
tails may still yield convex programming problems for functions $\phi$ which
may not be bilinear (ii) computational considerations, it will be of interest
to consider scaling regimes where also the dimension of the decision variable
$x$ and/or the dimension of $L$ may grow; we kept it fixed in this study;
(iii) the large-deviations insights we developed could be used to develop new
algorithms for chance-constrained optimization which are asymptotically
optimal in the rare event regime, and also produce close to optimal solutions
for fixed values of $\delta$, just like the notion of asymptotic efficiency
has led to novel algorithms in the area of rare-event simulation.

##  7\. Proofs

We prove our theorems using the following generic convergence result. There
are related results in the literature, cf. Bonnans and Shapiro (2000),
Rockafellar and Wets (1998), but our direct approach seems the most efficient
for our purposes.

######  Proposition 1.

Assume that $c\geq 0$, $g_{\ell}$ is non-decreasing in $y$ for each $\ell$ and
that $g_{\ell}\rightarrow g$ pointwise as $\ell\rightarrow\infty$, with $g$
continuous and strictly increasing. Let $v_{*}=\sup\\{c^{T}y:g(y)\leq 1,y\geq
0\\}$. Further, define $v_{\ell}=\sup\\{c^{T}y:g_{\ell}(y)\leq 1,y\geq 0$}.
Then, $v_{\ell}\rightarrow v_{*}$.

If $c=0$ the result is immediate. So, we assume that at least one entry of $c$
is positive. We first consider the case $v_{*}=\infty$. For every $a>0$ there
exists $y\geq 0$ such that $c^{T}x>a$. By continuity of $g$ there exists
$\gamma\in(0,1)$ such that $g(\gamma y)<1$ and $c^{T}y>a/2$. This implies that
$g_{\ell}(\gamma y)\leq 1$ and therefore $v_{\ell}\geq a/2$, for all $\ell$
large enough. Thus, we may assume that $v_{*}<\infty$. Next, observe that for
each $\eta<1$, $g_{\ell}(\eta y^{*})<1$ for $\ell$ sufficiently large.
Consequently, $c^{T}y_{\ell}\geq\eta c^{T}y^{*}$ for $\ell$ sufficiently
large, and by letting $\eta\uparrow 1$ we obtain that

$\liminf_{\ell\rightarrow\infty}c^{T}y_{\ell}\geq c^{T}y^{*}.$ (27)

To prove a corresponding asymptotic upper bound, we assume to the contrary
that for some $\varepsilon>0$ and a sub-sequence $(v_{\ell})$ (which we also
index by $\ell$) we have $v_{\ell}\geq v_{*}+2\varepsilon$.

We want to contradict the optimality of $y^{*}$. Let us partition
$c=(\bar{c},0)$, where the part $\bar{c}$ has strictly positive elements. We
then note that $v_{*}=\max\\{\bar{c}^{T}\bar{y}:g(\bar{y},0)\leq 1,\bar{y}\geq
1\\}$. Because $g_{\ell}$ is non-decreasing, we have that
$v_{\ell}=\max\\{\bar{c}^{T}\bar{y}:g_{\ell}(\bar{y},0)\leq 1,\bar{y}\geq
0\\}.$ So, we may assume without loss of generality that $c$ has strictly
positive elements. Now, let $y_{\ell}$ be an optimal solution corresponding to
$v_{\ell}$ (if the solution is not achieved, we can work with an
$\varepsilon/2$ solution). Since $g_{\ell}$ is non-decreasing, we can make the
coordinates of $y_{\ell}$ smaller, obtaining $y^{\prime}_{\ell}$ which remains
feasible and $c^{T}y^{\prime}_{\ell}=c^{T}y^{*}+\varepsilon$. Since the
elements of $c$ can be assumed to be strictly positive, the set
$C=\\{x:c^{T}y=c^{T}y^{*}+\varepsilon\\}$ is compact. Since the $g_{\ell}$’s
are non-decreasing, the convergence $g_{\ell}\rightarrow g$ occurs uniformly
on compact sets. By choosing $\ell$ sufficiently large we can construct
$y^{\prime\prime}_{\ell}$ (possibly by reducing the coordinates of
$y^{\prime}_{\ell}$) such that $g(y^{\prime\prime}_{\ell})\leq 1$,
contradicting the optimality of $y^{*}$. This completes the proof of this
proposition. $\square$

###  7.1. Proof of Theorem 1

We first show that the feasible region of ($\rm{LT}$) is convex and compact.
For this, it suffices to show that $\\{b:I(b)\geq 1\\}$ is convex and compact,
where $I$ is defined in (9). Let $b^{(1)}$ and $b^{(2)}$ be such that
$I(b^{(i)})\geq 1,i=1,2$. Let $a\in(0,1)$. For the desired convexity property,
we have to show that $I(ab^{(1)}+(1-a)b^{(2)})\geq 1$. Since $\lambda$,
defined in (7), is continuous, non-decreasing, and satisfies the relation
$\lambda(rx)=r^{\beta}\lambda(x)$, there exists an $x^{(a)}\geq 0$ such that
$I(ab^{(1)}+(1-a)b^{(2)})=\lambda(x^{(a)})$. Since
$(ab^{(1)}+(1-a)b^{(2)})x^{(a)}\geq 1$, we have that either
$b^{(1)}x^{(a)}\geq 1$ or $b^{(2)}x^{(a)}\geq 1$. Suppose for concreteness
that $b^{(1)}x^{(a)}\geq 1$. Then, $I(ab^{(1)}+(1-a)b^{(2)})\geq
I(b^{(1)})\geq 1$ implies the desired convexity property. In fact, a similar
argument as that provided in the proof of Proposition 1, shows that $I(\cdot)$
is continuous (using that the minimization defining $I(b)$ can be assumed to
be taken over a compact set); this implies that $\\{b:I(b)\geq 1\\}$ is
closed. The compactness of the set $\\{b:I(b)\geq 1\\}$ follows from the
relation $\lambda(rx)=r^{\beta}\lambda(x)$, $\beta>0$.

Next, we show that (10) holds, using standard textbook arguments from large-
deviation theory, see e.g. Dembo and Zeitouni (1998): note that, if
$B_{\delta}(x)$ is the $L_{\infty}$ ball of radius $\delta$ around $x$, the
assumption on $L$ implies that $-\log{\bf P}(L/r\in
B_{\delta}(x))\sim\inf_{u\in B_{\delta}(x)}\lambda(u)$. Therefore, using also
the continuity of $\lambda$, the same estimate holds for the closure of
$B_{\delta}(x)$, from which it follows that $L/r$ satisfies a weak large-
deviations principle with rate function $\lambda$ and speed $q(r)$. Since
$\lambda(ax)=a^{\beta}\lambda(x)$, $\lambda$ has compact level sets.
Therefore, $L/r$ satisfies a full large-deviations principle with rate
function $\lambda$ and speed $q(r)$. Since $b^{T}L$ is a continuous function
of $L$ and $\lambda$ is continuous, (10) follows from the contraction
principle.

We now turn to the proof of Theorem 1. We apply Proposition 1 with
$\ell=1/\delta$, $g_{\ell}(x)=(\log\delta)/\log{\bf P}(\phi(x,L)>q^{-1}(\log
1/\delta))$, and $g(x)=1/J(x)$. The pointwise convergence
$g_{\ell}(x)\rightarrow g(x)$ follows from (10). The fact that $g$ is strictly
increasing follows from the homogeneity property of $J$. $\square$

###  7.2. Proof of Theorem 2

Since $v_{{\bf CVaR},\delta}\leq v_{\delta}$, it suffices to prove an
asymptotic lower bound. For this, we utilize Theorem 1. Let $y^{*}$ be a
solution of ($\rm{LT}$) and for fixed $\gamma>0$, define $y:=(1-\gamma)y^{*}$.
Moreover, let $r_{\rm{LT},\delta}:=q^{-1}(\log 1/\delta)$. We will prove that
$x=y/r_{\rm{LT},\delta}$ will be feasible for ($\rm{CVaR}_{\delta}$), for
$\delta$ sufficiently small.

The optimal value $\tau=\tau_{\delta}$ in (5) with $x=y/r_{\rm{LT},\delta}$
solves a first-order optimality condition given by

$\delta={\bf P}\big{(}\phi(y,L)>r_{\rm{LT},\delta}(\tau_{\delta}+1)\big{)}.$
(28)

Since this equality is true for all $\delta>0$, this implies that
$r_{\rm{LT},\delta}(\tau_{\delta}+1)\to\infty$. Combining this with the tail
behavior of $\phi(y,L)$ given in Equation (11), we see that as
$\delta\downarrow 0$

$\log(1/\delta)=(1+o(1))J(y)q\big{(}r_{\rm{LT},\delta}(\tau_{\delta}+1)\big{)},$
(29)

and therefore

$r_{\rm{LT},\delta}(\tau_{\delta}+1)=q^{-1}\Big{(}(1+o(1))\frac{\log(1/\delta)}{J(y)}\Big{)}=J(y)^{-1/\beta}q^{-1}(\log
1/\delta)(1+o(1)).$ (30)

Recalling that $r_{\rm{LT},\delta}=q^{-1}(\log 1/\delta)$ by definition for
light-tailed distibutions, we see that

$\tau_{\delta}+1=J(y)^{-1/\beta}(1+o(1)).$ (31)

Now, recall $y=(1-\gamma)y^{*}$, so that
$J(y)=(1-\gamma)^{-\beta}J(y^{\ast})>1$ by (11), (9), as $J(y^{\ast})=1$. We
conclude that $\tau_{\delta}\rightarrow-\gamma$ as $\delta\downarrow 0$.

Next we use that the constraint in ($\rm{CVaR}_{\delta}$) with
$x=y/r_{\rm{LT},\delta}$ is equivalent to

${\bf
E}\big{[}\phi(y,L)/r_{\rm{LT},\delta}\mid\phi(y,L)>r_{\rm{LT},\delta}(\tau_{\delta}+1)\big{]}\leq
1.$ (32)

We compute that

$\displaystyle{\bf E}[\phi(y,L)/r_{\rm{LT},\delta}$
$\displaystyle\mid\phi(y,L)/r_{\rm{LT},\delta}>\tau_{\delta}+1]$
$\displaystyle\leq 1-\gamma/2+{\bf
E}[(\phi(y,L)/r_{\rm{LT},\delta}-(1-\gamma/2))^{+}\mid\phi(y,L)/r_{\rm{LT},\delta}>\tau_{\delta}+1]$
$\displaystyle=1-\gamma/2+\frac{\int_{1-\gamma/2}^{\infty}{\bf
P}(\phi(y,L)>r_{\rm{LT},\delta}u)\mathrm{d}u}{{\bf
P}(\phi(y,L)>(1+\tau_{\delta})r_{\rm{LT},\delta})}.$

For fixed $y$, we can write ${\bf P}(\phi(y,L)>z)=\exp\\{-\tilde{q}(z)\\}$,
with $\tilde{q}$ regularly varying of index $\beta>0$. Due to the Potter
bounds (see e.g. Bingham _et al._ (1987), Theorem 1.5.6), for every
$\varepsilon>0$ there exists a constant $z_{\varepsilon}$ such that, for
$u\geq s\geq 1/2$, and $r\geq z_{\varepsilon}$

$\tilde{q}(ru)/\tilde{q}(rs)\geq(1-\varepsilon)(u/s)^{\beta-\varepsilon}.$
(33)

Using this bound, for sufficiently small $\varepsilon$ (w.r.t. $\gamma$ and
$\beta$) and recalling that $\tau_{\delta}\rightarrow-\gamma$, it follows that

$\displaystyle\frac{\int_{1-\gamma/2}^{\infty}{\bf
P}(\phi(y,L)>r_{\rm{LT},\delta}u)\mathrm{d}u}{{\bf
P}(\phi(y,L)>(1+\tau_{\delta})r_{\rm{LT},\delta})}$
$\displaystyle\int_{1-\gamma/2}^{\infty}\exp\big{(}-\tilde{q}(r_{\rm{LT},\delta}u)+\tilde{q}((1+\tau_{\delta})r_{\rm{LT},\delta})\big{)}\mathrm{d}u$
$\displaystyle=\int_{1-\gamma/2}^{\infty}\exp\Big{(}-\tilde{q}((1+\tau_{\delta})r_{\rm{LT},\delta})\Big{(}\frac{\tilde{q}(r_{\rm{LT},\delta}u)}{\tilde{q}((1+\tau_{\delta})r_{\rm{LT},\delta})}-1\Big{)}\Big{)}\mathrm{d}u$
$\displaystyle\leq\int_{1-\gamma/2}^{\infty}\exp\Big{(}-\tilde{q}((1+\tau_{\delta})r_{\rm{LT},\delta})(1-\varepsilon)\Big{(}\big{(}\frac{u}{1+\tau_{\delta}}\big{)}^{\beta-\varepsilon}-1\Big{)}\Big{)}\mathrm{d}u$
$\displaystyle\rightarrow 0$

as $r_{\rm{LT},\delta}\to\infty$. We conclude that for every $\gamma>0$, the
choice $y=y^{*}(1-\gamma)$ satisfies (32) for sufficiently small $\delta$,
which implies that $\limsup_{\delta\downarrow 0}v_{{\bf
CVaR},\delta}/v_{\delta}\geq 1-\gamma$. Theorem 2 now follows by letting
$\gamma\downarrow 0$. $\square$

###  7.3. Proof of Theorem 3

Let $L^{(j)},j=1,...,k$ be $k$ i.i.d. samples of $L$. Combining (10) with
Theorem 1 of Resnick and Tomkins (1973), we have that

$\frac{1}{q^{-1}(\log k)}\max_{j=1,...,k}b^{T}L^{(j)}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\frac{1}{I(b)^{1/\beta}},\qquad k\to\infty.$ (34)

almost surely. We apply (34) to establish the following pointwise convergence
for $y\geq 0$:

$\displaystyle g_{k}(y)$ $\displaystyle:=\frac{1}{q^{-1}(\log
k)}\max_{j=1,...,k}\max_{i=1,...,d}y^{T}A^{(i)}L^{(j)}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\max_{i=1,...,d}1/I(y^{T}A^{(i)})^{1/\beta}=:g(y).$ (35)

Observe that $g(y)$ is strictly increasing by (9). The statement follows now
immediately from Proposition 1. $\square$

###  7.4. Proof of Theorem 4

We recall that $r_{\delta}=\bar{F}^{-1}(\delta)$. We first derive the
pointwise limit of ${\bf P}(\phi(y,L)>r_{\delta})/{\bf P}(|L|>r_{\delta})$ as
$\delta\downarrow 0$, using (12). Let
$\varepsilon\in(0,1/\max_{i}y^{T}A^{(i)}e)$, with $e$ a vector of ones. As we
assumed that the entries of all $A^{(i)}$ are non-negative,
$|L|\leq\varepsilon r_{\delta}$ implies that $\phi(y,L)\leq r_{\delta}$.
Consequently,

$\displaystyle\frac{{\bf P}(\phi(y,L)>r_{\delta})}{{\bf P}(|L|>r_{\delta})}$
$\displaystyle=\frac{{\bf P}(\phi(y,L)>r_{\delta};|L|>\varepsilon
r_{\delta})}{{\bf P}(|L|>r_{\delta})}$ (36) $\displaystyle={\bf
P}\big{(}\phi(y,L)>r_{\delta}\,\big{|}\,|L|>\varepsilon
r_{\delta}\big{)}\frac{{\bf P}(|L|>\varepsilon r_{\delta})}{{\bf
P}(|L|>r_{\delta})}$ (37) $\displaystyle={\bf
P}\Big{(}\phi\Big{(}y,\frac{L}{|L|}\Big{)}\frac{|L|}{r_{\delta}}>1\,\big{|}\,|L|>\varepsilon
r_{\delta}\Big{)}\frac{{\bf P}(|L|>\varepsilon r_{\delta})}{{\bf
P}(|L|>r_{\delta})}.$ (38)

Due to univariate regular variation, the ratio converges to
$\varepsilon^{-\alpha}$. Due to the multivariate regular variation property
(12), the conditional (on $|L|>\varepsilon r_{\delta}$) weak limit of
$(L/|L|,|L|/r_{\delta})$ is $(\Theta,\varepsilon T)$, with $T$ a random
variable with density $\alpha u^{-1-\alpha}$ independent of $\Theta$.
Consequently,

$\displaystyle\frac{{\bf P}(\phi(y,L)>r_{\delta})}{{\bf P}(|L|>r_{\delta})}$
$\displaystyle\rightarrow\varepsilon^{-\alpha}{\bf
P}(\phi(y,\Theta)\varepsilon T>1)$ (39)
$\displaystyle=\varepsilon^{-\alpha}\int_{1}^{\infty}\alpha u^{-1-\alpha}{\bf
P}(\phi(y,\Theta)\varepsilon>1/u)\mathrm{d}u$ (40)
$\displaystyle=\int_{0}^{1/\varepsilon}\alpha v^{\alpha-1}{\bf
P}(\phi(y,\Theta)>v)\mathrm{d}v$ (41) $\displaystyle=\int_{0}^{\infty}\alpha
v^{\alpha-1}{\bf P}(\phi(y,\Theta)>v)\mathrm{d}v$ (42) $\displaystyle={\bf
E}[\phi(y,\Theta)^{\alpha}].$ (43)

In the derivation above, the second equality follows from the change of
variable $v=1/\varepsilon u$, the third equality comes from the fact that
$\phi(y,\Theta)\leq\max_{i}y^{T}A^{(i)}e\leq 1/\varepsilon$, and the last
equality from a representation for the $\alpha$th moment of a non-negative
random variable.

To finish the proof, we define $k=1/\delta$, $g_{k}(y)={\bf
P}(\phi(y,L)>r_{\delta})/{\bf P}(|L|>r_{\delta})$ and $g(y)={\bf
E}[\phi(y,\Theta)^{\alpha}]$, and apply Proposition 1. $\square$

###  7.5. Proof of Theorem 5

The constraint in ($\rm{CVaR}_{\delta}$) with $x=y/r_{\delta}$ is equivalent
to

${\bf
E}\big{[}\phi(y,L)/r_{\delta}\mid\phi(y,L)>r_{\delta}(\tau_{\delta}+1)\big{]}\leq
1$ (44)

where $\tau_{\delta}$ is the optimizing value in (5). Using (12), it satisfies

$\displaystyle\delta={\bf P}(\phi(y/r_{\delta},L)-1>\tau_{\delta})$
$\displaystyle={\bf P}\big{(}\phi(y,L)>r_{\delta}(\tau_{\delta}+1)\big{)}$
$\displaystyle=(1+o(1)){\bf E}[\phi(y,\Theta)^{\alpha}]{\bf P}(|L|\geq
r_{\delta}(\tau_{\delta}+1)).$

As a result, we obtain that
$(\tau_{\delta}+1)r_{\delta}=\bar{F}^{-1}(\frac{(1+o(1))\delta}{{\bf
E}[\phi(y,\Theta)^{\alpha}]})\to\infty$ as $\delta\downarrow 0$. We substitute
this in the left-hand side in (44), and recall that
$\bar{F}^{-1}(\delta)=r_{\delta}$ by definition. Then, as $\delta\downarrow
0$,

$\displaystyle{\bf
E}\big{[}\phi(y,L)/r_{\delta}\mid\phi(y,L)>r_{\delta}(\tau_{\delta}+1)\big{]}$
$\displaystyle=\frac{\bar{F}^{-1}\Big{(}\frac{(1+o(1))\delta}{{\bf
E}[\phi(y,\Theta)^{\alpha}]}\Big{)}}{\bar{F}^{-1}(\delta)}\cdot{\bf
E}\left[\frac{\phi(y,L)}{r_{\delta}(\tau_{\delta}+1)}\,\Big{|}\,\frac{\phi(y,L)}{r_{\delta}(\tau_{\delta}+1)}>1\right]$
$\displaystyle\rightarrow({\bf
E}[\phi(y,\Theta)^{\alpha}])^{1/\alpha}\cdot\frac{\alpha}{\alpha-1}.$

The convergence of the first factor follows from regular variation. The
convergence of the second factor holds because
$\phi(y,L)/r_{\delta}(\tau_{\delta}+1)$, conditional upon exceeding 1,
converges to a Pareto distributed random variable on $[1,\infty)$ with
exponent $\alpha$. This convergence also holds in expectation; to show this we
can apply a Potter bound for the ratio of regularly varying functions as we
did in the proof of Theorem 2.

We now apply Proposition 1 with the values $k=1/\delta$, $g_{k}(y)={\bf
E}[\phi(y,L)/r_{\delta}\mid\phi(y,L)>r_{\delta}(\tau_{\delta}+1)]$, and
$g(y)=({\bf E}[\phi(y,\Theta)^{\alpha}])^{1/\alpha}\frac{\alpha}{\alpha-1}$.
We conclude that the normalized version of ($\rm{CVaR}_{\delta}$) converges to
the optimal value of

$\mbox{maximize }c^{T}y\hskip 28.45274pt{\mbox{ subject to }}y\geq 0\mbox{ and
}{\bf E}[\phi(y,\Theta)^{\alpha}]\leq(1-1/\alpha)^{\alpha}.$ (45)

Noting that ($\rm{HT}$) has a unique solution $y^{*}$, we conclude that the
optimal solution to (45) equals $y^{*}(1-1/\alpha)$ by the homogeneity of
$\phi$. Combining this with Theorem 4 completes the proof. $\square$

###  7.6. Proof of Theorem 6

Write $a_{k}:=\bar{F}^{-1}(1/k)$ for the normalizing constants in
(${\rm{SCP}\textnormal{-}\rm{HT}}_{k}$). To show convergence in distribution
of the optimal value $V_{SCP}(k)$ to $V^{\ast}$, it suffices to show that the
non-feasible set spanned by the constraints $\phi(y,L^{(i)}/a_{k})\leq 1$,
$i=1,\ldots,k$, converges weakly to the non-feasible set spanned by the
constraints $\phi(y,z)\leq 1$, $z\in N$, where $N$ is the Poisson point
process defined above (14). To do so, we rely on the theory of random sets, as
developed in Molchanov (2017).

To translate our setting to that of Section 4.2 in Molchanov (2017), define
the multi-function $M(L):=\\{y:\phi(y,L)\geq 1\\}$. By the homogeneity of
$\phi$, the closure of the set of non-feasible solutions for
(${\rm{SCP}\textnormal{-}\rm{HT}}_{k}$) can be written as

$U_{k}=\bigcup_{i=1}^{k}\Big{\\{}y:\phi\Big{(}y,\frac{L^{(i)}}{a_{k}}\Big{)}\geq
1\Big{\\}}=\frac{1}{a_{k}}\bigcup_{i=1}^{k}M(L^{(i)}),$ (46)

where we use the notation $bS=\\{x:bx\in S\\}$. Applying Theorem 4.2.9 and
4.2.10 of Molchanov (2017), we obtain that $U_{k}$ converges in distribution
to $U=\\{y:\phi(y,z)\geq 1\mbox{ for some }z\in N\\}$. We define for $t>0$ the
half-space $K_{t}=\\{y:c^{T}y>t\\}$. Therefore, we see that for $t>0$,

$\displaystyle{\bf P}(V_{SCP}(k)\leq t)$ $\displaystyle={\bf P}\big{(}\forall
y\in K_{t}\exists
i\in\\{1,\ldots,k\\}:\phi\big{(}y,L^{(i)}/a_{k}\big{)}>1\big{)}$ (47)
$\displaystyle={\bf P}\big{(}\forall y\in K_{t}\exists
i\in\\{1,\ldots,k\\}:\phi\big{(}y,L^{(i)}/a_{k}\big{)}\geq 1\big{)}$ (48)
$\displaystyle={\bf P}(U_{k}\cap K_{t}\neq\emptyset)\rightarrow{\bf P}(U\cap
K_{t}\neq\emptyset)={\bf P}(V^{*}\leq t).$ (49)

This completes the proof of Theorem 6. $\square$

## References

  * Alsenwi _et al._ (2019) M. Alsenwi, S. R. Pandey, Y. K. Tun, K. T. Kim, and C. S. Hong.  A chance constrained based formulation for dynamic multiplexing of eMBB-URLLC traffics in 5G new radio.  In _2019 International Conference on Information Networking (ICOIN)_ , pages 108–113. IEEE, 2019.
  * Barrera _et al._ (2016) J. Barrera, T. Homem-de Mello, E. Moreno, B. K. Pagnoncelli, and G. Canessa.  Chance-constrained problems and rare events: an importance sampling approach.  _Mathematical Programming_ , 157(1):153–189, 2016.
  * Basrak _et al._ (2002) B. Basrak, R. A. Davis, and T. Mikosch.  A characterization of multivariate regular variation.  _Annals of Applied Probability_ , pages 908–920, 2002.
  * Bingham _et al._ (1987) N. H. Bingham, C. M. Goldie, and J. L. Teugels.  _Regular Variation_.  Encyclopedia of Mathematics and its Applications. Cambridge University Press, 1987.
  * Blanchet _et al._ (2024) J. Blanchet, F. Zhang, and B. Zwart.  Efficient scenario generation for heavy-tailed chance constrained optimization.  _Stochastic Systems_ , 14(1):22–46, 2024.
  * Bonnans and Shapiro (2000) J. Bonnans and A. Shapiro.  _Perturbation analysis of optimization problems_.  Springer Series in Operations Research. Springer-Verlag, New York, 2000.
  * Calafiore and Campi (2006) G. Calafiore and M. C. Campi.  The scenario approach to robust control design.  _IEEE Transactions on Automatic Control_ , 51(5):742–753, 2006.
  * Dembo and Zeitouni (1998) A. Dembo and O. Zeitouni.  _Large deviations techniques and applications_ , volume 38 of _Applications of Mathematics (New York)_.  Springer-Verlag, New York, second edition, 1998.
  * Frank (2008) B. Frank.  Municipal bond fairness act.  In _110th Congress, 2d Session, House of Representatives, Report_ , pages 110–835, 2008.
  * Ganesh _et al._ (2004) A. Ganesh, N. O’Connell, and D. Wischik.  _Big queues_ , volume 1838 of _Lecture Notes in Mathematics_.  Springer-Verlag, Berlin, 2004.
  * Klein Haneveld _et al._ (2020) W. K. Klein Haneveld, M. H. van der Vlerk, and W. Romeijnders.  _Stochastic programming_.  Graduate Texts in Operations Research. Springer, Cham, 2020.  Modeling decision problems under uncertainty.
  * Luedtke _et al._ (2010) J. Luedtke, S. Ahmed, and G. L. Nemhauser.  An integer programming approach for linear programs with probabilistic constraints.  _Mathematical Programming_ , 122(2):247–272, 2010.
  * Molchanov (2017) I. Molchanov.  _Theory of random sets_ , volume 87 of _Probability Theory and Stochastic Modelling_.  Springer-Verlag, London, second edition, 2017.
  * Molchanov and Strokorb (2016) I. Molchanov and K. Strokorb.  Max-stable random sup-measures with comonotonic tail dependence.  _Stochastic Processes and their Applications_ , 126(9):2835–2859, 2016.
  * Mosca (2024) C. Mosca.  https://consultations.entsoe.eu/system-operations/methodology-for-performing-the-probabilistic-dimen/, 2024.
  * Nair _et al._ (2022) J. Nair, A. Wierman, and B. Zwart.  _The fundamentals of heavy tails—properties, emergence, and estimation_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 2022.
  * Nemirovski and Shapiro (2006a) A. Nemirovski and A. Shapiro.  Convex approximations of chance constrained programs.  _SIAM Journal on Optimization_ , 17(4):969–996, 2006a.
  * Nemirovski and Shapiro (2006b) A. Nemirovski and A. Shapiro.  Scenario approximations of chance constraints.  In _Probabilistic and Randomized Methods for Design under Uncertainty_ , pages 3–47. Springer, 2006b.
  * Nesti _et al._ (2019) T. Nesti, J. Nair, and B. Zwart.  Temperature overloads in power grids under uncertainty: A large deviations approach.  _IEEE Transactions on Control of Network Systems_ , 6(3):1161–1173, 2019.
  * Prekopa (1970) A. Prekopa.  On probabilistic constrained programming.  In _Proceedings of the Princeton Symposium on Mathematical Programming_ , volume 113, page 138. Princeton, NJ, 1970.
  * Prékopa (2003) A. Prékopa.  Probabilistic programming.  _Handbooks in Operations Research and Management Science_ , 10:267–351, 2003.
  * Resnick (2007) S. I. Resnick.  _Heavy-tail phenomena_.  Springer Series in Operations Research and Financial Engineering. Springer, New York, 2007.  Probabilistic and statistical modeling.
  * Resnick and Tomkins (1973) S. I. Resnick and R. J. Tomkins.  Almost sure stability of maxima.  _Journal of Applied Probability_ , 10(2):387–401, 1973.
  * Rockafellar and Wets (1998) R. Rockafellar and R. J.-B. Wets.  _Variational Analysis_.  Springer Verlag, Heidelberg, Berlin, New York, 1998.
  * Segers (2012) J. Segers.  Max-stable models for multivariate extremes.  _Revstat-statistical Journal_ , 10:61–82, 2012.
  * Shapiro _et al._ (2021) A. Shapiro, D. Dentcheva, and A. Ruszczynski.  _Lectures on Stochastic Programming: Modeling and Theory, Third Edition_.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2021.
  * Shi and Holzapfel (2021) D. Shi and F. Holzapfel.  Rare-event chance-constrained flight control optimization using surrogate-based subset simulation.  In _AIAA Scitech 2021 Forum_. American Institute of Aeronautics and Astronautics, Jan. 2021.
  * Sion (1958) M. Sion.  On general minimax theorems.  _Pacific Journal of Mathematics_ , 8(1):171 – 176, 1958.
  * Subramanyam (2020) A. Subramanyam.  _Chance-Constrained Programming: Rare Events_ , pages 1–6.  Springer International Publishing, Cham, 2020.
  * Tong _et al._ (2022a) S. Tong, A. Subramanyam, and V. Rao.  Optimization under rare chance constraints.  _SIAM Journal on Optimization_ , 32(2):930–958, 2022a.
  * Tong _et al._ (2022b) S. Tong, A. Subramanyam, and V. Rao.  Optimization under rare chance constraints.  _SIAM Journal on Optimization_ , 32(2):930–958, 2022b.