text:
# Effects of Pore Walls and Randomness on Phase Transitions in Porous Media

## I Introduction

Critical phenomena are generally well understood but the effects of randomness on the nature of the transitions are less well-studied. This is especially true in the case of phase transitions that take place in porous media where the effects of quenched randomness are provided by the pore walls. Among the best studied, are phase transitions in highly porous aerogels  — both the liquid - vapor transition  and the $`\lambda `$ \- transition of 4He  have been found to be remarkably sharp. Even more interestingly, the topology of the phase diagram for 3He-4He mixtures in aerogels has been found to be different from that in the bulk .

The simplest theoretical framework for studies of critical phenomena in non-random systems is the Ising model. The important role played by multiple length scales at a critical temperature $`T_c`$ leads to universality  — binary alloys which are about to order, binary liquids which are about to phase separate, certain kinds of magnets with uniaxial anisotropy which are about to become magnetized all exhibit the same critical behavior as the Ising model. Perhaps the simplest extension that incorporates randomness is the random field Ising model (RFIM) wherein a quenched random field is applied at each site . One example of the probability distribution is the symmetric bimodal distribution which corresponds to a situation in which half of the sites experience an up field and the other half, a down field of equal strength. Another example involves fields which are Gaussian-distributed. In both cases, the symmetry between the up and down directions is not broken by hand and thus provides a scope for a spontaneous symmetry breaking and a phase transition associated with it.

Recent research  has lead to the result that the two probability distributions may correspond, at least in dimensions 3 and larger, to distinct behaviors associated with the RFIM. The Gaussian case is governed by a $`T`$=0 fixed point while the bimodal model’s phase diagram is qualitatively different. The origins of the two distinct scenarios relate to their quite different $`T`$=0 phase diagrams.

Experimental realizations of the RFIM include dilute antiferromagnets in a uniform field  and binary liquid mixtures in a porous medium . In both these cases, many of the expected signatures associated with the RFIM with a Gaussian distribution of random fields were observed  (but significant deviations were also found when the disorder was correlated ). However, the sluggish dynamics and irreversibility predicted by the theory  precluded accurate measurements of the exponents for binary liquids in porous media. Nevertheless, the exponents determined for the dilute Ising antiferromagnet in a uniform field were in accord with the theory.

A major surprise in this field were the measurements by Chan and his collaborators  on the liquid-vapor transition of helium and hydrogen in a variety of porous media. While the liquid-vapor coexistence region was considerably shrunk compared to the bulk uniform case, the exponents were found to be much more akin to those of the uniform Ising model instead of the RFIM. It has been suggested  that these features are related to the properties of the RFIM with the fields being distributed bimodally (though not necessarily symmetrically).

Another surprise was the novel shape of the experimentally determined  phase diagram of the 3He-4He mixtures in porous media which allowed for superfluidity at large concentrations of 3He. The classical aspects of this phase separation are captured by the Blume-Emery-Griffiths (BEG)  spin-1 model with an anisotropy. The presence of a porous medium can be modeled by making this anisotropy random with a bimodal distribution.

In this paper, we focus on the role played by pore walls in liquid-vapor transitions in porous media, as studied in their corresponding spin $`\frac{1}{2}`$ Ising spin systems. There are two aspects to the role of walls in a porous medium. First, there is a preference for one of the phases over the other in the vicinity of the walls. This mechanism alone ought to lead to observable consequences even when the placement of the walls is substantially periodic, i.e. the different phases are connected but there is no inherent randomness. Second, the random placement of the walls in the porous medium provides quenched disorder and can induce further changes in the phase behavior. The principal result of our paper is that the former aspect is more crucial – indeed, we show that within the framework of simple models, the phase diagram does not change on incorporating randomness. This finding is consistent with the analysis by Galam and Aharony  indicating that the mean field results of a ferromagnet in a random longitudinal field are the same as a uniaxially anistropic antiferromagnet in a uniform field. Our results suggest that liquid-vapor transitions in designed porous media, with a periodic geometrical pattern , ought to exhibit behavior quite akin to that observed in random porous media. We demonstrate these findings in simple mean field Ising models and two distinct values of the local magnetic fields. We then generalize these studies to spin 1 Ising models with non-uniform anisotropy and show that such systems behave like the random anisotropy BEG systems with a bimodal distribution of the anisotropies . Our results suggest that the 3He-4He phase separation is also primarily governed by the mere presence of walls in the porous medium and not randomness.

It should be noted that there have been several recent mean field studies of phase transitions in random porous media . One would expect that fluctuations in lower dimensions  could play a key role in qualitatively modifying the mean field picture.

## II Effects of Confinement Random in the Field Ising Model – the Symmetric Case

We start by considering the simplest case – that of four Ising spins located at two kinds of sites, 0 and 1, as shown in Figure 1a. Periodic boundary conditions are adopted in the plane of the figure. Furthermore, it is assumed that above and below this plane, there are spins which sit in locations that repeat the pattern shown and allow for a connected string of nearest neighbor 0 and 1 sites in the direction perpendicular to the plane of the paper. Physically, this geometry corresponds to a periodic arrangements of one-dimensional strings of 0 and 1 arranged on two sublattices. All spins are coupled by a uniform exchange constant, $`J`$. The magnetic field on sites 0 is denoted by $`h_0`$ and it points up. On sites 1, on the other hand, the magnetic field is equal to $`h_1`$ and it points down. (The case of a simple ferromagnet with a staggered field is obtained when periodic boundary conditions are adopted in all directions. One then does not have connected strings of 0 or 1 sites which lead to an inability to sustain certain phases at non-zero temperatures.)

Our objective here is to determine the phase diagram of this non-random system within the mean field approximation and compare it to the corresponding mean field results  of the RFIM in which the probability distribution of the magnetic fields is bimodal: half of the randomly selected sites have an up-pointing field $`h_0`$ and the other half – a down-pointing field $`h_1`$. The RFIM may be thought of as modeling porous media with the sites with field $`h_0`$ corresponding to locations near pore walls and the sites with field $`h_1`$ to the interior locations. The underlying assumption here is that there is a different environment near the pore wall than in the interior.

The phase diagram is obtained in the three-dimensional space of $`h_0`$, $`h_1`$, and $`T`$ and is determined by solving the following equations for the magnetizations $`m_0`$ and $`m_1`$:

$$m_0=tanh[(h_0+2Jm_0+4Jm_1)/k_BT],$$
(1)
$$m_1=tanh[(h_1+4Jm_0+2Jm_1)/k_BT],$$
(2)
on sites with field $`h_0`$ and $`h_1`$ respectively. The solution is obtained in an iterative manner that leads to self-consistency. The form of equations (1) and (2) reflects the fact that each site of a given kind has four neighbors of the other kind and two neighbors of the same kind – the latter resulting from the out-of plane connectivity. Once the solutions for the local magnetizations are found, one can determine the free energy, $`F=UTS`$, by calculating the internal energy

$$U=2J(m_0m_0+m_1m_1+4m_0m_1)2h_0m_0+2h_1m_1$$
(3)
and the entropy

$$S=2s_0+2s_1,$$
(4)
where

$$s_i=k_B\frac{1}{2}[(1+m_i)ln(1+m_i)/2+(1m_i)ln(1m_i)/2].$$
(5)
A first order phase transition is identified by the presence of a cusp in the free energy.

It is easy to show that there are three possible phases at $`T`$=0 in this system. We shall denote them by +, $``$, and $`+`$ and their energies by $`E_+`$, $`E_{}`$, and $`E_+`$ respectively. In the first phase, all spins are up and in the second all spins are down. In the third phase, on the other hand, the spins point in the directions of the local magnetic fields.

At $`T`$=0, the + and $``$ phases coexist along the diagonal direction in the $`h_0h_1`$ plane (until $`h_0=h_1=4J`$), as shown in the top left panel of figure 2. The $`+`$ phase coexists with the + phase along $`h_1=4J`$ for $`h_0>4J`$ and with the $``$ phase along $`h_0=4J`$ for $`h_1>4J`$. All of the phase boundaries at $`T`$=0 are first order lines denoted by dashed lines. The solid lines denote lines of continuous transitions. Two of these lines occur close to $`T=2J/k_B`$ and separate the $`+`$ phase from the $`+`$ and $``$ phases respectively. The star in the main phase diagram, where three critical lines come together, is a tricritical point. The critical line corresponding to the transition between + and $``$ starts at 6$`J/k_B`$ when $`h_0`$ and $`h_1`$ tend to zero and then decreases steadily as the fields are increased. In the vicinity of $`h_0=h_1`$=$`\mathrm{\hspace{0.33em}4}J`$, the descent towards the tricritical point is almost vertical.

A particularly simple case is obtained on fixing $`h_0`$ at the value of $`6J`$ (corresponding to strong pinning at the pore wall) and varying $`h_1`$ and $`T`$ to map out the coexistence curve between the + and $`+`$ phases. In order to cancel the effective field introduced by the nearly fully aligned $`m_0`$ spins, one needs to impose a field $`h_1=4J`$ (note that each site ’1’ has four ’0’ neighbors) and effectively one is left with the bulk Ising model. Figure 3 shows that the coexistence manifests itself in the presence of two values of $`m_1`$ but a unique value of $`m_0`$. Of course, a similar scenario takes place when the boundary between $``$ and $`+`$ is crossed. What is quite remarkable is that the topology of the phase diagram does not change even when randomness is introduced  in such a way that the symmetry between the + and $``$ phases is maintained. Thus within mean field theory, for the symmetric case, the random placement of the walls plays no role at all. We will show in the rest of the paper that the same result holds for more complex situations.

## III Effects of Confinement in the Random Field Ising Model – the Asymmetric Case

In porous media, the volume of a fluid near the pore walls is usually much less than the volume of the fluid in the interior. In the Ising spin model, this translates into an unequal number of sites with fields $`h_0`$ and $`h_1`$. In fact, in a random version of the model, we  considered a situation in which a fraction $`p`$ of the sites has a field $`h_0`$ – the symmetric case is obtained when $`p=\frac{1}{2}`$.

In order to study the effects of walls under such asymmetric conditions, we consider, for simplicity, the model shown in figure 1b which is a generalization of figure 1a. The plane of the figure shows nine sites. The central site, denoted by 0, has a local field of $`h_0`$. The remaining eight sites have a field of $`h_1`$ and they are denoted by 1 and 2. Thus $`p=\frac{1}{9}`$ but there is no randomness. The distinction between the two classes of sites, 1 and 2, is that the former have site 0 as a neighbor and the latter do not. Again, it is assumed that above and below the plane shown there are other planes which repeat the pattern of the central plane so that each site has a coordination number of six. Recall that the boundary conditions along the two directions within the plane are periodic.

The mean field equations for the three magnetizations read

$$m_0=tanh[(h_0+2Jm_0+4Jm_1)/k_BT],$$
(6)
$$m_1=tanh[(h_1+Jm_0+3Jm_1+2Jm_2)/k_BT],$$
(7)
$$m_2=tanh[(h_1+2Jm_1+4Jm_2)/k_BT].$$
(8)
The internal energy of the system is given by

