Inspiration: Carvajal-Bohorquez et al. 2025

In[]:=
psdCB[ω_,σ_,τ_,α_]:=σ/(1+(ω*τ)^α)
In[]:=
psdCBmod[ω_,σ_,τ_,α_]:=σ/(1+RealAbs[(ω*τ)]^α)
The authors draw SFH samples from this PSD in the same way we draw modes for the Initial Conditions (Timmer & König 1995).​
​- Components are drawn from the unit Gaussian for the real and imaginary components at a series of positive ω​
​- These modes are multiplied by the square root of the PSD at that frequency​
​- The negative frequency components are set to the complex conjugate of these sampled modes f(-ω) = (f(ω))*, to ensure the signal is real​
​- Time-series data is obtained via the inverse Fourier transform​
​Note: for α == 2 (and pretty much only for α==2), There is a simple form for this correlation function
In[]:=
Simplify[InverseFourierTransform[psdCBmod[ω,σ,Subscript[τ,0],2],ω,t],{t>0,Subscript[τ,0]>0}]
Out[]=
-
t
τ
0

π
2
σ
τ
0
However, this would require either storing a high-resolution SFH per halo, or running a smoothing kernel (FFT) for every halo. Furthermore, since halos merge, we would need to be able to condition the SFH in chunks, requiring a large number of operations due to the high resolution.​
​​
​Instead, we can choose a few timescales which we care about (e.g 10Myr, 100Myr, snapshots) and *only* sample those, using the PSD to obtain autocorrelation functions, and filling a covariance matrix.
​

Defining the Covariance Matrix

In order to sample our star formation histories, first we need to define the covariance matrix containing all the variables we wish to sample (current snapshot) and condition (previous/descendant snapshot) on. This means we need the auto- and cross-correlation functions for each timescale we care about, for both correlations within the snapshot and correlations between snapshots. For example, if we have 3 timescales (10Myr, 100Myr, snapshot) corresponding to subscripts
(1,2,3)
, with current and previous snapshots with subscripts
(c,p)
this results in a 6x6 covariance matrix.
Σ
p,p
Σ
p,c
Σ
c,p
Σ
c,c

ρ
1p1p
ρ
1p2p
ρ
1p3p
ρ
1p1c
ρ
1p2c
ρ
1p3c
ρ
2p1p
ρ
2p2p
ρ
2p3p
ρ
2p1c
ρ
2p2c
ρ
2p3c
ρ
3p1p
ρ
3p2p
ρ
3p3p
ρ
3p1c
ρ
3p2c
ρ
3p3c
ρ
1c1p
ρ
1c2p
ρ
1c3p
ρ
1c1c
ρ
1c2c
ρ
1c3c
ρ
2c1p
ρ
2c2p
ρ
2c3p
ρ
2c1c
ρ
2c2c
ρ
2c3c
ρ
3c1p
ρ
3c2p
ρ
3c3p
ρ
3c1c
ρ
3c2c
ρ
3c3c
​

PSD (Fourier-space) Approach

Getting the auto- and cross-correlation functions from a PSD is straightforward numerically, since the cross-PSD between the SFH smoothed on two separate scales is simply:
P
ij
(k)
H
i
(k)
H
j
(k)P(k)
Now, all we need to do is multiply the PSD by the desired filters in fourier space, inverse Fourier transform to get the cross-correlation function
P
ij
(τ)
, and evaluate at
τ0
and
τ
t
snap
to fill the covariance matrix.​​Note: since the time between the current and previous snapshot can be different from the time between the previous and second previous snapshots, we effectively have four filter radii
(1,2,
3
p
,
3
c
)
P
11
(0)
P
12
(0)
P
1
3
p
(0)
P
11
(
t
snap
)
P
12
(
t
snap
)
P
13c
(
t
snap
)
P
21
(0)
P
22
(0)
P
2
3
p
(0)
P
12
(
t
snap
)
P
22
(
t
snap
)
P
23c
(
t
snap
)
P
3
p
1
(0)
P
3
p
2
(0)
P
3
p

