%run proof_setup
First, the real case is relatively simple. The integral that we want to do is:
$$ k_\Delta(\tau) = \frac{1}{\Delta^2}\int_{t_i-\Delta/2}^{t_i+\Delta/2} \mathrm{d}t \,\int_{t_j-\Delta/2}^{t_j+\Delta/2}\mathrm{d}t^\prime\,k(t - t^\prime) $$For celerite kernels it helps to make the assumtion that $t_j + \Delta/2 < t_i - \Delta/2$ (in other words, the exposure times do not overlap).
import sympy as sm
cr = sm.symbols("cr", positive=True)
ti, tj, dt, t, tp = sm.symbols("ti, tj, dt, t, tp", real=True)
k = sm.exp(-cr * (t - tp))
k0 = k.subs([(t, ti), (tp, tj)])
kint = sm.simplify(
sm.integrate(
sm.integrate(k, (t, ti - dt / 2, ti + dt / 2)), (tp, tj - dt / 2, tj + dt / 2)
)
/ dt ** 2
)
res = sm.simplify(kint / k0)
print(res)
(exp(2*cr*dt) - 2*exp(cr*dt) + 1)*exp(-cr*dt)/(cr**2*dt**2)
This is the factor that we want. Let's make sure that it is identical to what we have in the note.
kD = 2 * (sm.cosh(cr * dt) - 1) / (cr * dt) ** 2
sm.simplify(res.expand() - kD.expand())
0
Excellent.
Let's double check that this reduces to the original kernel in the limit $\Delta \to 0$:
sm.limit(kD, dt, 0)
1
The complex cases proceeds similarly, but it's a bit more involved. In this case,
$$ k(\tau) = (a + i\,b)\,\exp(-(c+i\,d)\,(t_i-t_j)) $$a, b, c, d = sm.symbols("a, b, c, d", real=True, positive=True)
k = sm.exp(-(c + sm.I * d) * (t - tp))
k0 = k.subs([(t, ti), (tp, tj)])
kint = sm.simplify(sm.integrate(k, (t, ti - dt / 2, ti + dt / 2)) / dt)
kint = sm.simplify(sm.integrate(kint.expand(), (tp, tj - dt / 2, tj + dt / 2)) / dt)
print(sm.simplify(kint / k0))
(2*cos(I*dt*(c + I*d)) - 2)/(dt**2*(c**2 + 2*I*c*d - d**2))
That doesn't look so bad!
But, I'm going to re-write it by hand and make sure that it's correct:
coeff = (c - sm.I * d) ** 2 / (dt * (c ** 2 + d ** 2)) ** 2
coeff *= sm.exp((c + sm.I * d) * dt) + sm.exp(-(c + sm.I * d) * dt) - 2
sm.simplify(coeff * k0 - kint)
0
Good.
Now we need to work out nice expressions for the real and imaginary parts of this. First, the real part. I found that it was easiest to look at the prefactors for the trig functions directly and simplify those. Here we go:
res = (a + sm.I * b) * coeff
A = sm.simplify((res.expand(complex=True) + sm.conjugate(res).expand(complex=True)) / 2)
sm.simplify(sm.poly(A, sm.cos(dt * d)).coeff_monomial(sm.cos(dt * d)))
(a*c**2*exp(2*c*dt) + a*c**2 - a*d**2*exp(2*c*dt) - a*d**2 + 2*b*c*d*exp(2*c*dt) + 2*b*c*d)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
sm.simplify(
sm.poly(
sm.poly(A, sm.cos(dt * d)).coeff_monomial(1), sm.sin(dt * d)
).coeff_monomial(sm.sin(dt * d))
)
(2*a*c*d*exp(2*c*dt) - 2*a*c*d - b*c**2*exp(2*c*dt) + b*c**2 + b*d**2*exp(2*c*dt) - b*d**2)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
sm.simplify(
sm.poly(
sm.poly(A, sm.cos(dt * d)).coeff_monomial(1), sm.sin(dt * d)
).coeff_monomial(1)
)
2*(-a*c**2 + a*d**2 - 2*b*c*d)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
Then, same thing for the imaginary part:
B = sm.simplify(
-sm.I * (res.expand(complex=True) - sm.conjugate(res).expand(complex=True)) / 2
)
sm.simplify(sm.poly(B, sm.cos(dt * d)).coeff_monomial(sm.cos(dt * d)))
(-2*a*c*d*exp(2*c*dt) - 2*a*c*d + b*c**2*exp(2*c*dt) + b*c**2 - b*d**2*exp(2*c*dt) - b*d**2)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
sm.simplify(
sm.poly(
sm.poly(B, sm.cos(dt * d)).coeff_monomial(1), sm.sin(dt * d)
).coeff_monomial(sm.sin(dt * d))
)
(a*c**2*exp(2*c*dt) - a*c**2 - a*d**2*exp(2*c*dt) + a*d**2 + 2*b*c*d*exp(2*c*dt) - 2*b*c*d)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
sm.simplify(
sm.poly(
sm.poly(B, sm.cos(dt * d)).coeff_monomial(1), sm.sin(dt * d)
).coeff_monomial(1)
)
2*(2*a*c*d - b*c**2 + b*d**2)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
Ok.
Now let's make sure that the simplified expressions are right.
C1 = a * c ** 2 - a * d ** 2 + 2 * b * c * d
C2 = b * c ** 2 - b * d ** 2 - 2 * a * c * d
cos_term = (sm.exp(c * dt) + sm.exp(-c * dt)) * sm.cos(d * dt) - 2
sin_term = (sm.exp(c * dt) - sm.exp(-c * dt)) * sm.sin(d * dt)
denom = dt ** 2 * (c ** 2 + d ** 2) ** 2
A0 = (C1 * cos_term - C2 * sin_term) / denom
B0 = (C2 * cos_term + C1 * sin_term) / denom
sm.simplify(A.expand() - A0.expand())
0
sm.simplify(B.expand() - B0.expand())
0
Finally let's rewrite things in terms of hyperbolic trig functions.
sm.simplify(2 * (sm.cosh(c * dt) * sm.cos(d * dt) - 1).expand() - cos_term.expand())
0
sm.simplify(2 * (sm.sinh(c * dt) * sm.sin(d * dt)).expand() - sin_term.expand())
0
Looks good!
Let's make sure that this actually reproduces the target integral:
sm.simplify(((a + sm.I * b) * kint / k0 - (A + sm.I * B)).expand(complex=True))
0
Finally, let's make sure that this reduces to the original kernel when $\Delta \to 0$:
sm.limit(A, dt, 0), sm.limit(B, dt, 0)
(a, b)
sm.simplify(A.subs([(d, 0), (b, 0)]))
a*((exp(c*dt) - 2)*exp(c*dt) + 1)*exp(-c*dt)/(c**2*dt**2)
sm.simplify(B.subs([(d, 0), (b, 0)]))
0
If we directly evaluate the power spectrum of this kernel, we'll have some issues because there will be power from lags where our assumption of non-overlapping exposures will break down. Instead, we can evaluate the correct power spectrum by realizing that the integrals that we're doing are convolutions. Therefore, the power spectrum of the integrated kernel will be product of the original power spectrum with the square of the Fourier transform of the top hat exposure function.
omega = sm.symbols("omega", real=True)
sm.simplify(sm.integrate(sm.exp(sm.I * t * omega) / dt, (t, -dt / 2, dt / 2)))
Piecewise((2*sin(dt*omega/2)/(dt*omega), Ne(dt*omega, 0)), (1, True))
Therefore, the integrated power spectrum is
$$ S_\Delta(\omega) = \frac{\sin^2(\Delta\,\omega/2)}{(\Delta\,\omega/2)^2}\,S(\omega) = \mathrm{sinc}^2(\Delta\,\omega/2)\,S(\omega) $$For overlapping exposures, some care must be taken when computing the autocorrelation because of the absolute value. This also means that celerite cannot be used (as far as I can tell) to evaluate exposure time integrated models with overlapping exposures. In this case, the integral we want to do is:
$$ k_\Delta(\tau) = \frac{1}{\Delta^2}\int_{t_i-\Delta/2}^{t_i+\Delta/2} \mathrm{d}t \,\int_{t_j-\Delta/2}^{t_j+\Delta/2}\mathrm{d}t^\prime\,k(|t - t^\prime|) $$which can be broken into three integrals when $\tau = |t_i - t_j| \le \Delta$ (assuming still that $t_i \ge t_j$):
$$ \Delta^2\,k_\Delta(\tau) = \int_{t_j+\Delta/2}^{t_i+\Delta/2} \mathrm{d}t \,\int_{t_j-\Delta/2}^{t_j+\Delta/2}\mathrm{d}t^\prime\,k(t - t^\prime) + \int_{t_i-\Delta/2}^{t_j+\Delta/2} \mathrm{d}t \,\int_{t_j-\Delta/2}^{t}\mathrm{d}t^\prime\,k(t - t^\prime) + \int_{t_i-\Delta/2}^{t_j+\Delta/2} \mathrm{d}t \,\int_{t}^{t_j+\Delta/2}\mathrm{d}t^\prime\,k(t^\prime - t) $$tau = sm.symbols("tau", real=True, positive=True)
kp = sm.exp(-cr * (t - tp))
km = sm.exp(-cr * (tp - t))
k1 = sm.simplify(
sm.integrate(
sm.integrate(kp, (tp, tj - dt / 2, tj + dt / 2)), (t, tj + dt / 2, ti + dt / 2)
)
/ dt ** 2
)
k2 = sm.simplify(
sm.integrate(sm.integrate(kp, (tp, tj - dt / 2, t)), (t, ti - dt / 2, tj + dt / 2))
/ dt ** 2
)
k3 = sm.simplify(
sm.integrate(sm.integrate(km, (tp, t, tj + dt / 2)), (t, ti - dt / 2, tj + dt / 2))
/ dt ** 2
)
kD = sm.simplify((k1 + k2 + k3).expand())
res = sm.simplify(kD.subs([(ti, tau + tj)]))
res
(2*cr*dt*exp(cr*tau) - 2*cr*tau*exp(cr*tau) + exp(cr*(-dt + 2*tau)) - 2 + exp(-cr*dt))*exp(-cr*tau)/(cr**2*dt**2)
kint = (
2 * cr * (dt - tau)
+ sm.exp(-cr * (dt - tau))
- 2 * sm.exp(-cr * tau)
+ sm.exp(-cr * (dt + tau))
) / (cr * dt) ** 2
sm.simplify(kint - res)
0
Ok. That's the result for the real case. Now let's work through the result for the complex case.
arg1 = ((a + sm.I * b) * kint.subs([(cr, c + sm.I * d)])).expand(complex=True)
arg2 = ((a - sm.I * b) * kint.subs([(cr, c - sm.I * d)])).expand(complex=True)
res = sm.simplify((arg1 + arg2) / 2)
res
(2*dt*(a*c**3 + a*c*d**2 + b*c**2*d + b*d**3)*exp(c*(2*dt + 2*tau)) - 2*tau*(a*c**3 + a*c*d**2 + b*c**2*d + b*d**3)*exp(c*(2*dt + 2*tau)) + 2*(-a*c**2*cos(d*tau) + 2*a*c*d*sin(d*tau) + a*d**2*cos(d*tau) - b*c**2*sin(d*tau) - 2*b*c*d*cos(d*tau) + b*d**2*sin(d*tau))*exp(c*(2*dt + tau)) + (a*c**2*cos(d*(dt - tau)) - 2*a*c*d*sin(d*(dt - tau)) - a*d**2*cos(d*(dt - tau)) + b*c**2*sin(d*(dt - tau)) + 2*b*c*d*cos(d*(dt - tau)) - b*d**2*sin(d*(dt - tau)))*exp(c*(dt + 3*tau)) + (a*c**2*cos(d*(dt + tau)) - 2*a*c*d*sin(d*(dt + tau)) - a*d**2*cos(d*(dt + tau)) + b*c**2*sin(d*(dt + tau)) + 2*b*c*d*cos(d*(dt + tau)) - b*d**2*sin(d*(dt + tau)))*exp(c*(dt + tau)))*exp(-2*c*(dt + tau))/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
C1 = a * c ** 2 - a * d ** 2 + 2 * b * c * d
C2 = b * c ** 2 - b * d ** 2 - 2 * a * c * d
denom = dt ** 2 * (c ** 2 + d ** 2) ** 2
dpt = dt + tau
dmt = dt - tau
cos_term = (
sm.exp(-c * dmt) * sm.cos(d * dmt)
+ sm.exp(-c * dpt) * sm.cos(d * dpt)
- 2 * sm.exp(-c * tau) * sm.cos(d * tau)
)
sin_term = (
sm.exp(-c * dmt) * sm.sin(d * dmt)
+ sm.exp(-c * dpt) * sm.sin(d * dpt)
- 2 * sm.exp(-c * tau) * sm.sin(d * tau)
)
ktest = 2 * (a * c + b * d) * (c ** 2 + d ** 2) * dmt
ktest += C1 * cos_term + C2 * sin_term
ktest /= denom
sm.simplify(ktest - res)
0
This corresponds to adding a function to the kernel:
delta = sm.simplify(
ktest.expand(trig=True)
- (A * sm.exp(-c * tau) * sm.cos(d * tau) + B * sm.exp(-c * tau) * sm.sin(d * tau))
)
C1 = a * c ** 2 - a * d ** 2 + 2 * b * c * d
C2 = b * c ** 2 - b * d ** 2 - 2 * a * c * d
sinh = (sm.exp(c * (dt - tau)) - sm.exp(-c * (dt - tau))) / 2
cosh = (sm.exp(c * (dt - tau)) + sm.exp(-c * (dt - tau))) / 2
norm = dt ** 2 * (c ** 2 + d ** 2) ** 2
delta_test = (
2
* (
C2 * cosh * sm.sin(d * (dt - tau))
- C1 * sinh * sm.cos(d * (dt - tau))
+ (a * c + b * d) * (dt - tau) * (c ** 2 + d ** 2)
)
/ norm
)
sm.simplify(delta.expand(trig=True) - delta_test.expand(trig=True))
0
And now, finally, I think we're done.
By default, the diagonal terms in a celerite model are $\sum a_j$. In this case, that would be $A$, but that would be wrong because $\tau = 0 < \Delta$. Let's work out what extra term should be added to the diagonal:
delta = sm.simplify(ktest.subs([(tau, 0)]) - A)
sm.simplify(sm.poly(delta, sm.cos(dt * d)).coeff_monomial(sm.cos(dt * d)))
(-a*c**2*exp(2*c*dt) + a*c**2 + a*d**2*exp(2*c*dt) - a*d**2 - 2*b*c*d*exp(2*c*dt) + 2*b*c*d)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
sm.simplify(
sm.poly(
sm.poly(delta, sm.cos(dt * d)).coeff_monomial(1), sm.sin(dt * d)
).coeff_monomial(sm.sin(dt * d))
)
(-2*a*c*d*exp(2*c*dt) - 2*a*c*d + b*c**2*exp(2*c*dt) + b*c**2 - b*d**2*exp(2*c*dt) - b*d**2)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))
sm.simplify(
sm.poly(
sm.poly(delta, sm.cos(dt * d)).coeff_monomial(1), sm.sin(dt * d)
).coeff_monomial(1)
)
2*(a*c + b*d)/(dt*(c**2 + d**2))
Let's check that our simplified expression is correct:
C1 = a * c ** 2 - a * d ** 2 + 2 * b * c * d
C2 = b * c ** 2 - b * d ** 2 - 2 * a * c * d
norm = dt ** 2 * (c ** 2 + d ** 2) ** 2
sinh = (sm.exp(c * dt) - sm.exp(-c * dt)) / 2
cosh = (sm.exp(c * dt) + sm.exp(-c * dt)) / 2
delta_test = (
2
* (
C2 * cosh * sm.sin(d * dt)
- C1 * sinh * sm.cos(d * dt)
+ (a * c + b * d) * dt * (c ** 2 + d ** 2)
)
/ norm
)
sm.simplify(
delta.expand(trig=True, complex=True) - delta_test.expand(trig=True, complex=True)
)
0
And check the result for the real term:
sm.simplify(delta_test.subs([(d, 0), (b, 0)]) - 2 * a * (c * dt - sinh) / (dt * c) ** 2)
0
Now let's check to make sure that the diagonal will always be positive:
diag_test = (
2
* (
(a * c + b * d) * (c ** 2 + d ** 2) * dt
- C1
+ sm.exp(-c * dt) * (C1 * sm.cos(d * dt) + C2 * sm.sin(d * dt))
)
/ (dt * (c ** 2 + d ** 2)) ** 2
)
sm.simplify(ktest.subs([(tau, 0)]) - diag_test)
0
sm.limit(diag_test, dt, 0)
a
sm.limit(diag_test, dt, sm.oo)
0