$$U=J(m_0m_0+6m_1m_1+8m_2m_2+4m_0m_1+8m_1m_2)h_0m_0+4h_1m_1+4h_1m_2$$
(9)
and the entropy by

$$S=s_0+4s_1+4s_2.$$
(10)
This system continues to have three phases at $`T`$=0 as indicated in the top left panel of figure 4. The boundaries between the phases, however, are shifted to new locations. For instance, the + and $``$ phases coexist along the line $`h_0=8h_1`$, from the origin until $`h_0=4J`$ The + and $`+`$ phases coexist along $`h_1=\frac{1}{2}J`$ for $`h_0>4J`$ whereas the $``$ and $`+`$ phases coexist along $`h_0=4J`$ for $`h_1>J/2`$. The triple point at which all the three phases coexist is at $`T=0`$, $`h_0=4J`$ and $`h_1`$=$`\frac{1}{2}J`$.

The emergence of the three phases is the only similarity that exists between the symmetric and the asymmetric model. The way they coexist at non-zero temperatures, for instance, is quite different. The biggest distinction, shown in the main panel of figure 4, is that now the two sheets separating the + phase from the $`+`$ and the $`+`$ and $``$ phases along the diagonal direction combine together to form one surface. This surface has a tilt that is clearly visible on the right top panel of figure 4 which shows a section of the phase diagram at $`k_BT=5.5J`$. The surface terminates at a critical line which falls very gently from $`k_BT=6J`$ at $`h_0=0`$ and $`h_1`$=0 to about $`5.5J`$ at $`h_0=8.5J`$ and $`h_1=0.42J`$.

The coexistence surface of the $``$ and $`+`$ phases continues to be substantially planar with a critical line close to 2$`J`$ reflecting the one-dimensional connectivity of the 0 sites. This line of critical points intersects the combined $`+,`$ and $`+,+`$ coexistence surface at a critical end point at $`k_BT=2.05J`$, $`h_0=3.9612J`$, and $`h_1`$=$`\frac{1}{2}J`$. Quite remarkably, the topology of this phase diagram is exactly as in the random case .

The coexistence curves for $`h_0`$=$`6J`$ (again mimicking a strong pore-wall interaction) are shown in Figure 5. Physically, the transition corresponds to crossing from a phase in which the interior of the pore space is filled by liquid to one in which the liquid coats the walls and the vapor occupies the interior. Note the unusual geometry of the coexistence curve. The magnetizations $`m_1`$ and $`m_2`$ have broad coexistence curves, similar to $`m_1`$ of Figure 3 for the symmetric case. On the other hand, the coexistence curve for $`m_0`$ is much narrower than for $`m_1`$ and $`m_2`$ and its non-zero width arises when the values of $`m_1`$ and $`m_2`$ are distinct. When $`m_2=m_1`$, $`m_0`$ has a unique value in analogy to the symmetric case. It should be noted that there are just two coexisting solutions – the larger value of $`m_0`$ selects positive values of $`m_1`$ and $`m_2`$, whereas the smaller values correspond to negative $`m_1`$ and $`m_2`$.

## IV Novel Superfluid Phases in 3He–4He Mixtures in Aerogel

We turn now to a discussion of spin systems modeling the phase separation of 3He–4He mixtures. Figure 6(a) shows a sketch of the experimental phase diagram (in the temperature ($`T`$) - concentration of 3He ($`x`$) plane) of bulk 3He-4He mixtures in the vicinity of the superfluid transition of 4He. In the temperature range of interest, the superfluid transition involving the pairing of 3He atoms is not a factor and indeed the 3He atoms can be thought of as inert, annealed (i.e. they are not stuck in space but can move around) entities. At low 3He concentrations, on cooling the system, a superfluid transition denoted by the solid line (AB) is observed. However at higher 3He concentrations, the system opts to phase separate into a 4He rich region which becomes superfluid. The coexistence curve of the 3He-4He phase separation is shown as a dashed curve (CBD). Two interesting features of the phase diagram are the tricritical point B, where the superfluid transition line collides with the coexistence curve at its critical point and the miscibility gap at C – small amounts of 3He added to 4He do not lead to phase separation: a feature exploited in dilution refrigerators.

Perhaps, the simplest classical model that captures the topology of this phase diagram is the Blume-Emery-Griffiths model  (BEG) which is a lattice model populated with spins, $`S_i`$, that can take on one of three values 0, -1, or 1. The inert 3He is represented by 0 spins and 4He is denoted by +1 or 1 spin values. An exchange coupling between nearest neighbor non-zero spins, favoring alignment causes the analog of the superfluidity transition, with the broken symmetry phase having a non-zero magnetization (i.e. a mismatch in the number of +1 and 1 spins). The Hamiltonian reads

$$H=J\underset{<ij>}{}S_iS_j+\underset{i}{}\mathrm{\Delta }S_i^2,$$
(11)
where $`\mathrm{\Delta }`$ is an anisotropy field which controls the relative concentrations of the two isotopes. The presence of 4He corresponds to $`S_i=\pm 1`$, superfluidity of $`{}_{}{}^{4}He`$ to the existence of non-zero magnetization and the 3He atoms are represented by $`S_i=0`$. The random anisotropy field here does not break the $`\pm `$ symmetry.

The resulting phase diagram (shown in Figure 6b) has all the correct qualitative features, except for the absence of the miscibility gap at C which is believed to arise from a purely quantum mechanical effect. Even thought the BEG model is purely classical and does not have the correct symmetry of the spins (the superfluid transition has the same characteristics of the transition in a $`xy`$ model in which spins lie in a plane rather than having up-down symmetry), it nevertheless reproduces almost all the qualitative features of the experiment correctly.

Recent experiments of Chan and coworkers  on the phase separation of 3He-4He mixtures in aerogel in the vicinity of the superfluid transition have yielded a phase diagram shown Figure 6(c). The key features of the phase diagram are i) the absence of the tricritical point (the superfluid transition line no longer intersects the coexistence curve); ii) an enhancement of the superfluid transition temperature compared to the bulk at large 3He concentration iii) at low temperatures (below the critical point associated with the phase separation) and for a range of values of $`x`$ within the coexistence curve, 4He rich and 3He rich regions coexist, both of which are superfluid and iv) the experimental data, while restricted to temperatures above 0.35 K, are suggestive that the aerogel causes a miscibility gap to open up at large value of $`x`$. This is of fundamental importance, if true, since the superfluid phase observed is one in which a small quantity of 4He in 3He does not phase separate (as observed in the bulk), but is yet superfluid and probably represents the long sought after dilute Bose gas superfluid phase. Even more exciting, such a miscibility gap would lead to the extremely novel situation of two distinct coexisting superfluid phases at low temperatures, the dilute Bose gas phase of 4He and the superfluid phase of 3He. Two factors in support of the novel dilute superfluid phase are: a) Adding a small amount of 4He to the aerogel (in the absence of 3He) leads to a superfluid phase whose density is enhanced by the addition of 3He. b) Because a coexistence curve for the phase separation is found in the phase diagram, it is plausible that there is no phase separation in the region between B and D (6(c)), since it is unlikely to expect phase separation of an already phase separated phase.

## V The Random Anisotropy BEG Model

In the BEG model the effect of the aerogel is assumed to be present on a fraction p of the sites, these sites are randomly chosen and fixed (unlike the mobile 3He atoms, the aerogel is a manifestation of quenched randomness). The distribution of the single site anisotropies is bimodal and given by

$$P(\mathrm{\Delta }_i)==p\delta (\mathrm{\Delta }_i\mathrm{\Delta }_0)+(1p)\delta (\mathrm{\Delta }_i\mathrm{\Delta }_1).$$
(12)
The sites with anisotropy $`\mathrm{\Delta }_0`$ correspond to vicinity of pore walls and for the situation in which 4He prefers to be near the wall, $`\mathrm{\Delta }_0<0`$. The $`\mathrm{\Delta }_1`$ anisotropy characterizes the pore interior and its value controls the total number of 3He atoms. In a mean field approach, one obtains the phase diagram shown in Figure 6(d). Note that this is in accord with the experimental observations (i) and (iii), but does not reproduce (ii) and (iv). The tricritical point (where three phases go critical simultaneously) requires a special symmetry, which is absent when one incorporates the random anisotropy to mimic the aerogel. The line of superfluid transitions is, however, virtually insensitive to the presence of the random anisotropy. The coexistence curve, however, is shifted to lower temperatures and higher effective 4He concentration due to the space taken up by the aerogel, thus leading to the topology shown in Figure (6(d)).

In order to investigate whether this lack of complete agreement arises due to quantum mechanical effects and their neglect by the BEG model or due to the inherent simplicity of the mean field approach, we have carefully studied the BEG model within an improved mean field theory which is a generalization of the approach presented in Ref. . The improved method captures features such as a percolation threshold and yield better estimates of the transition temperature than the mean field theory.

The phase diagrams 6(e) and (f) are obtained depending on whether the fraction $`p`$ of the sites at which the aerogel is present is less than or greater than the percolation threshold $`p_c`$. Unlike the porous medium aerogel, which has a strongly correlated, connected interface, our model, in its simplest form, consists of randomly chosen interface sites allowing for a percolation threshold. In the experiment, in spite of the large porosity, one is always in the fully connected regime. Note the presence of a miscibility gap at high 3He concentrations. Unlike mean field theory (6(d)), the point D is at a concentration less than $`(1p)`$. Indeed, our calculations suggest a second coexistence curve between D and the point ($`x=1p`$; $`T`$=0) analogous to the one between C and D except that the critical temperature is shifted down to zero. Thus for concentrations of 3He corresponding to points between D and ($`x=1p`$; $`T`$=0), and temperatures less than the superfluid transition line AB, the model predicts the analog of the dilute Bose gas superfluid phase.

The superfluid transition temperature smoothly extrapolates to the value of the transition temperature of a coated phase of helium atoms residing on the aerogel surface. This is dramatically seen in Figure 6(f) where the transition plunges to zero at $`x=1p`$ when $`p`$ is less than the percolation threshold $`p_c`$ and is unable to sustain a phase transition at non-zero temperatures. Within the context of the BEG model, the novel superfluid phase is found to be one in which the magnetization (superfluidity) arises from the aerogel sites and from the sites in their vicinity. Indeed, in any classical model with short range interactions, the spins yielding a non-zero magnetization must lie on a connected cluster and are thus in an essentially phase separated phase. This phase separation, however, does not preclude a further bulk-like phase separation on increasing the concentration of 4He atoms. In our BEG model, the minimum number of 4He atoms is equal to the number of aerogel surface sites. A further addition of 4He atoms (in the absence of 3He) causes an increase in the magnetization, corresponding to the attachment of some of these atoms to the already existing spanning cluster at the aerogel surface. Subsequent addition of 3He atoms results in more of the 4He atoms going in the cluster, thereby enhancing the magnetization, as in the experiment. We have also studied the effects of correlation in the selection of aerogel surface sites: the probability of a nearest neighbor site of an aerogel surface site to be another aerogel surface site is enhanced compared to a completely random selection. We find this correlation enhances the superfluid transition temperature compared to the bulk in accord with experiment.

