Quadratic genetic modifications: a streamlined route to cosmological simulations with controlled merger history

TL;DR

This paper introduces quadratic genetic modifications for cosmological simulations, enabling controlled variation of merger histories by adjusting local variance, which enhances the study of galaxy evolution.

Contribution

It extends existing genetic modification techniques to include quadratic features like local variance, providing a new method for controlling galaxy merger histories in simulations.

Findings

Developed an efficient algorithm for quadratic modifications.

Demonstrated control of variance in a toy model.

Clarified the formulation of genetic modifications.

Abstract

Recent work has studied the interplay between a galaxy's history and its observable properties using "genetically modified" cosmological zoom simulations. The approach systematically generates alternative histories for a halo, while keeping its cosmological environment fixed. Applications to date altered linear properties of the initial conditions such as the mean overdensity of specified regions; we extend the formulation to include quadratic features such as local variance, which determines the overall importance of smooth accretion relative to mergers in a galaxy's history. We introduce an efficient algorithm for this new class of modification and demonstrate its ability to control the variance of a region in a one-dimensional toy model. Outcomes of this work are two-fold: (i) a clarification of the formulation of genetic modifications and (ii) a proof of concept for quadratic modifications leading the way to a forthcoming implementation in cosmological simulations.

Equations

$$P(\mathbf{\delta }|\text{b})=\frac{P(\text{b}|\mathbf{\delta })P(\mathbf{\delta })}{P(\text{b})}.$$

$$P(\mathbf{\delta }|\text{b})\propto {\delta }_D(\text{A}\mathbf{\delta }-\text{b})\exp (-\frac{1}{2}{\mathbf{\delta }}^{†}{\text{C}}_0^{-1}\mathbf{\delta }),$$

$$P(\mathbf{\delta }|\text{b})\propto e^{-\frac{1}{2}(\mathbf{\delta }-\bar{\mathbf{\delta }}{)}^{†}{\text{C}}^{-1}(\mathbf{\delta }-\bar{\mathbf{\delta }})},\text{with}$$

$$\bar{\mathbf{\delta }}(\text{b})={\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}\text{b},\text{and}$$

$$\text{C}={\text{C}}_0-{\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}\text{A}{\text{C}}_0,$$

$${\mathbf{\delta }}_{\text{residual}}$$

$$={\mathbf{\delta }}_0-{\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}{\text{b}}_0.$$

$${\mathbf{\delta }}_1={\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}\text{b}+{\mathbf{\delta }}_{\text{residual}}.$$

$${\mathbf{\delta }}_1={\mathbf{\delta }}_0-{\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}(\text{A}{\mathbf{\delta }}_0-\text{b}).$$

$${\chi }^2\equiv ||\mathbf{\delta }|{|}_{{\text{C}}_0^{-1}}^2={\mathbf{\delta }}^{†}{\text{C}}_0^{-1}\mathbf{\delta }.$$

$$\begin{array}{rl}& \underset{\mathbf{\delta }}{\min }||\mathbf{\delta }-{\mathbf{\delta }}_0|{|}_{{\text{C}}_0^{-1}}^2,\\ & \text{subject to}\text{A}\mathbf{\delta }=\text{b}.\end{array}$$

$$\mathcal{L}\equiv (\mathbf{\delta }-{\mathbf{\delta }}_0{)}^{†}{\text{C}}_0^{-1}(\mathbf{\delta }-{\mathbf{\delta }}_0)+{\mathbf{\lambda }}^{†}(\text{A}\mathbf{\delta }-\text{b}),$$

$${\mathbf{\delta }}_1={\mathbf{\delta }}_0-{\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}(\text{A}{\mathbf{\delta }}_0-\text{b}),$$

$$\Delta {\chi }^2={\mathbf{\delta }}_1^{†}{\text{C}}_0^{-1}{\mathbf{\delta }}_1-{\mathbf{\delta }}_0^{†}{\text{C}}_0^{-1}{\mathbf{\delta }}_0.$$

$$\begin{array}{rl}& \underset{\mathbf{\delta }}{\min }||\mathbf{\delta }-{\mathbf{\delta }}_0|{|}_{{\text{C}}_0^{-1}}^2,\\ & \text{subject to}{\mathbf{\delta }}^{†}\text{Q}\mathbf{\delta }=q,\end{array}$$

$$\mathcal{L}=(\mathbf{\delta }-{\mathbf{\delta }}_0{)}^{†}{\text{C}}_0^{-1}(\mathbf{\delta }-{\mathbf{\delta }}_0)+\mu ({\mathbf{\delta }}^{†}\text{Q}\mathbf{\delta }-q),$$

$${\mathbf{\delta }}_1=(\text{I}+\mu {\text{C}}_0\text{Q}{)}^{-1}{\mathbf{\delta }}_0,\text{and}$$

$${\mathbf{\delta }}_0^{†}(\text{I}+\mu {\text{C}}_0\text{Q}{)}^{-1}\text{Q}(\text{I}+\mu {\text{C}}_0\text{Q}{)}^{-1}{\mathbf{\delta }}_0=q.$$

$${\mathbf{\delta }}^{\prime }={\left(\text{I}+(\mu +{\mu }^{\prime }){\text{C}}_0\text{Q}+\mu {\mu }^{\prime }{\left({\text{C}}_0\text{Q}\right)}^2\right)}^{-1}{\mathbf{\delta }}_0\text{,}$$

$${\mathbf{\delta }}^{†}\text{Q}\mathbf{\delta }={\mathbf{\delta }}_0^{†}\text{Q}{\mathbf{\delta }}_0+2{\mathbf{\delta }}_0^{†}\text{Q}\mathbf{ϵ}+\mathcal{O}({\mathbf{ϵ}}^2),$$

$$\mathbf{ϵ}=\mu {\text{C}}_0\text{Q}{\mathbf{\delta }}_0,\text{with}$$

$$\mu =\frac{1}{2}\frac{q-{\mathbf{\delta }}_0^{†}\text{Q}{\mathbf{\delta }}_0}{{\mathbf{\delta }}_0^{†}\text{Q}{\text{C}}_0\text{Q}{\mathbf{\delta }}_0}.$$

