text:
# 1 Introduction

## 1 Introduction

From olfaction to vision and audition, a growing number of experiments are examining the responses of sensory neurons to natural stimuli (Creutzfeldt & Northdurft, 1978; Rieke et al., 1995; Baddeley et al., 1997; Stanley et al., 1999; Theunissen et al., 2000; Vinje & Gallant, 2000; Lewen et al., 2001; Sen et al., 2001; Vickers et al., 2001; Vinje & Gallant, 2002; Ringach et al., 2002; Weliky et al., 2003; Rolls et al., 2003; Smyth et al., 2003). Observing the full dynamic range of neural responses may require using stimulus ensembles which approximate those occurring in nature (Rieke et al., 1997; Simoncelli & Olshausen, 2001), and it is an attractive hypothesis that the neural representation of these natural signals may be optimized in some way (Barlow, 1961, 2001; von der Twer & Macleod, 2001; Bialek, 2002). Many neurons exhibit strongly nonlinear and adaptive responses that are unlikely to be predicted from a combination of responses to simple stimuli; for example neurons have been shown to adapt to the distribution of sensory inputs, so that any characterization of these responses will depend on context (Smirnakis et al., 1996; Brenner et al., 2000a; Fairhall et al., 2001). Finally, the variability of a neuron’s responses decreases substantially when complex dynamical, rather than static, stimuli are used (Mainen & Sejnowski, 1995; de Ruyter van Steveninck et al., 1997; Kara et al., 2000; de Ruyter van Steveninck et al., 2001). All of these arguments point to the need for general tools to analyze the neural responses to complex, naturalistic inputs.

The stimuli analyzed by sensory neurons are intrinsically high dimensional, with dimensions $`D10^210^3`$. For example, in the case of visual neurons, the input is commonly specified as light intensity on a grid of at least $`10\times 10`$ pixels. Each of the presented stimuli can be described as a vector $`𝐬`$ in this high dimensional stimulus space, see Fig. 1. The dimensionality becomes even larger if stimulus history has to be considered as well. For example, if we are interested in how the past $`N`$ frames of the movie affect the probability of a spike, then the stimulus $`𝐬`$, being a concatenation of the past $`N`$ samples, will have dimensionality $`N`$ times that of a single frame. We also assume that the probability distribution $`P(𝐬)`$ is sampled during an experiment ergodically, so that we can exchange averages over time with averages over the true distribution as needed.

Even though direct exploration of a $`D10^210^3`$ dimensional stimulus space is beyond the constraints of experimental data collection, progress can be made provided we make certain assumptions about how the response has been generated. In the simplest model, the probability of response can be described by one receptive field (RF) or linear filter (Rieke et al., 1997). The receptive field can be thought of as a template or special direction $`\widehat{e}_1`$ in the stimulus space111The notation $`\widehat{e}`$ denotes a unit vector, since we are interested only in the direction the vector specifies and not in its length. such that the neuron’s response depends only on a projection of a given stimulus $`𝐬`$ onto $`\widehat{e}_1`$, although the dependence of the response on this projection can be strongly nonlinear, cf. Fig. 1. In this simple model, the reverse correlation method (de Boer & Kuyper, 1968; Rieke et al., 1997; Chichilnisky, 2001) can be used to recover the vector $`\widehat{e}_1`$ by analyzing the neuron’s responses to Gaussian white noise. In a more general case, the probability of the response depends on projections $`s_i=\widehat{e}_i𝐬`$ of the stimulus $`𝐬`$ on a set of $`K`$ vectors $`\{\widehat{e}_1,\widehat{e}_2,\mathrm{},\widehat{e}_K\}`$:

$$P(\mathrm{spike}|𝐬)=P(\mathrm{spike})f(s_1,s_2,\mathrm{},s_K),$$
(1)
where $`P(\mathrm{spike}|𝐬)`$ is the probability of a spike given a stimulus $`𝐬`$ and $`P(\mathrm{spike})`$ is the average firing rate. In what follows we will call the subspace spanned by the set of vectors $`\{\widehat{e}_1,\widehat{e}_2,\mathrm{},\widehat{e}_K\}`$ the relevant subspace (RS)222Since the analysis does not depend on a particular choice of a basis within the full D-dimensional stimulus space, for clarity we choose the basis in which the first $`K`$ basis vectors span the relevant subspace and the remaining $`DK`$ vectors span the irrelevant subspace.. We reiterate that vectors $`\{\widehat{e}_i\}`$, $`1iK`$ may also describe how the time dependence of the stimulus $`𝐬`$ affects the probability of a spike. An example of such a relevant dimension would be a spatiotemporal receptive field of a visual neuron. Even though the ideas developed below can be used to analyze input–output functions $`f`$ with respect to different neural responses, such as patterns of spikes in time (de Ruyter van Steveninck & Bialek, 1988; Brenner et al., 2000b; Reinagel & Reid, 2000), for illustration purposes we choose a single spike as the response of interest.333 We emphasize that our focus here on single spikes is not equivalent to assuming that the spike train is a Poisson process modulated by the stimulus. No matter what the statistical structure of the spike train is we always can ask what features of the stimulus are relevant for setting the probability of generating a single spike at one moment in time. From an information theoretic point of view, asking for stimulus features that capture the mutual information between the stimulus and the arrival times of single spikes is a well posed question even if successive spikes do not carry independent information; note also that spikes carrying independent information is not the same as spikes being generated as a Poisson process. On the other hand, if (for example) different temporal patterns of spikes carry information about different stimulus features, then analysis of single spikes will result in a relevant subspace of artefactually high dimensionality. Thus it is important that the approach discussed here carries over without modification to the analysis of relevant dimensions for the generation of any discrete event, such as a pattern of spikes across time in one cell or synchronous spikes across a population of cells. For a related discussion of relevant dimensions and spike patterns using covariance matrix methods see (de Ruyter van Steveninck & Bialek, 1988; Agüera y Arcas et al., 2003).

Equation (1) in itself is not yet a simplification if the dimensionality $`K`$ of the RS is equal to the dimensionality $`D`$ of the stimulus space. In this paper we will assume that the neuron’s firing is sensitive only to a small number of stimulus features, i.e. $`KD`$. While the general idea of searching for low dimensional structure in high dimensional data is very old, our motivation here comes from work on the fly visual system where it was shown explicitly that patterns of action potentials in identified motion sensitive neurons are correlated with low dimensional projections of the high dimensional visual input (de Ruyter van Steveninck & Bialek, 1988; Brenner et al., 2000a; Bialek & de Ruyter van Steveninck, 2003). The input–output function $`f`$ in Eq. (1) can be strongly nonlinear, but it is presumed to depend only on a small number of projections. This assumption appears to be less stringent than that of approximate linearity which one makes when characterizing neuron’s response in terms of Wiener kernels \[see for example the discussion in Section 2.1.3 of Rieke et al. (1997)\]. The most difficult part in reconstructing the input–output function is to find the RS. Note that for $`K>1`$, a description in terms of any linear combination of vectors $`\{\widehat{e}_1,\widehat{e}_2,\mathrm{},\widehat{e}_K\}`$ is just as valid, since we did not make any assumptions as to a particular form of the nonlinear function $`f`$.

Once the relevant subspace is known, the probability $`P(\mathrm{spike}|𝐬)`$ becomes a function of only a few parameters, and it becomes feasible to map this function experimentally, inverting the probability distributions according to Bayes’ rule:

$$f(s_1,s_2,\mathrm{},s_K)=\frac{P(s_1,s_2,\mathrm{},s_K|\mathrm{spike})}{P(s_1,s_2,\mathrm{},s_K)}.$$
(2)
If stimuli are chosen from a correlated Gaussian noise ensemble, then the neural response can be characterized by the spike–triggered covariance method (de Ruyter van Steveninck & Bialek, 1988; Brenner et al., 2000a; Schwartz et al., 2002; Touryan et al., 2002; Bialek & de Ruyter van Steveninck, 2003). It can be shown that the dimensionality of the RS is equal to the number of nonzero eigenvalues of a matrix given by a difference between covariance matrices of all presented stimuli and stimuli conditional on a spike. Moreover, the RS is spanned by the eigenvectors associated with the nonzero eigenvalues multiplied by the inverse of the a priori covariance matrix. Compared to the reverse correlation method, we are no longer limited to finding only one of the relevant dimensions $`\{\widehat{e}_i\}`$, $`1iK`$. Both the reverse correlation and the spike–triggered covariance method, however, give rigorously interpretable results only for Gaussian distributions of inputs.

In this paper we investigate whether it is possible to lift the requirement for stimuli to be Gaussian. When using natural stimuli, which certainly are non–Gaussian, the RS cannot be found by the spike–triggered covariance method. Similarly, the reverse correlation method does not give the correct RF, even in the simplest case where the input–output function in Eq. (1) depends only on one projection; see Appendix A for a discussion of this point. However, vectors that span the RS are clearly special directions in the stimulus space independent of assumptions about $`P(𝐬)`$. This notion can be quantified by the Shannon information. We note that the stimuli $`𝐬`$ do not have to lie on a low-dimensional manifold within the overall $`D`$ dimensional space.444 If one suspects that neurons are sensitive to low dimensional features of their input, one might be tempted to analyze neural responses to stimuli that explore only the (putative) relevant subspace, perhaps along the line of the subspace reverse correlation method (Ringach et al., 1997). Our approach (like the spike–triggered covariance approach) is different because it allows the analysis of responses to stimuli that live in the full space, and instead we let the neuron “tell us” which low dimensional subspace is relevant. However, since we assume that the neuron’s input-output function depends on a small number of relevant dimensions, the ensemble of stimuli conditional on a spike may exhibit clear clustering. This makes the proposed method of looking for the RS complimentary to the clustering of stimuli conditional on a spike done in the information bottleneck method (Tishby et al., 1999); see also (Dimitrov & Miller, 2001). Non–information based measures of similarity between probability distributions $`P(𝐬)`$ and $`P(𝐬|\mathrm{spike})`$ have also been proposed to find the RS (Paninski, 2003a).

To summarize our assumptions:

* The sampling of the probability distribution of stimuli $`P(𝐬)`$ is ergodic and stationary across repetitions. The probability distribution is not assumed to be Gaussian. The ensemble of stimuli described by $`P(𝐬)`$ does not have to lie on a low-dimensional manifold embedded in the overall $`D`$-dimensional space.
* We choose a single spike as the response of interest (for illustration purposes only). An identical scheme can be applied, for example, to particular interspike intervals or to synchronous spikes from a pair of neurons.
* The subspace relevant for generating a spike is low dimensional and Euclidean, cf. Eq. 1.
* The input–output function, which is defined within the low dimensional RS, can be arbitrarily nonlinear. It is obtained experimentally by sampling the probability distributions $`P(𝐬)`$ and $`P(𝐬|\mathrm{spike})`$ within the RS.

The paper is organized as follows. In Sec. 2 we discuss how an optimization problem can be formulated to find the RS. A particular algorithm used to implement the optimization scheme is described in Sec. 3. In Sec. 4 we illustrate how the optimization scheme works with natural stimuli for model orientation sensitive cells with one and two relevant dimensions, much like simple and complex cells found in primary visual cortex, as well as for a model auditory neuron responding to natural sounds. We also discuss the convergence of our estimates of the RS as a function of data set size. We emphasize that our optimization scheme does not rely on any specific statistical properties of the stimulus ensemble, and thus can be used with natural stimuli.

## 2 Information as an objective function

When analyzing neural responses, we compare the a priori probability distribution of all presented stimuli with the probability distribution of stimuli which lead to a spike (de Ruyter van Steveninck & Bialek, 1988). For Gaussian signals, the probability distribution can be characterized by its second moment, the covariance matrix. However, an ensemble of natural stimuli is not Gaussian, so that in a general case neither second nor any finite number of moments is sufficient to describe the probability distribution. In this situation, Shannon information provides a rigorous way of comparing two probability distributions. The average information carried by the arrival time of one spike is given by (Brenner et al., 2000b)

$$I_{\mathrm{spike}}=𝑑𝐬P(𝐬|\mathrm{spike})\mathrm{log}_2\left[\frac{P(𝐬|\mathrm{spike})}{P(𝐬)}\right],$$
(3)
where $`d𝐬`$ denotes integration over full $`D`$–dimensional stimulus space. The information per spike as written in (3) is difficult to estimate experimentally, since it requires either sampling of the high–dimensional probability distribution $`P(𝐬|\mathrm{spike})`$ or a model of how spikes were generated, i.e. the knowledge of low–dimensional RS. However it is possible to calculate $`I_{\mathrm{spike}}`$ in a model–independent way, if stimuli are presented multiple times to estimate the probability distribution $`P(\mathrm{spike}|𝐬)`$. Then,