In summary, a simplified model for 3He-4He mixtures in aerogel reproduces many (but the miscibility gap at low 3He concentrations) of the features observed in experiments and suggests the opening of a miscibility gap at low 4He concentration. The analog of the dilute Bose gas phase within the classical model is one in which the superfluidity arises from 4He adsorbed on the aerogel. Note that this does not however preclude further phase separation. Cooling the system to ultralow temperatures in the miscibility gap at high 3He concentrations should lead to two coexisting superfluid phases: 3He and the 4He near the aerogel. It would be very exciting if quantum mechanical effects (not considered here) delocalize the 4He atoms leading to a dilute Bose gas phase at higher temperatures and interpenetrating 4He and 3He superfluid phases at low temperatures.

It should be noted that much of the physics corresponding to the scenarios of Figure 6 has been captured by the renormalization group analysis of Berker and his collaborators . They considered random and non-random ”jungle-gym” models of the aerogel and explained the phase diagrams by the connectivity and tenuousness of the aerogel.

## VI Effects of Confinement in the BEG model

In order to study the separate roles of the presence of walls and randomness on 3He–4He separation, we again consider the geometry shown in Figure 1b and set up mean field equations that correspond to the spin-1 problem. These equations for the magnetizations, i.e. the expectation values of $`S_i`$, and for the three parameters $`q_0`$, $`q_1`$, and $`q_2`$ which are the expectation values of $`S_i^2`$, are:

$$m_0=q_0tanh[(2Jm_0+4Jm_1)/k_BT],$$
(13)
$$m_1=q_1tanh[(Jm_0+3Jm_1+2Jm_2)/k_BT],$$
(14)
$$m_2=q_2tanh[(2Jm_1+4Jm_2)/k_BT],$$
(15)
and

$$q_i=\frac{4\mathrm{exp}(2\mathrm{\Delta }_i/k_BT)\sqrt{m_i^2+4\mathrm{exp}(2\mathrm{\Delta }_i/k_BT)(1m_i^2)}}{4\mathrm{exp}(2\mathrm{\Delta }_i/k_BT)1}.$$
(16)
The internal energy is given by

$$U=J(m_0m_0+6m_1m_1+8m_2m_2+4m_0m_1+8m_1m_2)+\mathrm{\Delta }_0q_0+4\mathrm{\Delta }_1q_1+4\mathrm{\Delta }_1q_2$$
(17)
and the entropy by

$$S=\stackrel{~}{s}_0+4\stackrel{~}{s}_1+4\stackrel{~}{s}_2,$$
(18)
where

$$\stackrel{~}{s}_i=k_B[(1q_i)ln(1q_i)+\frac{1}{2}(q_i+m_i)ln(q_i+m_i)/2+(q_im_i)ln(q_im_i)/2].$$
(19)
In direct analogy to the random case , there are four posssible phases at $`T`$=0:

phase 1 in which all three $`m_i`$s are non-zero;

phase 2 in which $`m_0>0`$, $`m_1=m_2=0`$;

phase 3 in which $`m_0=0`$, $`m_1`$, $`m_2>0`$; and

phase 4 in which $`m_0=m_1=m_2=0`$.

In each phase, $`q_i`$=$`m_i`$ at $`T`$=0. The non-zero magnetization persists to higher temperatures and its disappearance corresponds to the $`\lambda `$ line for superfluid 4He. (Note that our analysis lumps in any inert or dead layer of 4He as belonging to the pore wall.) The analog of the 3He concentration is given by

$$x=\mathrm{\hspace{0.33em}1}(q_0+4q_1+4q_2)/9.$$
(20)
The overall topology of the phase diagram is shown on figure 7 and several isothermal slices through it are shown in figure 8. At low temperatures all four phases exist and their number goes down on increasing the temperature. The $`T`$=0 boundaries are given by

1 – 2 coexistence at $`\mathrm{\Delta }_1=\frac{13}{4}J`$, $`\mathrm{\Delta }_0<J`$;

1 – 3 coexistence at $`\mathrm{\Delta }_0=5J`$, $`\mathrm{\Delta }_1<\frac{11}{4}J`$;

1 – 4 coexistence $`\mathrm{\Delta }_1=\frac{27}{8}J\frac{1}{8}\mathrm{\Delta }_0`$, $`J<\mathrm{\Delta }_0<5J`$ ;

2 – 4 coexistence at $`\mathrm{\Delta }_0=J`$, $`\mathrm{\Delta }_1>\frac{13}{4}J`$;

3 – 4 coexistence at $`\mathrm{\Delta }_1=\frac{11}{4}J`$, $`\mathrm{\Delta }_0>5J`$.

These first-order boundaries become vertical surfaces on considering the $`T`$-axis. The top edges of these surfaces are critical lines. This has as its roof a critical surface at which the magnetization disappears. Of course phase 4, which is paramagnetic, is not covered by a roof. This roof corresponds to the superfluid transition of 4He or the $`\lambda `$ \- line. The shape pf the roof is illustrated in figure 9 for several values of $`\mathrm{\Delta }_0`$. In the figure, for $`\mathrm{\Delta }_0/J0.8`$, the roof continues indefinitely for large $`\mathrm{\Delta }_1`$, becaouse $`m_0`$ remians zero at sufficiently low temperatures. However, when $`\mathrm{\Delta }_0/J`$, the roof smoothly terminates at the top of the wall separating planes 1 and 4. This is necessitated by the fact that phase 4 has no roof.

Figure 10 illustrates the nature of the phase diagram for selected values of $`\mathrm{\Delta }_0`$. The insets show the transition lines as a function of $`\mathrm{\Delta }_1`$ and the main figures – as a function of $`x`$ – the analog of the 3He concentration. The top two panels of Figure 10 refer to the uniform anisotropy case - when $`\mathrm{\Delta }_0`$ is equal to $`\mathrm{\Delta }_1`$ and confirm that this simple nine-spin model captures the topology of the phase diagram of the uniform BEG model .

The physically interesting regime is that of negative $`\mathrm{\Delta }_0`$ which favors 4He near the pore walls. It is seen that, as $`\mathrm{\Delta }_1`$ is varied, the $`\lambda `$ line becomes disconnected from the phase separation coexistence lines (the middle panels). The details depend on how one moves on the $`\mathrm{\Delta }_0`$ \- $`\mathrm{\Delta }_1`$ plane. For instance, if one crosses the 1–2 and 2–4 boundaries at an angle, one gets the situation depicted in the bottom panels of Figure 10. All of these phase diagrams are in accord with the random anisotropy version of the model except that the $`\lambda `$ line of the bottom panels of figure 10 does not reemerge from the 3He-4He coexistence region because, this region ends at $`x`$=1 in this simple model, and not at an $`x`$ which is less than 1. Figure 11 shows the coexistence curves for the $`q`$ order parameters for $`\mathrm{\Delta }_0`$=-0.5$`J`$. They are remarkably similar to the magnetization coexistence curves of Figure 6. In particular, the non-zero width for the coexistence region in $`q_0`$ reflects inequality of $`q_1`$ and $`q_2`$.

The basic message of our analysis is that the topology of the phase diagram changes qualitatively (in accord with experiment) when one moves in the $`(\mathrm{\Delta }_0\mathrm{\Delta }_1)`$ plane such that instead of going directly from phase 1 to the paramagnetic phase 4, one moves through the intermediate phase 2. Somewhat surprisingly and at odds with expectations, one finds that the mean field topology of the phase diagram is more sensitive to the presence of walls in a porous medium than to the role played by their random placement. A similar conclusion has been reached by Pricaupenko and Treiner  within a nonlocal density functional analysis of 3He-4He mixtures in a channel geometry. In particular their analysis shows possibility of the detachment of the superfluid line from the coexistence region. It would be interesting to consider whether fluctuations make a qualitative difference in the conclusions reached in our simple mean field analysis.

We are indebted to Moses Chan for stimulating discussions. This work was supported by the Center of Collective Phenomena in Restricted Geometries, the Penn State MRSEC under NSF grant DMR-0080019, INFM, MURST, and NASA.

FIGURE CAPTIONS

Figure 1. a) The basic unit of the model used to study the symmetric case in which there are as many sites with the local field $`h_0`$, denoted by 0, as with the field $`h_1`$, denoted by 1. The 0s and 1s are placed on two sublattices in the plane as shown and is periodically repeated in both directions in the plane of the paper. The model has lines of 0 and 1 sites respectively perpendicular to the plane of the paper. In other words, the pattern shown is repeated in other parallel planes. b) The geometry of the model used to study the asymmetric case. The site denoted by 0 has a field $`h_0`$ and the remaining sites a field $`h_1`$. As before, there are periodic boundary conditioons in the plane and repeat boundary conditions in the direction perpendicular to the page.

Figure 2. The main panel shows the phase diagram corresponding to the symmetric case in a three dimensional representation $`h_0h_1T`$. The dashed lines correspond to first order transitions whereas the thick solid lines correspond to the continuous transitions. The star indicates a tricritical point. The top panels show constant temperature slices of the phase diagram for the temperatures indicated. The asterisk on the right hand panel indicates a critical point and the square on the top left panel is a triple point.

Figure 3. The temperature dependence of the local magnetizations, $`m_0`$ and $`m_1`$ on crossing the boundary between the + and $`+`$ phases at $`h_0=6J`$. For temperatures below $`T=2J`$, there is a coexistence of two values of $`m_1`$ but $`m_0`$ stays essentially fully magnetized up to $`T=\mathrm{\hspace{0.33em}2}J`$.

Figure 4. The main panel shows the phase diagram corresponding to the asymmetric case in a three dimensional representation $`h_0h_1T`$. The dashed lines correspond to first order transitions whereas the thick solid lines correspond to continuous transitions. The pentagon indicates a critical end point. The top panels show constant temperature slices of the phase diagram for the temperatures indicated. The asterisk on the right hand panel indicates a critical point, whereas the square on the top left hand panel is a triple point.

Figure 5. The temperature dependence of the local magnetizations, $`m_0`$, $`m_1`$, and $`m_2`$ on crossing the boundary between the + and $`+`$ phases at $`h_0=6J`$. The coexistence curve for $`m_0`$ is much narrower than for $`m_1`$ and $`m_2`$. There are two coexisting solutions and the larger value of $`m_0`$ selects positive values of $`m_1`$ and $`m_2`$, whereas the smaller values correspond to negative $`m_1`$ and $`m_2`$.

Figure 6. Schematic representations of the experimentally (panels a and c) and theoretically determined phase diagrams for 3He-4He separation as described in more detail in the text. The phase diagram is shown in the temperature ($`T`$) – concentration of 3He ($`x`$) plane. Panels a) and b) correspond to experimental and theoretical results for the bulk case. All other panels refer to random situations. Panel c) corresponds to the experiments in aerogels. Panel d) corresponds to the mean field analysis of the random anisotropy BEG model with a bimodal distribution of the anisotropies. Panels e) and f) correspond to theoretical results in which $`p`$, the fraction of randomly chosen sites corresponding to the pore walls, is smaller or larger than the percolation threshold respectively.

Figure 7. The phase diagram for the BEG model with a bimodal non-random distribution of the anisotropies in a three-dimensional representation $`\mathrm{\Delta }_0\mathrm{\Delta }_1T`$. The broken lines correspond to the first order transitions whereas the solid lines correspond to continuous transitions.