$${\mathbf{\delta }}_1=\prod _{j=0}^N\left(\text{I}+{\mu }_j{\text{C}}_0\text{Q}\right){\mathbf{\delta }}_0,$$

$${\mathbf{\delta }}_1$$

$$=\prod _{j=0}^{\infty }\exp \left({\mu }_j{\text{C}}_0\text{Q}\right)\equiv \exp \left(\alpha {\text{C}}_0\text{Q}\right){\mathbf{\delta }}_0,$$

$${\eta }_{\text{initial}}\equiv {\mathbf{\delta }}_{1,\text{initial}}^{†}\text{Q}{\mathbf{\delta }}_{1,\text{initial}}^{}-q.$$

$$N=N_{\text{initial}}\sqrt{\frac{{\eta }_{\text{initial}}}{{\eta }_{\text{target}}}},$$

$$\mathbf{ϵ}=-\mu {\text{C}}_0\text{Q}\mathbf{\delta }+\mu {\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}\text{A}{\text{C}}_0\text{Q}\mathbf{\delta },$$

$$\mu =\frac{1}{2}\frac{q-{\mathbf{\delta }}^{†}\text{Q}\mathbf{\delta }}{{\mathbf{\delta }}^{†}\text{Q}{\text{C}}_0{\text{A}}^{†}(\text{A}{\text{C}}_0{\text{A}}^{†}{)}^{-1}\text{A}{\text{C}}_0\text{Q}\mathbf{\delta }-{\mathbf{\delta }}^{†}\text{Q}{\text{C}}_0\text{Q}\mathbf{\delta }}.$$

$${\text{Q}}_{{\sigma }^2}=\frac{1}{R^2}{\text{W}}^{†}\left(R\text{I}-\text{1}\otimes \text{1}\right)\text{W}.$$

$${\tilde{F}}_{lm}={\delta }_{lm}\left(1-\exp \left[-\frac{1}{2}{\left(\frac{k_l}{k_f}\right)}^2\right]\right).$$

$$\equiv {\text{W}}^{†}\text{W}{\text{F}}^{†}{\text{Q}}_{{\sigma }^2}\text{F}{\text{W}}^{†}\text{W}$$

$$=\frac{1}{R^2}{\text{W}}^{†}\text{W}{\text{F}}^{†}{\text{W}}^{†}\left(R\text{I}-\text{1}\otimes \text{1}\right)\text{W}\text{F}{\text{W}}^{†}\text{W}\text{.}$$

$$\mathbf{ϵ}=\sum _i{\mu }_i{\text{C}}_0{\text{Q}}_i\mathbf{\delta },$$

Full text

††pubyear: 2017††pagerange: Quadratic genetic modifications: a streamlined route to cosmological simulations with controlled merger history–B

1 Introduction

Mergers and accretion are thought to play a key role in shaping the observed galaxy population; in the prevailing cosmological paradigm merger histories are in turn seeded from random inflationary perturbations. Numerical studies must make inferences about the galaxy population from a finite sample of such histories. Due to the limited computer time available, this generates a tension between resolution (for resolving the interstellar medium) and volume (for adequately sampling histories).

One attempt to sidestep this problem is to create and study a small number of carefully controlled tests of the relationship between a galaxy’s history and its observable properties. This has long been attempted in idealised, non-cosmological settings (e.g. Hernquist 1993; Di Matteo et al. 2005; Hopkins et al. 2012). More recently, Roth et al. (2016) proposed performing such tests within a fully cosmological environment by constructing a series of closely-related initial conditions with targeted “genetic modifications” (hereafter GMs). The formalism resembles that of constrained realisations (Bardeen et al. 1986; Bertschinger 1987; Hoffman & Ribak 1991, hereafter HR91) which generates realisations of Gaussian random fields satisfying user-defined constraints on initial densities, velocities or potentials (e.g. Bertschinger 2001). Simulations based on constrained realisations have been extensively applied to recreating the local universe using observed galaxy distributions as constraints (for recent examples see Heß et al. 2013; Wang et al. 2016; Sorce et al. 2016; Hoffman et al. 2017).

Despite a resemblance, genetically modified simulations are markedly different from constrained simulations. The process of GM involves creating multiple versions of the initial conditions, each with carefully selected small changes. By re-simulating each scenario it becomes possible to study how the changes affect the non-linear evolution of structure. For example, modifications can be chosen such that they enhance or suppress merger ratios in incremental steps and so vary a galaxy’s history in a systematic and controlled way. The first application of this technique in a hydrodynamic simulation was made by Pontzen et al. (2017); that work focuses on the response of a galaxy’s central black hole and its ability to quench star formation as the merger history is changed gradually. Unlike studies based on fully idealised merger simulations, the GM-based approach is able to capture the effects of gradual gas accretion from filaments which is essential when probing the balance between star formation and black hole feedback.

On a technical level, Pontzen et al. (2017) used multiple linear modifications to alter the merger history. Such a method requires human effort on two fronts: (i) to identify and track particles forming the merging substructures; and (ii) to tune the modifications and understand their effects on one another. For instance, GMs suppressing a merger tend to increase the mass of other nearby substructures, which complicates interpretation of the final results (see section 2.3 and figure 2 of Pontzen et al. 2017). Bypassing this behaviour would be possible by individually identifying all substructures and demanding the algorithm fix each one. However, the spiralling complexity of the setup makes this option unattractive.

Another possibility, which is the primary aim of the present paper, is to find a new type of modification which automatically suppresses the merger ratios of all large substructures in a target galaxy’s history. Such a modification would smooth the expected history while keeping its final mass and overall environment fixed. These modifications must be applicable to cosmological simulations, so our objective is an algorithm that remains tractable even working with fields on multidimensional grids. To achieve this goal, we start by clarifying the formulation of GMs (Section 2). We then expand the framework to quadratic modifications (Section 3), allowing control over the variance at different scales to tackle the problem of multiple mergers. We demonstrate the feasibility of our method on a one-dimensional model (Section 4); in forthcoming work we will demonstrate the implementation for a full 3D zoom simulation. Results are discussed in Section 5 and we conclude in Section 6.