$$I_{\mathrm{spike}}=\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}\mathrm{log}_2\left[\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}\right]_𝐬,$$
(4)
where the average is taken over all presented stimuli. This can be useful in practice (Brenner et al., 2000b), because we can replace the ensemble average $`_𝐬`$ with a time average, and $`P(\mathrm{spike}|𝐬)`$ with the time dependent spike rate $`r(t)`$. Note that for a finite dataset of $`N`$ repetitions, the obtained value $`I_{\mathrm{spike}}(N)`$ will be on average larger than $`I_{\mathrm{spike}}(\mathrm{})`$. The true value $`I_{\mathrm{spike}}`$ can be found either by subtracting an expected bias value, which is of the order of $`1/(P(\mathrm{spike})N\mathrm{\hspace{0.17em}2}\mathrm{ln}2)`$ (Treves & Panzeri, 1995; Panzeri & Treves, 1996; Pola et al., 2002; Paninski, 2003b), or by extrapolating to $`N\mathrm{}`$ (Brenner et al., 2000b; Strong et al., 1998). Measurement of $`I_{\mathrm{spike}}`$ in this way provides a model independent benchmark against which we can compare any description of the neuron’s input–output relation.

Our assumption is that spikes are generated according to a projection onto a low dimensional subspace. Therefore to characterize relevance of a particular direction $`𝐯`$ in the stimulus space, we project all of the presented stimuli onto $`𝐯`$ and form probability distributions $`P_𝐯(x)`$ and $`P_𝐯(x|\mathrm{spike})`$ of projection values $`x`$ for the a priori stimulus ensemble and that conditional on a spike, respectively:

$`P_𝐯(x)`$ $`=`$ $`\delta (x𝐬𝐯)_𝐬,`$ (5)
$`P_𝐯(x|\mathrm{spike})`$ $`=`$ $`\delta (x𝐬𝐯)|\mathrm{spike}_𝐬,`$ (6)
where $`\delta (x)`$ is a delta-function. In practice, both the average $`\mathrm{}_𝐬𝑑𝐬\mathrm{}P(𝐬)`$ over the a priori stimulus ensemble, and the average $`\mathrm{}|\mathrm{spike}_𝐬𝑑𝐬\mathrm{}P(𝐬|\mathrm{spike})`$ over the ensemble conditional on a spike, are calculated by binning the range of projections values $`x`$. The probability distributions are then obtained as histograms, normalized in a such a way that the sum over all bins gives 1. The mutual information between spike arrival times and the projection $`x`$, by analogy with Eq. (3), is

$$I(𝐯)=𝑑xP_𝐯(x|\mathrm{spike})\mathrm{log}_2\left[\frac{P_𝐯(x|\mathrm{spike})}{P_𝐯(x)}\right],$$
(7)
which is also the Kullback-Leibler divergence $`D[P_𝐯(x|\mathrm{spike})||P_𝐯(x)]`$; notice that this information is a function of the direction $`𝐯`$. The information $`I(𝐯)`$ provides an invariant measure of how much the occurrence of a spike is determined by projection on the direction $`𝐯`$. It is a function only of direction in the stimulus space and does not change when vector $`𝐯`$ is multiplied by a constant. This can be seen by noting that for any probability distribution and any constant $`c`$, $`P_{c𝐯}(x)=c^1P_𝐯(x/c)`$ \[see also Theorem 9.6.4 of Cover & Thomas (1991)\]. When evaluated along any vector $`𝐯`$, $`I(𝐯)I_{\mathrm{spike}}`$. The total information $`I_{\mathrm{spike}}`$ can be recovered along one particular direction only if $`𝐯=\widehat{e}_1`$, and only if the RS is one dimensional.

By analogy with (7), one could also calculate information $`I(𝐯_1,\mathrm{},𝐯_n)`$ along a set of several directions $`\{𝐯_1,\mathrm{},𝐯_n\}`$ based on the multi-point probability distributions of projection values $`x_1`$, $`x_2`$, … $`x_n`$ along vectors $`𝐯_1`$, $`𝐯_2`$, … $`𝐯_n`$ of interest:

$`P_{𝐯_1,\mathrm{},𝐯_n}(\{x_i\}|\mathrm{spike})`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}\delta (x_i𝐬𝐯_i)|\mathrm{spike}_𝐬,`$ (8)
$`P_{𝐯_1,\mathrm{},𝐯_n}(\{x_i\})`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}\delta (x_i𝐬𝐯_i)_𝐬.`$ (9)
If we are successful in finding all of the directions $`\{\widehat{e}_i\}`$, $`1iK`$ contributing to the input–output relation (1), then the information evaluated in this subspace will be equal to the total information $`I_{\mathrm{spike}}`$. When we calculate information along a set of $`K`$ vectors that are slightly off from the RS, the answer, of course, is smaller than $`I_{\mathrm{spike}}`$ and is initially quadratic in small deviations $`\delta 𝐯_i`$. One can therefore hope to find the RS by maximizing information with respect to $`K`$ vectors simultaneously. The information does not increase if more vectors outside the RS are included. For uncorrelated stimuli, any vector or a set of vectors that maximizes $`I(𝐯)`$ belongs to the RS. On the other hand, as discussed in Appendix B, the result of optimization with respect to a number of vectors $`k<K`$ may deviate from the RS if stimuli are correlated. To find the RS, we first maximize $`I(𝐯)`$, and compare this maximum with $`I_{\mathrm{spike}}`$, which is estimated according to (4). If the difference exceeds that expected from finite sampling corrections, we increment the number of directions with respect to which information is simultaneously maximized.

## 3 Optimization algorithm

In this section we describe a particular algorithm we used to look for the most informative dimensions in order to find the relevant subspace. We make no claim that our choice of the algorithm is most efficient. However, it does give reproducible results for different starting points and spike trains with differences taken to simulate neural noise. Overall, choices for an algorithm are broader because the information $`I(𝐯)`$ as defined by (7) is a continuous function, whose gradient can be computed. We find (see Appendix C for a derivation)

$$_𝐯I=𝑑xP_𝐯(x)\left[𝐬|x,\mathrm{spike}𝐬|x\right]\left[\frac{d}{dx}\frac{P_𝐯(x|\mathrm{spike})}{P_𝐯(x)}\right],$$
(10)
where

$$𝐬|x,\mathrm{spike}=\frac{1}{P(x|\mathrm{spike})}𝑑𝐬𝐬\delta (x𝐬𝐯)P(𝐬|\mathrm{spike}),$$
(11)
and similarly for $`𝐬|x`$. Since information does not change with the length of the vector, we have $`𝐯_𝐯I=0`$, as also can be seen directly from Eq. (10).

As an optimization algorithm, we have used a combination of gradient ascent and simulated annealing algorithms: successive line maximizations were done along the direction of the gradient (Press et al., 1992). During line maximizations, a point with a smaller value of information was accepted according to Boltzmann statistics, with probability $`\mathrm{exp}[(I(𝐯_{i+1})I(𝐯_i))/T]`$. The effective temperature $`T`$ is reduced by factor of $`1ϵ_T`$ upon completion of each line maximization. Parameters of the simulated annealing algorithm to be adjusted are the starting temperature $`T_0`$ and the cooling rate $`ϵ_T`$, $`\mathrm{\Delta }T=ϵ_TT`$. When maximizing with respect to one vector we used values $`T_0=1`$ and $`ϵ_T=0.05`$. When maximizing with respect to two vectors, we either used the cooling schedule with $`ϵ_T=0.005`$ and repeated it several times (4 times in our case) or allowed the effective temperature $`T`$ to increase by a factor of 10 upon convergence to a local maximum (keeping $`TT_0`$ always), while limiting the total number of line maximizations.

The problem of maximizing a function often is related to the problem of making a good initial guess. It turns out, however, that the choice of a starting point is much less crucial in cases where the stimuli are correlated. To illustrate this point we plot in Fig. 2 the probability distribution of information along random directions $`𝐯`$ both for white noise and for naturalistic stimuli in a model with one relevant dimension. For uncorrelated stimuli, not only is information equal to zero for a vector that is perpendicular to the relevant subspace, but in addition the derivative is equal to zero. Since a randomly chosen vector has on average a small projection on the relevant subspace (compared to its length) $`v_r/|𝐯|\sqrt{n/d}`$, the corresponding information can be found by expanding in $`v_r/|𝐯|`$:

$$I\frac{v_r^2}{2|𝐯|^2}𝑑xP_{\widehat{e}_{ir}}(x)\left(\frac{P_{\widehat{e}_{ir}}^{}(x)}{P_{\widehat{e}_{ir}}(x)}\right)^2[𝐬\widehat{e}_r|\mathrm{spike}𝐬\widehat{e}_r]^2$$
(12)
where vector $`v=v_r\widehat{e}_r+v_{ir}\widehat{e}_{ir}`$ is decomposed in its components inside and outside the RS, respectively. The average information for a random vector is, therefore, $`(v_r^2/|𝐯|^2)=K/D`$.

In cases where stimuli are drawn from a Gaussian ensemble with correlations, an expression for the information values has a similar structure to (12). To see this, we transform to Fourier space and normalize each Fourier component by the square root of the power spectrum $`S(𝐤)`$. In this new basis, both the vectors $`\{e_i\}`$, $`1iK`$, forming the RS and the randomly chosen vector $`𝐯`$ along which information is being evaluated are to be multiplied by $`\sqrt{S(𝐤)}`$. Thus, if we now substitute for the dot product $`v_r^2`$ the convolution weighted by the power spectrum, $`_i^K(𝐯\widehat{e}_i)^2`$, where

$$𝐯\widehat{e}_i=\frac{_𝐤v(𝐤)\widehat{e}_i(𝐤)S(𝐤)}{\sqrt{_𝐤v^2(𝐤)S(𝐤)}\sqrt{_𝐤\widehat{e}_i^2(𝐤)S(𝐤)}},$$
(13)
then Eq. (12) will describe information values along randomly chosen vectors $`𝐯`$ for correlated Gaussian stimuli with the power spectrum $`S(𝐤)`$. Even though both $`v_r`$ and $`v(𝐤)`$ are Gaussian variables with variance $`1/D`$, the weighted convolution has not only a much larger variance but is also strongly non–Gaussian \[the non–Gaussian character is due to the normalizing factor $`_𝐤v^2(𝐤)S(𝐤)`$ in the denominator of Eq. (13)\]. As for the variance, it can be estimated as $`<(𝐯\widehat{e}_1)^2>=4\pi /\mathrm{ln}^2D`$, in cases where stimuli are taken as patches of correlated Gaussian noise with the two-dimensional power spectrum $`S(𝐤)=A/k^2`$. The large values of the weighted dot product $`𝐯\widehat{e}_i`$, $`1iK`$ result not only in significant information values along a randomly chosen vector, but also in large magnitudes of the derivative $`I`$, which is no longer dominated by noise, contrary to the case of uncorrelated stimuli. In the end, we find that randomly choosing one of the presented frames as a starting guess is sufficient.

## 4 Results

We tested the scheme of looking for the most informative dimensions on model neurons that respond to stimuli derived from natural scenes and sounds. As visual stimuli we used scans across natural scenes, which were taken as black and white photos digitized to 8 bits with no corrections made for the camera’s light intensity transformation function. Some statistical properties of the stimulus set are shown in Fig. 3. Qualitatively, they reproduce the known results on the statistics of natural scenes (Ruderman & Bialek, 1994; Ruderman, 1994; Dong & Atick, 1995; Simoncelli & Olshausen, 2001). Most important properties for this study are strong spatial correlations, as evident from the power spectrum S(k) plotted in panel (b), and deviations of the probability distribution from a Gaussian one. The non–Gaussian character can be seen in panel (c), where the

probability distribution of intensities is shown, and in panel (d) which shows the distribution of projections on a Gabor filter \[in what follows the units of projections, such as $`s_1`$, will be given in units of the corresponding standard deviations\]. Our goal is to demonstrate that even though the correlations present in the ensemble are non–Gaussian, they can be removed successfully from the estimate of vectors defining the RS.

### 4.1 A model simple cell

Our first example is based on the properties of simple cells found in the primary visual cortex. A model phase and orientation sensitive cell has a single relevant dimension $`\widehat{e}_1`$ shown in Fig. 4(a). A given stimulus $`𝐬`$ leads to a spike if the projection $`s_1=𝐬\widehat{e}_1`$ reaches a threshold value $`\theta `$ in the presence of noise:

$$\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}f(s_1)=H(s_1\theta +\xi ),$$
(14)
where a Gaussian random variable $`\xi `$ of variance $`\sigma ^2`$ models additive noise, and the function $`H(x)=1`$ for $`x>0`$, and zero otherwise. Together with the RF $`\widehat{e}_1`$, the parameters $`\theta `$ for threshold and the noise variance $`\sigma ^2`$ determine the input–output function.

When the spike–triggered average (STA), or reverse correlation function, is computed from the responses to correlated stimuli, the resulting vector will be broadened due to spatial correlations present in the stimuli (see Fig. 4b). For stimuli that are drawn from a Gaussian probability distribution, the effects of correlations could be removed by multiplying $`𝐯_{\mathrm{sta}}`$ by the inverse of the a priori covariance matrix, according to the reverse correlation method, $`\widehat{v}_{Gaussianest}C_{apriori}^1𝐯_{\mathrm{sta}}`$, Eq. (20). However this procedure tends to amplify noise. To separate errors due to neural noise from those due to the non–Gaussian character of correlations, note that in a model the effect of neural noise on our estimate of the STA can be eliminated by averaging the presented stimuli weighted with the exact firing rate, as opposed to using a histogram of responses to estimate $`P(\mathrm{spike}|𝐬)`$ from a finite set of trials. We have used this “exact” STA,

$$𝐯_{\mathrm{sta}}=𝑑𝐬𝐬P(𝐬|\mathrm{spike})=\frac{1}{P(\mathrm{spike})}𝑑𝐬P(𝐬)𝐬P(\mathrm{spike}|𝐬),$$
(15)
in calculations presented in Figs. 4 (b) and (c). Even with this noiseless STA (the equivalent of collecting an infinite data set), the standard decorrelation procedure is not valid for non–Gaussian stimuli and nonlinear input–output functions, as discussed in detail in Appendix A. The result of such a decorrelation in our example is shown in Fig. 4(c). It clearly is missing some of the structure in the model filter, with projection $`\widehat{e}_1\widehat{v}_{Gaussianest}0.14`$. The discrepancy is not due to neural noise or finite sampling, since the “exact” STA was decorrelated; the absence of noise in the exact STA also means that there would be no justification for smoothing the results of the decorrelation. The discrepancy between the true receptive field and the decorrelated STA increases with the strength of nonlinearity in the input–output function.

In contrast, it is possible to obtain a good estimate of the relevant dimension $`\widehat{e}_1`$ by maximizing information directly, see panel (d). A typical progress of the simulated annealing algorithm with decreasing temperature $`T`$ is shown in Fig. 4(e). There we plot both the information along the vector, and its projection on $`\widehat{e}_1`$. We note that while information $`I`$ remains almost constant, the value of projection continues to improve. Qualitatively this is because the probability distributions depend exponentially on information. The final value of projection depends on the size of the data set, as discussed below. In the example shown in Fig. 4 there were $`50,000`$ spikes with average probability of spike $`0.05`$ per frame, and the reconstructed vector has a projection $`\widehat{v}_{max}\widehat{e}_1=0.920\pm 0.006`$. Having estimated the RF, one can proceed to sample the nonlinear input-output function. This is done by constructing histograms for $`P(𝐬\widehat{v}_{\mathrm{max}})`$ and $`P(𝐬\widehat{v}_{\mathrm{max}}|\mathrm{spike})`$ of projections onto vector $`\widehat{v}_{\mathrm{max}}`$ found by maximizing information, and taking their ratio, as in Eq. (2). In Fig. 4(f) we compare $`P(\mathrm{spike}|𝐬\widehat{v}_{\mathrm{max}})`$ (crosses) with the probability $`P(\mathrm{spike}|s_1)`$ used in the model (solid line).

### 4.2 Estimated deviation from the optimal dimension

When information is calculated from a finite data set, the (normalized) vector $`\widehat{v}`$ which maximizes $`I`$ will deviate from the true RF $`\widehat{e}_1`$. The deviation $`\delta 𝐯=\widehat{v}\widehat{e}_1`$ arises because the probability distributions are estimated from experimental histograms and differ from the distributions found in the limit of infinite data size. For a simple cell, the quality of reconstruction can be characterized by the projection $`\widehat{v}\widehat{e}_1=1\frac{1}{2}\delta 𝐯^2`$, where both $`\widehat{v}`$ and $`\widehat{e}_1`$ are normalized, and $`\delta 𝐯`$ is by definition orthogonal to $`\widehat{e}_1`$. The deviation $`\delta 𝐯A^1I`$, where $`A`$ is the Hessian of information. Its structure is similar to that of a covariance matrix:

$$A_{ij}=\frac{1}{\mathrm{ln}2}𝑑xP(x|\mathrm{spike})\left(\frac{d}{dx}\mathrm{ln}\frac{P(x|\mathrm{spike})}{P(x)}\right)^2(s_is_j|xs_i|xs_j|x).$$
(16)
When averaged over possible outcomes of $`N`$ trials, the gradient of information is zero for the optimal direction. Here in order to evaluate $`\delta 𝐯^2=\mathrm{Tr}[A^1II^TA^1]`$, we need to know the variance of the gradient of $`I`$. Assuming that the probability of generating a spike is independent for different bins, we can estimate $`I_iI_jA_{ij}/(N_{\mathrm{spike}}\mathrm{ln}2)`$. Therefore an expected error in the reconstruction of the optimal filter is inversely proportional to the number of spikes. The corresponding expected value of the projection between the reconstructed vector and the relevant direction $`\widehat{e}_1`$ is given by:

$$\widehat{v}\widehat{e}_11\frac{1}{2}\delta 𝐯^2=1\frac{\mathrm{Tr}^{}[A^1]}{2N_{\mathrm{spike}}\mathrm{ln}2},$$
(17)
where $`\mathrm{Tr}^{}`$ means that the trace is taken in the subspace orthogonal to the model filter555 By definition $`\delta v_1=\delta 𝐯\widehat{e}_1=0`$, and therefore $`\delta v_1^2A_{11}^1`$ is to be subtracted from $`\delta 𝐯^2\mathrm{Tr}[A^1]`$. Because $`\widehat{e}_1`$ is an eigenvector of $`A`$ with zero eigenvalue, $`A_{11}^1`$ is infinite. Therefore the proper treatment is to take the trace in the subspace orthogonal to $`\widehat{e}_1`$.. The estimate (17) can be calculated without knowledge of the underlying model, it is $`D/(2N_{\mathrm{spike}})`$. This behavior should also hold in cases where the stimulus dimensions are expanded to include time. The errors are expected to increase in proportion to the increased dimensionality. In the case of a complex cell with two relevant dimensions, the error is expected to be twice that for a cell with single relevant dimension, also discussed in section 4.3.

We emphasize that the error estimate according to Eq. (17) is of the same order as errors of the reverse correlation method when it is applied for Gaussian ensembles. The latter are given by $`(\mathrm{Tr}[C^1]C_{11}^1)/[2N_{\mathrm{spike}}f^{}_{}{}^{}2(s_1)]`$. Of course, if the reverse correlation method were to be applied to the non–Gaussian ensemble, the errors would be larger. In Fig. 5 we show the result of simulations for various numbers of trials, and therefore $`N_{\mathrm{spike}}`$. The average projection of the normalized reconstructed vector $`\widehat{v}`$ on the RF $`\widehat{e}_1`$ behaves initially as $`1/N_{\mathrm{spike}}`$ (dashed line). For smaller data sets, in this case $`N_{\mathrm{spikes}}\text{ }<30,000`$, corrections $`N_{\mathrm{spikes}}^2`$ become important for estimating the expected errors of the algorithm. Happily these corrections have a sign such that smaller data sets are more effective than one might have expected from the asymptotic calculation. This can be verified from the expansion $`\widehat{v}\widehat{e}_1=\left[1\delta 𝐯^2\right]^{1/2}1\frac{1}{2}\delta 𝐯^2+\frac{3}{8}\delta 𝐯^4`$, where only the first two terms where taken into account in Eq. 17.

### 4.3 A model complex cell

A sequence of spikes from a model cell with two relevant dimensions was simulated by projecting each of the stimuli on vectors that differ by $`\pi /2`$ in their spatial phase, taken to mimic properties of complex cells, as in Fig. 6. A particular frame leads to a spike according to a logical OR, that is if either $`|s_1|`$ or $`|s_2|`$ exceeds a threshold value $`\theta `$ in the presence of noise, where $`s_1=𝐬\widehat{e}_1`$, $`s_2=𝐬\widehat{e}_2`$. Similarly to (14),

$$\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}=f(s_1,s_2)=H(|s_1|\theta \xi _1)H(|s_2|\theta \xi _2),$$
(18)
where $`\xi _1`$ and $`\xi _2`$ are independent Gaussian variables. The sampling of this input–output function by our particular set of natural stimuli is shown in Fig. 6(c). As is well known, reverse correlation fails in this case because the spike–triggered average stimulus is zero, although with Gaussian stimuli the spike–triggered covariance method would recover the relevant dimensions (Touryan et al., 2002). Here we show that searching for maximally informative dimensions allows us to recover the relevant subspace even under more natural stimulus conditions.

We start by maximizing information with respect to one vector. Contrary to the result Fig. 4(e) for a simple cell, one optimal dimension recovers only about 60% of the total information per spike \[Eq. (4)\]. Perhaps surprisingly, because of the strong correlations in natural scenes, even a projection onto a random vector in the $`D10^3`$ dimensional stimulus space has a high probability of explaining 60% of total information per spike, as can be seen in Fig. 2. We therefore go on to maximize information with respect to two vectors. As a result of maximization, we obtain two vectors $`𝐯_1`$ and $`𝐯_2`$ shown in Fig. 6. The information along them is $`I(𝐯_1,𝐯_2)0.90`$ which is within the range of information values obtained along different linear combinations of the two model vectors $`I(\widehat{e}_1,\widehat{e}_2)/I_{\mathrm{spike}}=0.89\pm 0.11`$. Therefore, the description of neuron’s firing in terms of vectors $`𝐯_1`$ and $`𝐯_2`$ is complete up to the noise level, and we do not have to look for extra relevant dimensions. Practically, the number of relevant dimensions can be determined by comparing $`I(𝐯_1,𝐯_2)`$ to either $`I_{\mathrm{spike}}`$ or $`I(𝐯_1,𝐯_2,𝐯_3)`$, the later being the result of maximization with respect to three vectors simultaneously. As mentioned in the Introduction, information along set a of vectors does not increase when extra dimensions are added to the relevant subspace. Therefore, if $`I(𝐯_1,𝐯_2)=I(𝐯_1,𝐯_2,𝐯_3)`$ (again, up to the noise level), then this means that there are only 2 relevant dimensions. Using $`I_{\mathrm{spike}}`$ for comparison with $`I(𝐯_1,𝐯_2)`$ has the advantage of not having to look for an extra dimension, which can be computationally intensive. However $`I_{\mathrm{spike}}`$ might be subject to larger systematic bias errors than $`I(𝐯_1,𝐯_2,𝐯_3)`$.

Vectors $`𝐯_1`$ and $`𝐯_2`$ obtained by maximizing $`I(𝐯_1,𝐯_2)`$ are not exactly orthogonal, and are also rotated with respect to $`\widehat{e}_1`$ and $`\widehat{e}_2`$. However, the quality of reconstruction, as well as the value of information $`I(𝐯_1,𝐯_2)`$, is independent of a particular choice of basis with the RS. The appropriate measure of similarity between the two planes is the dot product of their normals. In the example of Fig. 6, $`\widehat{n}_{(\widehat{e}_1,\widehat{e}_2)}\widehat{n}_{(𝐯_1,𝐯_2)}=0.82\pm 0.07`$, where $`\widehat{n}_{(\widehat{e}_1,\widehat{e}_2)}`$ is a normal to the plane passing through vectors $`\widehat{e}_1`$ and $`\widehat{e}_2`$. Maximizing information with respect to two dimensions requires a significantly slower cooling rate, and consequently longer computational times. However, the expected error in the reconstruction, $`1\widehat{n}_{(\widehat{e}_1,\widehat{e}_2)}\widehat{n}_{(𝐯_1,𝐯_2)}`$, scales as $`1/N_{\mathrm{spike}}`$ behavior, similarly to (17), and is roughly twice that for a simple cell given the same number of spikes. We make vectors $`𝐯_1`$ and $`𝐯_2`$ orthogonal to each others upon completion of the algorithm.

### 4.4 A model auditory neuron with one relevant dimension