Figure 8. Constant temperature slices of the phase diagram shown in Figure 7 at the temperatures indicated.

Figure 9. Plot of the critical lines (the temperature at which the magnetization goes to zero) as a function of $`\mathrm{\Delta }_1`$ for the indicated values of $`\mathrm{\Delta }_0`$.

Figure 10. Top left panel: The phase diagram for the uniform (bulk) anisotropy BEG model in the $`Tx`$ plane, where $`x`$ is the spin analog of the concentration of He3. Top right panel: The same phase diagram but in the $`T\mathrm{\Delta }_1`$ plane. $`\lambda `$ indicates the critical line for the magnetization and I the first order transition in the $`q`$-order parameter. Middle left panel: The phase diagram for the non-uniform (periodic) anisotropy BEG model in the $`Tx`$ plane. The data are for $`\mathrm{\Delta }_0=0.5J`$. Middle right panel: The corresponding phase diagram in the $`T\mathrm{\Delta }_1`$ plane. The bottom panels: The phase diagrams for for $`\mathrm{\Delta }_0=\mathrm{\Delta }_1/5`$.

Figure 11. The coexistence curves for the order parameter $`q`$ on the three sites when $`\mathrm{\Delta }_0=0.5J`$. There are two coexisting solutions and the larger value of $`q_0`$ selects larger values of $`q_1`$ and $`q_2`$.


format: text

text:
# A learning approach to the detection of gravitational wave transients

## 1 Motivations

A number of large scale gravitational wave interferometric detectors such as LIGO, TAMA, GEO600 and VIRGO  are taking scientific data or will reach this goal soon. The objective of this paper is to contribute to the arsenal of detection algorithms able to locate the weak and short signature of a gravitational wave of astrophysical origin in the long data stream produced by the detector.

In the list of candidates having a good chance for the first detection, there are sources for which we can only make a rough guess or simulate the highly non-linear Physics which is involved. This causes the expected gravitational waveforms to be poorly modeled. Most of these are collapses of very massive astrophysical objects in the final stage, and this yields the resulting gravitational wave to be a burst.

In such cases, computer simulations may give good indications of what can be the waveform for some choices of the parameter values describing the physical phenomenon. However, it is generally not possible to have a tight sampling of the parameter space i.e., to scan a large range of physical configurations. We only have at our disposal a catalog of waveforms, whose members are selected representatives of the large set of possibilities.

Two examples of such sources are supernovae core collapse and binary black hole merger. Despite some recent progresses, work is still needed in the latter case to produce reliable waveforms in a realistic set up (including spins for instance). Concerning the former case, hydrodynamic simulations of relativistic supernovae have been recently computed  and the expected waveforms for different parameter configurations were made available.

In this paper, we propose a method for designing systematically a decision statistic for the detection of gravitational transients by extracting the necessary information from a catalog of test waveforms emitted from a targeted source. We use the supernovae waveforms as one possible application. Although not considered in this paper, the presented approach may also apply to other problems encountered when analyzing the output of a gravitational wave interferometer, such as the classification and the characterization of the noisy transients. (Because they worsen the detector sensitivity, such interferences deserve the development of algorithms to determine their actual origin.) In this context, the initial database could be a collection of characteristic individuals extracted “by hand” or with some other simple algorithm. In any case, we refer the (gravitational wave or noise) transient(s) we want to detect to as signal.

Some attention has been paid to the choice of an adequate vector formalism to treat the problem and make the implementation on computers easier. The resulting notations are defined in Sect. 2. This section also includes the formulation, within the chosen framework, of classical results such as the Plancherel formula which will be of use further.

In Sect. 3, we describe the detection problem we consider with the accompanying hypothesis. We assume the signal to be random and of unknown probability density function (PDF). This assumption translates explicitly the lack of knowledge about the signal. The only piece of information at our disposal is its first and second order statistical moments (i.e., its mean and covariance). In the situation of interest here, these two quantities are not known, but they can be estimated with sufficient accuracy from available sample sets.

We consider that the noise is Gaussian and stationary and that we know its correlation function (or equivalently its power spectral density). Analogously to the signal, an extension to the case where there is no reasonable noise model, is possible through the use of estimates done with “noise only” data streams.

Since we do not have the signal PDF, it is impossible to write the exact form of the likelihood ratio, and thus to obtain the optimal statistic. However, a satisfactory solution can be obtained by first imposing the mathematical structure of the statistic and second look in the selected set of functions for the best element by maximizing a contrast criterion.

The difficulty then lies in making the choice of a sufficiently general class of statistics and a sensible criterion for the considered problem. In Sect. 4, it is explained why the family of quadratic detectors (i.e., the detection statistic is a quadratic form of the observed data) and a measurement of the deflection (a quantity homogeneous to the signal to noise ratio) are good candidates. Coming back to our initial detection problem, we individuate in the selected class the statistic which performs best according to the chosen criterion and show that it can be expressed easily with the signal and noise covariance. Our proof was inspired by the work presented in  and  which we adapt to the case of interest (finite vector spaces and non-centered signals). We finish Sect. 4 showing that the proposed approach is not unfamiliar since it can be related to the well-known method of matched filtering.

In Sect. 5, we put the quadratic detector of best deflection into practice with the problem of detecting gravitational wave transients from supernovae core collapse. In this specific case, we show how the signal covariance matrix can be estimated from the catalog of simulated waveforms. This yields a simplification of the detector and an efficient implementation for the online detection. We give also some results about the determination of the decision threshold required to get a chosen false alarm rate.

The vector formalism introduced here is a general framework in which all quadratic detectors can be easily related and compared. In Sect. 6, we do this comparative study between the solution with optimal deflection and other detection techniques  proposed in the literature that are also belonging to the quadratic class.

## 2 Notations and basic algebra

We will denote scalar quantities and scalar-valued functions with plain italics, e.g., $`x`$; vectors by boldface letters, e.g. $`𝐱`$; and matrices by boldface capitals, e.g., $`𝐗`$. We will represent the components of vectors and matrices with superscripts within brackets, e.g. $`𝐗^{(m,n)}`$ designates the element located in $`m`$th row and $`n`$th column of the matrix $`𝐗`$. Finally, $`𝐱^t`$ denotes the transpose of the vector $`𝐱`$. The symbol $``$ will be used in the following to define our variables and therefore stands for “equal by definition.”

We use curve brackets for denoting continuous (random and/or deterministic) time (or frequency) series (e.g. $`x(t)`$), whereas squared brackets are employed for discrete time (sampled) processes (e.g., $`x[k]`$). The samples of a discrete time signal are collected in a single column vector of $`^N`$, e.g.

$$𝐱(x[k]=x(t_sk),k=0\mathrm{}N1)^t$$
(1)
where $`t_s1/f_s`$ is the sampling period and $`f_s`$ the sampling rate.

We define the Fourier transform of $`x[k]`$ by :

$$X(f)t_s\underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}x[k]e^{2i\pi kf/f_s}.$$
(2)
As a general rule, Fourier transforms are denoted with the same (capital) letter used for its associated time sequence. The function $`X(f)`$ is $`f_s`$–periodic (i.e., $`X(f)=X(f+f_s)`$ for all $`f`$) and its inverse may be calculated with the following inversion equation :

$$x[k]=_{f_s/2}^{f_s/2}X(f)e^{2i\pi kf/f_s}𝑑f.$$
(3)
We recall that the Plancherel formula relates scalar products expressed in the time and frequency domains :

$$t_s\underset{j=\mathrm{}}{\overset{+\mathrm{}}{}}x[j]y[j]=_{f_s/2}^{+f_s/2}X(f)\overline{Y(f)}𝑑f.$$
(4)
This equation may also be expressed using vector scalar product, provided the assumption that one of the two signals $`x[k]`$ or $`y[k]`$ has a finite support (denoted with $`\mathrm{supp}\{\}`$), as specified in the following lemma :

###### Lemma 1.

Assuming that $`\mathrm{supp}\{y[]\}\{0\mathrm{}N1\}`$, the Plancherel formula writes :

$$t_s𝐱^t𝐲=_{f_s/2}^{+f_s/2}X(f)\overline{Y(f)}𝑑f$$
(5)
Continuing along the same idea, the convolution of two signals $`x[k]`$ and $`y[k]`$ defined by :

$$z[k]=t_s\underset{j=\mathrm{}}{\overset{+\mathrm{}}{}}x[kj]y[j]$$
(6)
or equivalently with $`Z(f)=X(f)Y(f)`$, may also be rewritten using vectors under some support constraints as in the next lemma.

###### Lemma 2.

Let $`N_y<(N1)/2`$ be a positive integer and suppose that $`\mathrm{supp}\{y[]\}=\{N_y\mathrm{}N_y\}`$, the collection of samples $`𝐳=(z[N_y]\mathrm{}z[N1N_y])^t`$ where $`z[k]`$ is the convolution of $`x[k]`$ and $`y[k]`$ as defined in eq. (6), can be expressed as :

$$𝐳=t_s\mathrm{𝐘𝐱}$$
(7)
where $`𝐘^{N2N_y\times N}`$ is a matrix of the form

$$𝐘=\left(\begin{array}{cccccc}y[N_y]& \mathrm{}& y[N_y]& 0& 2& 0\\ 0& y[N_y]& \mathrm{}& y[N_y]& 0& \mathrm{}& 0\\ 7\\ 0& 2& 0& y[N_y]& \mathrm{}& y[N_y]\end{array}\right).$$
(8)
Proofs of both Lemma 1 and 2 are simple and left to the reader.

## 3 Problem statement

The question we consider here is to detect a (possibly non-stationary) random signal in a stationary Gaussian noise (the signal and noise covariance function being known or could be estimated with accuracy). Using the notations defined in the previous section, the problem is to distinguish between two statistical hypothesis $`(H_0)`$ and $`(H_1)`$ :

$`(H_0):𝐱`$ $`=𝐧`$ (9a)
$`(H_1):𝐱`$ $`=𝐬+𝐧`$ (9b)
with the following assumptions

1. the signal $`𝐬`$ is a random vector of mean $`𝐬_m𝖨𝖤[𝐬]`$ (where $`𝖨𝖤[]`$ denotes the expectation operator) and correlation matrix $`𝐑_s𝖨𝖤[(𝐬𝐬_m)(𝐬𝐬_m)^t]`$,
2. the noise $`𝐧`$ is a zero-mean, stationary Gaussian vector with correlation matrix $`𝐑_n`$. Note that, since $`𝐧`$ is stationary, $`𝐑_n`$ is a Toeplitz symmetric matrix (the terms of $`𝐑_n`$ are given by the autocorrelation function $`𝐑_n^{(j,j+k)}𝖨𝖤[n[j]n[j+k]]`$),
3. the signal and the noise are decorrelated, meaning that $`𝖨𝖤[𝐧^t(𝐬𝐬_m)]=0`$,

This set of assumptions correspond to several different practical situations. A first one is when the signal or noise models are good enough to get reliable close form expressions of $`𝐬_m`$, $`𝐑_s`$ and $`𝐑_n`$. The second situation is when a sufficiently large set of “signal only” and/or “noise only” realizations is available and can be used to obtain a good estimate of the first and second order moments of the signal and the noise. Note that, except for its first and second order moments, we made no hypothesis about the PDF of $`𝐬`$.