2 Linear constraints and modified fields

In this Section, we contrast the method of constraints (HR91) against that of linear genetic modifications. The aim is to clarify the status of the latter as a building block for non-linear GMs, which are introduced in Section 3.

2.1 Constrained ensemble

We start by reviewing the construction of constrained ensembles (see bottom panel of Figure 1). In this case, constraints must be known a priori, i.e. independently of any specific realisation. Constrained ensembles are therefore particularly useful when using observations as external inputs to constrain numerical simulations.

Consider a Gaussian random field sampled at $n$ points to create a vector $\bm{\delta}$ with covariance matrix $\textbf{{C}}_{0}=\langle\,\bm{\delta}\,\bm{\delta}^{\dagger}\,\rangle$. The HR91 algorithm allows for an arbitrary number (denoted $p$) of linear constraints to be placed on $\bm{\delta}$; these can be expressed as $\textbf{{A}}\,\bm{\delta}=\textbf{{b}}$ where A is a $p\times n$ matrix and b is a length-$p$ vector.

We start by constructing the ensemble of all fields $\bm{\delta}$ satisfying the constraint for a chosen b, i.e. $P(\bm{\delta}|\textbf{{b}})$. Applying Bayes’ theorem, the probability reads

[TABLE]

Using the fact that $P(\bm{\delta})$ is Gaussian and disregarding normalization, this relation becomes

[TABLE]

where $\delta_{D}$ is the ($p$-dimensional) Dirac delta function.

This expression suggests a brute force sampling solution: we could draw many trial $\bm{\delta}$s from the original ensemble and keep only the ones satisfying the constraints (within some tolerance). This solution is, however, computationally inefficient. Making use of the fact that the Dirac delta function can be represented as the zero-variance limit of a Gaussian, we can instead derive the following results (Bertschinger 1987):

[TABLE]

where $\bar{\bm{\delta}}$ and C are the expectation and the covariance of the Gaussian distribution $P(\bm{\delta}|\textbf{{b}})$. By construction, all fields drawn from this distribution will satisfy the constraints ($\textbf{{A}}\bm{\delta}=\textbf{{b}}$).

Hoffman & Ribak (1991) pointed out a convenient shortcut for efficiently sampling from the distribution specified by Equation (2.1). Starting from a draw of the unconstrained ensemble, $\bm{\delta}_{0}$, we calculate $\textbf{{b}}_{0}=\textbf{{A}}\,\bm{\delta}_{0}$. One can then rewrite $\bm{\delta}_{0}$ as the sum of the mean field $\bar{\bm{\delta}}(\textbf{{b}}_{0})$ from Equation (2.1) and a residual term $\bm{\delta}_{\mathrm{residual}}$, defined by:

[TABLE]

From here, a draw from the constrained ensemble $\bm{\delta}_{1}$ can be generated by recombining the residuals with the corrected mean $\bar{\bm{\delta}}(\textbf{{b}})$:

[TABLE]

To verify this procedure draws samples $\bm{\delta}_{1}$ from the constrained distribution, one first writes the mapping from $\bm{\delta}_{0}$ to $\bm{\delta}_{1}$ in a single step:

[TABLE]

Then, by calculating $\langle\bm{\delta}_{1}\rangle$ and $\langle\bm{\delta}_{1}\bm{\delta}_{1}^{\dagger}\rangle$, it is possible to check that the ensemble has the correct mean and covariance from Equation (2.1). The fact that $\bm{\delta}_{1}$ is Gaussian follows from its construction as a linear transformation of $\bm{\delta}_{0}$. The underlying efficiency of this method is that the covariance matrix in Equation (2.1), does not depend on the value of b, allowing the $\bm{\delta}_{\mathrm{residual}}$ term to be the same for both expressions.

In summary, the HR91 algorithm creates a draw from the constrained ensemble in two steps, using the realisation $\bm{\delta}_{0}$ as an intermediate construction tool. It provides a computationally efficient way of generating Gaussian constrained fields.

2.2 Genetic modifications

We now turn to GMs (see upper panel of Figure 1) to constrast their formulation with that of constrained fields. The GM procedure can be summarized as follows:

Draw the unmodified realisation $\bm{\delta}_{0}$. 2. 2.

Define the modifications by choosing which properties of $\bm{\delta}_{0}$ are to be modified. Unlike in the constrained field case, this is accomplished with reference to specific features of the $\bm{\delta}_{0}$ realisation (e.g. the location and properties of particular haloes). This reflects how GMs are intended for constructing numerical experiments rather than for recreating observationally motivated scenarios. We focus first on linear modifications, i.e. of the form $\textbf{{A}}\,\bm{\delta}=\textbf{{b}}$. 3. 3.

Create the modified field (or multiple modified fields with different values of b). We require changes between fields to be as small as possible, which relies on the definition of a distance in field space. In the context of Gaussian fields, the only available metric is defined by the $\chi^{2}$ distance,

[TABLE]

Consequently, GMs can be formulated as finding the modified field solution of the following optimization problem:

[TABLE]

The problem is solved by minimising the Lagrangian

[TABLE]

where $\bm{\lambda}$ is a vector of size $p$ containing the Lagrange multipliers for each modification.

By differentiating to find critical points with respect to $\bm{\delta}$ and $\bm{\lambda}$, we obtain a system of two vector equations with the solution

[TABLE]

where $\bm{\delta}_{1}$ is the modified field.

Equation (10) has regenerated Equation (6) using a different motivation and derivation. To summarise:

• •

In the case of (6), $\bm{\delta}_{0}$ is an intermediate construct that is never used in a simulation; it only exists to aid finding $\bm{\delta}_{1}$, which is a sample from the distribution (2.1).

• •

In the case of (10), $\bm{\delta}_{0}$ and $\bm{\delta}_{1}$ are put on equal footing. They are both initial condition fields drawn from the original, underlying ensemble $P(\bm{\delta})$. The fact that the modifications (choice of A and b) depend on $\bm{\delta}_{0}$, as emphasised by Porciani (2016), does not impact this interpretation.