Because stimuli $`𝐬`$ are treated as vectors in an abstract space, the method of looking for the most informative dimensions can be applied equally well to auditory as well as to visual neurons. Here we illustrate the method by considering a model auditory neuron with one relevant dimension, which is shown in Fig. 7(c) and is taken to mimic the properties of cochlear neurons. The model neuron is probed by two ensembles of naturalistic stimuli: one is a recording of a native Russian speaker reading a piece of Russian prose, and the other one is a recording of a piece of English prose read by a native English speaker. Both of the ensembles are non–Gaussian and exhibit amplitude distributions with long, nearly exponential tails, cf. Fig. 7(a), which are qualitatively similar to those of light intensities in natural scenes (Voss & Clarke, 1975; Ruderman, 1994). However, the power spectrum is different in the two cases, as can be seen in Fig. 7(b). The differences in the correlation structure in particular lead to different STAs across the two ensembles, cf. panel (d). Both of the STAs also deviate from the model filter shown in panel (c).

Despite differences in the probability distributions $`P(𝐬)`$, it is possible to recover the relevant dimension of the model neuron by maximizing information. In panel (e) we show the two most informative vectors found by running the algorithm for the two ensembles, and replot the model filter from (c) to show that the three vectors overlap almost perfectly. Thus different non–Gaussian correlations can be successfully removed to obtain an estimate of the relevant dimension. If the most informative vector changes with the stimulus ensemble, this can be interpreted as caused by adaptation to the probability distribution.

## 5 Summary

Features of the stimulus that are most relevant for generating the response of a neuron can be found by maximizing information between the sequence of responses and the projection of stimuli on trial vectors within the stimulus space. Calculated in this manner, information becomes a function of direction in stimulus space. Those vectors that maximize the information and account for the total information per response of interest span the relevant subspace. The method allows multiple dimensions to be found. The reconstruction of the relevant subspace is done without assuming a particular form of the input–output function. It can be strongly nonlinear within the relevant subspace, and is estimated from experimental histograms for each trial direction independently. Most importantly, this method can be used with any stimulus ensemble, even those that are strongly non–Gaussian as in the case of natural signals. We have illustrated the method on model neurons responding to natural scenes and sounds. We expect the current implementation of the method to be most useful in cases where several most informative vectors ($`10`$, depending on their dimensionality) are to be analyzed for neurons probed by natural scenes. This technique could be particularly useful in describing sensory processing in poorly understood regions of higher level sensory cortex (such as visual areas V2, V4 and IT and auditory cortex beyond A1) where white noise stimulation is known to be less effective than naturalistic stimuli.

## Acknowledgments

We thank K. D. Miller for many helpful discussions. Work at UCSF was supported in part by the Sloan and Swartz Foundations and by a training grant from the NIH. Our collaboration began at the Marine Biological Laboratory in a course supported by grants from NIMH and the Howard Hughes Medical Institute.

## Appendix A Limitations of the reverse correlation method

Here we examine what sort of deviations one can expect when applying the reverse correlation method to natural stimuli even in the model with just one relevant dimension. There are two factors that, when combined, invalidate the reverse correlation method: the non–Gaussian character of correlations and the nonlinearity of the input/output function (Ringach et al., 1997). In its original formulation (de Boer & Kuyper, 1968), the neuron is probed by white noise and the relevant dimension $`\widehat{e}_1`$ is given by the STA $`\widehat{e}_1𝐬r(𝐬)`$. If the signals are not white, i.e. the covariance matrix $`C_{ij}=s_is_j`$ is not a unit matrix, then the STA is a broadened version of the original filter $`\widehat{e}_1`$. This can be seen by noting that for any function $`F(𝐬)`$ of Gaussian variables $`\{s_i\}`$ the identity holds:

$$s_iF(𝐬)=s_is_j_{s_j}F(𝐬),_j_{s_j}.$$
(19)
When property (19) is applied to the vector components of the STA, $`s_ir(𝐬)=C_{ij}_jr(𝐬)`$. Since we work within the assumption that the firing rate is a (nonlinear) function of projection onto one filter $`\widehat{e}_1`$, $`r(𝐬)=r(s_1)`$, the later average is proportional to the model filter itself, $`_jr=\widehat{e}_{1j}r^{}(s_1)`$. Therefore, we arrive at the prescription of the reverse correlation method

$$\widehat{e}_{1i}[C^1]_{ij}s_jr(𝐬).$$
(20)
The Gaussian property is necessary in order to represent the STA as a convolution of the covariance matrix $`C_{ij}`$ of the stimulus ensemble and the model filter. To understand how the reconstructed vector obtained according to Eq. (20) deviates from the relevant one, we consider weakly non–Gaussian stimuli, with the probability distribution

$$P_{nG}(𝐬)=\frac{1}{Z}P_0(𝐬)e^{ϵH_1(𝐬)},$$
(21)
where $`P_0(𝐬)`$ is the Gaussian probability distribution with covariance matrix $`C`$, and the normalization factor $`Z=e^{ϵH_1(𝐬)}`$. The function $`H_1`$ describes deviations of the probability distribution from Gaussian, and therefore we will set $`s_iH_1=0`$ and $`s_is_jH_1=0`$, since these averages can be accounted for in the Gaussian ensemble. In what follows we will keep only the first order terms in perturbation parameter $`ϵ`$. Using the property (19), we find the STA to be given by

$$s_ir_{nG}=s_is_j\left[_jr+ϵr_j(H_1)\right],$$
(22)
where averages are taken with respect to the Gaussian distribution. Similarly, the covariance matrix $`C_{ij}`$ evaluated with respect to the non–Gaussian ensemble is given by:

$$C_{ij}=\frac{1}{Z}s_is_je^{ϵH_1}=s_is_j+ϵs_is_ks_j_k(H_1)$$
(23)
so that to the first order in $`ϵ`$, $`s_is_j=C_{ij}ϵC_{ik}s_j_k(H_1)`$. Combining this with Eq. (22), we get

$$s_ir_{nG}=\mathrm{const}\times C_{ij}\widehat{e}_{1j}+ϵC_{ij}\left(rs_1r^{}\right)_j(H_1).$$
(24)
The second term in (24) prevents the application of the reverse correlation method for non–Gaussian signals. Indeed, if we multiply the STA (24) with the inverse of the a priori covariance matrix $`C_{ij}`$ according to the reverse correlation method (20), we no longer obtain the RF $`\widehat{e}_1`$. The deviation of the obtained answer from the true RF increases with $`ϵ`$, which measures the deviation of the probability distribution from Gaussian. Since natural stimuli are known to be strongly non–Gaussian, this makes the use of the reverse correlation problematic when analyzing neural responses to natural stimuli.

The difference in applying the reverse correlation to stimuli drawn from a correlated Gaussian ensemble vs. a non–Gaussian one is illustrated in Figs. 8 (b) and (c). In the first case, shown in (b), stimuli are drawn from a correlated Gaussian ensemble with the covariance matrix equal to that of natural images. In the second case, shown in (c), the patches of photos are taken as stimuli. The STA is broadened in both cases. Even though the two-point correlations are just as strong in the case of Gaussian stimuli as they are in the natural stimuli ensemble, Gaussian correlations can be successfully removed from the STA according to Eq. (20) to obtain the model filter. On the contrary, an attempt to use reverse correlation with natural stimuli results in an altered version of the model filter. We reiterate that for this example the apparent noise in the decorrelated vector is not due to neural noise or finite datasets, since the “exact” STA has been used (15) in all calculations presented in Figs. 8 and 9.

The reverse correlation method gives the correct answer for any distribution of signals if the probability of generating a spike is a linear function of $`s_i`$, since then the second term in Eq. (24) is zero. In particular, a linear input-output relation could arise due to a neural noise whose variance is much larger than the variance of the signal itself. This point is illustrated in Figs. 9 (a), (b), and (c), where the reverse correlation method is applied to a threshold input–output function at low, moderate, and high signal-to-noise ratios. For small signal-to-noise ratios where the noise standard deviation is similar to that of projections $`s_1`$, the threshold nonlinearity in the input-output function is masked by noise, and is effectively linear. In this limit, the reverse correlation can be applied with the exact STA. However, for experimentally calculated STA at low signal-to-noise ratios the decorrelation procedure results in strong noise amplification. On the other hand, at higher signal-to-noise ratios decorrelation fails due to the nonlinearity of the input-output function in accordance with (24).

## Appendix B Maxima of $`I(v)`$: what do they mean?

The relevant subspace of dimensionality $`K`$ can be found by maximizing information simultaneously with respect to $`K`$ vectors. The result of maximization with respect to a number of vectors that is less than the true dimensionality of the relevant subspace may produce vectors which have components in the irrelevant subspace. This happens only in the presence of correlations in stimuli. As an illustration, we consider the situation where the dimensionality of the relevant subspace $`K=2`$, and vector $`\widehat{e}_1`$ describes the most informative direction within the relative subspace. We show here that even though the gradient of information is perpendicular to both $`\widehat{e}_1`$ and $`\widehat{e}_2`$, it may have components outside the relevant subspace. Therefore the vector $`v_{\mathrm{max}}`$ that corresponds to the maximum of $`I(v)`$ will then lie outside the relevant subspace. We recall from Eq. (10) that

$$I(\widehat{e}_1)=𝑑s_1P(s_1)\frac{d}{ds_1}\frac{P(s_1|\mathrm{spike})}{P(s_1)}(𝐬|s_1,\mathrm{spike}𝐬|s_1),$$
(25)
We can rewrite the conditional averages $`𝐬|s_1=𝑑s_2P(s_1,s_2)𝐬|s_1,s_2/P(s_1)`$ and
$`𝐬|s_1,\mathrm{spike}=𝑑s_2f(s_1,s_2)P(s_1,s_2)𝐬|s_1,s_2/P(s_1|\mathrm{spike})`$, so that

$`I(\widehat{e}_1)={\displaystyle 𝑑s_1𝑑s_2P(s_1,s_2)𝐬|s_1,s_2\frac{P(\mathrm{spike}|s_1,s_2)P(\mathrm{spike}|s_1)}{P(\mathrm{spike})}\frac{d}{ds_1}\mathrm{ln}\frac{P(s_1|\mathrm{spike})}{P(s_1)}}.`$ (26)
Because we assume that the vector $`\widehat{e}_1`$ is the most informative within the relevant subspace, $`\widehat{e}_1I=\widehat{e}_2I=0`$, so that the integral in (26) is zero for those directions in which the component of the vector $`𝐬|s_1,s_2`$ changes linearly with $`s_1`$ and $`s_2`$. For uncorrelated stimuli this is true for all directions, so that the most informative vector within the relevant subspace is also the most informative in the overall stimulus space. In the presence of correlations, the gradient may have non-zero components along some irrelevant directions if projection of the vector $`𝐬|s_1,s_2`$ on those directions is not a linear function of $`s_1`$ and $`s_2`$. By looking for a maximum of information we will therefore be driven outside the relevant subspace. The deviation of $`v_{\mathrm{max}}`$ from the relevant subspace is also proportional to the strength of the dependence on the second parameter $`s_2`$, because of the factor $`[P(s_1,s_2|\mathrm{spike})/P(s_1,s_2)P(s_1|\mathrm{spike})/P(s_1)]`$ in the integrand.

## Appendix C The gradient of information

According to the expression (7), the information $`I(𝐯)`$ depends on the vector $`𝐯`$ only through the probability distributions $`P_𝐯(x)`$ and $`P_𝐯(x|\mathrm{spike})`$. Therefore we can express the gradient of information in terms of gradients of those probability distributions:

$`_𝐯I={\displaystyle \frac{1}{\mathrm{ln}2}}{\displaystyle 𝑑x\left[\mathrm{ln}\frac{P_𝐯(x|\mathrm{spike})}{P_𝐯(x)}_𝐯(P_𝐯(x|\mathrm{spike}))\frac{P_𝐯(x|\mathrm{spike})}{P_𝐯(x)}_𝐯(P_𝐯(x))\right]},`$ (27)
where we took into account that $`𝑑xP_𝐯(x|\mathrm{spike})=1`$ and does not change with $`𝐯`$. To find gradients of the probability distributions, we note that

$`_𝐯P_𝐯(x)=_𝐯\left[{\displaystyle 𝑑𝐬P(𝐬)\delta (x𝐬𝐯)}\right]={\displaystyle 𝑑𝐬P(𝐬)𝐬\delta ^{}(x𝐬𝐯)}={\displaystyle \frac{d}{dx}}\left[p(x)𝐬|x\right],`$ (28)
and analogously for $`P_𝐯(x|\mathrm{spike})`$:

$$_𝐯P_𝐯(x|\mathrm{spike})=\frac{d}{dx}\left[p(x|\mathrm{spike})𝐬|x,\mathrm{spike}\right].$$
(29)
Substituting expressions (28) and (29) into Eq. (27) and integrating once by parts we obtain:

$$_𝐯I=𝑑xP_𝐯(x)\left[𝐬|x,\mathrm{spike}𝐬|x\right]\left[\frac{d}{dx}\frac{P_𝐯(x|\mathrm{spike})}{P_𝐯(x)}\right],$$

which is the expression (10) of the main text.


format: text

text:
# 1. Introduction

## 1. Introduction

According to the Standard cosmological model , the Universe has been expanding and cooling from a split second after the Big Bang to the present moment and remained uniform and isotropic overall in doing so. In the process of its evolution, the Universe has experienced a chain of phase transitions with spontaneous symmetry breaking, including Grand Unification, electroweak phase transition, formation of neutrons and protons from quarks, recombination, and so forth. Regions with spontaneously broken symmetry, which are more than the correlation length apart, are statistically independent. At interfaces between these regions, so-called topological defects necessarily arise. A systematic exposition of the potential role of topological defects in our Universe is provided by Vilenkin and Shellard . The particular types of defects: domain walls, strings, monopoles, or textures are determined by the topological properties of vacuum . If the vacuum manifold after the breakdown is not shrinkable to a point, then the Polyakov-t’Hooft monopole-type solutions  appear in quantum field theory.

Spontaneous symmetry breaking (SSB) plays a fundamental role in modern attempts to construct particle theories. A symmetry in this context is not necessarily associated with space-time transformations. It can be a kind of “internal” symmetry as well, such as the Grand Unification symmetry, the electroweak and the isotopic symmetry, or even supersymmetry, whose transformations mix bosons and fermions. Topological defects, caused by spontaneous breaking of internal symmetries (independent of space-time coordinates), are called global.

A fundamental property of global symmetry violation is the Goldstone degree of freedom. In the monopole case, the term related to the Goldstone boson in the energy-momentum tensor decreases rather slowly away from the center. As a result, the total energy of a global monopole grows linearly with distance, in other words, diverges. Without gravity such a divergence is a general property of spontaneously broken global symmetries. In his pioneering paper  Polyakov mentioned two possibilities of avoiding this difficulty. The first one was to combine a monopole with a Yang-Mills field. This idea was independently considered by t’Hooft . This, among other reasons, gave rise to numerous papers on gauge (magnetic) monopoles. The second possibility was to consider a bound monopole-antimonopole system, whose total energy would be large (proportional to the distance between the components) but finite.

One more opportunity is to take into account the self-gravity of global monopoles, which can in principle remove the above self-energy problem. This is also necessary for potential astrophysical applications. Such a study was first performed by Barriola and Vilenkin  who found that the gravitational field outside a monopole is characterized by a solid angle deficit proportional to the SSB energy scale. Harari and Lousto showed that the gravitational mass of a global monopole, calculated using the Tolman integral, is negative . Solutions with a horizon for supermassive global monopoles were found by Liebling , who also confirmed the estimate of Ref.  for the upper value of the symmetry breaking energy compatible with a static configuration. The existence of de Sitter cores inside global monopoles and other topological defects gave rise to the idea of “topological inflation” .

For global strings in flat space, the energy per unit length (without gravitation) also diverges with growing distance from the axis, but only logarithmically. However, in general relativity integration over the cross-section yields a finite result . The gravitational interaction thus leads to self-localization of a global string. Will a similar effect take place for a global monopole? An attempt to answer this question, which appears not to be answered in the existing papers, was one of motivations for reconsidering the gravitational properties of a global monopole.

The previous studies have used the boundary condition according to which the symmetry-breaking potential should vanish at spatial infinity. Our approach is different: we do not even assume the existence of a spatial asymptotic but require regularity at the center and try to observe the properties of the whole set of global monopole solutions. In doind so, among other quantities, we discuss the behavior of the total scalar field energy, which turns out to be finite in static regions of supermassive global monopoles.

In Sec. 2 we present the complete sets of equations for a static spherically symmetric gravitating global monopole in two most convenient coordinate systems, namely, with quasiglobal and harmonic radial coordinates. The general properties of static global monopoles are summarized in Sec. 3. In Sec. 4 we analyze analytically and numerically the specific features of a global monopole in the particular case of the “Mexican hat” potential. Sec. 5 contains a general discussion of our results, including their possible cosmological interpretation.

## 2. Equations and boundary conditions

### 2.1. General problem setting

We begin with the most general form of a static, spherically symmetric metric, without specifying the radial coordinate $`x^1=u`$:

$$ds^2=g_{\mu \nu }dx^\mu dx^\nu =\mathrm{e}^{2F_0}dt^2\mathrm{e}^{2F_1}du^2\mathrm{e}^{2F_\mathrm{\Omega }}d\mathrm{\Omega }^2.$$
(1)
Here $`d\mathrm{\Omega }^2=d\theta ^2+\mathrm{sin}^2\theta d\phi ^2`$ is the linear element on a unit sphere and $`F_0`$, $`F_1`$, and $`F_\mathrm{\Omega }`$ are functions of $`u`$.

The nonzero components of the Ricci tensor are (the prime denotes $`d/du`$)

$`R_0^{0}`$ $`=`$ $`\mathrm{e}^{2F_1}[F_0^{\prime \prime }+F_0^{}(F_1^{}+2F_\mathrm{\Omega }^{}+F_0^{})];`$
$`R_1^{1}`$ $`=`$ $`\mathrm{e}^{2F_1}[F_0^{\prime \prime }+2F_\mathrm{\Omega }^{\prime \prime }+2F_\mathrm{\Omega }^{}{}_{}{}^{2}+F_0^{}{}_{}{}^{2}F_1^{}(2F_\mathrm{\Omega }^{}+F_0^{})];`$
$`R_2^{2}`$ $`=`$ $`R_3^{3}=\mathrm{e}^{2F_\mathrm{\Omega }}+\mathrm{e}^{2F_1}[F_\mathrm{\Omega }^{\prime \prime }+F_\mathrm{\Omega }^{}(F_1^{}+2F_\mathrm{\Omega }^{}+F_0^{})].`$ (2)
Consider the Lagrangian describing a triplet of real scalar fields $`\varphi ^a`$ ($`a=1,2,3`$) in general relativity:

$$L=\frac{R}{16\pi G}+\frac{1}{2}g^{\mu \nu }_\mu \varphi ^a_\nu \varphi ^aV(\varphi ),$$
(3)
where $`R`$ is the scalar curvature, $`V(\varphi )`$ is a potential depending on $`\varphi =\pm \sqrt{\varphi ^a\varphi ^a}`$ and $`G`$ is the gravitational constant. We use the natural units such that

$$\mathrm{}=c=1,$$
(4)
so that $`G=m_{\mathrm{pl}}^2`$, where $`m_{\mathrm{pl}}=1.22\times 10^{19}`$ GeV is the Planck mass.

To obtain a global monopole with unit topological charge , let us assume that the metric has the form (1) while $`\varphi ^a`$ comprise the following “hedgehog” configuration:

$$\varphi ^a=\varphi (u)n^a,n^a=\{\mathrm{sin}\theta \mathrm{cos}\phi ,\mathrm{sin}\theta \mathrm{sin}\phi ,\mathrm{cos}\theta \}.$$
(5)
The Einstein equations can be written in the form

$$R_\mu ^\nu =8\pi G\stackrel{~}{T}_\mu ^\nu =8\pi G(T_\mu ^\nu \frac{1}{2}\delta _\mu ^\nu T_\alpha ^\alpha )$$
(6)
where $`T_\mu ^\nu `$ is the energy-momentum tensor and the nonzero components of $`\stackrel{~}{T}_\mu ^\nu `$ are

$$\stackrel{~}{T}_0^0=V,\stackrel{~}{T}_1^1=V\mathrm{e}^{2F_1}\varphi ^{}{}_{}{}^{2},\stackrel{~}{T}_2^2=\stackrel{~}{T}_3^3=V\mathrm{e}^{2F_\mathrm{\Omega }}\varphi ^2$$
(7)
The regular center conditions for the metric (1) are that at the corresponding value $`u_c`$ of the coordinate $`x^1=u`$

$$\mathrm{e}^{F_\mathrm{\Omega }}0;F_0=F_{0c}+O(\mathrm{e}^{2F_\mathrm{\Omega }}),\mathrm{e}^{F_1+F_\mathrm{\Omega }}|F_\mathrm{\Omega }^{}|1.$$
(8)
The last condition is necessary for local flatness and provides the correct circumference to radius ratio for coordinate circles at small $`r=\mathrm{e}^{F_\mathrm{\Omega }}`$.

The scalar field energy, defined as the partial time derivative of the scalar field action, $`E=S/t`$, is a conserved quantity for our static system:

$$E=\sqrt{g}T_0^0d^3x=4\pi \mathrm{e}^{F_0+F_1+2F_\mathrm{\Omega }}[\frac{1}{2}\mathrm{e}^{2F_1}\varphi ^{}{}_{}{}^{2}+\mathrm{e}^{2F_\mathrm{\Omega }}\varphi ^2+V]du,$$
(9)
where $`g`$ is the metric tensor determinant.

In what follows we will make some general inferences without specifying the potential $`V(\varphi )`$ and then perform a more detailed study for the simplest and most frequently used symmetry-breaking potential

$$V(\varphi )=\frac{1}{4}\lambda (\varphi ^a\varphi ^a\eta ^2)^2=\frac{1}{4}\eta ^4\lambda (f^21)^2,$$
(10)
where $`\eta >0`$ characterizes the energy of symmetry breaking, $`\lambda `$ is a dimensionless constant and $`f(u)=\varphi (u)/\eta `$ is the normalized field magnitude playing the role of an order parameter. The model has a global $`SO(3)`$ symmetry, which can be spontaneously broken to $`SO(2)`$ due to the potential wells ($`V=0`$) at $`f=\pm 1`$.

Let us now write down the Einstein equations and the boundary conditions explicitly in two coordinate frames to be used.

### 2.2. The quasiglobal coordinate $`\rho `$

The first choice is the coordinate $`u=\rho `$ specified by the condition $`F_0+F_1=0`$. Denoting $`\mathrm{e}^{2F_0}=\mathrm{e}^{2F_1}=A(\rho )`$ and $`\mathrm{e}^{F_\mathrm{\Omega }}=r(\rho )`$, we obtain the metric in the form

$$ds^2=A(\rho )dt^2\frac{d\rho ^2}{A(\rho )}r^2(\rho )d\mathrm{\Omega }^2.$$
(11)
The scalar field equation

$$\varphi ^a+V/\varphi ^a=0,$$
(12)
where $`=^\alpha _\alpha `$ is the d’Alembert operator, and certain combinations of the Einstein equations have the form

$`(Ar^2\varphi ^{})^{}2\varphi `$ $`=`$ $`r^2dV/d\varphi ;`$ (13)
$`(A^{}r^2)^{}`$ $`=`$ $`16\pi Gr^2V;`$ (14)
$`2r^{\prime \prime }/r`$ $`=`$ $`8\pi G\varphi ^{}{}_{}{}^{2};`$ (15)
$`A(r^2)^{\prime \prime }r^2A^{\prime \prime }`$ $`=`$ $`2(18\pi G\varphi ^2);`$ (16)
$`A^{}rr^{}+Ar_{}^{}{}_{}{}^{2}1`$ $`=`$ $`8\pi G[\frac{1}{2}Ar^2\varphi ^{}{}_{}{}^{2}\varphi ^2r^2V],`$ (17)
where the prime denotes $`d/d\rho `$. Only three of these five equations are independent: the scalar field equation (13) follows from the Einstein equations, while Eq. (17) is a first integral of the others. Given a potential $`V(\phi )`$, this is a determined set of equations for the unknowns $`r,A,\varphi `$.

This choice of the coordinates is preferable for considering Killing horizons, which correspond to zeros of the function $`A(\rho )`$, since such zeros are regular points of Eqs. (13)–(17); moreover, in a close neighborhood of a horizon, the coordinate $`\rho `$ defined in this manner varies (up to a positive constant factor) like manifestly well-behaved Kruskal-like coordinates used for analytic continuation of the metric . Therefore one can jointly consider regions at both sides of a horizon in terms of $`\rho `$, and, in general, the whole range of $`\rho `$ can contain several horizons. For this reason the coordinate $`\rho `$ can be called quasiglobal.

The regular center conditions (8) are fulfilled if, near some value $`\rho _c`$ of the coordinate $`\rho `$,

