Sensitivity quantile distributions

Michele Vallisneri and Chad R. Galley, Jan 17 2012

Jet Propulsion Laboratory, California Insitute of Technology
Copyright 2012. All rights reserved.

This notebook may be used to compute the distribution of squared-SNRs and sensitivities to monochromatic gravitational wave sources distributed either isotropically or in a Galactic disk (for several models). The last section below gives an example of how to run the code. The previous sections are mostly function definitions for the motion of the spacecraft, the single-link response, various TDI observables and noises, and the code for Monte Carlo simulations of squared-SNRs and sensitivities for a given distribution of monochromatic sources. See Vallisneri & Galley (2012) for further exposition. 

The code may be extended to include chirp waveforms of inspiraling binaries, horizon distance computations, other source distributions and polarizations (e.g., Poincare), and alternative spacecraft configurations and orbits.

Preliminaries

Off[General::spell1];

Define some numerical constants.

G=6.67428 10^(-11); (* in m^3/kg/s^2 *)
c=2.99792458 10^8; (* in m/s *)
Msun=1.98892 10^30; (* in kg *)
year=365.25*24*3600; (* in sec *)

lisaL=16.67; (* LISA arm lengths in sec *)

LISA geometry

Define a stationary LISA. n vectors (i=1,...3) are light-propagation directions, p vectors are spacecraft positions. They agree with the conventions used below in "moving LISA geometry".

armn[i_]:={Cos[σ[i]],Sin[σ[i]],0};

armp[i_]:={-Cos[2σ[i]],Sin[2σ[i]],0}/Sqrt[3];

σ[i_]:=3Pi/2-2(i-1)Pi/3;

n[1]=armn[1];n[2]=armn[2];n[3]=armn[3];
n[-1]=-armn[1];n[-2]=-armn[2];n[-3]=-armn[3];

Sending spacecraft

ps[1]=armp[3];ps[2]=armp[1];ps[3]=armp[2];
ps[-1]=armp[2];ps[-2]=armp[3];ps[-3]=armp[1];

Receiving spacecraft

pr[1]=armp[2];pr[2]=armp[3];pr[3]=armp[1];
pr[-1]=armp[3];pr[-2]=armp[1];pr[-3]=armp[2];

Wave

Wave propagation vector (lat and lon are ecliptic longitude and latitude)

myk[lat_,lon_]:=-{Cos[lon]Cos[lat],Sin[lon]Cos[lat],Sin[lat]};

GW rotation matrix

Amat={
{Sin[lon]Cos[psi]-Cos[lon]Sin[lat]Sin[psi],-Sin[lon]Sin[psi]-Cos[lon]Sin[lat]Cos[psi],-Cos[lon]Cos[lat]},
{-Cos[lon]Cos[psi]-Sin[lon]Sin[lat]Sin[psi],Cos[lon]Sin[psi]-Sin[lon]Sin[lat]Cos[psi],-Sin[lon]Cos[lat]},
{Cos[lat]Sin[psi],Cos[lat]Cos[psi],-Sin[lat]}
};

Polarization tensors

hp[β_,λ_,ψ_]:=Amat.{{1,0,0},{0,-1,0},{0,0,0}}.Transpose[Amat]/.{lat->β,lon->λ,psi->ψ};

hc[β_,λ_,ψ_]:=Amat.{{0,1,0},{1,0,0},{0,0,0}}.Transpose[Amat]/.{lat->β,lon->λ,psi->ψ};

Sanity check: transversality...

Simplify[hp[β,λ,ψ].myk[β,λ]]

Moving LISA geometry

Following Królak A, Tinto M and Vallisneri M 2004 Phys. Rev. D 70 022003

Cartwheeling motion and inclination:
eta = Omega t + eta0 is the phase of the guiding center
xi = -Omega t + xi0 is the phase of each spacecraft
zeta = -pi/6 is the inclination angle of plane of orbit to the ecliptic