• •

We show in Appendix A that in the case of non-linear constraints, there is no joint expression for GMs and HR91, formalising their intrinsic difference.

GMs should therefore be seen as a mapping between fields of the same ensemble. A family of modified fields is generated by choosing multiple values for b; the resulting mapping between members of the family is continuous and invertible. These properties are highly valuable for providing controlled tests, allowing for systematic exploration of the effects of formation history on a galaxy.

While the algorithm makes the minimal changes to the field, $\bm{\delta}_{1}$ may still not be a particularly likely draw from $P(\bm{\delta})$ if the modifications are too extreme. To quantify the level of alteration, the relative likelihood of the two fields is given by $\exp{(-\Delta\chi^{2}/2)}$ with

[TABLE]

As long as $\Delta\chi^{2}$ stays small, we can regard the modified and unmodified fields as similarly likely draws from $\Lambda\text{CDM}$ initial conditions.

Turning $\Delta\chi^{2}$ into a precise quantitative statement about the relative abundance of a particular galactic history remains a topic for future research. It relies on knowing the detailed Jacobian relating the initial conditions to properties of the final galaxy. This can so far only be estimated, and only in simple scenarios such as small modifications to the halo mass (Roth et al., 2016). There are multiple possible modifications (i.e. choices of A and b) leading to a given effect in the target galaxy history (Porciani, 2016); some will carry a smaller $\Delta\chi^{2}$ cost than others. Finding the minimum-cost route to a given change in the non-linear universe is not the aim of GMs; to perform galaxy formation experiments, we only need to find one choice of modification with an acceptably small $\Delta\chi^{2}$ penalty.

3 Extension to quadratic modifications

The main aim of this paper is to formulate modifications that control the variance of a field. The variance on scales smaller than the parent halo scale relates to the number of substructures in haloes (Press & Schechter 1974; Bond et al. 1991), and is therefore a proxy for the overall importance of mergers. It is important to distinguish variance modifications of a region from alterations to the power spectrum. The power spectrum defines only the average variance over the entire box, and over all possible realisations. We propose on the other hand to modify the local variance, targeting one region of interest and making minimal changes to the remaining structures. Another way to picture this goal is as follows. In any one stochastic ensemble, two realisations might by chance have enhanced or reduced variance in an area. Our procedure aims to map between such realisations rather than to modify the underlying power spectrum.

Variance is quadratic in the field value and therefore the approach in Section 2.2 cannot be applied directly. One natural formulation of the problem is through a new minimisation problem (analogous to the original linear case):

[TABLE]

where Q is a $n\times n$ matrix and $q$ is a scalar. We can assume without loss of generality that Q is Hermitian. For a suitable choice of Q (see Section 4), $q$ specifies the variance of a chosen region.

Following a similar approach to the linear modifications, we introduce the Lagrangian

[TABLE]

where $\mu$ is a scalar Lagrange multiplier associated with the quadratic modification. Searching for critical points, we obtain two equations relating the modified field $\bm{\delta}_{1}$ and the multiplier:

[TABLE]

Equation (14) and (15) provide a closed system for $\mu$ and $\bm{\delta}_{1}$ given a target $q$. Unlike the linear case, the system can not be solved analytically. A possibility would be to solve Equation (15) numerically for $\mu$ but direct matrix inversions are prohibited due to their computational cost. One would therefore need to perform approximate matrix inversion at each step of a root-finding scheme for $\mu$, making the worst-case complexity of such method infeasible.

There are moreover deeper reasons why such procedures can not be straightforwardly adapted to GMs. In the linear case discussed above, we defined GMs as a continuous and invertible mapping. Both of these properties are lost when looking at Equations (14) and (15). First, it is not clear that Equation (15) has a real solution for $\mu$. Consequently a real-valued $\bm{\delta}_{1}$ may not exist111We note in passing that, since variance is a positive quantity, Q is a positive semi-definite matrix. By definition, $\textbf{{C}}_{0}$ is positive definite. These conditions ensure that $\bm{\delta}_{1}$ is unique if it exists – but they do not guarantee existence. for any chosen value of $q$.

Second, the relationship between $\bm{\delta}_{0}$ and $\bm{\delta}_{1}$ is asymmetric: if a new field $\bm{\delta}^{\prime}$ is constructed by taking $q$ back to its original value $q_{0}$, we will have

[TABLE]

for suitable choices of $\mu$ and $\mu^{\prime}$. To obtain a solution $\mu^{\prime}$ allowing recovery of the initial field ($\bm{\delta}_{0}=\bm{\delta}^{\prime}$), it must hold that $\textbf{{C}}_{0}\textbf{{Q}}\propto(\textbf{{C}}_{0}\textbf{{Q}})^{2}$. This will not generally be the case for our applications, and so we conclude that in general $\bm{\delta}^{\prime}\neq\bm{\delta}_{0}$. Such asymmetry would be problematic for GM; the sense of a unique ‘family’ of fields is lost.

The combination of computational intractability and loss of key properties for GMs lead us to focus on an alternate method. We describe next a Newton-like method which efficiently approximates a solution to the optimization problem, Equation (12), while reinstating the desired properties of the GM mapping.

3.1 Linearised solution

In this section, we restate the quadratic problem in a way that has a guaranteed solution and that generates a single family as a function of $q$. The trick is to make only infinitesimally small changes to the value of $q$, building up finite changes by following a path through field space that is locally minimal. This leads to an iterative procedure for quadratic genetic modifications, which we will demonstrate is both unique and computationally tractable.

3.1.1 One infinitesimal step

We start by defining the displacement $\bm{\epsilon}$ from the unmodified field $\bm{\delta}=\bm{\delta}_{0}+\bm{\epsilon}$; for sufficiently small changes we may then neglect $\mathcal{O}(\bm{\epsilon}^{2})$ terms. We will discuss in Section 3.1.3 how to practically decompose a macroscopic change into a series of such minor modifications.

At first order, the updated variance (or other quadratic property) is given by

[TABLE]