$$A(\rho )=A_c+O((\rho \rho _c)^2),r(\rho )(\rho \rho _c)/\sqrt{A_c}.$$
(18)
In regions where $`A<0`$ (sometimes called T-regions), if any, the coordinate $`\rho `$ is timelike and $`t`$ is space-like. Changing the notations, $`tx`$, and introducing the proper time of a comoving observer in a T-region,

$$\tau =𝑑\rho /\sqrt{|A(\rho )|},$$
(19)
we can rewrite the metric in the form

$$ds^2=d\tau ^2|A(\tau )|dx^2r^2(\tau )d\mathrm{\Omega }^2.$$
(20)
The space-time geometry then corresponds to a homogeneous anisotropic cosmological model of Kantowski-Sachs (KS) type , where spatial sections have the topology $`\times 𝕊^2`$.

### 2.3. The harmonic coordinate $`u`$

Another convenient variable, well simplifying the form of the equations, is the harmonic coordinate $`u`$ specified by the condition 333A cylindrical version of the harmonic radial coordinate has been used previously in the analysis of gravitational properties of current-conducting filaments  and cosmic strings .

$$F_1=2F_\mathrm{\Omega }+F_0,$$
(21)
so that $`u=0`$. Then the field equations may be written as

$`\varphi ^{\prime \prime }2\mathrm{e}^{F_0+F_1}\varphi `$ $`=`$ $`\mathrm{e}^{2F_1}dV/d\varphi ,`$ (22)
$`F_0^{\prime \prime }`$ $`=`$ $`8\pi G\mathrm{e}^{2F_1}V,`$ (23)
$`F_1^{\prime \prime }2F_\mathrm{\Omega }^{}(F_\mathrm{\Omega }^{}+2F_0^{})`$ $`=`$ $`8\pi G(\varphi ^{}{}_{}{}^{2}+\mathrm{e}^{2F_1}V),`$ (24)
$`F_\mathrm{\Omega }^{\prime \prime }\mathrm{e}^{2(F_0+F_\mathrm{\Omega })}`$ $`=`$ $`8\pi G(\varphi ^2\mathrm{e}^{2(F_0+F_\mathrm{\Omega })}+\mathrm{e}^{2F_1}V),`$ (25)
$`\mathrm{e}^{2F_\mathrm{\Omega }}+\mathrm{e}^{2F_1}(F_\mathrm{\Omega }^{}{}_{}{}^{2}+2F_\mathrm{\Omega }^{}F_0^{})`$ $`=`$ $`8\pi G(\frac{1}{2}\mathrm{e}^{2F_1}\varphi ^{}{}_{}{}^{2}\mathrm{e}^{2F_\mathrm{\Omega }}\varphi ^2V)`$ (26)
where the prime denotes $`d/du`$.

It is straightforward to obtain that a regular center can only correspond to $`u\pm \mathrm{}`$; we choose $`u\mathrm{}`$, where one should have

$$\mathrm{e}^{F_\mathrm{\Omega }}1/|u|,\mathrm{e}^{F_0}=\sqrt{A_c}(1+O(u^2)),\mathrm{e}^{F_1}1/u^2,$$
(27)
and $`A_c`$ is the same as in (18).

## 3. General properties of global monopoles

### 3.1. Monopoles in Minkowski space-time

The Minkowski metric in the usual spherical coordinates

$$ds^2=dt^2dr^2r^2d\mathrm{\Omega }^2$$
(28)
is a special case of (11) with $`r\rho `$ and $`A1`$. The only unknown in flat space-time is $`\varphi (r)`$, and the only field equation is (13) which takes the form

$$(r^2\varphi ^{})^{}2\varphi =r^2dV/d\varphi ,$$
(29)
where, in particular, for the potential (10) $`dV/d\varphi =\lambda \varphi (\varphi ^2\eta ^2)`$. In this case the scalar field equation can be written in terms of $`f=\varphi /\eta `$ as

$$r^2(r^2f^{})^{}2fr^2+\lambda \eta ^2f(1f^2)=0.$$
(30)
The energy integral (9) takes the form

$$E=4\pi r^2[\frac{1}{2}\varphi ^{}{}_{}{}^{2}+\frac{\varphi ^2}{r^2}+V]dr,$$
(31)
and for its convergence, in case $`V(\varphi )0`$, all the three terms should vanish quickly enough at infinity:

$$\varphi =o(r^{1/2}),\varphi ^{}=o(r^{3/2}),V=o(r^3)$$
(32)
as $`r\mathrm{}`$. This actually means that a configuration with finite energy is only possible with $`V(0)=0`$, contrary to the symmetry breaking assumption according to which $`V`$ has minima in nonsymmetric states, $`\varphi 0`$. In particular, the potential (10) does not give rise to global monopoles with finite energy. A consideration of self-gravity of the field triad $`\varphi ^a`$ is one of the ways to overcome this difficulty.

The harmonic coordinate $`u`$ in flat space-time is connected with $`r`$ by the relation $`uu_0=\pm 1/r`$, where $`u_0`$ is an arbitrary constant; choosing the minus sign, we find that $`u`$ ranges from $`\mathrm{}`$, which corresponds to the center $`r=0`$, to $`u_0`$ corresponding to spatial infinity.

### 3.2. Solutions with constant $`\varphi `$

Under the assumption $`\varphi =\varphi _0=\mathrm{const}`$, the corresponding value of the potential $`V(\varphi _0)=V_0`$ (times $`8\pi G`$) plays the role of a cosmological constant, and the Einstein equations can be integrated explicitly.

Indeed, in a region where $`\varphi =\mathrm{const}`$, Eq. (15) reduces to $`r^{\prime \prime }=0`$, whence $`r=\alpha \rho +r_0`$, $`\alpha ,r_0=\mathrm{const}`$. It remains to find $`A(r)`$, and this is immediately done by integrating Eq. (16):

$$A(r)=\frac{1\mathrm{\Delta }}{\alpha ^2}\frac{2GM}{r}+Cr^2,\mathrm{\Delta }=8\pi G\varphi _0^2,$$
(33)
where $`M`$ and $`C`$ are integration constants. Substituting (33) into (14), we find

$$C=8\pi GV_0/3.$$
(34)
Thus the solution is essentially determined by the values of $`\varphi _0`$, $`V_0`$ and $`M`$. One more constant, $`\alpha `$, reflects the freedom in choosing the unit of time. It should be noted that this is not a monopole solution. Even if we put $`M=0`$, which is evidently necessary for regularity at $`r=0`$, this solution with constant $`\varphi 0`$ is singular at the center: for (33) with $`M=0`$ the Kretschmann scalar at small $`r`$ is $`𝒦=R_{\alpha \beta \gamma \delta }R^{\alpha \beta \gamma \delta }4\mathrm{\Delta }^2/r^4`$.

With respect to a global monopole, two cases of the solution (33) are of interest. The case $`\varphi _00`$ describes the symmetric state, while the case $`V_0=0`$ gives a possible asymptotic behavior at spatial or temporal infinity.

In case $`\varphi _0=0`$ (the symmetric state), putting $`M=0`$ (which is necessary for a regular center), we arrive at the de Sitter metric

$$ds^2=\left(1\frac{r^2}{r_h^2}\right)dt^2\left(1\frac{r^2}{r_h^2}\right)^1dr^2r^2d\mathrm{\Omega }^2,r_h^2=\frac{8\pi GV_0}{3}.$$
(35)
This metric has a horizon at $`r=r_h`$. At $`r>r_h`$, outside the horizon, $`r`$ becomes a timelike coordinate, and $`t`$ is a spacelike one. Changing the notations as in (19), (20), we obtain the metric

$$ds^2=d\tau ^2\mathrm{sinh}^2(\tau /r_h)dx^2r_h^2\mathrm{cosh}^2(\tau /r_h)d\mathrm{\Omega }^2.$$
(36)
This is a KS cosmology with isotropic inflationary expansion at late times ($`\tau \mathrm{}`$).

In the other case, $`\varphi _00`$ but $`V_0=0`$ (the case of broken symmetry, such as $`\varphi =\eta `$ in the potential (10)), the metric takes the form

$$ds^2=\left(\frac{1\mathrm{\Delta }}{\alpha ^2}\frac{2GM}{r}\right)dt^2\left(\frac{1\mathrm{\Delta }}{\alpha ^2}\frac{2GM}{r}\right)^1dr^2r^2d\mathrm{\Omega }^2,$$
(37)
where the constant $`M`$ has the meaning of mass in the sense that test particles at rest at large $`r`$ experience an acceleration equal to $`GM/r^2`$ in the gravitational field (37). Furthermore, a nonzero value of $`\varphi _0`$ leads to a solid angle deficit $`\mathrm{\Delta }`$ defined in (33) at the asymptotic $`r\mathrm{}`$ (see  for more detail) and to a linear divergence of the integral (9) at large $`r`$.

The general case of Eq. (33) describes the large $`r`$ asymptotic of any solution to Eqs. (13)–(17), provided that such an asymptotic exists and $`\varphi `$ tends quickly enough to a constant value.

For the monopoles to be studied, the metric (37) gives a large $`r`$ asymptotic in the case $`\mathrm{\Delta }<1`$. We will also consider solutions with $`\mathrm{\Delta }>1`$, for which a static asymptotic is absent. Then the metric (37) describes cosmological evolution at late times.

### 3.3. General properties of solutions with varying $`\varphi `$

Consider now the general form of Eqs. (13)–(17) with varying $`\varphi `$, without specifying the potential $`V(\varphi )`$.

First of all, we note that, due to (15), $`r^{\prime \prime }0`$, which rules out any nonsingular configurations without a center such as wormholes and horns (see Theorem 1 of Ref.  for more detail).

Second, Eq. (16) can be rewritten in the form

$$(r^4B^{})^{}=2(18\pi G\varphi ^2),B\stackrel{\mathrm{def}}{=}A/r^2,$$
(38)
and at a point where $`B^{}=0`$ we have $`r^4B^{\prime \prime }=2(18\pi G\varphi ^2)`$. Hence it follows that, as long as $`\varphi ^2<1/(8\pi G)`$ (the $`\varphi `$ field does not reach trans-Planckian values), $`B^{\prime \prime }<0`$ at possible extrema of the function $`B`$. In other words, $`B`$ cannot have a regular minimum.

Our interest is in systems with a regular center satisfying the conditions (18), so that both $`A(\rho )>0`$ and $`B(\rho )>0`$ near $`\rho =\rho _c`$. At a possible horizon $`\rho =h`$, both $`A`$ and $`B`$ vanish, and, since it cannot be a minimum of $`B`$, $`B<0`$ at $`\rho >h`$ near the horizon. At greater $`\rho `$ the function $`B(\rho )`$, having no minima, can only decrease and will never return to zero; hence $`A=Br^2<0`$ at $`\rho >h`$. We conclude that there can be no more than one horizon, and, if it exists, it is simple (corresponds to a simple zero of $`A(\rho )`$). Since the global causal structure of space-time is determined (up to possible identifications of isometric hypersurfaces) by the number and disposition of Killing horizons , we have the following result:

Statement 1. Under the assumption that $`\varphi ^2<1/(8\pi G)`$ in the whole space, our system with a regular center can have either no horizon, or one simple horizon, and in the latter case its global structure is the same as that of de Sitter space-time.

The above reasoning is in essence the same as in the proof of Theorem 2 of Ref.  on the disposition of horizons in scalar-vacuum space-times. It uses only Eq. (16) where the potential $`V`$ does not enter. Therefore the conclusion is valid for systems with any potentials, positive or negative.

Let us return to Eq. (15), according to which $`r^{\prime \prime }0`$. Since at a regular center $`r^{}>0`$, this leaves three possibilities for the function $`r(\rho )`$ (see Fig. 1):

monotonic growth with a decreasing slope, but $`r\mathrm{}`$ as $`\rho \mathrm{}`$,

monotonic growth with $`rr_{\mathrm{max}}<\mathrm{}`$ as $`\rho \mathrm{}`$, and

growth up to $`r_{\mathrm{max}}`$ at some $`\rho _1<\mathrm{}`$ and further decrease, reaching $`r=0`$ at some finite $`\rho _2>\rho _1`$.

In each case, according to Statement 1, a horizon can occur at some $`\rho =h`$ within the range of $`\rho `$, so that at $`\rho >h`$ we have a T-region with the geometry of a KS cosmological model.