## 4 Quadratic detectors

Deciding $`(H_1)`$ or $`(H_0)`$ is classically done by finding a partition function $`\mathrm{\Lambda }()`$ dividing the observation space (here, $`^N`$) into two disjoint subsets :

$`\mathrm{\Lambda }(𝐱)`$ $`\eta `$ decide $`(H_1)`$ (10a)
$`\mathrm{\Lambda }(𝐱)`$ $`<\eta `$ decide $`(H_0)`$ (10b)
where the detection threshold $`\eta `$ defines the border between the two decision area. Its value is given by fixing to a reasonable value the probability of deciding upon hypothesis $`(H_1)`$ although no signal $`𝐬`$ is present, which we refer to “false alarm probability”. The function $`\mathrm{\Lambda }()`$ is referred to as detection statistic or simply detector.

### 4.1 Intuitive background

It is intuitive to search for some unknown signal by looking for abnormal excess of power in one or several frequency bands of the observed signal spectrum. To implement this idea, we define the power of a signal $`x`$ using a $`l^2`$ measure :

$$E_x1/N\underset{k=0}{\overset{N1}{}}x[k]^2=𝐱^t𝐱/N,$$
(11)
and a bank of filters which select adequately the signal in the frequency bands of interest. Let $`\{g_m[k],k=0\mathrm{}N1\text{ and }m=0\mathrm{}M1\}`$ be the impulse responses (assumed to be of finite support) of the chosen bandpass filters. We get the signal $`y_m[k]`$ at the output of each filter by convolving the observed signal $`x[k]`$ to the corresponding impulse response. With the constraint that $`\mathrm{supp}\{g_m\}\{N_g\mathrm{}N_g\}`$ for all $`m`$ where $`N_g<(N1)/2`$, we can apply the Lemma 2 and express the output signal in a vector form as :

$$𝐲_m=t_s𝐆_m𝐱,$$
(12)
where $`𝐲_m`$ and $`𝐆_m`$ are as defined in Lemma 2.

Hence, we can write down the detection statistic corresponding to the basic idea mentioned above by summing up the power in all $`M`$ bands which yields :

$$\mathrm{\Lambda }_{intuitive}(𝐱)\underset{m=0}{\overset{M1}{}}E_{y_m}^2=𝐱^t\left(\frac{t_s^2}{N2N_g}\underset{m=0}{\overset{M1}{}}𝐆_m^t𝐆_m\right)𝐱.$$
(13)
We conclude that the heuristic principle we chose, is practically implemented with a detection statistic which is a quadratic form of the data. Extending this result to cases where the kernel of the form is an arbitrary symmetric matrix, this leads us to define the following family of detectors :

###### Definition 1.

Let $`𝐀^{N\times N}`$ be a symmetric real matrix, the function

$$\mathrm{\Lambda }_𝐀(𝐱)=𝐱^t\mathrm{𝐀𝐱},$$
(14)
defines a quadratic detector of kernel $`𝐀`$.

Quadratic detectors will be a central ingredient in this paper. It is worth noting that they have been extensively used in many applications (see e.g.  for a list of examples). In the case of gravitational wave detection, the specific area of interest here, several works  are based on such detector structure.

Beyond qualitative arguments, the following theoretical result is an important motivation for the use of quadratic detectors  : with the additional assumption of a zero-mean Gaussian signal (i.e., both signal and noise are Gaussian), the optimal solution in the Neymann-Pearson sense of the problem (9) belongs to the family defined in Def. 1. Although the problem considered here excludes the signal to be Gaussian, we expect quadratic detectors to retain their good performances, when the signal PDF is close to Gaussian.

### 4.2 Optimal quadratic detectors

For our problem (9), we don’t have a complete knowledge of the statistics of the input signal (the PDF of $`𝐬`$ is unknown). In consequence, we cannot write out the likelihood ratio which gives the best (in the Neymann-Pearson sense) detector among all possibilities.

We propose to overcome this difficulty by first, reducing the class of possible solutions to the family of the quadratic detectors defined in Def. 1 and second, extracting from this smaller set the statistic which will best perform for our problem. More precisely, our objective is to get the quadratic detector (i.e., get the kernel matrix $`𝐀`$ which identifies this quadratic detector in the whole family) which maximizes the following contrast criterion based on the first and second order moments :

$$d^2(\mathrm{\Lambda }_𝐀)=\frac{(𝖨𝖤[\mathrm{\Lambda }_𝐀(𝐱)|H_1]𝖨𝖤[\mathrm{\Lambda }_𝐀(𝐱)|H_0])^2}{\mathrm{var}\{\mathrm{\Lambda }_𝐀(𝐱)|H_0\}}.$$
(15)
This criterion is generally referred to as signal-to-noise ratio (statistics) or deflection (signal processing). The terminology “signal-to-noise ratio” is generally associated in most of the literature about the gravitational wave data analysis to the quantity in eq. (15) where $`\mathrm{\Lambda }()`$ is set to the matched filter statistic. In consequence, we adopt the term “deflection” to avoid confusion. The deflection may be viewed as a contrast measurement between the two statistical hypothesis $`H_0`$ and $`H_1`$ in the sense it measures the distance between the centers of the PDF of the statistic in the two hypothesis relatively to the PDF width in the null hypothesis $`H_0`$.

In the context of the problem (9), we apply now this approach for selecting the best statistic among all quadratic detectors.

###### Lemma 3.

In the situation described in eq. (9), the deflection of a quadratic detector as defined in Def. 1 is

$$d^2(\mathrm{\Lambda }_𝐀)=\frac{\mathrm{tr}^2\{\mathrm{𝐀𝐂}_s\}}{2\mathrm{t}\mathrm{r}\{(\mathrm{𝐀𝐑}_n)^2\}},$$
(16)
where $`𝐂_s𝖨𝖤[𝐬^t𝐬]`$ defines the (non-central) covariance matrix and $`\mathrm{tr}\{\}`$ is the trace 111Let $`𝐀^{N\times N}`$, the operator $`\mathrm{tr}\{𝐀\}_{n=0}^N𝐀^{(n,n)}`$ defines the trace of $`𝐀`$. operator.

###### Proof.

The proof of this lemma may be found in other articles  for zero mean signals and infinite vector spaces. Here, we give the proof for non-central signals (i.e, $`𝐬_m\mathrm{𝟎}`$) and in the case of discrete signals forming vectors of finite size.

We compute the first two statistical moments of $`\mathrm{\Lambda }_𝐀(𝐱)`$ under hypothesis $`H_0`$ and $`H_1`$. Starting conditionally to $`H_0`$, we have

$$𝖨𝖤[\mathrm{\Lambda }_𝐀(𝐱)|H_0]=𝖨𝖤[𝐧^t\mathrm{𝐀𝐧}].$$
(17)
Using the identity $`𝐱^t𝐱=\mathrm{tr}\{\mathrm{𝐱𝐱}^t\}`$, this can be reduced in

$$𝖨𝖤[\mathrm{\Lambda }_𝐀(𝐱)|H_0]=\mathrm{tr}\{𝐀𝖨𝖤[\mathrm{𝐧𝐧}^t]\}=\mathrm{tr}\{\mathrm{𝐀𝐑}_n\}.$$
(18)
Under the “signal+noise” $`H_1`$ hypothesis, we expand the quadratic form

$$𝖨𝖤[\mathrm{\Lambda }_𝐀(𝐱)|H_1]=𝖨𝖤[𝐬^t\mathrm{𝐀𝐬}+𝐬^t\mathrm{𝐀𝐧}+𝐧^t\mathrm{𝐀𝐬}+𝐧^t\mathrm{𝐀𝐧}]$$
(19)
Then, with the identity mentioned previously, we can simplify the auto-terms in the expansion $`𝖨𝖤[𝐬^t\mathrm{𝐀𝐬}]=\mathrm{tr}\{\mathrm{𝐀𝐂}_s\}`$ with $`𝐂_s𝖨𝖤[𝐬^t𝐬]`$ while the cross terms vanish : $`𝖨𝖤[𝐧^t𝐬]=𝖨𝖤[𝐧^t]𝐬_m=0`$ and $`𝖨𝖤[𝐧^t\mathrm{𝐀𝐬}]=𝖨𝖤[𝐬^t\mathrm{𝐀𝐧}]=0`$, which yields

$$𝖨𝖤[\mathrm{\Lambda }_𝐀(𝐱)|H_1]=\mathrm{tr}\{𝐀(𝐂_s+𝐑_n)\}.$$
(20)
For the general expressions of higher order moments of Gaussian quadratic forms in , we get the variance under $`H_0`$ :

$$\mathrm{var}(\mathrm{\Lambda }_𝐀(𝐱)|H_0)=\mathrm{var}(𝐧^t\mathrm{𝐀𝐧})=2\mathrm{tr}\{(\mathrm{𝐀𝐑}_n)^2\}.$$
(21)
The combination of all these ingredients leads to the result. ∎

In order to find the best detector in the quadratic family, we need now to find which kernel matrix maximizes the deflection. This is stated by the following proposition.

###### Proposition 1.

There exists a unique symmetric matrix $`𝐇`$ such that $`d^2(\mathrm{\Lambda }_𝐇)`$ obtained in Lemma 3 is maximum, and this matrix is

$$𝐇=\mathrm{argmax}_𝐀\{d^2(\mathrm{\Lambda }_𝐀)\}=𝐑_n^1𝐂_s𝐑_n^1,$$
(22)
###### Proof.

We first note that the noise autocorrelation $`𝐑_n`$ is symmetric positive definite matrix. Therefore, there exists a triangular matrix $`𝐓_n`$ with positive diagonal such that $`𝐑_n=𝐓_n𝐓_n^t`$. This factorization method is referred to as Cholesky factorization .

It is useful to introduce the two following matrices $`𝐆𝐓_n\mathrm{𝐀𝐓}_n^t`$ and $`𝐂\left(𝐓_n^1\right)^t𝐂_s𝐓_n^1`$, which we make appear in the expression of the deflection we got in Lemma 3, then reducing to :

$$d^2(\mathrm{\Lambda }_𝐀)=\frac{\mathrm{tr}^2\{\mathrm{𝐆𝐂}\}}{2\mathrm{tr}\{𝐆^2\}}.$$
(23)
Let $`𝒮_N()`$ be the vector space of real symmetric matrices of $`^{N\times N}`$. It is easily shown that $`𝐀,𝐁_{𝒮_N()}\mathrm{tr}\{\mathrm{𝐀𝐁}\}`$ defines a scalar product on $`𝒮_N()`$. Since the matrices $`𝐆`$ and $`𝐂`$ belong to $`𝒮_N()`$, we can rewrite the deflection as a ratio of scalar products, also referred to as Rayleigh quotient , namely :

$$d^2(\mathrm{\Lambda }_𝐀)=\frac{𝐆,𝐂_{𝒮_N()}^2}{2𝐆,𝐆_{𝒮_N()}}.$$
(24)
Using the Cauchy-Schwarz inequality, we conclude that $`d^2(\mathrm{\Lambda }_𝐀)`$ is maximum if and only if $`𝐆𝐂`$. Setting $`𝐆=𝐂`$ without loss of generality and replacing $`𝐆`$ and $`𝐂`$ by their definition directly yields eq. (22). ∎