where we have assumed $\bm{\delta}$ is real and made use of the previously stated Hermitian assumption, $\textbf{{Q}}^{\dagger}=\textbf{{Q}}$. Having linearised the modification, we can now find an analytic solution for the displacement and the multiplier $\mu$:

[TABLE]

Equation (19) does not involve matrix inversions and can therefore be efficiently evaluated, even in a 3D cosmological simulation context.

3.1.2 Building finite changes by successive infinitesimal updates

We now want to construct a macroscopic change in the field by iterating the infinitesimal steps of Equation (18). Performing a finite number of steps $N$, the modified field reads:

[TABLE]

where $\mu_{j}$ is the Lagrange multiplier at step $j$. The value of each $\mu_{j}$ depends on how the fixed interval is divided, i.e. implicitly on $N$. In the limit of increasing number of steps, each individual $\mu_{j}$ becomes infinitesimally small and the final solution is

[TABLE]

where $\alpha=\sum_{j=0}^{\infty}\mu_{j}$ is the overall displacement and is finite. The right-hand side of Equation (21) defines the matrix exponential operator, which is guaranteed to exist and is invertible.

The matrix exponential is a useful formal expression to show that there is a unique result, but does not help computationally since the required value of $\alpha$ to reach the objective $\bm{\delta}_{1}^{\dagger}\textbf{{Q}}\,\bm{\delta}_{1}=q$ is unknown. In practice, we use the finite approximation, Equation (20). The $\mu_{j}$ at each step are chosen by targeting $N$ intermediate modifications linearly spaced between the starting value $q_{0}\equiv\bm{\delta}_{0}^{\dagger}\textbf{{Q}}\,\bm{\delta}_{0}$ and the target $q$. At each step, $\mu_{j}$ is calculated using Equation (19); $\bm{\epsilon}_{j}$ is deduced with Equation (18); and the field is updated, $\bm{\delta}\to\bm{\delta}+\bm{\epsilon}_{j}$.

3.1.3 Step choice for a practical

algorithm

When calculating Equation (20) as an approximation to Equation (21), the accuracy will increase with the number of steps $N$. One must choose a minimal $N$ (for computational efficiency) while ensuring that linearly approximating the modification at each step is sufficiently accurate.

We first perform the calculation with a fixed number of steps $N_{\text{initial}}$. This gives rise to an initial estimate for the modified field that we denote $\bm{\delta}_{1,\text{initial}}$. The error on the resulting modification can be characterised by the magnitude of $\eta_{\text{initial}}$, where

[TABLE]

Because second-order terms are neglected in the modification, the error term $\eta_{\text{initial}}$ should scale inverse-quadratically with the number of steps $N_{\text{initial}}$. We verified this behaviour numerically for a variety of fields and modifications. If $\eta_{\text{initial}}$ is smaller than a desired precision, $\eta_{\text{target}}$, we retain the initial estimate as our final output field. Otherwise, the calculation must be repeated; the required number of steps to achieve the target precision is inferred from the quadratic scaling as

[TABLE]

Note that $N_{\text{initial}}$ should be kept small to avoid unnecessary iterations; $N_{\text{initial}}=10$ has been chosen for our test scenarios below.

The final algorithm has a worst-case complexity of $\mathcal{O}(\eta_{\mathrm{target}}^{-1/2}\,n^{3})$, where $n$ is the number of elements in the field $\bm{\delta}$. The $n^{3}$ arises from matrix multiplications required to compute each step; in practice the matrices will be sparse either in Fourier space (for the covariance matrix) or in real space (for the variance Q matrix). Therefore, one can speed up the matrix multiplications by transforming back and forth from real to Fourier space, improving the complexity to $\mathcal{O}(\eta_{\mathrm{target}}^{-1/2}\,n\,\log n)$.

The final procedure shares numerous similarities with Newton methods, used in large-scale optimization (see Nocedal & Wright 2006 for a comprehensive review). It retains quadratic information in the objective and linear information in the modification at each step and has a quadratic rate of convergence to the solution.

3.2 Joint quadratic and linear

modifications

The algorithm above can be generalised to the case where we have both a quadratic modification and $p$ linear modifications of the form $\textbf{{A}}\,\bm{\delta}=\textbf{{b}}$. We first apply the linear modifications using Equation (8), then turn to the iterative quadratic modifications. However if Equation (20) is applied directly, the linear objective will no longer be satisfied; in other words we need to enforce $\textbf{{A}}\,\epsilon=\textbf{{0}}$ at each step. Constructing and solving the appropriate minimisation, expression (18) is replaced by

[TABLE]

where

[TABLE]

These results can be iterated to achieve the final modified field, in exactly the same way as for the pure-quadratic modification.

Despite the complexity of these expressions, the evaluation will remain $\mathcal{O}(\eta_{\mathrm{target}}^{-1/2}\,n\,\log n)$ for reasons discussed previously. To help interpret the method, there is a clear geometric meaning for each term, which we present in Appendix B.

4 Demonstration

In this Section we demonstrate our algorithm in a $n$-pixel, one-dimensional setting as a proof of concept and as a reference for future implementation on cosmological simulations. We choose an example red power spectrum, as typically encountered on the scales from which galaxies collapse. Specifically, we adopt $P(k)=P_{0}\,(k_{0}+k)^{-2}$, where $P_{0}$ is an arbitrary normalisation and $k_{0}=2\pi/n$, an offset that prevents divergence of $P(k)$ at $k=0$.

4.1 Defining an example modification

The framework developed in Section 3 can alter any property that is quadratic in the field by suitable choice of Q. We now specialise to the case that Q corresponds to the variance of a length-$R$ region of the field. We start by defining the windowing operator W as a rectangular matrix picking out the desired $R$ entries from the $n$ pixels in $\bm{\delta}$. To calculate the variance of the region, one then calculates $\bm{\delta}^{\dagger}\textbf{{Q}}_{\sigma^{2}}\,\bm{\delta}$ where $\textbf{{Q}}_{\sigma^{2}}$ can be written

[TABLE]

Here, I is the $R\times R$ identity matrix and 1 is a length-$R$ vector of ones. Expression (26) is readily verified by constructing $\bm{\delta}^{\dagger}\textbf{{Q}}_{\sigma^{2}}\,\bm{\delta}$ and seeing that it does boil down to the variance of the chosen region.

