avalable Rubin-LSST France Nov 2021 talk

## Weak Gravitational lensing

### How do we make the most of the available data?

HSC cosmic shear power spectrum
HSC Y1 constraints on $(S_8, \Omega_m)$
(Hikage,..., Lanusse, et al. 2018)
• Data reduction

• Compute summary statistics based on 2pt functions,
e.g. the power spectrum

• Run an MCMC to recover a posterior on model parameters, using an analytic likelihood $$p(\theta | x ) \propto \underbrace{p(x | \theta)}_{\mathrm{likelihood}} \ \underbrace{p(\theta)}_{\mathrm{prior}}$$
Main limitation: The lensing convergence field is inherently and significantly non-Gaussian
The two-point statistics do not fully capture the non-Gaussian information encoded in the peaks of the matter distribution

## How to maximize the information gain?

(Ajani, et al. 2020)
• Approaches based on measuring high-order correlations to access the non-Gaussian information
e.g. the Peaks count

• Simulation-based approaches

• Costly as they require a large number of simulations
• Intractable for more than 3 or 4 cosmological parameters.
• ## Numerical Differentiation

$$$\left.\frac{df(x)}{dx}\right|_{x_1} \approx \frac{f(x_1+h)-f(x_1)}{h}$$$
Flaws :
• It’s numerically very unstable

• It’s very expensive in term of simulation time

• Computes an approximation

Different approach :

### Automatic Differentiation and Gradients in TensorFlow

• Automatic differentiation allows you to compute analytic derivatives of arbitraty expressions:
If I form the expression $y = a * x + b$, it is separated in fundamental ops: $$y = u + b \qquad u = a * x$$ then gradients can be obtained by the chain rule: $$\frac{\partial y}{\partial x} = \frac{\partial y}{\partial u} \frac{ \partial u}{\partial x} = 1 \times a = a$$

• To differentiate automatically, TensorFlow remember what operations happen in what order during the forward pass and traverses this list of operations in reverse order to compute gradients.

• Derivative as functions

• No numerical approximation

• High speed

## How do we simulate the Universe in a fast and differentiable way?

### The Particle-Mesh scheme for N-body simulations

The idea: approximate gravitational forces by estimating densities on a grid.
• The numerical scheme:

• Estimate the density of particles on a mesh
$\Longrightarrow$ compute gravitational forces by FFT

• Interpolate forces at particle positions

• Update particle velocity and positions, and iterate

• Fast and simple, at the cost of approximating short range interactions.
$\Longrightarrow$ Only a series of FFTs and interpolations.

### Introducing FlowPM: Particle-Mesh Simulations in TensorFlow


import tensorflow as tf
import flowpm
# Defines integration steps
stages = np.linspace(0.1, 1.0, 10, endpoint=True)

initial_conds = flowpm.linear_field(32,        # size of the cube
100,       # Physical size
ipklin,    # Initial powerspectrum
batch_size=16)

# Sample particles and displace them by LPT
state = flowpm.lpt_init(initial_conds, a0=0.1)

# Evolve particles down to z=0
final_state = flowpm.nbody(state, stages, 32)

# Retrieve final density field
final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions),
final_state[0])

with tf.Session() as sess:
sim = sess.run(final_field)

• Seamless interfacing with deep learning components

Flow N-body PM simulation:

• Fast (limited number of time steps while enforcing the correct linear growth)

• Cannot give accurate halo matter profiles or matter power spectrum

The PGD idea : mimics the physics that is missing
$\Longrightarrow$ Halo virialization
(Biwei Dai et al. 2018)

### Mocking the weak lensing universe: The Born Approximation

• Numerical simulation of WL features rely on ray-tracing through the output of N-body simulations​

• Knowledge of the Gravitational potential and accurate solvers for light ray trajectories is computationally expensive

• Born approximation , only requiring knowledge of the density field, can be implemented more efficiently and at a lower computational cost

$$$\kappa_{born}(\boldsymbol{\theta},\chi_s)= \frac{3H_0^2 \Omega_m}{2c^2} \int_0^{\chi_s} d\chi \frac{\chi}{a(\chi)} W(\chi,\chi_s) \delta(\chi \boldsymbol{\theta},\chi).$$$

### Proof of Concept

• Convergence map at z=1.0, based on a 3D simulation of $128^3$ particles for side. The 2D lensing map has an angular extent of $5^{\circ}.$
• Angular Power Spectrum $C_\ell$

### Proof of Concept: Peak counts

• Local maxima in the map

• Non-Gaussian information

• Trace overdense regions
• ### Proof of Concept: $l_1$norm

• We take the sum of all wavelet coefficients of the original image $\kappa$ map in a given bin $i$ defined by two values, $B_i$ and $B_{i+1}$,

• Non-Gaussian information

• Information encoded in all pixels

• Avoids the problem of defining peaks and voids

## A first application: Fisher forecasts

We are testing the simulations to reproduce a LSST Y1-like setting
(Preliminary results, presented on behalf of LSST DESC)

$$$F_{\alpha, \beta} =\sum_{i,j} \frac{d\mu_i}{d\theta_{\alpha}} C_{i,j}^{-1} \frac{d\mu_j}{d\theta_{\beta}}$$$

• Use Fisher matrix to estimate the information content extracted with a given statistic

• Derivative of summary statistics respect to the cosmological parameters.
e.g. the $\Omega_c$, $\sigma_8$

• Fisher matrices are notoriously unstable
$\Longrightarrow$ they rely on evaluating gradients by finite differences.

• They do not scale well to large number of parameters.

## A first application: Fisher forecasts

• Use Fisher matrix to investigate the constraining power of different statistics
• ## A first application: Fisher forecasts

• Use Fisher matrix to estimate the information content extracted with a given statistic for more cosmological parameters simultaneously
• Fisher constraints from the full peak count statistics

# Conclusion

### Conclusion

What can be gained by simulating the Universe in a fast and differentiable way?

• Investigate the constraining power of various map-based higher order weak lensing statistics and control the systematics