We now have the final expression of the quadratic detector reaching the deflection optimum, namely

$$\mathrm{\Lambda }_𝐇(𝐱)=𝐱^t\mathrm{𝐇𝐱}\text{given }𝐇=𝐑_n^1𝐂_s𝐑_n^1.$$
(25)
Before looking how the approach proposed here may be practically implemented, we first give some interpretations of the detection statistic we obtained.

### 4.3 Interpretation and relation to matched filtering

With a direct calculation from its definition, we can separate $`𝐂_s𝖨𝖤[𝐬^t𝐬]`$ into two terms $`𝐂_s=𝐬_m𝐬_m^t+𝐑_s`$. The mean $`𝐬_m`$ can be viewed as a trend of the signal which happens systematically. In this sense, we refer the first term to as “deterministic”. The correlation matrix is related to the typical amplitude of the random fluctuations superimposed to the mean. For this reason, we refer the second term to as “random”.

Injecting this expansion in (25), we obtain a similar separation of $`\mathrm{\Lambda }_𝐇(𝐱)`$

$$\mathrm{\Lambda }_𝐇(𝐱)=\mathrm{\Lambda }_𝐇^{det}(𝐱)+\mathrm{\Lambda }_𝐇^{rand}(𝐱),$$
(26)
where $`\mathrm{\Lambda }_𝐇^{det}(𝐱)=(𝐬_m^t𝐑_n^1𝐱)^2`$ is related to the deterministic part of the signal model and $`\mathrm{\Lambda }_𝐇^{rand}(𝐱)=𝐱^t𝐑_n^1𝐑_s𝐑_n^1𝐱`$ to its random part.

The two contributions of the detection statistic are worth to be investigated further. An interesting interpretation and a link to matched filtering  results from the reformulation in the frequency domain of $`𝐲^t𝐑_n^1𝐱`$ where $`𝐱`$ and $`𝐲^N`$. This is the objective of the Proposition whose proof is detailed in the next section.

#### 4.3.1 Towards matched filtering

As the preamble, we define formally the whitening operation (i.e., filtering the signal with the inverse of the square root of the noise power spectral density) in the vector formalism used here. Let $`\mathrm{\Gamma }_n(f)𝖨𝖤[|N(f)|^2]`$ be the power spectral density of the noise $`n[k]`$ and $`\stackrel{ˇ}{x}[k]`$ be the result of the whitening of a given signal $`x[k]`$, we define

$$\stackrel{ˇ}{X}(f)\frac{X(f)}{\sqrt{\mathrm{\Gamma }_n(f)}}.$$
(27)
Clearly, $`\stackrel{ˇ}{x}[k]`$ is the result of the convolution of the signal $`x[k]`$ with the whitening filter of impulse response $`w[k]`$,

$$w[k]=_{f_s/2}^{+f_s/2}\frac{1}{\sqrt{\mathrm{\Gamma }_n(f)}}e^{2\pi ikf/f_s}𝑑f.$$
(28)
Using the Lemma 2, we deduce that the whitening filter defined in (27) can be written in the vector formalism with

$$\stackrel{ˇ}{𝐱}=t_s\mathrm{𝐖𝐱},$$
(29)
where $`\stackrel{ˇ}{𝐱}`$ and $`𝐖`$ are as given in Lemma 2.

This expression of the whitening operation is needed in the proof of the following Proposition, in which we get the vector form of the matched filtering statistic provided some care for the cancellation of finite size effects.

###### Proposition 2.

Let $`N_w<(N1)/2`$ be the half-size of $`\mathrm{supp}\{w\}=\{N_w\mathrm{}N_w\}`$ i.e., the support of the impulse response of the whitening filter defined in eq. (28), and let $`y[k]`$, a signal whose whiten version has a finite size support, $`\mathrm{supp}\{\stackrel{ˇ}{y}\}\{0\mathrm{}N1\}`$, then

$$𝐲^t𝐑_n^1𝐱=_{f_s/2}^{+f_s/2}\frac{X(f)\overline{Y(f)}}{\mathrm{\Gamma }_n(f)}𝑑f.$$
(30)
###### Proof.

We first compute $`𝖨𝖤[\stackrel{ˇ}{𝐧}\stackrel{ˇ}{𝐧}^t]`$ and get a first expression using eq. (29) :

$$𝖨𝖤[\stackrel{ˇ}{𝐧}\stackrel{ˇ}{𝐧}^t]=t_s^2\mathrm{𝐖𝐑}_n𝐖^t.$$
(31)
A second expression is obtained by writing each terms of the considered matrix in the Fourier domain. The component located at the $`j`$th row and $`k`$th column may be expressed as

$$𝖨𝖤[\stackrel{ˇ}{n}[j]\stackrel{ˇ}{n}[k]]=_{f_s/2}^{+f_s/2}\frac{𝖨𝖤[N(f)\overline{N(f^{})}]}{\sqrt{\mathrm{\Gamma }_n(f)\mathrm{\Gamma }_n(f^{})}}e^{2\pi i(jfkf^{})/f_s}𝑑f𝑑f^{}.$$
(32)
The integrand can be work out applying the Wiener-Khinchine theorem  :

$$𝖨𝖤[N(f)\overline{N(f^{})}]=\delta (ff^{})\mathrm{\Gamma }_n(f)$$
(33)
where $`\delta (f)=t_s_{k=\mathrm{}}^+\mathrm{}e^{2\pi ikf/f_s}`$.

The function $`\delta `$ acts on the elements of the set $`𝒫(f_s)`$ (containing the periodic functions of period $`f_s`$) the same way the Dirac operator acts on functions of $`L^2()`$. Let $`\mathrm{\Phi }`$ be a test function of $`𝒫(f_s)`$ and $`\varphi `$ its inverse Fourier transform, we have

$$_{f_s/2}^{+f_s/2}\mathrm{\Phi }(f)\delta (f)𝑑f=t_s\underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}\varphi [k]=\mathrm{\Phi }(0).$$
(34)
Note also that $`\delta `$ is $`f_s`$–periodic i.e., $`\delta (f+f_s)=\delta (f)`$, for all $`f`$. The proofs of these two properties of $`\delta `$ are left to the reader.

With (33), we have

$$𝖨𝖤[\stackrel{ˇ}{n}[j]\stackrel{ˇ}{n}[k]]=_{f_s/2}^{+f_s/2}\frac{e^{2\pi ijf/f_s}}{\sqrt{\mathrm{\Gamma }_n(f)}}\left(_{f_s/2}^{+f_s/2}\sqrt{\mathrm{\Gamma }_n(f^{})}\delta (ff^{})e^{2\pi ikf^{}/f_s}𝑑f^{}\right)𝑑f,$$
(35)
which simplifies with the change of variables $`u=ff^{}`$ and using the periodicity of the functions $`\delta `$ and $`\mathrm{\Gamma }_n`$,

$$𝖨𝖤[\stackrel{ˇ}{n}[j]\stackrel{ˇ}{n}[k]]=_{f_s/2}^{+f_s/2}\frac{e^{2\pi i(jk)f/f_s}}{\sqrt{\mathrm{\Gamma }_n(f)}}\left(_{f_s/2}^{+f_s/2}\sqrt{\mathrm{\Gamma }_n(fu)}\delta (u)e^{2\pi iku/f_s}𝑑u\right)𝑑f,$$
(36)
leading finally using (34) to

$$𝖨𝖤[\stackrel{ˇ}{n}[j]\stackrel{ˇ}{n}[k]]=_{f_s/2}^{+f_s/2}e^{2\pi i(jk)f/f_s}𝑑f,$$
(37)
or equivalently to

$$𝖨𝖤[\stackrel{ˇ}{n}[j]\stackrel{ˇ}{n}[k]]=f_s\delta _{jk}$$
(38)
where $`\delta _{jk}`$ is the Kronecker symbol (by definition, $`\delta _{jk}=1`$ if $`j=k`$, 0 otherwise).

We conclude that

$$𝖨𝖤[\stackrel{ˇ}{𝐧}\stackrel{ˇ}{𝐧}^t]=f_s𝐈.$$
(39)
Assuming that $`𝐑_n`$ is invertible and combining both eqs. (31) and (39), we get the relation between the whitening operator and the noise correlation matrix, namely

$$𝐑_n^1=t_s^3𝐖^t𝐖.$$
(40)
With eq. (29) and the Plancherel formula in eq. (1)

$$𝐲^t𝐑_n^1𝐱=t_s\stackrel{ˇ}{𝐲}^t\stackrel{ˇ}{𝐱}=_{f_s/2}^{+f_s/2}\stackrel{ˇ}{X}(f)\overline{\stackrel{ˇ}{Y}(f)}𝑑f,$$
(41)
provided that $`\stackrel{ˇ}{y}[k]`$ has support in $`\{0\mathrm{}N1\}`$. Replacing $`\stackrel{ˇ}{Y}(f)`$ and $`\stackrel{ˇ}{X}(f)`$ by their definitions completes the proof. ∎

#### 4.3.2 Deterministic and random components

If we choose $`N`$ large enough so that no finite size effect appears (i.e., the supports of all required signals respect the condition of Prop. 2), we can rewrite the deterministic term of the detection statistic as

$$\mathrm{\Lambda }_𝐇^{det}(𝐱)=\left(_{f_s/2}^{+f_s/2}\frac{X(f)\overline{S_m(f)}}{\mathrm{\Gamma }_n(f)}𝑑f\right)^2.$$
(42)
From the following eigen expansion of the signal correlation matrix

$$𝐑_s=\underset{k=0}{\overset{N1}{}}\sigma _k𝐯_k𝐯_k^t,$$
(43)
the random component may be expressed as :

$$\mathrm{\Lambda }_𝐇^{rand}(𝐱)=\underset{k=0}{\overset{N1}{}}\sigma _k(𝐯_k^t𝐑_n^1𝐱)^2.$$
(44)
Similarly to the deterministic component, assuming again that $`N`$ is sufficiently large, it follows that :

$$\mathrm{\Lambda }_𝐇^{rand}(𝐱)=\underset{k=0}{\overset{N1}{}}\sigma _k\left(_{f_s/2}^{+f_s/2}\frac{X(f)\overline{V_k(f)}}{\mathrm{\Gamma }_n(f)}𝑑f\right)^2$$
(45)
Eqs. (42) and (45) show that the quadratic detector with optimal deflection $`\mathrm{\Lambda }_𝐇(𝐱)`$ is closely related to the well-known technique of matched filtering . The complete statistic can be equivalently implemented as a bank of $`N+1`$ matched filters (using the templates given by $`𝐬_m`$ and $`𝐯_k`$, $`k=0\mathrm{}N1`$) whose output energies are combined with a weighted sum. Being a covariance matrix, $`𝐑_s`$ is positive definite. All eigenvalues $`\sigma _k`$ are then real and positive numbers. In consequence, the weighs (equaled to the eigenvalues) put favor or attenuate the contribution of a corresponding term in the summation.