We wish to consider the field variance only on scales smaller than the region size (corresponding to substructures with mass lower than that of the parent halo). To achieve this, $\textbf{{Q}}_{\sigma^{2}}$ can be high-pass filtered; we use a standard Gaussian high-pass filter F where in Fourier space the elements of $\tilde{\textbf{{F}}}$ are given by

[TABLE]

Here, $k_{l}=2\pi l/n$ is the wavenumber of the $l$-th Fourier series element and $k_{f}$, the filtering scale, is defined in our case by $k_{f}=2\pi/R$. The most appropriate choice of filtering scales and shapes in the context of cosmological simulations will be discussed in a forthcoming paper.

In real space the matrix F is defined by $\textbf{{F}}=\textbf{{U}}^{\dagger}\tilde{\textbf{{F}}}\textbf{{U}}$ where U is the unitary Fourier transform matrix. Finally, to localise the target modification fully, we can re-window the matrix after smoothing. The operator $\textbf{{W}}^{\dagger}\textbf{{W}}$ achieves this by setting pixels outside the target window to zero. With this set of choices, the final quadratic objective is set by

[TABLE]

In practice, we never calculate the matrix Q explicitly but rather implement a routine to efficiently calculate $\textbf{{Q}}\,\bm{\delta}$ for any field $\bm{\delta}$, which is then used by the algorithm described in Section 3. The ability to bypass storing or manipulating Q is essential to permit the computation to operate on a 3D cosmological simulation.

4.2 Results

Figure 2 shows examples of modified fields obtained with our algorithm. We alter the variance of a region of width $R=100$ pixels enclosed by vertical lines, showing two quadratic modifications with the variance reduced by a factor 3 (light grey) and a factor 10 (dark blue). In both cases, the mean of the field is held fixed at the unmodified value (horizontal line). In the setting of a cosmological simulation, we expect to be able to fix the parent halo mass (through the mean value) while modifying the smoothness of accretion (through the variance).

We verified that these fields achieve the linear modification $\textbf{{A}}\,\bm{\delta}_{1}-\textbf{{b}}$ to within numerical accuracy and the quadratic modification $\bm{\delta}_{1}^{\dagger}\textbf{{Q}}\bm{\delta}_{1}-q$ to $\eta_{\text{target}}=10^{-6}$ accuracy. The heights of small-scale peaks inside the enclosed region are successfully reduced and brought closer to the mean value for the modified fields. Visually, it can be seen that the changes to the field are minimal, maintaining as much as possible of the structure of the unmodified field in the modified versions. This underlines how the analytic minimisation, Equation (12), and its refinement to a linearised procedure (Section 3.1) agrees with the intuitive sense of making minimal changes. The different versions of the field form a continuous family as illustrated by the smooth deformation when reducing the variance by different factors.

Despite the modification objective Q being strictly confined to the target region, modifications can be seen to “leak” outside (beyond the vertical lines). This effect, which is also seen in linear GMs, is an intentional aspect of the minimisation construction – any sharp discontinuities in the field value or its gradients would give rise to a power spectrum inconsistent with the ensemble. In this specific example, the leakage appears more significant to the left than to the right of the target region. In general, the algorithm is spatially symmetric but its effect in any given case is not.

5 Discussion

5.1 The advantage of quadratic over linear modifications

Pontzen et al. (2017) showed that using multiple linear modifications was sufficient to change the merger ratios in the history of a galaxy; substructures can be diminished or enhanced by manually modifying individual peak heights.

Nonetheless, we expect the new quadratic approach to bring considerable benefits when making such manipulations; the advantages are illustrated in Figure 3. The top panel shows a field representing the density in initial conditions expected for a halo. The field has a broad overdensity enclosed by vertical lines and two narrower peaks labelled (1) and (2). According to the excursion-set formalism (Bond et al., 1991), (1) and (2) will collapse to form two separate haloes that later merge. This, together with smooth accretion, will form the final halo.

Suppose we wish to generate a smoother accretion history by reducing the heights of peaks (1) and (2) while maintaining the large-scale overdensity. In the original approach, we use linear GMs to set the mean values of the peaks to the mean value of the broad overdensity (horizontal line). The middle panel of Figure 3 presents the resulting field. However, a number of problems arise when performing the alteration using this approach.

We had to identify (1) and (2) as the most interesting substructures and define specific modifications for each. In the context of N-body simulations, this requires manually identifying which particles of the initial conditions constitute each individual subhalo. 2. 2.

More importantly, spatially neighbouring modifications interact and create new substructures (peak labelled (3) in our example). One solution to prevent the appearance of new substructures could be to add a new linear objective forcing problematic regions such as (3) to remain unmodified. Identifying and mitigating side effects in this way adds a layer of complexity to the linear GM procedure. Depending on the specific problem and the number of modifications at play, time spent at this tuning phase can rise steeply.

On the other hand, a single quadratic modification can avoid these problems by defining a variance target across the region. The third panel of Figure 3 shows the same field with variance reduced by a factor 10 (using the method from Section 4). The two local peaks are successfully reduced in amplitude while conserving the remaining small-scale structure of the parent halo. By construction, the variance modification naturally avoids compensation problems inherent to linear GMs. For this reason, quadratic GMs provide a cleaner, streamlined way to control merger histories.

5.2 Multiple quadratic modifications

The formalism discussed so far applies a single quadratic modification to a field (possibly in combination with linear objectives). Simultaneously applying multiple quadratic modifications would allow one to act concurrently on two separate haloes, or to further fine-tune the merger history of a single object. For instance, decreasing the variance on intermediate scales while increasing on small scales should increase the frequency of minor mergers.

To study this generalisation, we introduce $i=1,\cdots,P$ quadratic modifications, each with matrix $\textbf{{Q}}_{i}$. For an infinitesimal update, the change in the field $\bm{\epsilon}$ is then given by

[TABLE]