We conclude that there are six classes of qualitative behaviors of the solutions, i.e., (a), (b), (c), each with or without a horizon, the latter circumstance to be labelled with the symbol 1 or 0, respectively. Thus, all solutions with a spatial asymptotic belong to class (a0). Class (b0) includes space-times ending with a “tube” consisting of two-dimensional spheres of equal radius. Class (c0) solutions contain a second center at $`\rho =\rho _2`$, and this center can a priori be regular or singular. We thus obtain a static analogue of closed cosmologies. Classes (a1), (b1), (c1) describe different late-time cosmological behaviors in the two directions corresponding to $`𝕊^2`$, whereas the fate of the third spatial direction ($``$) is determined by the function $`A(\rho )`$. In particular, the possible de Sitter asymptotic (36) belongs to class (a1) solutions, and in this case the expansion is isotropic at late times. On the other hand, class (c1) contains models which behave at late times like the Schwarzschild space-time inside the horizon, contracting to $`r=0`$.

This classification is obtained without any assumptions about $`V(\varphi )`$. Solutions with given $`V(\varphi )`$ will contain some of these classes, not necessarily all of them.

In case $`V0`$, Eq. (14) leads to one more important observation: since at a regular center $`A^{}r^2=0`$, we can write (14) in the integral form

$$A^{}r^2=16\pi G_0^\rho V(\overline{\rho })r^2(\overline{\rho })𝑑\overline{\rho }.$$
(39)
Thus, unless $`V0`$, $`A(\rho )`$ is a decreasing function. Eq. (39) leads to the following conclusions.

Statement 1a. If $`V(\varphi )0`$, our system with a regular center can have either no horizon, or one simple horizon, and in the latter case its global structure is the same as that of de Sitter space-time.

Statement 2. If $`V(\varphi )0`$, the second center in class (c0) solutions is singular.

Statement 3. If $`V(\varphi )0`$ and the solution is asymptotically flat, the mass $`M`$ of the global monopole is negative.

Statement 1a shows that, for nonnegative potentials, the assumption $`\varphi ^2<1/(8\pi G)`$ in Statement 1 is unnecessary, and the causal structure types are known for any magnitudes of $`\varphi `$.

Statement 2 follows from $`A^{}(\rho _2)<0`$, whereas at a regular center it should be $`A^{}=0`$, see (18). $`A^{}(\rho _2)=0`$ could only be possible with $`V0`$, but in this case the only solution with a regular center is trivial (flat space, $`\varphi =0`$).

In Statement 3, asymptotic flatness is understood up to the solid angle deficit, i.e., $`r=\rho `$ and $`A`$ is given by (33) with $`C=0`$ at large $`\rho `$. Then, (39) for $`\rho \mathrm{}`$ gives $`2GM`$ on the left-hand side and a negative quantity on the right.

To our knowledge, this simple conclusion, valid for all nonnegative potentials, has been so far obtained only numerically for the particular potential (10) . Note that Statement 3 is an extension to global monopoles of the so-called generalized Rosen theorem , previously known for scalar-vacuum configurations.

Thus, even before studying particular solutions with the potential (10), we have a more or less complete knowledge of what can be expected from such global monopole systems.

## 4. Mexican hat potential

### 4.1. Equations and boundary conditions

Further analysis is performed for the particular “Mexican hat” potential (10). For numerical integration we prefer to use the harmonic coordinate $`u`$ and to work with Eqs. (22)–(25). This variable enters into the equations only via derivatives and is thus invariant under translations $`uu+\mathrm{const}`$.

Introducing the dimensionless quantities

$$\stackrel{~}{u}=u/(\sqrt{\lambda }\eta ),\mathrm{e}^{\stackrel{~}{F}_\mathrm{\Omega }}=\sqrt{\lambda }\eta \mathrm{e}^{F_\mathrm{\Omega }},\mathrm{e}^{\stackrel{~}{F}_1}=\lambda \eta ^2\mathrm{e}^{F_1}$$
(40)
we exclude the parameter $`\lambda `$ from the equations. Indeed, omitting the tildes, we obtain

$`f^{\prime \prime }`$ $`=`$ $`\mathrm{e}^{2(F_0+F_\mathrm{\Omega })}[2\mathrm{e}^{2F_\mathrm{\Omega }}(1f^2)]f,`$ (41)
$`F_0^{\prime \prime }`$ $`=`$ $`{\displaystyle \frac{\gamma }{4}}\mathrm{e}^{2(F_0+2F_\mathrm{\Omega })}(f^21)^2,`$ (42)
$`F_\mathrm{\Omega }^{\prime \prime }`$ $`=`$ $`\mathrm{e}^{2(F_0+F_\mathrm{\Omega })}[1\gamma f^2{\displaystyle \frac{\gamma }{4}}\mathrm{e}^{2F_\mathrm{\Omega }}(1f^2)^2]`$ (43)
The condition (21) is preserved for the newly defined quantities, but the metric now reads

$$ds^2=\mathrm{e}^{2F_0}dt^2\frac{\mathrm{e}^{2F_1}du^2+\mathrm{e}^{2F_\mathrm{\Omega }}d\mathrm{\Omega }^2}{\lambda \eta ^2}.$$
(44)
The boundary conditions at $`u\mathrm{}`$ are

$$f=0,F_0=0,F_0^{}=0,F_\mathrm{\Omega }=\mathrm{ln}(u)+o(1/|u|),\mathrm{as}u\mathrm{},$$
(45)
They follow from the requirement of regularity at the center and a particular choice of the time unit ($`F_0=0`$) and the origin of the $`u`$ coordinate (the fourth condition).

There remains only one dimensionless parameter in Eqs. (41)–(43),

$$\gamma =8\pi G\eta ^2,$$
(46)
it is the squared energy of symmetry breaking in Planck units.

It is easy to obtain that $`\gamma =1`$ is a critical value of this parameter. Indeed, if we suppose the existence of a large $`r`$ asymptotic at which $`f1`$, i.e., the field tends to the minimum of the potential (10), then the asymptotic form of the metric at large $`r`$ is (37) with $`\mathrm{\Delta }=\gamma `$. Consequently, the asymptotic can be static only if $`\gamma 1`$, whereas for $`\gamma >1`$ the large $`r`$ asymptotic can be only cosmological (KS type), and there is a horizon separating such an outer region from the static monopole core.

On the other hand, if a configuration with $`\gamma <1`$ possesses a horizon, there is again a KS cosmology outside it, but there cannot be a large $`r`$ asymptotic and, according to Sec. 3, the solutions belong to classes (b1) or (c1).

Now, leaving aside the rather well studied case of solutions with a static asymptotic , belonging to class (a0) according to Sec. 3, let us suppose that there is a horizon and return to Eqs. (41)–(43). The horizon corresponds to $`u+\mathrm{}`$. For such cases, in addition to (45), we impose the boundary condition

$$f(u)f_h,|f_h|<\mathrm{}\mathrm{as}u+\mathrm{}.$$
(47)
This condition is necessary for regularity of a solution on the horizon and is applicable to classes (a1), (b1), (c1).

For class (a0) solutions, having a spatial asymptotic and no horizon, the condition (47) is meaningless. Moreover, the coordinate $`u`$ then ranges from $`\mathrm{}`$ to some $`u_0<\mathrm{}`$ such that $`r(u_0)=\mathrm{}`$.

For configurations of classes (a0) and (a1), the commonly used boundary condition at large $`r`$ is

$$f1\mathrm{as}r\mathrm{}.$$
(48)
It is of interest that in the case (a1), to which both condition are applicable, the condition (47), being less restrictive, still leads to solutions satisfying (48) due to the properties of the physical system itself.

The set of equations (41)–(43) with the boundary conditions (45) and (47) comprise a well-posed nonlinear eigenvalue problem. Its trivial solution, with $`f=0`$ and the de Sitter metric (35), describes the symmetric state (with unbroken symmetry). Nontrivial solutions, describing hedgehog configurations with spontaneously broken symmetry, can be found numerically and yield a sequence of eigenvalues $`\gamma _n`$, $`n=0,1,\mathrm{}`$ and the corresponding values of the horizon radius $`r_{h,n}`$ for each given value of $`f_h`$. Conversely, for given (admissible) value of $`\gamma `$ one obtains a sequence of values of $`f_h`$ and $`r_h`$.

### 4.2. Linear eigenvalue problem

Liebling has found empirically the upper critical value $`\gamma _03`$ for the existence of static solutions 444In the notations of Ref.  $`\eta ^{}\sqrt{3/(8\pi )}.`$. In this section we find a theoretical ground for this limit. Actually, we find analytically a sequence of critical values $`\gamma _n,`$ $`n=0,1,\mathrm{}`$, such that for $`\gamma <\gamma _n`$ there exist static configurations with the field magnitude $`f(u)`$ changing its sign $`n`$ times.

In the literature one can find only an analysis for $`f(u)>0`$. Our numerical integration of Eqs. (41)–(43) shows that, in addition to solutions with monotonically growing $`f(u)`$ (which exist for $`\gamma <\gamma _0=3`$, see Fig. 2a), for $`\gamma <\gamma _1=2/3`$ there are also regular solutions with $`f(u)`$ changing sign once, see Fig. 2b. For $`\gamma <\gamma _2=0.3`$ there are solutions with two zeros of $`f(u)`$ (Fig. 2c), etc. All these solutions have a horizon, and the absolute value of $`f`$ on the horizon $`|f_{h,n}|=|f_n(\mathrm{})|`$ is a decreasing function of $`\gamma `$, vanishing as $`\gamma \gamma _n0`$, see Fig. 3.

As $`\gamma \gamma _n`$, the function $`f(u)`$ vanishes in the whole range of $`u`$, and it is this circumstance that allows us to find the critical values $`\gamma _n`$ analytically. In a close neighborhood of $`\gamma _n`$ the field $`f(u)`$ within the horizon is small, $`f^21`$, so that Eq. (41) reduces to a linear equation with given background functions $`F_0`$ and $`F_\mathrm{\Omega }`$, corresponding to the de Sitter metric (35). In terms of the dimensionless spherical radius $`r`$, Eq. (41) takes the form

$$\frac{d}{dr}\left[r^2\left(1\frac{r^2}{r_h^2}\right)\frac{df}{dr}\right](2r^2)f=0,$$
(49)
where $`r_h=\sqrt{12/\gamma }`$ is the value of $`r`$ on the horizon. The boundary conditions are

$$f|_{r=0}=0,|f(r_h)|<\mathrm{}.$$
(50)
Nontrivial solutions of (49) with these boundary conditions exist for a sequence of eigenvalues $`\gamma =\gamma _n`$, $`n=0,1,2,\mathrm{},`$ and the corresponding eigenfunctions $`f_n(r)`$, regular in the interval $`0rr_h`$, are simple polynomials:

$$f_n(r)=\underset{k=0}{\overset{n}{}}a_k\left(\frac{r}{r_h}\right)^{2k+1}.$$
(51)
Substituting (51) into (49), we find the eigenvalues

$$r_{h,n}^2=2(2n+1)(n+2),\gamma _n=\frac{3}{(n+1/2)(n+2)},$$
(52)
and the recurrent relation

$$a_k=a_{k1}\frac{(2k1)(2k+2)r_{h,n}^2}{(2k+1)(2k+2)2},k=1,2,\mathrm{},$$
(53)
allowing one to express all $`a_k,`$ $`k=1,2,\mathrm{},n,`$ in terms of $`a_0`$. Eq. (49) is linear and homogeneous, so $`a_0`$ is an arbitrary constant555At $`\gamma \gamma _n0`$, the general equation (41) has the same solution as (49), with $`a_01`$. To find the dependence $`a_0(\gamma )`$ one has to take into account the next terms nonlinear in $`f`$.. For fixed $`n`$ the coefficients $`a_k`$ in (51) are

$$a_k=a_0\underset{i=1}{\overset{k}{}}\frac{(2i1)(2i+2)r_{h,n}^2}{(2i+1)(2i+2)2},n>0,1kn.$$
(54)
The case

$$n=0,r_{h,0}=2,f_0(r)=a_0r/r_{h,0}$$

gives a monotonically growing function $`f(u)`$ in a close vicinity of $`\gamma =\gamma _0=3`$, see Fig. 2a. Thus the upper limit $`\gamma _0=3`$ for the existence of static monopole solutions, previously found numerically by Liebling , is now obtained analytically.

The case

$$n=1,r_{h,1}=3\sqrt{2},f_1(r)=a_0\frac{r}{r_{h,1}}[1\frac{7}{5}(\frac{r}{r_{h,1}})^2]$$

describes the function $`f(u)`$, changing its sign once, at $`\gamma `$ close to $`\gamma _1=2/3`$, see Fig. 2b. The case $`n=2,`$ $`\gamma _2=3/10`$, $`r_{h,2}=2\sqrt{10}`$,