## 5 Application to the detection of gravitational wave transients

Because the physics driving supernovae explosions is highly non-linear, the expected gravitational radiation is difficult to obtain through analytical means. However, numerical simulations are available  and a catalog of the reference waveforms associated to typical situations is accessible on the Internet. The waveforms of this catalog, we refer to as DFM (i.e., the initial letters of authors’ names, Dimmelmeier, Font and Müller) present an intrinsic diversity which motivates to look at them as if they were produced by a single random mechanism.

Consequently, the detection problem we face is similar to the one exposed in eqs. (9) provided that the second order statistics of both signal $`s(t)`$ and noise $`n(t)`$ are known. Strictly speaking, the covariance matrix $`𝐂_s`$ of the signal is not available but if we assume that the DFM gravitational waveforms are noise-free and independent realizations of the random process $`s(t)`$, they can be used to get a sufficiently accurate estimate.

We can then apply the method proposed in Sect. 4.2 to optimally detect $`s(t)`$. From the signal covariance estimate and a realistic noise model, we can calculate the quadratic detector with best deflection.

### 5.1 Finding the best quadratic detector

From the database publicly available on the Internet , we have extracted the $`N_z=25`$ waveforms, which we have resampled at the constant rate of $`f_s=20`$ kHz. Actually three sets of waveforms can be used : the one drawn from a Newtonian simulation set up and the other from a fully relativistic code. We preferred to use the latter since the result is likely to be closer to reality. The waveforms are stored in column vectors $`\{𝐳_k\}_{k=0\mathrm{}N_z1}`$ of size $`N=2251`$ rows. This corresponds to a maximum burst duration of about $`112.5`$ ms.

For each signal, we fixed the time axis origin to be located at the minimum of the largest negative bump of the whiten signal (in most cases, this is synchronized with the core bounce ), precisely : $`t_0=\text{cst}=\mathrm{argmin}_j(\stackrel{ˇ}{z}_k[j]),k`$. Figure 1 presents three individuals of each type of supernovae taken from the processed DFM catalog .

Assuming that the DFM gravitational waveforms are noise-free and independent realizations of the random process $`s(t)`$, we use these waveforms to estimate the covariance matrix of $`s(t)`$. This is done with the following empirical unbiased estimator :

$$\widehat{𝐂}_s=\frac{1}{N_z}\underset{k=0}{\overset{N_z1}{}}𝐳_k𝐳_k^t.$$
(46)
For simplicity, we consider that the noise power spectrum is known a priori and is given by the expected sensitivity curve for the planned detectors. (Note that the noise correlation matrix may also be estimated from “noise only” data streams.) We restrict our study to the case of the Virgo detector using the noise model available at the following address . Extensions to other large scale interferometers are straightforward. We can get the noise correlation matrix from the power spectrum applying the inverse Fourier transformation. From the obtained values of $`\widehat{𝐂}_s`$ and $`𝐑_n`$, we deduce the kernel of the best quadratic detector as given in eq. (22).

This computation requires $`O(N^3)`$ operations to get the inverse of $`𝐑_n`$ and we need roughly the same number of operations to make the two matrix multiplications in eq. (22). A more computationally efficient algorithm may be used. With this aim in view, we process the whole catalog of waveforms with the following operation $`\stackrel{˘}{𝐳}_k𝐑_n^1𝐳_k`$. Roughly speaking, the multiplication by $`𝐑_n^1`$ is equivalent to whitening the signal twice as suggested by the factorization of $`𝐑_n^1`$ in eq. (40). The more precise relation $`\stackrel{˘}{Z}_k(f)=t_sZ_k(f)/\mathrm{\Gamma }_n(f)`$ can be stated by applying the result in Prop. 2 with $`y[j]=\delta _{jm}`$ for any $`m\{0\mathrm{}N1\}`$ and $`x[j]=z_k[j]`$. This double whitening operation amounts to filtering the signal in the frequency bandwidth of interest where the noise is low and removing the remaining part where the noise is large. From eq. (25) and (46), it is easily shown that the objective kernel can be computed directly using the modified catalog with the relation :

$$\widehat{𝐇}=\frac{1}{N_z}\underset{k=0}{\overset{N_z1}{}}\stackrel{˘}{𝐳}_k\stackrel{˘}{𝐳}_k^t.$$
(47)
Since the correlation matrix $`𝐑_n`$ is Toeplitz symmetric, the computation of $`\stackrel{˘}{𝐳}_k`$ is equivalent to solving a $`N_z\times N_z`$ Toeplitz linear system. This can be done efficiently with a variety of fast $`O(N^2)`$ algorithms. We selected and applied the Levinson algorithm .

The total gravitational energy radiated during the collapse varies according to the selected models. The peak amplitudes of the waveforms of the DFM catalog have values ranging in an interval as large as one order of magnitude. To ensure that all types of supernovae are treated equitably in the sum of (47), we scale all $`\stackrel{˘}{𝐳}_k`$ by dividing by the expected signal-to-noise ratio defined as :

$$\mathrm{SNR}_k\left(_{f_s/2}^{f_s/2}|Z_k(f)|^2/\mathrm{\Gamma }_n(f)𝑑f\right)^{1/2}.$$
(48)
A practical expression of this quantity can be obtained by first deducing from Prop. 2 that $`\mathrm{SNR}_k=(𝐳_k^t𝐑_n^1𝐳_k)^{1/2}`$ and using the definition of the whitening operator in eq. (29) and its relation to $`𝐑_n`$ in eq. (40) thus leading to $`\mathrm{SNR}_k=\stackrel{ˇ}{𝐳}_k_2/\sqrt{f_s}`$ where $`𝐱_2^2𝐱^t𝐱`$ defines the $`l_2`$ norm.

At this point, it is worth noting that eq. (47) implements a learning scheme which extracts systematically the necessary information from a (possibly large and heterogeneous) database of a reference waveforms in order to find the best detector among a common used class of possibilities.

### 5.2 Approximated detector

A close look to the detector kernel $`\widehat{𝐇}`$ indicates that it is degenerated (its rank is much smaller than $`N`$). There are two reasons for that : first, as a result of the linear combining of $`N_zN`$ rank-1 matrices (see eq. (47)), the rank of the kernel cannot exceed $`N_z`$. The second reason is the fact that the waveforms of the DFM catalog have common features in their shapes (e.g., fundamental oscillation frequency, time duration, …). This causes the matrices $`\stackrel{˘}{𝐳}_k\stackrel{˘}{𝐳}_k^t`$ to share some linear dependency.

Precisely, it means that the kernel may be decomposed along a small number of preferred directions of the measurement space. The most adequate basis to check this is formed by the generalized eigen-vectors of $`\widehat{𝐂}_s`$ and $`𝐑_n`$ as explained in the following Section. We show that the kernel degeneracy may be used to simplify the detector and reduce its computational complexity.

#### 5.2.1 Truncating to principal directions

The vector $`𝐮`$ and scalar $`\gamma `$ are respectively the generalized eigen-vector and value of $`\widehat{𝐂}_s`$ and $`𝐑_n`$ if the following equation is satisfied  :

$$\widehat{𝐂}_s𝐮=\gamma 𝐑_n𝐮.$$
(49)
Since $`𝐑_n`$ is a definite positive matrix, it can be decomposed using the Cholesky factorization  as the product of invertible and triangular matrices, namely $`𝐑_n=𝐓_n𝐓_n^t`$. Multiplying to the left both sides of (49) by $`𝐓_n^1`$, the generalized eigen problem above turns out to be equivalent to the standard one given by :

$$𝚪𝐯=\gamma 𝐯,$$
(50)
provided that $`𝚪=𝐓_n^1\widehat{𝐂}_s𝐓_n^t`$ and $`𝐯=𝐓_n^t𝐮`$. Consequently, the matrix $`𝚪`$ may be expanded along its eigen-directions $`\{𝐯_k\}_{k=0\mathrm{}N1}`$, namely :

$$𝚪=\underset{k=0}{\overset{N1}{}}\gamma _k𝐯_k𝐯_k^t.$$
(51)
Since we have $`\widehat{𝐇}=𝐓_n^1𝚪𝐓_n^t`$, the previous expansion yields the one of the detector kernel along the generalized eigen-basis defined in eq. (49)

$$\widehat{𝐇}=\underset{k=0}{\overset{N1}{}}\gamma _k𝐮_k𝐮_k^t.$$
(52)
Combining adequately a Cholesky and a Schur decomposition , we computed the solutions of eq. (49). The eigen-values are sorted in decreasing order $`\gamma _0>\gamma _1>\mathrm{}>\gamma _{N1}`$ and presented in Fig. 2. It appears clearly that the resulting spectrum is essentially dominated by a first few eigenvalues.

A consequence is that the sum in eq. (52) can be fairly approximated by the summation truncated to the first terms. Let $`n<N`$ be the truncation limit, we get the following kernel

$$\stackrel{~}{𝐇}_n\underset{k=0}{\overset{n1}{}}\gamma _k𝐮_k𝐮_k^t$$
(53)
which we use to compute the approximated detection statistic :

$$\mathrm{\Lambda }_{\stackrel{~}{𝐇}_n}(𝐱)=\underset{k=0}{\overset{n1}{}}\gamma _k(\widehat{𝐮}_k^t𝐱)^2\mathrm{\Lambda }_{\widehat{𝐇}}(𝐱).$$
(54)
The value of the truncation index $`n`$ is essentially related to the intrinsic complexity of the initial waveform database. In the case of interest, $`n`$ is much smaller than $`N`$ by several order of magnitude (a non-empirical choice of $`n`$ is described in the next section), and its value remains stable when $`N`$ increases. In consequence, the approximated statistic (54) is computed with $`O(N^2)`$ floating point operations versus a total cost of $`O(N^3)`$ in the non approximated case.

#### 5.2.2 Loss in deflection due to approximation

The truncation to a few eigen directions causes $`\mathrm{\Lambda }_{\stackrel{~}{𝐇}_n}(𝐱)`$ to be sub-optimal i.e., the resulting deflection is smaller than the one obtained with $`\mathrm{\Lambda }_{\stackrel{ˇ}{𝐇}_n}(𝐱)`$. Two interesting questions are (i) how much deflection do we lose? and (ii) can we adjust $`n`$ so that the loss is acceptable? To address these questions, it is convenient to define the loss in deflection : $`l_nd^2(\mathrm{\Lambda }_{\stackrel{~}{𝐇}_n})/d^2(\mathrm{\Lambda }_{\widehat{𝐇}})`$. This index whose values are between 0 and 1 measures the degree of “sub-optimality” of the truncated detector.

Replacing $`𝐀`$ in the expression of the deflection obtained in Lemma 3 with the truncated sum in eq. (53), and using the fact that $`\{𝐮_k\}_{k=0\mathrm{}N1}`$ form a basis which diagonalizes simultaneously $`𝐑_n`$ and $`\widehat{𝐂}_s`$ i.e., more precisely $`𝐮_k^t𝐑_n𝐮_j=\delta _{jk}`$ and $`𝐮_k^t\widehat{𝐂}_s𝐮_j=\gamma _k\delta _{jk}`$ (see  for details), a straightforward calculation leads to