where $\mu_{i}$ are the Lagrange multipliers associated with each modification. Equation (30) defines a system of $P$ equations to be solved. The resulting value of a specific $\mu_{i}$ depends on the whole set of $q_{i}$ and $\textbf{{Q}}_{i}$, i.e. modifications are interdependent.

In the same way as Section 3.1, the update (29) can be iterated to create finite changes. Performing $N$ steps, the modified field reads

[TABLE]

where $\mu_{ij}$ is the multiplier $\mu_{i}$ at step $j$. However in the limit that the number of steps $N\to\infty$, convergence to the matrix exponential,

[TABLE]

is only guaranteed if either the $\textbf{{Q}}_{i}$ commute with respect to $\textbf{{C}}_{0}$ (i.e. $\textbf{{Q}}_{i}\textbf{{C}}_{0}\textbf{{Q}}_{j}=\textbf{{Q}}_{j}\textbf{{C}}_{0}\textbf{{Q}}_{i}$) or each $\mu_{ij}$ is directly proportional to $\alpha_{i}$. Because $\alpha_{i}$s are not known in advance, the latter option is hard to arrange; the previously noted interdependence of the $\mu_{i}$s on all $q_{i}$ and $\textbf{{Q}}_{i}$ exacerbates the difficulty.

With our current algorithms, convergence to the matrix exponential is therefore only assured when the $\textbf{{Q}}_{i}$ matrices commute. The easiest way to arrange for the commutation is to use orthogonal modifications, i.e.

[TABLE]

Physically, this requirement can be achieved by imposing modifications that are spatially separated by a sufficient number of correlation lengths or address distinct Fourier modes. This condition even allows one to apply the formalism of Section 3.1 to each modification one-by-one and still converge to the correct overall matrix exponential of Equation (32). We leave the case of non-orthogonal multiple quadratic modifications to further work.

6 Conclusions

We have presented an efficient algorithm to modify the variance in a particular region of a Gaussian random field realisation, with the aim of manipulating initial conditions for cosmological simulations. The modification produces a field that is as close as possible to the original realisation. In this way it provides a route for controlled tests of galaxy formation where multiple versions of the same galaxy are simulated within a fixed cosmological environment, but with altered accretion history.

We argued that quadratic controls, as developed here, offer a useful complement to the existing linear technique (Roth et al., 2016). In particular, variance on different filtering scales relates to dark matter halo substructure and merging history (Press & Schechter 1974; Bond et al. 1991). The new algorithm can construct GM fields with simultaneous control on the mean value and filtered variance of a region (Figure 2). This provides a route to altering merger history and accretion over the lifetime of a given halo in a way that is more streamlined than modifying individual substructures (see Figure 3).

In both linear and quadratic GM, the algorithm searches for fields which are nearby in the sense of the $\chi^{2}$ distance measure. In the quadratic case, this definition is further refined: for large shifts in the control parameter $q$ (which represents the variance in our test cases), the path through field space is defined by following a series of small shifts. Each of these individually minimize the $\chi^{2}$ distance travelled. We demonstrated a formal convergence property for this series and argued that the approach is desirable for (a) returning a continuously-deforming field $\bm{\delta}$ as a function of the changing target variance $q$; (b) being reversible, so that returning the variance to its initial value also returns the field to its initial state; (c) being numerically tractable even for 3D zoom simulations.

In the process, we clarified the mathematical formulation of GM, carefully distinguishing it from the constrained ensemble of HR91 (see Figure 1 for an overview). The status of fields constructed in the two approaches is distinct – unlike constrained realisations, GMs should be seen as a mapping between two fields from the same ensemble. In the case of quadratic objectives such as variance, even the cosmetic similarities between constraints and modifications are lost (Appendix A).

The next step is an implementation of the new algorithm in a full $N$-body initial conditions generator, including on varying-resolution grids appropriate to zoom simulations. This will be presented in a forthcoming paper where we will evaluate the effectiveness of quadratic GM (alongside the existing linear technique) for constructing controlled tests of galaxy formation.

Acknowledgements

We thank Amélie Saintonge for helpful discussions and Hiranya Peiris for useful comments on a draft manuscript. MR acknowledges supports from the Perren Fund and the IMPACT fund. AP is supported by the Royal Society.

Appendix A Constructing constrained ensembles for quadratic

constraints

In Section 2, we contrasted the notion of a linearly-constrained ensemble (Section 2.1) against that of genetic modifications (Section 2.2). While conceptually different, constrained ensembles can be sampled using the HR91 procedure which can in turn be seen as applying suitable modifications to realisations from the unconstrained pdf.

In this Appendix we show that there is no such similarity between quadratically-constrained ensembles and quadratically modified fields. To put it another way, there is no HR91-like method for generating samples from a quadratically-constrained ensemble.

Following the same Bayesian argument as in Section 2.1, we define the quadratically-constrained ensemble for a fixed Q and $q$ by

[TABLE]

We will show that the modification procedure does not generate samples from the ensemble (34), even when Q and $q$ are known and fixed in advance.

We start by defining the alternative ensemble $P(\bm{\delta}_{1}|q)$ to be that sampled by drawing an unconstrained field from $P(\bm{\delta})$ and using the GM procedure to enforce the constraint $\bm{\delta}^{\dagger}\textbf{{Q}}\bm{\delta}=q$. In Section 3.1, the mapping $\bm{\delta}_{0}\to\bm{\delta}_{1}$ was given by

[TABLE]

where $\textbf{{C}}_{0}$ is the covariance matrix of the Gaussian distribution $P(\bm{\delta})$ and the value of $\alpha(\bm{\delta}_{0},q)$ is implicitly defined by the need to satisfy the quadratic constraint $\bm{\delta}_{1}^{\dagger}\,\textbf{{Q}}\,\bm{\delta}_{1}=q$.

To incorporate this implicit requirement to choose the correct value of $\alpha$ into an expression for the ensemble, we make use of Bayes’ theorem:

[TABLE]

Note that the constraint demands $P(q|\alpha,\bm{\delta}_{0})=\delta_{D}(\bm{\delta}_{1}^{\dagger}\textbf{{Q}}\bm{\delta}_{1}-q)$, where $\bm{\delta}_{1}$ and $\bm{\delta}_{0}$ are related by the condition (35). Writing out the normalisation condition for $P(\alpha|q,\bm{\delta})$ then gives