3
p
(0)
P
1
3
p
(
t
snap
)
P
2
3
p
(
t
snap
)
P
3
p

3
c
(
t
snap
)
P
11
(
t
snap
)
P
12
(
t
snap
)
P
1
3
p
(
t
snap
)
P
11
(0)
P
12
(0)
P
1
3
c
(0)
P
12
(
t
snap
)
P
22
(
t
snap
)
P
2
3
p
(
t
snap
)
P
12
(0)
P
22
(0)
P
2
3
c
(0)
P
1
3
c
(
t
snap
)
P
2
3
c
(
t
snap
)
P
3
p

3
c
(
t
snap
)
P
1
3
c
(0)
P
2
3
c
(0)
P
3
c

3
c
(0)

Direct (Real-space) Approach

Hang on a second, If I need to do the inverse Fourier transform to get the correlation function, why not just choose a form for the correlation function directly, and use that? Especially since we only need to evaluate it at two delays, it seems wasteful to build the entire function via the IFFT. The downside here is that the formula for the cross-correlation involves a double integral:
P
ij
(τ)
∞
∫
-∞
du
∞
∫
-∞
dvH(u)H(v)P(τ+u-v)
For top-hat filter functions this becomes
P
ij
(τ)
τ
1
∫
0
du
τ
2
∫
0
dvP(τ+u-v)
If the function P is integrable (and positive?) over the rectangle
[[0,
τ
1
],[0,
τ
2
]]
then we can substitute
sτ+u-v
and swap the integral order (Jacobian == -1)
(Mathematica doesn't want to do this substitution, double check the Fubini Theorem)​​
P
ij
(τ)
τ
1
∫
0
du
τ+u
∫
τ+u-
τ
2
dsP(s)
τ+
τ
1
∫
τ-
τ
2
dsP(s)l(s)
where
l(s)
is the span of
u
values at a given
s
which lay inside the rectangle, given by the piecewise function (for t1 > t2), this is also equal to the trapezoid resulting from the convolution of the two top-hat filters:
l(s)
0,s<τ-
τ
2
||s>τ+
τ
1
s-(τ-
τ
2
),τ-
τ
2
<s<τ
τ
2
,τ<s<τ+
τ
1
-
τ
2
τ+
τ
1
-s,τ+
τ
1
-
τ
2
<s<τ+
τ
1
​​The function will evaluate the same for
t1<t2
in the end, but the general version is messier with a lot of min/max.
This approach will work well when the function
(A+Bs)P(s)
is analytically integrable. using the Autocorrelation function from BC25 (assuming
α2
):
In[]:=
abintpos[A_,B_,s0_,s1_]:=Integrate[(A+B*s)*Exp[-s/t0],{s,s0,s1},Assumptions{t0>0,t1>0,t2>0,τ>0}]
In[]:=
abintneg[A_,B_,s0_,s1_]:=Integrate[(A+B*s)*Exp[s/t0],{s,s0,s1},Assumptions{t0>0,t1>0,t2>0,τ>0}]
In[]:=
abintpos[A,B,s0,s1]
Out[]=
-
s0
t0

t0(A+B(s0+t0))-
-
s1
t0

t0(A+B(s1+t0))
In[]:=
abintneg[A,B,s0,s1]
Out[]=
A(-
s0/t0

+
s1/t0

)t0+Bt0(
s1/t0

(s1-t0)+
s0/t0

(-s0+t0))
Where one final piecewise split may be necessary to account for the modulus in this function (tau is always positive but s can be negative). So the function to calculate the integral will do the following for each segment of
l(s)
:​
​- determine if this interval is within our range of
s
​
​- determine if there there is a negative
s
portion of the interval​
​- evaluate the analytic integrals of
(A+Bs)P(s)
or
(A+Bs)P(-s)
as appropriate​
​- add the result to the total.​
​​
​As an example where all components are required (
τ<
τ
2
). Since we chose
t
1
>
t
2
, s can only be negative in the first segment of
l(s)
:
In[]:=
pwint=abintneg[t2-τ,1,τ-t2,0]+abintpos[t2-τ,1,0,τ]+abintpos[t2,0,τ,τ+t1-t2]+abintpos[t1+τ,-1,τ+t1-t2,τ+t1]
Out[]=
-
t1+τ
t0


t1
t0

-
t2
t0

t0t2+
-
t1+τ
t0

t0t0+
t2
t0

(-t0+t2)+t0-1+
-t2+τ
t0

t0+t2-τ+t0t0+t2-
-
τ
t0

(t0+t2)-τ
In[]:=
Simplify[pwint]
Out[]=
-
t1+t2+τ
t0

t0
t2
t0

t0-
2t2
t0

t0-
t1+t2
t0

t0+
t1+2τ
t0

t0+2
t1+t2+τ
t0

(t2-τ)
Let's verify that the piecewise integral is equal to the double integral:
In[]:=
dbint=Integrate[Exp[-Abs[τ+u-v]/t0],{u,0,t1},{v,0,t2},Assumptions{t0>0,t1>0,t2>0,τ>0,τ<t2,t2<t1}]
Out[]=
-
t1
t0
-
t2
t0
-
τ
t0

t0
t2
t0

t0-
2t2
t0

t0-
t1
t0
+
t2
t0

t0+
t1
t0
+
2τ
t0

t0+2
t1
t0
+
t2
t0
+
τ
t0

t2-2
t1
t0
+
t2
t0
+
τ
t0

τ
In[]:=
Simplify[dbint]
Out[]=
-
t1+t2+τ
t0

t0
t2
t0

t0-
2t2
t0

t0-
t1+t2
t0

t0+
t1+2τ
t0

t0+2
t1+t2+τ
t0

(t2-τ)
In[]:=
Simplify[pwint-dbint,{τ<t2,t2<t1}]
Out[]=
0
Once we find the cross-correlation function for each filter delay, we fill the covariance matrix in the same manner as the PSD method.

Summary

If we can define the autocorrelation such that (A + Bs)P(s) is analytically integrable, We can precompute the covariance matrix without worrying about any numerical error. Otherwise, we should define a PSD and sample it well enough in order to minimise aliasing in construction of the covariance matrix. In my testing I noticed non-neglible differences in the covariance matrices between sampling a 500/1000 Myr range at spacings of 1/2 Myr, however it's also possible my implementation contains bugs at this stage.

Sampling Correlated and Conditional Gaussian Variables

When we have the covariance matrix fully defined, we use the following formula to compute the covariance of each quantity of the current snapshot, *conditioned* on the quantities at the previous snapshots:
Σ
c|p

Σ
cc
-
Σ
cp
-1
Σ
pp
Σ
pc
μ
c|p

m
c
+
Σ
cp
-1
Σ
pp

f
p
-
m
p

Where
Σ
ij
are the blocks of the covariance matrix defined above (2 = current, 1 = previous),
m
i
are the means of the samples at that snapshot, and
f
p
is the actual previous sample that we are conditioning on. We are calculating the stochastic component of the SFH here, so we take
m
p

m
c
0
and calculate the mean adjustment and conditional covariance matrix from the full covariance matrix we computed earlier (we use the GSL BLAS libraries for this). Then we draw samples from this as follows:
f
c

m
c
+
Le
Where
L
T
L

Σ
c|p
and e are uncorrelated samples from the unit gaussian. Since
Σ
c|p
is a 3x3 matrix, this only involves sampling 3 random numbers from the unit gaussian and a 3x3 matrix multiplication per halo if the covariance is only dependent on halo properties via a constant multiplication (e.g. the σ parameter in the PSD above can depend on the properties, but not α or τ0).​​If the covariance timescale depends on the galaxy properties, I can't currently see a way around evaluating the entire covariance matrix per halo, but I can't rule out there existing some clever transformation which can get around it. The full calculation will be ~4 BLAS operations, plus the evaluation of the correlation functions, which is thousands of `exp` calls and an FFT in the PSD case, or dozens of exp calls in the direct case.​​For now I will go on assuming that the halo properties don't change the auto- and cross- correlation timescales of the SFH for simplicity.