$$d^2(\mathrm{\Lambda }_{\stackrel{~}{𝐇}_n})=\underset{k=0}{\overset{n1}{}}\gamma _k^2.$$
(55)
This results holds also for $`n=N`$ yielding the maximum value of the deflection, which we denote by $`d_{max}d^2(\mathrm{\Lambda }_{\widehat{𝐇}})=_{k=0}^{N1}\gamma _k^2`$. The loss in deflection can then be expressed as :

$$l_n=\underset{k=0}{\overset{n1}{}}\frac{\gamma _k^2}{d_{max}}$$
(56)
and is presented in Figure 2 (b). We conclude that, with $`n=2`$ i.e., keeping the first two leading eigen-directions, the truncated detector has a performance index of about 99% (1% from optimum). Figure 2 (c) details the waveforms of the two leading eigenvectors.

Provided that $`u_0[k]`$ and $`u_1[k]`$ have support in $`\{0\mathrm{}N1\}`$ (this is the case in the example presented here, see Fig. 2 (c)), we can apply Lemma 1 and get the truncated detector (54) expressed in the frequency domain :

$$\mathrm{\Lambda }_{\stackrel{~}{𝐇}}(𝐱)=\gamma _0f_s^2\left(_{f_s/2}^{f_s/2}X(f)\overline{U_0(f)}𝑑f\right)^2+\gamma _1f_s^2\left(_{f_s/2}^{f_s/2}X(f)\overline{U_1(f)}𝑑f\right)^2.$$
(57)
We conclude that the detection statistic is computed by first selecting the interesting frequency contents of the spectrum of the observed data with the two (bandpass, as shown by Figure 2 (d)) filters $`U_0(f)`$ and $`U_1(f)`$ and then combining the energy of the filter outputs with a weighted sum. The weight parameters can be interpreted as “confidence coefficients” in finding a (supernovae) signal in the corresponding frequency bands.

In other words, the proposed method extracts systematically from a database of reference signals the frequency bands which need to be considered in order to maximize the deflection.

#### 5.2.3 Detection threshold and false alarm probability

Under noise only $`(H_0)`$ assumption, the detector (54) is a finite sum of the squares of the random variables defined by $`n_k𝐮_k^t𝐧`$

$$\mathrm{\Lambda }_{\stackrel{~}{𝐇}}(𝐧)=\underset{k=0}{\overset{n1}{}}\gamma _kn_k^2.$$
(58)
These variables can be easily shown to be Gaussian and zero-mean. Furthermore, since $`\{𝐮_k\}_{k=0\mathrm{}N1}`$ diagonalizes the noise correlation $`𝐑_n`$, we have $`𝖨𝖤[n_jn_k]=\delta _{jk}`$, from which we conclude that $`\{n_k\}_{k=0\mathrm{}n1}`$ is a sequence of independent and identically distributed Gaussian variables of PDF $`𝒩(0,1)`$.

Let $`f(\lambda )`$ be the PDF of $`\mathrm{\Lambda }_{\stackrel{~}{𝐇}}(𝐧)`$ when there is only noise. In the case where $`n=2`$ eigen-vectors are sufficient to get a good approximation of the optimal detector, we have for $`\lambda >0`$

$$f(\lambda )=\frac{1}{2\sqrt{\gamma _0\gamma _1}}\mathrm{exp}\left(\frac{1}{4}\left(\frac{1}{\gamma _1}+\frac{1}{\gamma _0}\right)\lambda \right)I_0\left(\frac{1}{4}\left(\frac{1}{\gamma _1}\frac{1}{\gamma _0}\right)\lambda \right)$$
(59)
where $`I_0()`$ is the modified Bessel function of first kind  and $`f(\lambda )=0`$ if $`\lambda 0`$.

Integrating the PDF in eq. (59), we obtained the cumulative probability function $`F(\lambda )_0^\lambda f(\nu )𝑑\nu =(\mathrm{\Lambda }_{\stackrel{~}{𝐇}}(𝐱)<\lambda |H_0)`$. The threshold ensuring a given false alarm probability $`p_0`$ is thus given by $`\eta =F^1(1p_0)`$. Such function is presented in Fig. 3 in the range of useful values of $`p_0`$ and using the first two leading eigenvalues of $`\widehat{𝐇}`$, namely $`\gamma _00.627`$ and $`\gamma _10.181`$.

#### 5.2.4 Time running implementation

Until now, we have only considered the statistical test of the presence of a supernovae transient at a given time. The date of arrival of the gravitational wave being unknown, we must apply the detection procedure at any given time instants. To do this, we select the $`N`$ data samples starting from a given time index $`m`$, namely $`𝐱_m(x[m+k],k=0\mathrm{}N1)^t`$. We compute the detection statistic $`\mathrm{\Lambda }_𝐇(𝐱_m)=𝐱_m^t\mathrm{𝐇𝐱}_m`$ for increasing and equally spaced values of $`m=0,\delta _m,\mathrm{}`$. This is similar to select the data with a time sliding window.

Using the approximated statistic expressed as in (57) and noting that $`X_m(f)=X(f)e^{2\pi imf/f_s}`$, we get

$$\mathrm{\Lambda }_{\stackrel{~}{𝐇}}(𝐱_m)=\gamma _0f_s^2(y_0[m])^2+\gamma _1f_s^2(y_1[m])^2,$$
(60)
where $`y_{(0,1)}[m]=_{f_s/2}^{+f_s/2}X(f)\overline{U_{(0,1)}(f)}e^{2\pi imf/f_s}𝑑f`$ are obtained by passing the signal through a time-invariant linear filter. Assuming $`U_0(f)`$ and $`U_1(f)`$ are stored in memory, the computation of $`y_0[m]`$ and $`y_1[m]`$ can be efficiently computed with the FFT (and inverse) algorithm.

## 6 Relation to other detection techniques

We have shown in Sect. 4.3 that the quadratic detector with optimal deflection can be related to matched filtering. In fact, many of the methods for transient detection available in the literature, e.g.  belong to the class of quadratic detectors defined in Def. 1. The vector formalism used here constitutes a general framework in which all these methods can be reformulated and easily compared. Considering that the noise model remains the same than the one we used in the previous section, we get the shape of the kernel used by the two contributions described in  and  to which we limit the investigation. With this “back-engineering” approach, we can retrieve the a priori assumption on the signal covariance needed for the considered detector to have optimal deflection. We make this comparison by looking to the generalized eigen-basis of the obtained kernel.

### 6.1 Excess power statistic

The basic idea is to monitor the power into one (or several) given frequency band $`f_0\pm \delta f/2`$ (similarly to Sect. 4.1). Let $`X[j]_{k=0}^{N1}x[k]e^{2\pi ijk/N}`$ be the discrete Fourier transform (DFT). The excess power statistic presented in  reads :

$$\mathrm{\Lambda }_{eps}(𝐱)=\underset{j=j_{}}{\overset{j_+}{}}\frac{|X[j]|^2}{\mathrm{\Gamma }_n(jf_s/N)},$$
(61)
where the limit indices are defined as $`j_\pm =N(f_0\pm \delta f/2)/f_s`$.

The DFT can be re-expressed as a scalar product $`X[j]=\overline{𝐟}_j^t𝐱`$ with the Fourier exponentials $`𝐟_j(e^{2\pi ijk/N},k=0\mathrm{}N1)^t`$. It is straightforward to show that the above statistic is a quadratic detector as in Def. 1 with the kernel

$$𝐇_{eps}=\underset{j=j_{}}{\overset{j_+}{}}\frac{𝐟_j\overline{𝐟_j^t}}{\mathrm{\Gamma }_n(jf_s/N)}.$$
(62)
Roughly speaking, assuming the noise power spectral density is “flat” in the selected frequency bandwidth, $`\{𝐟_j\}_{j=j_{}\mathrm{}j_+}`$ diagonalizes $`𝐑_n`$. In this case, the kernel of the excess power statistic has a number of $`(j_+j_{}+1)=V/2`$ generalized eigenvectors. Fig. 4 presents the some of these eigenvectors and their Fourier transforms.

### 6.2 Linear fit filter (ALF)

The detection statistic $`\mathrm{\Lambda }_{alf}(𝐱)`$ is obtained from a local linear fit of the whiten signal $`\stackrel{ˇ}{x}`$. A mean square rule yields the two parameters of the fit :

$`a`$ $`={\displaystyle \frac{t\stackrel{ˇ}{x}t^2}{t^2t^2}}`$ (63)
$`b`$ $`=\stackrel{ˇ}{x}at`$ (64)
which are orthonormalized, squared and combined to get $`\mathrm{\Lambda }_{alf}(𝐱)`$.

It turns out to be convenient to set the time origin at the center of data chunk, which we assume to have an odd number $`N`$ of samples. We can do this with no lost of generality. In this set up, the fit parameters are given by the following scalar products $`a=𝐭^t\stackrel{ˇ}{𝐱}/𝐭_2^2`$ and $`b=\mathrm{𝟏}^t\stackrel{ˇ}{𝐱}`$ where we defined $`𝐭=(L\mathrm{}L)^t`$, $`L(N1)/2`$ being the half-size of the data chunk and $`\mathrm{𝟏}=(1\mathrm{}\mathrm{\hspace{0.17em}1})^t`$. After the orthonormalization and combining, the detection statistic appears to be a quadratic detector as in Def. 1 of kernel

$$𝐇_{alf}=t_s^2𝐖^t\left(\frac{\mathrm{𝐭𝐭}^t}{𝐭_2^2}+\frac{\mathrm{𝟏𝟏}^t}{\mathrm{𝟏}_2^2}\right)𝐖$$
(65)
where $`𝐖`$ is the whitening matrix defined in eq. (29).

It can be easily shown that $`𝐖^t\mathrm{𝟏}`$ and $`𝐖^t𝐭`$ are the two generalized eigenvectors of $`𝐇_{alf}`$ associated to the eigenvalue 1 (this is the only non-zero eigenvalue). They are presented in Fig. 4. In this degenerated case, similar calculations as the ones done in Sect. 5.2.3 yields the PDF of $`\mathrm{\Lambda }_{alf}(𝐱)`$ in the noise only case , namely :

$$f_{alf}(\lambda )=\frac{1}{2}e^{\lambda /2},$$
(66)
from which the threshold can be obtained for a given false alarm probability. This is a complementary contribution to the analysis made in  about the local linear fit method.

## 7 Conclusions

Quadratic detectors (i.e., statistics which are bilinear functions of the data) can be essentially viewed as a filtering of the data through a selection of frequency bands, the power of the filtered data being further linearly combined. We introduced a method which systematically extracts from a complicated and possibly large database of target signals, the important features which need to be considered in order to design this filter bank and choose the parameters for the energy combining. In the context of the detection of supernovae core collapses, we show that the method gives an intuitively appealing result of a filter bank composed of two elements (the one selecting the bounce pulse and the other, the few oscillations of the ringdown phase) whose output powers are combined in such a way to favor the bounce (which is the most energetic part of the signal).

The scope of the approach presented here is general. The algorithm can be adapted to other problems of the same type (for instance the detection of the final merger part of a binary black coalescence or of non-stationary noise interferences) provided that a sufficient number of training waveforms are available.