B={
{Sin[η]Cos[ξ]-Cos[η]Sin[ζ]Sin[ξ],-Sin[η]Sin[ξ]-Cos[η]Sin[ζ]Cos[ξ],-Cos[η]Cos[ζ]},
{-Cos[η]Cos[ξ]-Sin[η]Sin[ζ]Sin[ξ],Cos[η]Sin[ξ]-Sin[η]Sin[ζ]Cos[ξ],-Sin[η]Cos[ζ]},
{Cos[ζ]Sin[ξ],Cos[ζ]Cos[ξ],-Sin[ζ]}
};

ζ=-Pi/6;
movelisa={ξ->ξ0-Ω t,η->η0+Ω t};

The LISA rotation is modeled by applying the matrix B(t) to the n and p.

Looking at the single-link response (below), we see we obtain the same result by applying the transpose of B to k.

(1/2)(n[arm].H.n[arm])/(1-k.n[arm])Exp[I x k.ps[arm]](Exp[I x]-Exp[I x k.n[arm]])

So we compare the next two expressions:

myk[β,λ]

newk=Transpose[B].myk[β0,λ0]

Therefore we find β as

ruleβ=β->ArcSin[-newk[[3]]]

and also λ

cosβ=Cos[β]/.ruleβ;
ruleλ=λ->ArcTan[-newk[[1]]/cosβ,-newk[[2]]/cosβ];

As far as polarizations go, we must have

newAmat=(Transpose[B].(Amat/.{lat->β0,lon->λ0,psi->ψ0}));

which gives the rule for ψ

ruleψ=ψ->ArcTan[newAmat[[3,2]]/cosβ,newAmat[[3,1]]/cosβ];

Let's plot an apparent sky movement over a year...