$$f_2(r)=a_0\frac{r}{r_{h,2}}\left[1\frac{18}{5}\left(\frac{r}{r_{h,2}}\right)^2+\frac{99}{35}\left(\frac{r}{r_{h,2}}\right)^4\right]$$

gives the field function $`f(u)`$ changing its sign twice (Fig. 2c).

For $`n1`$ the function $`f_n(r)`$ rapidly oscillates:

$$f_n(r)=a_0\frac{\mathrm{cos}[r_{h,n}\mathrm{arcsin}\sqrt{{\displaystyle \frac{r^22}{r_{h,n}^22}}}\sqrt{2}\mathrm{arcsin}({\displaystyle \frac{\sqrt{2}}{r}}\sqrt{{\displaystyle \frac{1(r/r_{h,n})^2}{12/r_{h,n}^2}}}\left)\right]}{\sqrt[4]{r^2(r^22)\left[1(r/r_{h,n})^2\right]}}.$$

However, this semiclassical formula is not valid near the left turning point666Recall that in view of the substitution (40) the distances are measured in the units $`(\sqrt{\lambda }\eta )^1`$. $`r=\sqrt{2}`$, see dashed curve on Fig. 4. Its applicability range is $`1r<`$ $`r_{h,n}2n`$, $`n1`$.

We have not met so far regular monopole configurations with the field function $`f(u)`$ changing its sign. It seems that this is their first presentation.

### 4.3. Solutions with monotonically growing $`f(u)`$

As is clear from the aforesaid, the interval $`0<\gamma <3`$ of existence of nontrivial solutions with monotonically growing $`f(u)`$ splits into two qualitatively different regions, separated by $`\gamma =1`$.

In the interval $`0<\gamma <1`$ the solutions have the spatial asymptotic (37) and, according to our general classification, belong to class (a0). The spherical radius $`r(u)=\mathrm{e}^{F_\mathrm{\Omega }(u)}`$ varies from zero to infinity, $`f(u)`$ grows from zero to unity, $`A(u)`$ decreases from unity to its limiting positive value (compare with (37))

$$A|_r\mathrm{}=\frac{1\gamma }{\alpha ^2},\alpha =\frac{dr}{d\rho }|_\rho \mathrm{}=1\frac{\gamma }{2}_0^{\mathrm{}}f^{}{}_{}{}^{2}(\rho )r(\rho )d\rho ,$$
(55)
and the energy integral (9) diverges.

In the interval $`1<\gamma <3`$, the solutions with monotonically growing $`f(u)`$ belong to class (a1). Instead of a spatial asymptotic, there is a horizon and a KS cosmology outside it. The functions $`A`$ and $`r`$ inside and outside the horizon are presented in Fig. 5 for $`\gamma =2`$.

In the presence of a global monopole the cosmological expansion is slower than the de Sitter one (36). As $`\tau \mathrm{}`$, the radius $`r(\tau )`$ grows linearly, while $`A`$ tends to the negative constant value $`(\gamma 1)/\alpha ^2`$.

Within the horizon $`f(u)`$ monotonically grows from zero at $`u=\mathrm{}`$ to a value $`f_h=f_{h,0}(\gamma )`$ on the horizon, $`u\mathrm{}`$, see Fig. 2a. The $`f`$ value on the horizon $`f_{h,0}`$ as a function of $`\gamma `$ decreases from unity at $`\gamma =1`$ to zero at $`\gamma =\gamma _0=3`$, see Fig. 3. The integral (9) taken over the static region converges, and we can conclude that at $`1<\gamma <3`$ the gravitational field is strong enough to suppress the Goldstone divergence and to localize the monopole. At $`\gamma >3`$ gravity is probably so strong that it restores the high symmetry of the system.

Outside the horizon the field $`f`$ as a function of the proper time $`\tau `$ grows from $`f_{h,0}`$ on the horizon to unity at $`\tau \mathrm{}`$. Introducing the proper radial length $`l`$ inside the horizon by the relation $`dl=d\rho /\sqrt{A}`$, one can ascertain that the functions $`f(l(\rho ))`$ at $`\rho <h`$ and $`f(\tau (\rho ))`$ at $`\rho >h`$ are two parts of a single smooth curve, see Fig. 6.

When the parameter $`\gamma `$ is close to its critical value ($`\gamma =1`$), separating the (a0) and (a1) branches of the solution, i.e., for

$$0<\gamma 11.$$
(56)
one can find analytically the horizon radius $`r_{h,0}`$ and the scalar field value on the horizon $`f_{h,0}`$ under certain additional assumptions on the system behavior which follow from the results of numerical analysis. In particular, there is an “intermediate” region of the $`u`$ range, $`1uu_0=\mathrm{const}`$, where the first term $`f^{\prime \prime }`$ in the scalar field equation (41) is very small whereas the function $`\mathrm{e}^{2(F_0+F_\mathrm{\Omega })}`$ is quite large (despite the fact that this function eventually vanishes as $`u\mathrm{}`$). Therefore in this region the expression in square brackets should be small, i.e.,

$$\mathrm{e}^{2F_\mathrm{\Omega }}(1f^2)2,$$

and this relation can be used for further estimates.

The results are

$`\mathrm{ln}r_{h,0}\mathrm{ln}[1/(\gamma 1)]1,f_{h,0}1C(\gamma 1)^2,`$ (57)
where the constant $`C`$ can be found by comparison with the numerical results; our estimate is $`C0.3`$

The solution behavior in the critical regime, $`\gamma =1`$, can be characterized as a globally static model with a “horizon at infinity”  since $`A0`$ as $`r\mathrm{}`$.

The fact that monotonic solutions with horizons are absent for $`\gamma <1`$ becomes clear from an analysis of the inflection point $`u=u_{\mathrm{inf}}`$ of the function $`F_\mathrm{\Omega }(u)`$. A horizon, if any, corresponds to $`u\mathrm{}`$ where $`F_\mathrm{\Omega }`$ remains finite. Since it behaves logarithmically as $`u\mathrm{}`$, there is (at least one) inflection point, where the second-order derivative is zero, and from Eq. (43) we have

$$1\gamma f^2\frac{\gamma }{4}\mathrm{e}^{2F_\mathrm{\Omega }}(f^21)^2=0,u=u_{\mathrm{inf}}.$$

This is a quadratic equation with respect to $`1f^2`$, whence

$$1f^2=2\mathrm{e}^{2F_\mathrm{\Omega }}\left(1\pm \sqrt{1\frac{\gamma 1}{\gamma }\mathrm{e}^{2F_\mathrm{\Omega }}}\right).$$
(58)
A monotonically growing function $`f(u)`$ corresponds to greater values of $`f`$, i.e., to the “minus” branch of (58) (as is confirmed by numerical results). But then the r.h.s. of (58) is negative for $`\gamma <1`$, leading to $`f^2>1`$, which cannot happen since $`f^2=1`$ is the maximum attainable value for the solutions under study. Thus, for $`\gamma <1`$, all solutions with monotonically growing $`f(u)`$ belong to class (a0), possess a spatial asymptotic with a solid angle deficit and a divergent field energy.

Numerical integration confirms these conclusions. The different behavior of $`F_\mathrm{\Omega }(u)`$ for $`\gamma <1`$ and $`\gamma >1`$ is shown in Fig. 7.

### 4.4. Solutions with $`f(u)`$ changing its sign

For $`\gamma <\gamma _1=2/3`$ there are solutions with the function $`f(\varphi )`$ changing its sign once, see Fig. 2b. For $`\gamma <\gamma _2=0.3`$ there are solutions where $`f(\varphi )`$ changes its sign twice, see Fig. 2c, etc. Unlike the monotonic solutions discussed in Sec. 4.3, all of them possess a horizon, and, in agreement with the general inferences of Sec. 4.1, they belong to class (c1). This means that, beginning with a regular center, the spherical radius $`r(\rho )`$ first grows, then passes its maximum $`r_{\mathrm{max}}`$ at some $`\rho _1`$ and then decreases to zero at finite $`\rho =\rho _2`$ which is a singularity. The horizon occurs at some $`\rho =h<\rho _2`$, which can be greater or smaller than $`\rho _1`$, but in any case the singularity takes place in a T-region and is of cosmological nature. The dependence $`r(\rho )`$ before and after the horizon is a single smooth curve, see Fig. 8a.

Beyond the horizon, $`|A(\tau )|`$ as a function of the proper time $`\tau `$ of a comoving observer grows from zero at $`\tau =0`$ (the horizon) to infinity at $`\tau \tau _s=\tau (\rho _2)`$ (the singularity), see Fig. 8b. The scalar field magnitude $`|f(\tau )|`$ beyond the horizon first grows, then slightly varies around unity. Approaching the singularity, $`f(\tau )`$ changes its sign and finally $`|f(\tau )|\mathrm{}`$ as $`\tau \tau _s`$, see Fig. 8c.

## 5. Conclusion and discussion

We have performed a general study of the properties of static global monopoles in general relativity. We have shown that, independently of the shape of the symmetry breaking potential, the metric can contain either no horizon, or one simple horizon, and in the latter case the space-time global structure is the same as that of de Sitter space-time. Outside the horizon the geometry corresponds to homogeneous anisotropic cosmological models of KS type, where spatial sections have the topology $`\times 𝕊^2`$. In general, all possible solutions can be divided into six classes with different qualitative behavior. This classification is obtained without any assumptions about $`V(\varphi )`$. Solutions with given $`V(\varphi )`$ contain some of these classes, not necessarily all of them. This qualitative analysis gives a complete picture of what can be expected for global monopole systems with particular symmetry breaking potentials.

Our analytical and numerical analysis for the particular case of “Mexican hat” potential confirmed the previous results of other authors concerning the configurations with monotonically growing order parameter. Among other things, we have obtained analytically the upper limit $`\gamma _0=3`$ for the existence of static monopole solutions, previously found numerically by Liebling . We have also found and analyzed a new family of solutions with the field function $`f`$ changing its sign, which we have not met in the existing literature.

Of particular interest can be the class (a1) solutions with a static nonsingular monopole core and a KS cosmological model outside the horizon. Its anisotropic evolution is determined by the functions of the proper time $`|A(\tau )|`$ (the squared scale factor in the $``$ direction, $`r(\tau )`$ (the scale factor in the two $`𝕊^2`$ directions) and the field magnitude $`f(\tau )`$. For a comoving observer in the T-region, the expansion starts with a rapid growth of $`|A(\tau )|`$ from zero to finite values, resembling inflation, and ending with $`A\mathrm{const}`$ as $`\tau \mathrm{}`$. The expansion in the $`𝕊^2`$ directions, described by $`r(\tau )`$, is comparatively uniform and linear at late times, i.e., much slower than the de Sitter’s $`\mathrm{cosh}^2(\tau /r_h)`$, see (35). It should be stressed that all such models with de Sitter-like causal structure, i.e., a static core and expansion beyond a horizon, drastically differ from standard Big Bang models in that the expansion starts from a nonsingular surface, and cosmological comoving observers can receive information in the form of particles and light quanta from the static region, situated in the absolute past with respect to them. Moreover, in our case the static core is nonsingular, and it is thus an example of an entirely nonsingular cosmology in the spirit of papers by Gliner and Dymnikova .

The nonzero symmetry-breaking potential plays the role of a time-dependent cosmological constant, a kind of hidden vacuum matter. Since the field function $`f`$ tends to unity as $`\tau \mathrm{}`$, the potential vanishes, and the “hidden vacuum matter” disappears.

The lack of isotropization at late times does not seem to be a fatal shortcoming of the model for two reasons. First, if the model is applied for describing the near-Planck epoch of the Universe evolution, then, on the next stage, the anisotropy can probably be damped by diverse particle creation, and the further stages with lower energy densities may conform to the standard picture (with possible further phase transitions). Second, if we add a comparatively small positive quantity $`\mathrm{\Lambda }`$ to the potential (10) (“slightly raise the Mexican hat”), this must change nothing but the late-time asymptotic which will become de Sitter, corresponding to the cosmological constant $`\mathrm{\Lambda }`$. In our view, these ideas deserve a further study.

Evidently, the present simple model cannot be directly applied to our Universe. It would be too naive to expect that a macroscopic description based on a simple toy model of a global monopole with only one dimensionless parameter $`\gamma `$ can explain the whole variety of early-Universe phenomena. Nevertheless, it may be considered as an argument in favor of the idea that the standard Big Bang might be replaced with a nonsingular static core and a horizon appearing as a result of some symmetry-breaking phase transition on the Planck energy scale.

## 6. Acknowledgement

The authors are grateful to Acad. A.F. Andreev for a useful discussion at the seminar at P.L. Kapitza Institute of Physical Problems.