[TABLE]

Because Q and $\textbf{{C}}_{0}\textbf{{Q}}$ are positive semi-definite, $q$ is a monotonically increasing function of $\alpha$; there is only one value of $\alpha$ which satisfies the Dirac delta function on the right-hand-side. Consequently, we can perform the integration by a change of variables to yield

[TABLE]

Substituting this result back into Equation (36) and performing the integral over $\bm{\delta}_{0}$ using $P(\bm{\delta}_{1}|\alpha,\bm{\delta}_{0})=\delta_{D}(\bm{\delta}_{1}-\exp{(\alpha\,\textbf{{C}}_{0}\,\textbf{{Q}})}\bm{\delta}_{0})$, one obtains

[TABLE]

where normalisation factors depending only on $\textbf{{C}}_{0}$ have been dropped. This expression no longer has any explicit reference to $\bm{\delta}_{0}$, which was our primary aim. It can now be compared with the distribution for a true constrained ensemble, Equation (34). The two distributions appear different (as expected given our claim of inequivalence), but the comparison is complicated by the unsolved integral over $\alpha$ which obscures the content of the expression.

We can see that this integral will never regenerate the true quadratic constrained ensemble by taking an illustrative example. Consider a three-dimensional field $\bm{\delta}_{1}=(x,y,z)$ with unit power spectrum ($\textbf{{C}}_{0}=\textbf{{I}}$). Let us further choose an explicit form for Q such that

[TABLE]

Inserting these results into Equation (39) gives

[TABLE]

where we have made the substitution $\beta=e^{-\alpha}$. The integral over $\beta$ has an analytical solution using the further substitution $t=(\beta z)^{2}/2$ and introducing the modified Bessel function of the second kind

[TABLE]

Equation (39) can then be evaluated explicitly to obtain

[TABLE]

For comparison, the quadratic constrained ensemble in this example is given by

[TABLE]

The distributions defined by (43) and (44) have identical $x$-dependence. This is a general property: degrees of freedom for which Q has a null direction are unconstrained and, similarly, left unchanged by our GM transformation. The distribution generated by these degrees of freedom will therefore coincide at all times with the constrained ensemble. However, the $y$ and $z$ dependences differ between Equations (43) and (44). In general, non-null directions in field space will behave differently between the GM and constrained ensemble cases.

The result establishes that our formulation of quadratic GM as a matrix exponential mapping does not reproduce a quadratically-constrained ensemble when used analogously to the HR91 algorithm. A similar argument allows one to verify that applying the alternative non-linear modification specified by Equation (14) also fails to regenerate the constrained result. In fact, one can go even further and write a general power series expansion for the mapping between $\bm{\delta}_{0}$ and $\bm{\delta}_{1}$, writing

[TABLE]

without further specifying the power series coefficients $A_{i}$. Even in this case, which generalises away from a specific mapping, it is not possible to generate a constrained ensemble from the modification procedure. This underlines that modifications and constraints need to be regarded as entirely separate procedures. Only in the linear case do they appear to be cosmetically related.

Appendix B Geometrical interpretation

Throughout the main text, we used fields sampled at a finite number of points $n$; the resulting algorithms can therefore be interpreted geometrically as acting on vectors in an $n$-dimensional space. For instance, Roth et al. (2016) noted that the linear GM procedure is equivalent to an orthonormal projection of the unmodified field onto a subspace defining the modification objective (see their appendix A). In this Appendix we provide the geometric interpretation for our extended formulation of GM.

For the purposes of visualising the connection, we use fields with only two samples, $\bm{\delta}=(x,y)$. The arguments of this Appendix generalise to higher dimensions but are easiest to visualise with $n=2$. Figure 4 shows the resulting two-dimensional geometry in terms of the displacements $\Delta x$ and $\Delta y$ from the unmodified field. By construction, the unmodified field is at the origin.

The left panel shows the elliptical geometry generated by the covariance matrix in the $\Delta x$ – $\Delta y$ plane; specifically, the ellipses are of constant distance $\Delta s^{2}$ from the origin, where

[TABLE]

The linear objective $\textbf{{A}}\bm{\delta}=\textbf{{b}}$ defines a line in two dimensions. The modification consists of finding the value of $(\Delta x,\Delta y)$ lying on the line, while minimising $\Delta s^{2}$. Since $\Delta s^{2}$ is measured in terms of $\textbf{{C}}_{0}^{-1}$, the solution does not correspond to the closest point on the page but to the point at which a covariance ellipse is tangent to the modification line.

Similarly, the quadratic modifications (right panel of Figure 4) are associated with ellipses of constant $q=(x,y)\,\textbf{{Q}}\,(x,y)^{\top}$. These targets are shown as dotted lines; note that they are centred on $(x,y)=(0,0)$ and therefore appear offset from the origin in the $\Delta x$ – $\Delta y$ plane.

The right panel of Figure 4 also illustrates the algorithm for finding the modified field with a simultaneous quadratic and linear objective. For visual clarity, an unrealistically small ($N=3$) number of steps are taken. We start by defining three intermediate ellipses (red-dotted) between the value of the modification at the unmodified field and the target. As explained in Section 3.2, we first apply the global linear modifications from Equation (8)

[TABLE]

The algorithm then iterates the step $\bm{\epsilon}$ defined by Equation (24)

[TABLE]

These operations can be understood geometrically as:

(a) A projection of the current field on the linear modification. This term is similar to the case with linear modifications only.

(b) A displacement along the normal of the ellipse at the current field value. This term is towards the next intermediate ellipse.

(c) The projection of the previous term back onto the linear modification to ensure that both are always satisfied.

Term (c) ensures that the current field at the end of each step always lies on the linear constraint. Term (a) therefore vanishes after the first step; it is an overall offset that needs to be applied only once. Together, (b) and (c) are locally orthogonalizing the quadratic modification with respect to the global linear modification. The orthogonalisation must be repeated at each step since the local linearisation changes as we progress towards the final value of $q$.