movesrc3[βinit_,λinit_,ψinit_,Ωlisa_:1,ξinit_:0,ηinit_:0]=(({β,λ,ψ}//.{ruleβ,ruleλ,ruleψ,β0->βinit,λ0->λinit,ψ0->ψinit})/.movelisa)/.{ξ0->ξinit,η0->ηinit,Ω->Ωlisa}//Simplify;

Show[{
ParametricPlot[Reverse[movesrc3[0,0,0][[1;;2]]],{t,0,2Pi},PlotRange->{{-Pi,Pi},{-Pi/2,Pi/2}},Axes->False,Frame->True],
ParametricPlot[Reverse[movesrc3[Pi/3,Pi/6,0][[1;;2]]],{t,0,2Pi},PlotRange->{{-Pi,Pi},{-Pi/2,Pi/2}},PlotStyle->{RGBColor[1,0,0]},Axes->False,Frame->True],
ParametricPlot[Reverse[movesrc3[-Pi/4,3Pi/2,0][[1;;2]]],{t,0,2Pi},PlotRange->{{-Pi,Pi},{-Pi/2,Pi/2}},PlotStyle->{RGBColor[0,1,0]},Axes->False,Frame->True]
}]

Try it out for a random and explicit sky position. That is, compute the response due to the spacecraft motion by evolving the coordinates in time (blue curve) and also by keeping the coordinates fixed but multiplying the response by the appropriate time-dependent rotation matrix (red curve). If everything works out then the two curves should lie on top of each other.

βr=Pi/2-ArcCos[2*Random[]-1];
λr=2*Pi*Random[];
ψr=Pi*Random[];

{βm,λm,ψm}=movesrc3[βr,λr,ψr];

moveB=B/.movelisa/.{ξ0->0,η0->0,Ω->1};

GraphicsGrid[{{
Plot[{n[1].hp[βm,λm,ψm].n[1],(moveB.n[1]).hp[βr,λr,ψr].(moveB.n[1])},{t,0,2Pi}],
Plot[{n[2].hp[βm,λm,ψm].n[2],(moveB.n[2]).hp[βr,λr,ψr].(moveB.n[2])},{t,0,2Pi}],
Plot[{n[3].hp[βm,λm,ψm].n[3],(moveB.n[3]).hp[βr,λr,ψr].(moveB.n[3])},{t,0,2Pi}]
},{
Plot[{n[1].hc[βm,λm,ψm].n[1],(moveB.n[1]).hc[βr,λr,ψr].(moveB.n[1])},{t,0,2Pi}],
Plot[{n[2].hc[βm,λm,ψm].n[2],(moveB.n[2]).hc[βr,λr,ψr].(moveB.n[2])},{t,0,2Pi}],
Plot[{n[3].hc[βm,λm,ψm].n[3],(moveB.n[3]).hc[βr,λr,ψr].(moveB.n[3])},{t,0,2Pi}]
}}]

Single link response

Estabrook–Wahlquist response; note that delays are exp 2π i f L, so the complex factor is consistent with send minus receive: ψ((Subscript[t, send])-k·(Subscript[p, send]))-ψ((Subscript[t, recv])-k.(Subscript[p, recv]))=ψ(t-L-k·(Subscript[p, send]))-ψ((Subscript[t, recv])-k·((Subscript[p, send])+(Subscript[n, link])))=Δ[k·(Subscript[p, send])](Δ[L]-Δ[k·(Subscript[n, link])]). See Eq. (1) of the synthLISA paper (Vallisneri M 2005 Phys. Rev. D 71 022001)

ygw[arm_?NumberQ,H_,k_]:=If[1-k.n[arm]==0,0,(1/2)(n[arm].H.n[arm])/(1-k.n[arm])Exp[I x k.ps[arm]](Exp[I x]-Exp[I x k.n[arm]])];

Build table of frequency abscissae

Logarithmic between lowf and midf, linear between midf and highf. tlen and tlen2 set how many points are in each frequency interval.

tlen=101;
tlen2=251;

lowf=10^-4;
midf=0.01;
highf=0.1;

fs=Join[
Table[10^(ex),{ex,Log[10,lowf],Log[10,midf],(Log[10,midf]-Log[10,lowf])/(tlen/2)}],
Table[f,{f,midf+(highf-midf)/(tlen2/2),highf,(highf-midf)/(tlen2/2)}]
];

xs=N[(2Pi lisaL)fs];

Plot the points. More points are included at higher frequency since the sensitivity exhibits oscillations that need to be resolved.

ListLogLogPlot[Transpose[{fs,fs}],PlotRange->{{lowf,highf},{lowf,highf}}]

Observable and noise definitions

See Sec. 3 and Table 1 in paper

Spm1=2.53654 10^-48*(1.0+(f/10^-4)^-1)*(1.0+(f/0.008)^4)*f^(-2)/.f->x/(2 Pi lisaL);
Sop1=1.42319 10^-37*(1.0+(f/0.002)^-4)*f^2/.f->x/(2 Pi lisaL);

noiserules={Spm->Spm1,Sop->Sop1};

noise[X]=(16 Sop ((Sin[x])^2)+16 Spm (3+Cos[2 x]) ((Sin[x])^2)/.noiserules);

obs[X]=(-l[3]+l[-2] l[2] l[3]) Y[-3]+(1-l[-3] l[3]) Y[-2]+(l[-2]-l[-3] l[-2] l[3]) Y[2]+(-1+l[-2] l[2]) Y[3];
obs[Y]=(1-l[-1] l[1]) Y[-3]+(-l[1]+l[-3] l[1] l[3]) Y[-1]+(-1+l[-3] l[3]) Y[1]+(l[-3]-l[-3] l[-1] l[1]) Y[3];
obs[Z]=(-l[2]+l[-1] l[1] l[2]) Y[-2]+(1-l[-2] l[2]) Y[-1]+(l[-1]-l[-2] l[-1] l[2]) Y[1]+(-1+l[-1] l[1]) Y[2];

See Vallisneri M, Crowder J and Tinto M 2008 Classical and Quantum Gravity 25 065005 for a larger bestiary of TDI observables.

noise[α]=(6 Sop+8 Spm (5+4 Cos[x]+2 Cos[2 x]) ((Sin[(x/2)])^2)/.noiserules);

noise[ζ]=(6 Sop-12 Spm (-1+Cos[x])/.noiserules);

The 5-link eigenvectors are [(Subscript[P, 2])-(Subscript[E, 3])-(Subscript[X, 1])]/(Sqrt[2]) and [(Subscript[P, 2])+(Subscript[E, 3])]/(Sqrt[2]).

noise[U1]=noise[V1]=(8 Sop (4+4 Cos[x]+Cos[2 x]) ((Sin[(x/2)])^2)+16 Spm (5+5 Cos[x]+2 Cos[2 x]) ((Sin[(x/2)])^2)/.noiserules);
noise[E1]=noise[P1]=(16 Spm (3+Cos[x]) ((Sin[(x/2)])^2)+8 Sop (3+2 Cos[x]) ((Sin[(x/2)])^2)/.noiserules);

obs[U1]=(1-l[1] l[2] l[3]) Y[-1]+(l[-1]-l[2] l[3]) Y[1]+(-1+l[-1] l[1]) Y[2]+(-l[2]+l[-1] l[1] l[2]) Y[3];
obs[V1]=(-1+l[-1] l[1]) Y[-3]+(-l[-3]+l[-3] l[-1] l[1]) Y[-2]+(-l[-3] l[-2]+l[1]) Y[-1]+(1-l[-3] l[-2] l[-1]) Y[1];

obs[E1]=(-1+l[-1] l[1]) Y[-2]+(-l[-2]+l[1] l[3]) Y[-1]+(-l[-2] l[-1]+l[3]) Y[1]+(1-l[-1] l[1]) Y[3];
obs[E2]=(-1+l[-2] l[2]) Y[-3]+(-l[-3]+l[1] l[2]) Y[-2]+(1-l[-2] l[2]) Y[1]+(-l[-3] l[-2]+l[1]) Y[2];
obs[E3]=(-l[-1]+l[2] l[3]) Y[-3]+(-1+l[-3] l[3]) Y[-1]+(1-l[-3] l[3]) Y[2]+(-l[-3] l[-1]+l[2]) Y[3];

obs[P1]=(-l[2]+l[-1] l[1] l[2]) Y[-3]+(-l[-3]+l[1] l[2]) Y[-1]+(-l[-3] l[-1]+l[2]) Y[1]+(l[-3]-l[-3] l[-1] l[1]) Y[2];
obs[P2]=(-l[-1]+l[2] l[3]) Y[-2]+(-l[3]+l[-2] l[2] l[3]) Y[-1]+(-l[-2] l[-1]+l[3]) Y[2]+(l[-1]-l[-2] l[-1] l[2]) Y[3];
obs[P3]=(-l[-2]+l[1] l[3]) Y[-3]+(-l[1]+l[-3] l[1] l[3]) Y[-2]+(l[-2]-l[-3] l[-2] l[3]) Y[1]+(-l[-3] l[-2]+l[1]) Y[3];

obs[α]=Y[-2]+l[-2]Y[-1]+l[-1]l[-2]Y[-3]-(Y[3]+l[3]Y[1]+l[1]l[3]Y[2]);
obs[β]=Y[-3]+l[-3]Y[-2]+l[-2]l[-3]Y[-1]-(Y[1]+l[1]Y[2]+l[2]l[1]Y[3]);
obs[γ]=Y[-1]+l[-1]Y[-3]+l[-3]l[-1]Y[-2]-(Y[2]+l[2]Y[3]+l[3]l[2]Y[1]);

obs[ζ]=l[2](Y[1]-Y[-3])+l[3](Y[2]-Y[-1])+l[1](Y[3]-Y[-2]);

The 5-link eigenvectors are [(Subscript[P, 2])-(Subscript[E, 3])-(Subscript[X, 1])]/(Sqrt[2]) and [(Subscript[P, 2])+(Subscript[E, 3])]/(Sqrt[2]).

obs[M]=(obs[P2]-obs[E3]-obs[X])/Sqrt[2];
obs[J]=(obs[P2]+obs[E3])/Sqrt[2];

noise[M]=(8 Sop (((2+Cos[x]))^2) ((Sin[(x/2)])^2)+8 Spm (14+15 Cos[x]+6 Cos[2 x]+Cos[3 x]) ((Sin[(x/2)])^2)/.noiserules);
noise[J]=(32 Spm (2+Cos[x]) ((Sin[(x/2)])^2)+8 Sop (((2+Cos[x]))^2) ((Sin[(x/2)])^2)/.noiserules);

The 6-link eigenvectors are (Z - X)/(Sqrt[2]), (X - 2 Y + Z)/(Sqrt[6]), and (X + Y + Z)/(Sqrt[3]).

obs[A]=(obs[Z]-obs[X])/Sqrt[2];
obs[Eo]=(obs[X]-2obs[Y]+obs[Z])/Sqrt[6];
obs[T]=(obs[X]+obs[Y]+obs[Z])/Sqrt[3];

noise[A]=noise[Eo]=(8 Sop (2+Cos[x]) ((Sin[x])^2)+16 Spm (3+2 Cos[x]+Cos[2 x]) ((Sin[x])^2)/.noiserules);
noise[T]=(-16 Sop (-1+Cos[x]) ((Sin[x])^2)+128 Spm ((Sin[(x/2)])^4) ((Sin[x])^2)/.noiserules);

Functions for computing squared-SNR and sensitivity

Random sample generation

Make a sample of randomly oriented and placed GW binaries.

Isotropic distribution

IsotropicSample[num_,pol_:binary]:=Module[{i,sample},
Switch[pol,
binary,
For[i=1,i≤num,i++,
(* {β (Ecliptic latitude), λ, ψ, ι (Inclination)} *)
sample[i]={π/2-ArcCos[2 RandomReal[]-1],2 π RandomReal[],2 π RandomReal[],ArcCos[2 RandomReal[]-1]};
]
];

Array[sample,num]
];

Galactic distribution

The location of and distance to a source drawn from a Galactic-disk distribution (see Eq. (32) of paper).

(*Working in kpc*)
scaleheight=2.5;
thinheight=0.3;
centerdistance=8;

{deltaNGP,alphaNGP,lascend,eps}={27.1282511,192.8594813,33.0,23.4392911}*Pi/180;

galactic2ecliptic[l_,b_]:=Module[{delta,alpha,beta,lam},delta=ArcSin[Cos[b]*Cos[deltaNGP]*Sin[l-lascend]+Sin[b]*Sin[deltaNGP]];
alpha=ArcTan[Sin[b]*Cos[deltaNGP]-Cos[b]*Sin[l-lascend]*Sin[deltaNGP],Cos[b]*Cos[l-lascend]]+alphaNGP;
beta=ArcSin[Sin[delta]*Cos[eps]-Cos[delta]*Sin[alpha]*Sin[eps]];
lam=ArcTan[Cos[delta]*Cos[alpha],Cos[delta]*Sin[alpha]*Cos[eps]+Sin[delta]*Sin[eps]];
{beta,lam}];

diskbinaryd[sh_: scaleheight,th_: thinheight,dist_: "Sech"]:=Module[{R,z,phi,x,y,r,l,b},
(* This works only in Mathematica v8: R = RandomVariate[ExponentialDistribution[1/sh]]; *)
R=sh (-1-ProductLog[-1,(-1+Random[])/E]);  (* Variate with probability R Exp[-R/scaleheight] *)
phi=RandomVariate[UniformDistribution[{0,2 Pi}]];
z=If[dist=="Sech",RandomVariate[SechDistribution[0,th/2]],Sign[Random[]-1/2] th ArcTanh[Random[]]];
x=R Cos[phi]-centerdistance;
y=R Sin[phi];
r=Sqrt[x^2+y^2+z^2];
b=ArcSin[z/r];
l=ArcTan[x/r,y/r];
Flatten[{galactic2ecliptic[l,b],r}]
];

The location, distance, polarization, and inclination of a source drawn from a Galactic-disk sample. However, the distance is not used in computing the sensitivity.

DiskSample[num_,pol_:binary]:=Module[{i,sample},
Switch[pol,
binary,
For[i=1,i≤num,i++,
(* {β (Ecliptic latitude), λ, distance, ψ, ι (Inclination)} *)
sample[i]={diskbinaryd[],2 π RandomReal[],ArcCos[2 RandomReal[]-1]}//Flatten;
]
];

Array[sample,num]
];

Instrumental and Galactic white dwarf binary noise models

For comparison, simple fit to LISA's instrumental noise.

lisanoise[f_]:=Module[{Sa,Sac,So,ft,T2,L},
L=lisaL*c;

Sa=3 10^-15*Sqrt[1.0+(f/10^-4)^-1]*Sqrt[1.0+(f/0.008)^4];
Sac=Sa*2.0/(2.0*Pi*f)^2;

So=18 10^-12*Sqrt[1+(f/0.002)^-4];

ft=c/(2.0*L);
T2=1.0+(f/(0.41*ft))^2;

(20.0/3.0)*T2*(Sac^2+So^2)/L^2
];

simplesnr[f_]:=5/(Sqrt[16.0/5.0]*Sqrt[year]/Sqrt[lisanoise[f]]);

Partially subtracted confusion noise from Galactic white dwarf binaries.

ConfusionNoise[model_String]:=Module[{logf,p1,p2,p3,p4,p5,q1,q2,c2,α,β,dNdf,SGWDB},
logf=Log[10,fs];

Switch[model,
"NVNP",(* From Nissanke, Vallisneri, Nelemans, and Prince (2011). *)
α=2;β=2/year;
dNdf=5 10^(-3) fs^(-11/3);
SGWDB=1.4 10^(-45) fs^(-8/3);
3/20. 16 xs^2Sin[xs]^2/(1+0.6 xs^2)SGWDB Tanh[1/2 β dNdf]^α,

(* feel free to roll your own... *)

"None",
0
]
];

Check the instrumental noise (noise[X]), the confusion noise, and their sum to make sure it looks right for the given model.

ListLogLogPlot[{{fs,ConfusionNoise["NVNP"]}//Transpose,{fs,noise[X]/.x->xs}//Transpose,{fs,ConfusionNoise["NVNP"]+(noise[X]/.x->xs)}//Transpose},Joined->True,PlotRange->{{10^-4,2 10^-2},{5 10^-45,10^-39}},PlotStyle->{Blue,Red,Black},Frame->True,FrameLabel->{"f (Hz)","Noise (1/Hz)"}]

Modulation of the above confusion noise as LISA orbits the Sun.

Modulation[pos_,model_String]:=Module[{a1,b1,a2,b2,ft},
Switch[model,
"NVNP",(* pos in months *)
a1=-0.21;b1=-0.05;
a2=-0.45;b2=0.01;
1+a1 Cos[π pos/6]+b1 Sin[π pos/6]+a2 Cos[2π pos/6]+b2 Sin[2π pos/6],

"None",
0
]
];

Plot the modulation of the confusion noise to make sure it looks right.

ListPlot[Table[{i,Modulation[i,"NVNP"]},{i,0,12,0.1}],Joined->True,Mesh->All,PlotStyle->Red]

Response functions

Response function for the stationary spacecraft configuration.

ResponseFn[stationary,obs_,onesample_List,pol_:binary]:=
Module[{obslist,ranβ,ranλ,rand,ranψ,ranι,twoϵ,twoτ,ranΓ,ranϕ,Ap,Ax,phasor,phases},

Switch[pol,
binary,
If[Length[onesample]==4,{ranβ,ranλ,ranψ,ranι}=onesample];  (* For member of an isotropic population *)
If[Length[onesample]==5,{ranβ,ranλ,rand,ranψ,ranι}=onesample];  (* For member of an Galactic disk population, onesample[[3]] is distance to Galactic source. *)
Ap=(1+Cos[ranι]^2)/2; (* According to MLDC definitions, with 1/2 for reality constraint... *)
Ax=I(-2Cos[ranι])/2, (* ...which may however need further justification than given above. *)

(* in case you want a random distribution on the Poincaré sphere *)
poincare,
twoϵ=ArcSin[2 RandomReal[]-1];
twoτ=2 π RandomReal[];
ranΓ=ArcCos[Cos[twoτ] Cos[twoϵ]]/2;
ranϕ=ArcTan[Sin[twoτ],Tan[twoϵ]];
Ap=Sin[ranΓ]/2; 
Ax=Exp[I ranϕ]Cos[ranΓ]/2
];

phasor=obs/.{l[link_]:>Exp[I x],Y[link_]:>ygw[link,Ap hp[ranβ,ranλ,ranψ]+Ax hc[ranβ,ranλ,ranψ],myk[ranβ,ranλ]]};
phases=phasor/.x->xs;

4Abs[phases]^2//Re
];

Response function for the moving spacecraft configuration. The default for measuring the position of the spacecraft throughout one year is 26 times, or once every two weeks.

ResponseFn[moving,obs_,onesample_List,currentpos_,positions_:26]:=
Module[{ranβ,ranλ,rand,ranψ,ranι,Ap,Ax,movβ,movλ,movψ,phasor,phases,res},

If[Length[onesample]==4,{ranβ,ranλ,ranψ,ranι}=onesample];  (* For member of an isotropic population *)
If[Length[onesample]==5,{ranβ,ranλ,rand,ranψ,ranι}=onesample];  (* For member of an Galactic disk population. onesample[[3]] is distance to Galactic source. *)
Ap=(1+((Cos[ranι])^2))/2; (* According to MLDC definitions, with 1/2 for reality constraint... *)
Ax=I(-2Cos[ranι])/2;  (* ...which may however need further justification than given above. *)

{movβ,movλ,movψ}=movesrc3[ranβ,ranλ,ranψ]/.t->2Pi currentpos/positions;

phasor=obs/.{l[link_]->Exp[I x],Y[link_]->ygw[link,Ap hp[movβ,movλ,movψ]+Ax hc[movβ,movλ,movψ],myk[movβ,movλ]]};
phases=phasor/.x->xs;

4 Abs[phases]^2//Re
];

Squared (optimal) signal-to-noise ratios (SNRs)

Optimal squared-SNRs for a single monochromatic source and for either the stationary (default) or moving spacecraft configurations.

OneBinarySNR2[monochromatic,observable_,noise_,onesample_List,GWDBmodel_String,response_:stationary,positions_:26]:=
Module[{instnoises,obslist,noiselist,weight,snr2},

obslist=If[Head[observable]=!=List,{observable},observable];
noiselist=If[Head[noise]=!=List,{noise},noise];

Switch[response,
stationary,(* Average the SNR^2, sum over observables... *)
snr2= year  ResponseFn[stationary,observable,onesample]/(ConfusionNoise[GWDBmodel]+(noise/.x->xs)),

moving,(* Average the SNR^2 with moving response... Note pol can only be binary... *)
snr2= year/positions Sum[ResponseFn[moving,observable,onesample,pp,positions]/(ConfusionNoise[GWDBmodel] Modulation[12pp/positions,GWDBmodel]+(noise/.x->xs)),{pp,0,positions-1}]
];

snr2
];

Optimal squared-SNRs for a given number of monochromatic sources and for either the stationary (default) or moving spacecraft configurations.

ManyBinarySNR2[monochromatic,observable_,noise_,iter_,samples_List,GWDBmodel_String,response_:stationary,positions_:26]:=
Module[{instnoises,obslist,noiselist,normweight,snr2,i},

obslist=If[Head[observable]=!=List,{observable},observable];
noiselist=If[Head[noise]=!=List,{noise},noise];

Switch[response,
stationary,(* Average the SNR^2, sum over observables... *)
snr2= Table[OneBinarySNR2[monochromatic,observable,noise,samples[[i]],GWDBmodel,stationary,positions],{i,1,iter}],

moving,(* Average the SNR^2 with moving response... Note pol can only be binary... *)
snr2= Table[OneBinarySNR2[monochromatic,observable,noise,samples[[i]],GWDBmodel,moving,positions],{i,1,iter}]
];

(* Rows correspond to frequencies *)
(* Columns correspond to the SNR^2 for each sampled GW source. *)
snr2//Transpose
];

Sensitivity

Sensitivity[snr2_,thresholdSNR_]:=thresholdSNR/Sqrt[#]&/@snr2

Running the code:  An example for stationary LISA and a Galactic population of monochromatic sources

Draw 1000 samples from a population of monochromatic sources from the Galactic-disk. (We use 10^5 sources in the paper so, here, the statistics will not be as good.)

sample=DiskSample[1000];

For the sensitivity of a single, monochromatic source (say, the first one of the above 1000 sources) one would execute the following:

oneSNR2=OneBinarySNR2[monochromatic,obs[X],noise[X],sample[[1]],"NVNP",stationary]

More specifically, this would be the optimal squared-SNR for the observable X, including the partially subtracted Galactic WD noise, for the case when LISA is moving. The corresponding sensitivity plotted as a function of the frequency (fs, defined in section `Building table of frequency abscissae’) relative to a threshold squared-SNR of 5 is then

oneSensitivity=Sensitivity[oneSNR2,5];

ListLogLogPlot[{fs,oneSensitivity}//Transpose,Joined->True]

Computing the squared-SNRs of all the 1000 samples generated above can be done by executing the following [11 seconds on our machine]:

manySNR2=ManyBinarySNR2[monochromatic,obs[X],noise[X],1000,sample,"NVNP",stationary]//Transpose;

The corresponding sensitivities are, for a threshold SNR of 5 (not squared-SNR):

manySensitivities=Sensitivity[manySNR2,5];

Plot the sensitivities of the first 3 sources:

GraphicsGrid[{
{ListLogLogPlot[{fs,manySensitivities[[1]]}//Transpose,Joined->True],ListLogLogPlot[{fs,manySensitivities[[2]]}//Transpose,Joined->True],ListLogLogPlot[{fs,manySensitivities[[3]]}//Transpose,Joined->True]}},ImageSize->900]

There is a lot of information contained in manySNR2 but a convenient representation for all 1000 sources is via quantiles. For example, let’s plot the 5% and 95% quantiles (upper and lower red curves, respectively) along with the median (dashed black) and mean (solid black) curves using the inverse-RMS, as discussed in the paper. The plot below thus shows the spread in sensitivity of 90% of the sampled sources.

ListLogLogPlot[{
Transpose[{fs,5/Sqrt[Quantile[manySNR2,0.05]]}],
Transpose[{fs,5/Sqrt[Quantile[manySNR2,0.95]]}],
Transpose[{fs,5/Sqrt[Quantile[manySNR2,0.50]]}],
Transpose[{fs,5/Sqrt[Mean[manySNR2]]}]
},Joined->True,PlotStyle->{Red,Red,{Black,Dashed},Black},PlotRange->{{10^-4,10^-1},Automatic},Frame->True,FrameLabel->{"Frequency [Hz]","Sensitivity"}]

Let's redo this for a moving detector [less than 5 mins on our machine]:

manySNR2Moving=ManyBinarySNR2[monochromatic,obs[X],noise[X],1000,sample,"NVNP",moving]//Transpose;

manySensitivitiesMoving=Sensitivity[manySNR2Moving,5];

ListLogLogPlot[{
Transpose[{fs,5/Sqrt[Quantile[manySNR2Moving,0.05]]}],
Transpose[{fs,5/Sqrt[Quantile[manySNR2Moving,0.95]]}],
Transpose[{fs,5/Sqrt[Quantile[manySNR2Moving,0.50]]}],
Transpose[{fs,5/Sqrt[Mean[manySNR2Moving]]}]
},Joined->True,PlotStyle->{Red,Red,{Black,Dashed},Black},PlotRange->{{10^-4,10^-1},Automatic},Frame->True,FrameLabel->{"Frequency [Hz]","Sensitivity"}]