text:
# How many cosmological parameters?

## 1 Introduction

Since the release of microwave anisotropy data from the Wilkinson Microwave Anisotropy Probe (WMAP, Bennett et al. 2003), it has been widely acknowledged that cosmology has entered a precision era, with many of the key cosmological parameters being determined at the ten percent level or better. By now, a wide range of analyses have been published, uniting this dataset with other cosmological datasets such as galaxy power spectrum information from the Two degree field (2dF) survey or the Sloan Digital Sky Survey (SDSS).

While the various analyses are in broad agreement with one another, typically some differences do arise in the precise constraints, for two reasons. One is that separate analyses often use slightly different data compilations, which of course should lead to differing results, hopefully consistent within the uncertainties. However, further differences arise due to the choice of cosmological model made, usually meaning the number of cosmological parameters allowed to vary. The standard approach thus far has been to first choose the set of parameters to be varied on a fairly ad hoc basis, and then use a likelihood method to find the best-fit model and confidence ranges for those parameters. Some papers analyze several combinations of parameters, primarily with the aim of investigating how the parameter confidence ranges are affected by modifying these assumptions.

So far, however, there have been few attempts to allow the data to determine which combination of parameters gives the preferred fit to the data. This is the statistical problem of model selection, which arises across many branches of science; for example, in studies of medical pathologies, one wishes to know which set of indicators, out of many potential factors, are best suited to predicting patient susceptibility. The emphasis is usually on ensuring the elimination of parameters which play an insufficient role in improving the fit to the data available. A key tool is this area is information criteria, specifically the Akaike information criterion (Akaike 1974) and the Bayesian information criterion (Schwarz 1978). These have led to considerable advances in understanding of statistical inference and its relation to information theory; Akaike’s 1974 paper now has over 3000 citations and is the subject of a complete textbook (Sakamoto, Ishiguro & Kitagawa 1986). However, so far they seem to have had minimal application in astronomy — keyword search on the abstracts of the entire astro-ph archive yields only four journal papers (Mukherjee et al. 1998; Takeuchi 2000; Connolly et al. 2000; Nakamichi & Morikawa 2003). In this paper I will apply the information criteria to the problem of selection of cosmological parameters.

## 2 The information criteria

The information criteria have a deep underpinning in the theory of statistical inference, but fortunately have a very simple expression. The key aim is to make an objective comparison of different models (here interpretted as different selections of cosmological parameters to vary) which may feature different numbers of parameters. Usually in cosmology a basic selection of ‘essential’ parameters is considered, to which additional parameters might be added to make a more general model. It is assumed that the models will be compared to a fixed dataset using a likelihood method.

Typically, the introduction of extra parameters will allow an improved fit to the dataset, regardless of whether or not those new parameters are actually relevant.111In cosmology, a new parameter will usually be a quantity set to zero in the simpler base model, and as the likelihood is a continuous function of the parameters, it will increase as the parameter varies in either the positive or negative direction. However some parameters are restricted to positive values (e.g. the amplitude of tensor perturbations), and in that case it may be that the new parameter does not improve the maximum likelihood. A simple comparison of the maximum likelihood of different models will therefore always favour the model with the most parameters. The information criteria compensate for this by penalizing models which have more parameters, offsetting any improvement in the maximum likelihood that the extra parameters might allow.

The simplest procedure to compare models is the likelihood ratio test (Kendall & Stuart 1979, ch. 24), which can be applied when the simple model is nested within a more complex model. The quantity $`2\mathrm{ln}_{\mathrm{simple}}/_{\mathrm{complex}}`$, where $``$ is the maximum likelihood of the model under consideration, is approximately chi-squared distributed and standard statistical tables can be used to look up the significance of any increase in likelihood against the number of extra parameters introduced. However the assumptions underlying the test are often violated in astrophysical situations (Protassov et al. 2002). Further, one is commonly interested in comparing models which are not nested.

The Akaike information criterion (AIC) is defined as

$$\mathrm{AIC}=2\mathrm{ln}+2k,$$
(1)
where $``$ is the maximum likelihood and $`k`$ the number of parameters of the model (Akaike 1974). The best model is the model which minimizes the AIC, and there is no requirement for the models to be nested. Typically, models with too few parameters give a poor fit to the data and hence have a low log-likelihood, while those with too many are penalized by the second term. The form of the AIC comes from minimizing the Kullback–Leibler information entropy, which measures the difference between the true distribution and the model distribution. The AIC arises from an approximate minimization of this entropy; an explanation geared to astronomers can be found in Takeuchi (2000), while the full statistical justification can be found in Sakamoto et al. (1986) and Burnham & Anderson (2002).

The Bayesian information criterion (BIC) was introduced by Schwarz (1978), and can be defined as

$$\mathrm{BIC}=2\mathrm{ln}+k\mathrm{ln}N,$$
(2)
where $`N`$ is the number of datapoints used in the fit (in current cosmological applications, this will be of order one thousand). It comes from approximating the Bayes factor (Jeffreys 1961; Kass & Raftery 1995), which gives the posterior odds of one model against another presuming that the models are equally favoured prior to the data fitting. Although expressed in terms of the maximum likelihood, it is therefore related to the integrated likelihood.

It is unfortunate that there are different information criteria in the literature, which forces one to ask which is better. Extensive Monte Carlo testing has indicated that the AIC tends to favour models which have more parameters than the true model (see e.g. Harvey 1993; Kass & Raftery 1995). Formally, this was recognized in a proof that the AIC is ‘dimensionally inconsistent’ (Kashyap 1980), meaning that even as the size of the dataset tends to infinity, the probability of the AIC incorrectly picking an overparametrized model does not tend to zero. By contrast, the BIC is dimensionally consistent, as the second term in its definition ever more harshly penalizes overparametrized models as the dataset increases in size, and hence the BIC does always pick the correct model for large datasets. Burnham & Anderson (2002) generally favour the AIC, but note that the BIC is well justified whenever the complexity of the true model does not increase with the size of the dataset and provided that the true model can be expected to be amongst the models considered, which one can hope is the case in cosmology. Accordingly, it seems that that BIC should ordinarily be preferred. Note though that for any likely dataset $`\mathrm{ln}N>2`$, and hence the AIC is always more generous towards extra parameters than the BIC. Hence the AIC remains useful as it gives an upper limit to the number of parameters which should be included.

In either case, the absolute value of the criterion is not of interest, only the relative value between different models. A difference of 2 for the BIC is regarded as positive evidence, and of 6 or more as strong evidence, against the model with the larger value (Jeffreys 1961; Mukherjee et al. 1998).

The rather limited literature on cosmological model selection has thus far not used the information criteria, but has instead used the more sophisticated idea of Bayesian evidence (see e.g. Jaynes 2003; MacKay 2003). This compares the total posterior likelihoods of the models, obtained as a product of the Bayes factor and the prior relative likelihood. This requires an integral of the likelihood over the whole model parameter space, which may be lengthy to calculate, but avoids the approximations used in the information criteria and also permits the use of prior information if required. It has been used in a variety of cosmological contexts by Jaffe (1996), Drell, Loredo & Wasserman (2000), John & Narlikar (2002), Hobson, Bridle & Lahav (2002), Slosar et al. (2003), Saini, Weller & Bridle (2004), and Niarchou, Jaffe & Pogosian (2004).

## 3 Application to present cosmological data

### 3.1 Choice of parameters

Most of the recent work on cosmological parameters has chosen a particular parameter set or sets, and investigated parameter constraints when faced with different observational datasets. However, the information criteria ask how well different models fit the same dataset. First we need to decide which models to consider.

A useful division of parameters is into those which are definitely needed to give a reliable fit to the data, which I will call the base parameter set, and those which have proved irrelevant, or of marginal significance, in fits to the present data. The base parameter set is actually extraordinarily small, and given in Table 1. At present it seems that a scale-invariant spectrum of adiabatic gaussian density perturbations, requiring specification of just a single parameter (the amplitude), is enough to give a good fit to the data. The Universe can be taken as spatially-flat, with the dark matter, baryon, and radiation densities requiring to be specified as independent parameters. The base model includes a cosmological constant/dark energy, whose density is fixed by the spatial flatness condition. To complete the parameter set, we need the Hubble constant. Accordingly, a minimal description of the Universe requires just five fundamental parameters.222To be more precise, this base model assumes all the parameters to be listed in Table 2 are zero. Analyses may use different parameter definitions equivalent to those given here, for instance using the physical densities $`\mathrm{\Omega }h^2`$ in place of the density parameters. Further, the radiation density $`\mathrm{\Omega }_\mathrm{r}`$ is directly measured at high accuracy from the cosmic microwave background temperature and is not normally varied in fits to other data.

In addition to these fundamental parameters, comparisons with microwave anisotropy and galaxy power spectrum data require knowledge of the reionization optical depth $`\tau `$ and the galaxy bias parameter $`b`$ respectively. These are not fundamental parameters, as they are in principle computable from the above, but present understanding does not allow an accurate first-principles derivation and instead typically they are taken as additional phenomenological parameters to be fit from the data.

Complementary to this base parameter set is what I will call the list of candidate parameters. These are parameters which are not convincingly measured with present data, but some of which might be required by future data. Many of them are available in model prediction codes such as cmbfast (Seljak & Zaldarriaga 1996). Cosmological observations seek to improve the measurement of the base parameters, and also to investigate whether better data requires the promotion of any parameters from the candidate set into the standard cosmological model. Table 2 shows a list of parameters which have already been discussed in the literature, and although already rather long is likely to be incomplete.

The upper portion of Table 2 lists possible additional parameters associated with the background space-time, while the lower part contains those specifying the initial perturbations. The base cosmological model assumes these are all zero (as defined in the table), and indeed it is a perfectly plausible cosmological model that they are indeed all zero, with the sole exception of the neutrino masses, for which there is good non-cosmological evidence that they are non-zero. One should be fairly optimistic about learning something about neutrino masses from cosmology, which is why they are included as cosmological parameters. It is also possible that one day they might be pinned down accurately enough by other measurements that cosmologists no longer need to worry about varying them, and then neutrino masses will not be cosmological parameters any more than the electron or proton mass are.

It is of course highly unlikely that all the parameters on the candidate list will be relevant (if they were, observational data would have little chance of constraining anything), and on theoretical grounds some are thought much more likely than others. In most cases parameters can be added individually to the base model, but there are some dependences; for example, it doesn’t make much sense to include spectral index running as a parameter unless the spectral index itself is included. Quite a lot of the parameters in Table 2 have now been added to a base parameter set (usually not the one I have adopted here, however) and compared to observational data. There is also the possibility that the simultaneous inclusion of two extra parameters, which are unrelated, might significantly improve the fit where neither parameter separately did. This is hard to fully test as there are so many possible combinations.

### 3.2 Application to WMAP+SDSS data

I will use the results from comparison of models to WMAP plus SDSS data given in Tegmark et al. (2004, henceforth T04). Much of the analysis in that paper focusses on a simple parameter set called the ‘vanilla’ model or sometimes the ‘six parameter’ model. Confusingly, it actually features seven parameters (they do not count the bias parameter, although it is an independent fit parameter). They are not quite the set given in Table 1; the radiation density parameter is omitted for reasons I explained above, while the spectral index $`n`$ is included as an independent parameter. However $`n1`$ is not actually detected to be non-zero; its 1-sigma confidence range (table 4, column 6 of T04) is $`0.952<n<1.016`$. In light of the above discussion, we might expect that the information criteria reject the inclusion of $`n1`$ as a useful parameter, and indeed that is the case.

The $`\chi ^2`$ values quoted by T04 are derived using the WMAP likelihood code \[see Verde et al. (2003) and Spergel et al. (2003) for details\] combined with a calculation of the likelihood from the SDSS data, and are defined as $`2\mathrm{ln}`$. The total number of datapoints $`N`$ (not corrected for the number of parameters in the fit) is $`N=1367`$ (899 WMAP temperature spectrum, 449 WMAP polarization cross-correlation, 19 SDSS). I note that their Markov chains were designed to estimate confidence intervals rather than to accurately determine the precise maximum likelihood, and a modest bias might occur from the maximum being less well pinpointed the greater the model dimensionality. Once the approximate locations of the maxima are determined via a Markov Chain Monte Carlo procedure, a variation on that method could be used to determine the maximum likelihoods accurately as an additional part of the data analysis process.

As seen in the upper two rows of Table 3, both information criteria prefer the base model, with $`n`$ fixed at one, as opposed to letting $`n`$ vary. As has been remarked before, there is presently no evidence that the parameter $`n1`$ is needed to fit present data. T04 draw the same conclusion on subjective grounds, and refer to the base model as ‘vanilla lite’.

A similar argument applies to other cosmological parameters. Unfortunately the other models analyzed by T04 include variation of $`n`$ (their table 3) and so other parameters are not directly compared with the base model, but anyway the trend seen in Table 3 is clear — the more parameters included the higher the AIC and BIC as compared to the base model. The need for these additional parameters is strongly rejected by the information criteria, particularly the BIC which strongly penalizes additional parameters for a dataset of this size.333It is interesting to note that recent applications of the Bayesian evidence to cosmological model selection have also found no significant evidence against the simplest model considered (Slosar et al. 2003; Saini et al. 2004; Niarchou et al. 2004). For example, the information criteria reject the need for $`\mathrm{\Omega }_k`$ as an independent parameter, instead identifying spatially-flat models as the preferred description of the data.

## 4 Candidate parameters and statistical significance

The information criteria are clearly a powerful tool for establishing the appropriate set of cosmological parameters. How do they relate to the standard approach in cosmology of looking at confidence levels of parameter detection?

Use of fairly low confidence levels, such as 95%, to identify new parameters is inherently very risky because of the large number of candidate parameters. If there were only one candidate parameter and it were detected at 95% confidence, that certainly be interesting. However there are many possible parameters, and if one analyzes a several of them and finds one at 95% confidence, then one can no longer say that the base model is ruled out at that level, because there were several different parameters any of which might, by chance, have been at its 95% limit. As an extreme example, if one considered 20 parameters it would be no surprise at all to find one at 95% confidence level, and that certainly wouldn’t mean the base model was excluded at that confidence. Consequently the true statistical significance of a parameter detection is always likely to be less than indicated by its confidence levels (e.g. Bromley & Tegmark 2000). This issue can arise both within a single paper which explores many parameters, and in a broader sense because the community as a whole investigates many different parameters.

This is a form of publication bias — the tendency for authors to preferentially submit, and editors to preferentially accept, papers showing positive statistical evidence. This bias is well recognized in the field of medical trials (see e.g. Sterne, Gavaghan & Egger 2000), where it can literally be a matter of life and death and tends to lead to the introduction of treatments which are at best ineffectual and may even be harmful. The stakes are not so high in cosmology, but one should be aware of its possible effects. Publication bias comes in several forms, for example if a single paper analyzes several parameters, but then focusses attention on the most discrepant, that in itself is a form of bias. The more subtle form is where many different researchers examine different parameters for a possible effect, but only those who, by chance, found a significant effect for their parameter, decided to publicize it strongly.

Publication bias is notoriously difficult to allow for, as it mainly arises due to unpublished analyses of null results. However a useful guide comes from considering the number of parameters which have been under discussion in the literature. Given the list in Table 2, it is clear that, even if the base cosmological model is correct, there are enough parameters to be investigated that one should not be surprised to find one or two at the 95% confidence level.

I conclude that when considering whether a new parameter should be transferred from the candidate parameter list to the base parameter list, a 95% confidence detection should not be taken as persuasive evidence that the new parameter is needed. Because there are so many candidate parameters, a more powerful threshold is needed. The BIC provides a suitably stringent criterion, whereas this line of argument supports the view that the AIC is too weak a criterion for cosmological model selection.

Another subtle point relating to cosmological data is the inability to fully repeat an experiment. Conventionally in statistics, once a dataset has identified an effect which looks interesting (e.g. spectral index running at 95% confidence), one is expected to throw away all that data and seek confirmation from a completely new dataset. This procedure is necessary to minimize publication bias effects, and failure to follow it is regarded as poor practice. Unfortunately, for the microwave anisotropies much of the noise comes from cosmic variance rather than instrumental effects, and so remeasuring does not give an independent realization of statistical noise. For example, if one analyzes the second-year WMAP data (once it becomes available) separately from the first-year data, there will be a tendency for the same cosmological parameter values to be obtained. Finding the same outlying parameter values therefore will have less statistical significance than were the datasets genuinely independent. Even Planck data will have noise significantly correlated to WMAP data in this sense, and properly allowing for that in determining statistical significance of parameter detections would be tricky. This supports the use of information criteria for model selection, rather than parameter confidence levels.

## 5 Conclusions

Various conclusions can be drawn from the information criterion approach. Most importantly, they provide a simple objective criterion for the inclusion of new parameters into the standard cosmological model. For example, it is sometimes said that the WMAP analysis actually mildly favours a closed cosmological model, as their best-fit value is $`\mathrm{\Omega }_k=1.02\pm 0.02`$ (at 1-sigma). However, the information criteria lead to the opposite conclusion: they say that the most appropriate conclusion to draw is that the spatial curvature is not needed as a parameter, and hence it is more likely that the observations were generated in a spatially-flat Universe. That’s not to say that future observations might not show that the Universe is closed, but a much higher significance level than 1-sigma is needed before it becomes the best description of the data in hand. Similar arguments can be applied also to parameters such as running of the spectral index; even in the absence of controversy over the use of lyman-alpha forest data, it seems likely that the information criteria would reject the running as a useful parameter (I can’t test it, as the WMAP team were unable to quote a maximum likelihood due to unknown error covariances). In general, a 95% ‘detection’ of a particular new parameter cannot be taken to imply that the base model, without that parameter, is ruled out at anything like that significance.

According to the information criteria, the best current cosmological model features only five fundamental parameters and two phenomenological ones, as listed in Table 1. While there is an elegant simplicity to this model which is satisfying, such simplicity does come at a cost, because the cosmological parameters are what tells about the physical processes relevant to the evolution of the Universe. That there are so few parameters is telling us that there is very little physics that we are currently able to probe observationally. Accordingly, we should be hoping that new observational data is powerful enough to promote parameters from the candidate list to the base list; for example, we won’t be able to say anything quantitative about how cosmological inflation might have taken place unless $`n1`$, and ideally $`r`$ as well, make their way into the standard cosmological model.

The information criteria appear well suited to providing an objective criterion for the incorporation of new parameters, and have had considerable testing across many scientific disciplines. The BIC appears to be preferred to the AIC for cosmological applications. For the size of the current dataset the BIC penalizes extra parameters very strongly, indicating that a very high-significance detection is needed to justify adoption of a new parameter.

## Acknowledgments

This research was supported in part by PPARC. I thank Charles Goldie for directing me to the literature on the information criteria, and thank Sarah Bridle, Martin Kunz, Sam Leach, Max Tegmark, and the referee Håvard Sandvik for helpful discussions and comments.


format: text

text:
# Unsteady Crack Motion and Branching in a Phase-Field Model of Brittle Fracture

## Abstract

Crack propagation is studied numerically using a continuum phase-field approach to mode III brittle fracture. The results shed light on the physics that controls the speed of accelerating cracks and the characteristic branching instability at a fraction of the wave speed.

The quest for a fundamental understanding of dynamic brittle fracture has been an ongoing challenge. The traditional continuum approach to computing the speed of brittle cracks consists of solving the equations of linear elasticity with boundary conditions on the moving fracture surfaces up to the crack tip Freund . The solutions have stress fields that diverge near the tip representing a finite energy flow rate to the tip. The crack speed $`v`$ is then assumed to be uniquely determined by this energy flow rate. In this theory, all the nonlinear physics of failure inside a microscopic region around the tip $``$ the so-called process zone $``$ is buried in a phenomenological function that relates the fracture energy $`\mathrm{\Gamma }`$ of the material (energy needed to advance the tip per unit length of crack front and per unit of crack extension) and $`v`$.

Detailed experiments have been carried out to test this theory by measuring the speed of accelerating cracks with different initial lengths and different loads exp3 . The theory should be valid if a unique $`\mathrm{\Gamma }(v)`$ curve applies to the different cracks for the same material. The main difficulty in any experimental test of this theory is that the energy flux to the tip cannot be measured directly. Therefore, this flux needs to be inferred from a time-dependent solution of linear elasticity for accelerating cracks. Using an approximate solution, Sharon and Fineberg have concluded that the data for different cracks collapse on a unique $`\mathrm{\Gamma }(v)`$ curve, whereas, using an exact solution of Eshelby Esh69 , Kessler and Levine have concluded that this collapse does not occur KesLev02 .

Also of long-standing interest in brittle fracture is the existence of a dynamic instability that limits the speed of fast moving cracks. This generic instability has been seen in both experiments exp3 ; cryst ; exp1  and molecular dynamic simulations md . The fact that different amorphous materials with entirely different microscopic details (such as glass and PMMA) exhibit strikingly similar branching instabilities exp3  strongly suggests that a continuum theory may be appropriate for understanding this phenomenon. Devising such a theory, however, has proven to be difficult. Cohesive zone theories modify the boundary conditions on the stress field near the tip to take into account the short scale force between crack surfaces. Langer and Lobkovsky LanLob98 , however, have shown that these models are unsuitable for stability calculations since the results depend singularly on the details of the cohesive zone. In addition, intrinsically discrete branching instabilities in lattice models MarLiu93 ; KesLev01  seem unlikely to bear relevance to experiments in amorphous materials.

In this letter, we study dynamic brittle fracture using a recently developed continuum phase-field approach for mode III cracks Karetal01 . This approach has the chief advantage that it incorporates both the short-scale physics of failure and macroscopic linear elasticity within a self-consistent set of coupled nonlinear partial differential equations that can be solved numerically. Moreover, solutions free of discretization artifacts can be obtained. Earlier simulations of this model were limited to small system sizes owing to the extreme stiffness of the equations, and no branching instability was found even for cracks approaching the wave speed Karetal01 . Here, we briefly discuss how to overcome this stiffness and carry out simulations in systems large enough to make contact with results from the fracture community. In addition, these simulations enable us to study the onset of the characteristic branching instability at a fraction of the wave speed that turns out to be a robust feature of this phase-field model.

The basic variables of the model Karetal01  are the scalar displacement $`u(x,y)`$ perpendicular to the $`x`$-$`y`$ plane of mass points from their original positions, and the phase-field, $`\varphi (x,y)`$, which describes the state of the material. The unbroken solid, which behaves purely elastically, corresponds to $`\varphi =1`$, whereas the fully broken material that cannot support stress corresponds to $`\varphi =0`$. The total energy (kinetic plus elastic) of the system per unit length of the crack front is

$$E=𝑑x𝑑y\left[\frac{\rho }{2}\left(\frac{u}{t}\right)^2+\frac{\kappa }{2}|\stackrel{}{}\varphi |^2+hf(\varphi )+\frac{\mu }{2}g(\varphi )\left(|\stackrel{}{ϵ}|^2ϵ_c^2\right)\right],$$
(1)
where $`\rho `$ is the density, $`\stackrel{}{ϵ}\stackrel{}{}u`$ is the strain, $`f(\varphi )`$ is a double-well potential with minima at $`\varphi =1`$ and $`\varphi =0`$, $`\mu `$ is the elastic shear modulus, and $`ϵ_c`$ is the critical strain magnitude such that the unbroken (broken) state is energetically favored for $`|\stackrel{}{ϵ}|<ϵ_c`$ ($`|\stackrel{}{ϵ}|>ϵ_c`$). The function $`g(\varphi )`$ is a monotonously increasing function of $`\varphi `$ with limits $`g(0)=0`$ and $`g(0)=1`$, which controls the softening of the elastic energy as more bonds are broken.

Taking the first variations of the energy with respect to the strain and to $`\varphi `$, we obtain the stress, $`\stackrel{}{\sigma }=\delta E/\delta \stackrel{}{ϵ}`$, and the equations of motion

$`{\displaystyle \frac{\varphi }{t}}`$ $`=`$ $`\chi {\displaystyle \frac{\delta E}{\delta \varphi }},`$ (2a)
$`\rho {\displaystyle \frac{^2u}{t^2}}`$ $`=`$ $`\stackrel{}{}\stackrel{}{\sigma }.`$ (2b)
Energy is dissipated in the process zone around the crack tip where $`\varphi `$ varies rapidly in space and time. Eq. (2a) implies that the size of the process zone is $`\xi =\sqrt{\kappa /\mu ϵ_c^2}`$ and the characteristic time of energy dissipation in this zone is $`\tau =1/(\chi \mu ϵ_c^2)`$. Furthermore, by rescaling lengths by $`\xi `$, time by $`\xi /c`$, where $`c\sqrt{\mu /\rho }`$ is the shear wave speed, and $`u`$ by $`\xi ϵ_c`$ in Eqs. (2), we find that crack propagation is controlled by two dimensionless parameters $`\delta =h/(\mu ϵ_c^2)`$ and $`\beta =c\tau /\xi `$. The first determines the surface energy normalized by $`\mu ϵ_c^2\xi `$ Karetal01

$$\stackrel{~}{\gamma }\frac{\gamma }{\mu ϵ_c^2\xi }=_0^1𝑑\varphi \sqrt{1g(\varphi )+2\delta f(\varphi )}.$$
(3)
For the choices $`f(\varphi )=16\varphi ^2(1\varphi )^2`$ and $`g(\varphi )=4\varphi ^33\varphi ^4`$ here, $`\stackrel{~}{\gamma }`$ approximately doubles when $`\delta `$ increases from $`0`$ to $`2`$.

Parameter $`\beta `$ controls the importance of inertia relative to dissipation in the process zone. When $`\beta 1`$ (inertia-dominated regime), failure is rapid and crack propagation is governed by elastic energy flow. When $`\beta 1`$ (dissipation-dominated regime), failure is sluggish. The elastic displacements are quasistatic and propagation is governed by the kinetics of the failure process.

We study fracture numerically in a strip of width $`2W`$ with a fixed displacement $`u(x,\pm W)=\pm \mathrm{\Delta }`$ at the strip edges. The stored elastic energy per unit area ahead of the crack tip is $`G=\mu \mathrm{\Delta }^2/W`$ and the Griffith’s threshold load for a semi-infinite crack is $`G_c=2\gamma `$.

The equations of the model are inherently stiff due to strain localization behind the crack tip Karetal01 . The strain is localized within a narrow region whose width vanishes in the large strip limit $`W/\xi \mathrm{}`$, while the width of the $`\varphi `$ variation remains of order $`\xi `$ in this limit. Therefore, while it is possible to resolve the spatial variation of $`\varphi `$ everywhere, resolving the strain concentration behind the crack tip is computationally impractical for large strip widths. The key to overcoming this difficulty is to recognize that the model remains well-defined when the strain variation is resolved only in the process zone, but is allowed to be discontinuous on the lattice scale far behind the tip. Therefore, it is possible to simply discretize the energy $`E`$ on a square lattice and to formulate the dynamics of the model from partial derivatives of $`E`$ with respect to the variables $`\varphi `$ and $`u`$ at lattice points. One subtlety of the present simulations is that high-frequency waves are radiated behind the crack tip above a small threshold tip speed. These waves increase slightly the fracture toughness and originate from the region behind the tip where the displacement becomes discontinuous on the lattice scale. In addition, multiple reflections of these waves can cause tip oscillations at long times. This time, however, is longer than the time for propagation to reach a steady-state or branching to occur such that this feature is not problematic. Furthermore, this radiation is easily suppressed by adding a small Kelvin viscosity to Eq. (2b) with no qualitative change of the results presented here.

We have verified that the important observables converge quickly in the limit of the vanishing lattice spacing, and that we operate in the regime where these observables differ by at most 15% from their continuum limits for the lattice spacing $`dx`$ in the range $`0.3\xi `$ to $`0.4\xi `$.

To study unsteady crack propagation, we use stationary solutions of Eqs. (2) as initial conditions that correspond to cracks of different initial length at rest. For a given load $`G`$ and corresponding initial crack length $`\mathrm{}`$, these stationary solutions are found by relaxing $`u`$ with the Gauss-Seidel iteration scheme, which amounts to solving Eq. (2b) without inertia, while $`\varphi `$ is relaxed using Eq. (2a). This procedure yields stationary cracks whose initial length $`\mathrm{}`$ decreases with increasing $`G`$. We then simulate the full equations of motion with inertia to study accelerating cracks with zero initial velocity. These simulations parallels previous experiments exp3  and lattice simulations KesLev02 . To study steady-state features of crack propagation independent of initial conditions, we run long simulations in strips that are effectively infinite along the propagation direction. To keep the computations tractable, we periodically translate the fields $`u`$ and $`\varphi `$ by one lattice spacing such that the crack tip remains in the middle of a strip of length much larger than $`W`$. We have checked that the results of these “tread mill” simulations are independent of boundary effects.

To compute the fracture energy during unsteady motion, we equate $`v\mathrm{\Gamma }`$ with the expression for the energy flow rate to the tip of mode III cracks Esh69

$$v\mathrm{\Gamma }=_C𝑑C\left[\mu \dot{u}_nu+vn_x_t\left(\rho \dot{u}^2/2+\mu |\stackrel{}{ϵ}|^2/2\right)\right],$$
(4)
where $`C`$ is a closed circuit around the moving tip and $`\widehat{n}`$ is the outward normal to $`C`$. Since energy is conserved in the model where $`\varphi =1`$, we can equivalently obtain the energy flow rate to the tip by calculating the time rate of change of the total energy energy (elastic plus kinetic) in the region of the system where $`\varphi `$ is larger than a threshold value $`\varphi _c`$ arbitrarily close to unity. This quantity is precisely the energy flow rate into the process zone defined as the region where $`\varphi <\varphi _c`$.

Plots of tip speed versus tip position for cracks accelerating from rest are shown in Fig. 1 for the inertia-dominated regime $`\beta =1`$. Corresponding plots of $`\mathrm{\Gamma }`$ versus $`v`$ are shown in Fig. 2. The initial crack acceleration increases with load and the crack tip splits into two symmetric branches above a critical onset load $`G_{\mathrm{onset}}`$. Cracks for smaller loads which do not split, or split only transiently for $`GG_{\mathrm{onset}}`$, reach a steady-state propagation velocity that coincides with the steady-state velocity calculated on the tread mill (solid line in Fig. 2). The fact that $`\mathrm{\Gamma }=G`$ in steady-state is a self-consistency check of our method of calculating $`\mathrm{\Gamma }`$ for unsteady cracks.

From the long simulations on the tread-mill, three basic regimes of crack propagation can be distinguished: a stable regime for small loads, where $`v`$ is constant in time and the crack is rectilinear, an asymmetric tip-splitting (“snake”) mode for intermediate loads, where $`v`$ fluctuates in time around some average value while the crack follows a sinusoidal trajectory, and a chaotic tip-splitting regime for large loads. Examples of these regimes are shown in Fig. 3. The simulations also reveal a strong coupling between the crack tip dynamics and elastic waves that we have not explored in detail. A clear signature of this coupling is the fact that the temporal period of branching in the snake mode is very close to the period $`4W/c`$ of the first harmonic standing wave of the strip. The transient coupling of the tip to higher frequency waves is also seen for higher loads in the chaotic branching regime.

The steady state velocity $`v`$ saturates at some value $`v_c`$. Rectilinear cracks cannot propagate faster than $`v_c`$. We conclude that off-axis branching in the present model is due to the absence of steady-state crack solutions above a critical speed, rather than a linear instability of solutions that exist up to the wave speed, as in the lattice models KesLev01 . The dependence of $`v_c`$ on the parameters of the model is shown in Fig. 4. The maximum crack speed is well defined in the inertia dominated limit $`\beta 0`$. It grows monotonically with the scaled surface energy $`\stackrel{~}{\gamma }`$ (Eq. (3)) and has a minimum at $`\delta =0`$ of $`v_c0.41c`$. At this speed the linear elastic field in the unbroken material around the crack tip is quasi-isotropic Yoffe . We therefore conjecture, along the lines of Gao Gao96 , that tip blunting which leads to velocity saturation and ultimately to tip-splitting is due to the relativistic contraction of stress fields in the nonlinear process zone where the sound speed is small due to the softening of the material.

The smallest value of $`v_c0.41c`$ in the $`\beta 0`$ limit is consistent with the calculation of Adda Bedia Add04  which shows that tip splitting is energetically possible for mode III fracture for speeds above $`0.39c`$ with a maximum angle between symmetric branches of 80 degrees, which is about 10-15 percent larger than the angle in our simulations. It should be noted, however, that this calculation Add04  considers a steady-state crack that stops abruptly and splits into two daughter cracks, whereas cracks decelerate before splitting in our simulations. In addition, the onset of branching depends here on the short-scale parameters $`\beta `$ and $`\stackrel{~}{\gamma }`$. Therefore, it cannot be predicted from purely energetic considerations.

To interpret our results for accelerating cracks, it is useful to first derive an analytic expression for the speed of steady-state cracks close to the Griffith threshold. Energy balance for steady-state propagation implies that the stored elastic energy ahead of the crack in excess of twice the surface energy must be dissipated in the process zone. This implies that

$$v(GG_c)=\frac{dE}{dt}=\chi ^1𝑑\stackrel{}{x}\left(\frac{\varphi }{t}\right)^2,$$
(5)
where the second equality follows from Eqs. (2) with no energy flux through the boundaries. In a coordinate system ($`x^{}=xvt,y^{}=y`$) translating with the tip, $`\varphi /t=v\varphi /x^{}`$. Thus the steady-state velocity of rectilinear crack propagation is

$$\frac{v}{c}\frac{2\stackrel{~}{\gamma }(G/G_c1)}{\beta I},$$
(6)
where $`I𝑑x^{}𝑑y^{}\left(\varphi /x^{}\right)^2`$ is a dimensionless integral factor of order unity that depends on the profile of $`\varphi `$ for a stationary crack at $`GG_c`$. This prediction is consistent with the simulation results that show a steeper increase of velocity with load for smaller $`\beta `$.

The fact that the $`\mathrm{\Gamma }v`$ plots for different accelerating cracks in Fig. 2 do not superimpose on the steady-state curve clearly shows that the speed of these cracks is not uniquely determined by the energy flow rate to the tip. In our model, accelerating cracks require more energy per unit length of advance than steadily propagating cracks. This extra energy is reversibly stored in the process zone to be either consumed by branching or radiated away later. This effect is exacerbated for large loads and causes the crack to decelerate before branching. This deceleration is marked by the overshoot of the velocity plotted as a function of tip position (Fig. 1) and the fact that the corresponding $`\mathrm{\Gamma }v`$ plot becomes double-valued (Fig. 2). A qualitatively similar overshoot of the velocity was observed experimentally for very large accelerations in brittle fracture of glass exp3 , albeit for much larger systems than in the present simulations.

To estimate when crack acceleration should be important in our model, we compute the change of the velocity over a distance comparable to the process zone. Crack’s acceleration should be irrelevant when this change is small compared to the wave speed, or $`\xi d(v/c)/dx1`$. For cracks accelerating from rest, let us assume that Eq. (6) gives the instantaneous speed of the accelerating crack if $`G`$ is taken to mean the instantaneous elastic energy release rate $`\mathrm{\Gamma }`$ which varies as $`d\mathrm{\Gamma }/dx(GG_c)/W`$ as the crack extends into the stressed strip. This implies that the acceleration is small when

$$\frac{W}{\xi }\frac{GG_c}{G_c\beta }.$$
(7)
We find that acceleration remains important in our simulations even for $`W/\xi =100`$ when $`\beta `$ is of order unity. Since in experiments $`G`$ can be up to a hundred times larger than $`G_c`$ and $`\beta =c\tau /\xi `$ can be much smaller than unity, acceleration may be important even in macroscopic strips several orders of magnitude larger than the process zone size. We conclude that the validity of the continuum theory of brittle fracture depends not only on the ratio of the system size to the process zone size, but also on the importance of inertia relative to bond breaking dissipation measured here by $`\beta `$.

This conclusion does not conflict with the fact that this theory was recently validated by lattice simulations in the limit of zero dissipation KesLev02 . In these simulations, the velocity jumps discontinuously to a fraction of the wave speed and the slope of the velocity-load relation is finite for velocities larger than this fraction. Therefore, $`1/\beta `$ should be replaced by $`d(v/c)/d(G/G_c)1`$ in the above estimate leading to Eq. (7). In contrast, for the velocity range from zero to a fraction of the wave speed pertinent to experiments, the continuum theory will break down for accelerating cracks in the $`\beta 0`$ limit of the present model for any finite system size. The crack dynamics is independent of $`\beta `$ and well-defined in this limit.

Clearly, the incorporation of more realistic microscopic mechanisms of failure in the present phase-field approach and the extension to mode I remain needed to make contact quantitatively with experiments.

This research is supported by U.S. DOE Grant No. DE-FG02-92ER45471.


text:
# Spin Polarization of Quasi Two-Dimensional Hole Systems

(January 6, 2004)

## Abstract

In quasi two-dimensional (2D) hole systems with an effective spin $`j=3/2`$, heavy hole–light hole splitting results in a quantization of angular momentum perpendicular to the 2D plane. The spin polarization of quasi 2D hole systems due to an in-plane magnetic field $`B_{}`$ thus competes with the heavy hole–light hole splitting. As a result, the spin polarization of hole systems is very different from the spin polarization of $`j=1/2`$ electron systems. In particular, it is shown that the spin polarization of quasi 2D heavy hole systems can change its sign at a finite value of $`B_{}`$.

A spin splitting of the energy levels of the electrons in a solid can be caused by an external magnetic field $`B`$, the Zeeman splitting, and we can also have a spin splitting due to the inversion asymmetry (IA) of the crystal structure Dresselhaus (1955) or the confining potential Bychkov and Rashba (1984). The IA-induced spin splitting can be ascribed to an effective magnetic field win . Intimately related with the spin splitting of the energy levels is the spin orientation or spin polarization of the wave functions. These quantities are of great importance for our understanding of the fundamental properties of electron systems. Recently, the interest in this subject has been renewed because it has been noted that the spin polarization of quasi two-dimensional (2D) systems due to a magnetic field $`B_{}`$ parallel to the 2D plane affects the apparent metallic behavior of these systems Kravchenko and Sarachik (2004). Furthermore, the spin polarization is an important parameter for possible applications in the field of spintronics Wolf et al. (2001). We distinguish here between the *spin orientation* which is a property of the individual (occupied or empty) states, and the *spin polarization* which is the average spin orientation of all occupied states. In quasi 2D systems the spin polarization has been studied by applying a magnetic field $`B_{}`$ parallel to the 2D plane Smith and Stiles (1972); Okamoto et al. (1999); Tutuc et al. (2002); Pudalov et al. (2002); Zhu et al. (2003); Tutuc et al. (2003); Papadakis et al. (2000); Tutuc et al. (2001); Proskuryakov et al. (2002); Noh et al. (2003). If the Zeeman energy splitting due to $`B_{}`$ becomes sufficiently large, the minority spin subband is completely depopulated. In such a situation, electron systems are fully spin-polarized Tutuc et al. (2002); Pudalov et al. (2002); Zhu et al. (2003); Tutuc et al. (2003).

Recently, there has been considerable interest in the spin polarization of quasi 2D hole systems Papadakis et al. (2000); Tutuc et al. (2001); Proskuryakov et al. (2002); Noh et al. (2003). Hole systems in the uppermost valence band of many common semiconductors such as Ge and GaAs are different from electron systems due to the fact that the hole states have an effective spin $`j=3/2`$ Luttinger (1956). Subband quantization in 2D hole systems yields a quantization of angular momentum with $`z`$ component of angular momentum $`m=\pm 3/2`$ for the so-called heavy hole (HH) states and $`m=\pm 1/2`$ for the light hole (LH) states. Often it is assumed that the spin degree of freedom of HH and LH states in quasi 2D hole systems behaves similar to the spin of $`j=1/2`$ electron states. This assumption appears to be justified when the HH-LH splitting is large. However, we will show here that such an approach can be misleading. The important point is that the quantization axis of angular momentum that is enforced by HH-LH splitting is pointing perpendicular to the plane of the quasi 2D system. On the other hand, in general the effective Hamiltonians for IA-induced spin splitting and Zeeman splitting at $`B_{}>0`$ tend to orient the spin vector parallel to the plane of the quasi 2D system win . The energy scale for IA-induced spin splitting and Zeeman splitting are typically much smaller than the HH-LH splitting. Thus it has been noted that for HH states the IA-induced spin splitting Winkler (2000) and the Zeeman splitting at $`B_{}>0`$ van Kesteren et al. (1990) are suppressed, i.e., in a model with axial symmetry, the IA-induced spin splitting and the Zeeman splitting are higher-order effects proportional to $`k_{}^3`$ and $`B_{}^3`$, respectively alm . Here $`k_{}`$ is the in-plane wave vector. This result reflects the fact that we can have only one quantization axis of angular momentum. A detailed discussion of HH-LH splitting can be found in Ref. Winkler (2003).

In the present work we thus want to investigate the spin orientation and spin polarization that can be achieved in quasi 2D hole systems. We will show that in quasi 2D hole systems these quantities behave fundamentally different from what is well-known about these quantities in spin-$`1/2`$ electron systems. Throughout, we will assume that only one orbital subband is occupied.

First, we briefly review for comparison the main features of spin orientation and spin polarization in electron systems. In $`j=1/2`$ electron systems at $`B=0`$ we have two degenerate states $`|`$ and $`|`$ that correspond to spin up and spin down. A finite magnetic field $`𝐁`$ gives rise to a Zeeman term $`𝝈𝐁`$ in the Hamiltonian which results in a spin splitting $`\pm B`$ of the energy levels. Here $`B`$ can be either an external magnetic field or the $`𝐤`$-dependent effective magnetic field due to a spin-orbit interaction such as the Dresselhaus Dresselhaus (1955) or Rashba term Bychkov and Rashba (1984). The symbol $`𝝈`$ denotes the vector of Pauli spin matrices. The eigenstates $`|\pm `$ of the Hamiltonian at $`B>0`$ are orthonormal linear combinations of the states $`\{|,|\}`$ with spin orientation $`𝐒_++|𝝈|+`$ parallel to $`𝐁`$ and $`𝐒_{}|𝝈|`$ antiparallel to $`𝐁`$. Furthermore, it is easy to show that for normalized states $`|\pm `$ the spin orientations $`𝐒_\pm `$ are always unit vectors win . We see here that the spin orientation vectors $`𝐒_\pm `$ of electron states always follow the magnetic field $`𝐁`$. The spin polarization at $`B>0`$ is obtained by averaging over the spin orientation of all occupied states. In the presence of an external magnetic field $`𝐁_{}`$ parallel to the 2D plane, the system will be fully spin polarized if the Zeeman energy splitting is large enough to completely depopulate the minority spin subband Tutuc et al. (2002). If $`𝐁_{}`$ is an effective magnetic field describing the spin splitting in inversion-asymmetric systems, the $`𝐤_{}`$ dependence of $`𝐁_{}`$ is such that the spin polarization is zero win .

Hole states with an effective spin $`j=3/2`$ are characterized by the $`4\times 4`$ Luttinger Hamiltonian Luttinger (1956). The corresponding wave functions are four-component spinors. In quasi 2D systems we get HH and LH states as discussed above. At $`B=0`$ the HH and LH states are each two-fold degenerate. We can choose the eigenstates

$`|\mathrm{\Psi }_{\alpha \pm }^{\mathrm{HH}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{e}^{i𝐤_{}𝐫_{}}}{2\pi }}\xi _\alpha ^{\mathrm{HH}}(z){\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}1\\ 0\\ 0\\ \pm 1\end{array}\right)`$ (1e)
$`|\mathrm{\Psi }_{\alpha \pm }^{\mathrm{LH}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{e}^{i𝐤_{}𝐫_{}}}{2\pi }}\xi _\alpha ^{\mathrm{LH}}(z){\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}0\\ 1\\ \pm 1\\ 0\end{array}\right).`$ (1j)
Here we have assumed that the Luttinger Hamiltonian is expressed in a basis of $`j=3/2`$ angular momentum eigenfunctions in the order $`m=3/2`$, $`1/2`$, $`1/2`$, and $`3/2`$; and the quantization axis of these basis functions is perpendicular to the 2D plane. In Eq. (1) the symbol $`\alpha `$ denotes the subband index. For simplicity we have neglected HH-LH mixing due to a nonzero in-plane wave vector $`𝐤_{}`$ and we also neglected the $`𝐤_{}`$ dependence of the envelope functions $`\xi _\alpha (z)`$. These aspects are fully taken into account in the numerical calculations discussed below. One can easily convince oneself that, in contrast to $`j=1/2`$ electron systems, one can only find a linear combination of the degenerate HH states $`|\mathrm{\Psi }_{\alpha +}^{\mathrm{HH}}`$ and $`|\mathrm{\Psi }_{\alpha +}^{\mathrm{HH}}`$ that are eigenstates of $`J_z`$, but there are no such linear combinations that represent eigenstates of $`J_x`$ or $`J_y`$. Here $`J_x`$, $`J_y`$, and $`J_z`$ are angular momentum matrices for $`j=3/2`$. This fact clearly illustrates the frozen angular momentum of HH states.

A finite in-plane magnetic field $`B_{}`$ gives rise to a mixing of HH and LH states. In first order degenerate perturbation theory we obtain the following expressions for the eigenstates of quasi 2D holes in an in-plane field $`B_{}`$ in $`x`$ direction

$`|\mathrm{\Psi }_{\alpha \pm }^{\mathrm{HH}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{e}^{i𝐤_{}𝐫_{}}}{2\pi }}\xi _\alpha ^{\mathrm{HH}}(z){\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}1\\ \sqrt{3}𝒦\sqrt{3}𝒢\\ \sqrt{3}𝒦\sqrt{3}𝒢\\ \pm 1\end{array}\right)`$ (2e)
$`|\mathrm{\Psi }_{\alpha \pm }^{\mathrm{LH}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{e}^{i𝐤_{}𝐫_{}}}{2\pi }}\xi _\alpha ^{\mathrm{LH}}(z){\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}\sqrt{3}𝒦\pm \sqrt{3}𝒢\\ 1\\ \pm 1\\ \pm \sqrt{3}𝒦+\sqrt{3}𝒢\end{array}\right),`$ (2j)
where

$`𝒦`$ $``$ $`{\displaystyle \frac{\kappa \mu _\mathrm{B}B_{}}{E_\alpha ^hE_\alpha ^l}},`$ (3a)
$`𝒢`$ $``$ $`\gamma _2{\displaystyle \frac{2m_0\xi _\alpha ^{\mathrm{HH}}|z^2|\xi _\alpha ^{\mathrm{LH}}}{\mathrm{}^2}}{\displaystyle \frac{(\mu _\mathrm{B}B_{})^2}{E_\alpha ^hE_\alpha ^l}}`$ (3b)
are positive constants, $`\kappa `$ is the isotropic $`g`$ factor and $`\gamma _2`$ is the second Luttinger parameter in the Luttinger Hamiltonian Luttinger (1956), $`\mu _\mathrm{B}`$ is the Bohr magneton, $`m_0`$ is the mass of free electrons, $`E_\alpha ^h`$ ($`E_\alpha ^l`$) is the energy of the $`\alpha `$th HH (LH) subband at $`B_{}=0`$, and we have used the gauge $`𝐀=(0,zB_{},0)`$. In our notation, the HH “+” states ($`\alpha =1`$) belong to the majority spin subband that remains occupied for large $`B_{}`$. For simplicity, we have neglected the small terms proportional to the anisotropic $`g`$ factor $`q`$. We also neglected the coupling to neighboring subbands $`\alpha ^{}\alpha `$. We note that $`𝒦w^2`$ and $`𝒢w^4`$, where $`w`$ is the width of the quasi-2D layer.

Similar to electron states the spin orientation of hole states is given by win

$$𝐒_{\alpha \pm }=\frac{2}{3}\mathrm{\Psi }_{\alpha \pm }|𝐉|\mathrm{\Psi }_{\alpha \pm }.$$
(4)
Here $`𝐉=(J_x,J_y,J_z)`$ is the vector of angular momentum matrices for $`j=3/2`$ and the prefactor $`2/3`$ ensures that $`|𝐒_{\alpha \pm }|1`$. For a magnetic field in $`x`$ direction only the $`x`$ components of the vectors $`𝐒_{\alpha \pm }`$ are nonzero. (As discussed below we also get a $`y`$ component if the $`𝐤_{}`$ dependence of $`|\mathrm{\Psi }_{\alpha \pm }`$ is taken into account.) Using Eq. (2) we obtain in first order of the HH-LH splitting $`E_\alpha ^hE_\alpha ^l`$

$`S_{\alpha \pm }^{\mathrm{HH}}`$ $`=`$ $`2𝒦2𝒢`$ (5a)
$`S_{\alpha \pm }^{\mathrm{LH}}`$ $`=`$ $`\pm 2/3+2𝒦\pm 2𝒢.`$ (5b)
For typical values of $`B_{}`$ the second and third term in Eq. (5b) are much smaller than $`2/3`$. Thus, similar to electron states at $`B_{}>0`$, the spin orientations $`S_{\alpha +}^{\mathrm{LH}}`$ and $`S_\alpha ^{\mathrm{LH}}`$ of the Zeeman-split LH states are essentially antiparallel to each other and independent of $`B_{}`$. For small fields $`B_{}`$ the Zeeman-split HH states have the *same* spin orientation antiparallel to $`𝐁_{}`$ which is “compensated” by a term of opposite sign in Eq. (5b). The spin orientations $`S_{\alpha \pm }^{\mathrm{HH}}`$ vanish for $`B_{}0`$. We see here that HH states are fundamentally different from $`j=1/2`$ electron states where the spin orientation vectors of the spin-split states are always antiparallel to each other, and they have always length one, independent of the magnitude of $`B_{}`$ win . According to Eq. (5a) the spin orientation vectors of HH states are the smaller in magnitude the larger the HH-LH splitting $`E_\alpha ^hE_\alpha ^l`$, which is opposite to what one expects according to the picture that the spin of hole systems with a large HH-LH splitting behaved similar to the spin of $`j=1/2`$ systems. Once again, these results reflect the frozen angular momentum of hole states due to the HH-LH splitting. For larger fields $`B_{}`$ the term $`𝒢`$ quadratic in $`B_{}`$ dominates over the linear term $`𝒦`$ in Eq. (5a). For the HH “+” states, we thus obtain a sign reversal of the spin orientation $`S_{\alpha +}^{\mathrm{HH}}`$ at a finite value of $`B_{}`$. The term $`𝒦`$ is proportional to the hole $`g`$ factor $`\kappa `$ whereas $`𝒢`$ is proportional to the Luttinger parameter $`\gamma _2`$ characterizing the orbital motion of the holes. We can thus attribute the sign reversal of $`S_{\alpha +}^{\mathrm{HH}}`$ to a competition between the effect of the motion of the orbital and the spin degree of freedom. We remark here that usually only the first HH subband is occupied in 2D hole systems. The sign reversal of $`S_{\alpha +}^{\mathrm{HH}}`$ as a function of $`B_{}`$ is thus an important feature of these systems. It depends sensitively on the width $`w`$ of the quasi-2D layer.

The functional form of the Zeeman energy splitting in an in-plane magnetic field $`B_{}`$ is very different from the expressions for the spin orientation in Eq. (5). We emphasize here that the spin orientation depends on the eigenfunctions of the Schrödinger equation whereas the Zeeman splitting depends on the eigenvalues. In second order perturbation theory the Zeeman energy splitting due to an in-plane field $`B_{}`$ in $`x`$ direction reads Winkler (2003)

$`\mathrm{\Delta }E_\alpha ^{\mathrm{HH}}`$ $`=`$ $`𝒵`$ (6a)
$`\mathrm{\Delta }E_\alpha ^{\mathrm{LH}}`$ $`=`$ $`4\kappa \mu _\mathrm{B}B_{}+𝒵,`$ (6b)
where

$$𝒵\kappa \gamma _2\frac{24m_0\xi _\alpha ^{\mathrm{HH}}|z^2|\xi _\alpha ^{\mathrm{LH}}}{\mathrm{}^2}\frac{(\mu _\mathrm{B}B_{})^3}{E_\alpha ^hE_\alpha ^l}$$
(7)
The product of $`\kappa `$ and $`\gamma _2`$ in Eq. (7) as well as the cubic dependence on $`B_{}`$ indicate that in contrast to the spin orientation of HH states the Zeeman energy splitting is due to the *combined* effect of the motion of the orbital and the spin degree of freedom. The expressions for the energy eigenvalues also contain terms of even order in $`B_{}`$. However, these terms do not contribute to the Zeeman splitting. We note that $`𝒵w^4`$.

If the spin splitting is caused by the $`𝐤_{}`$-dependent effective magnetic field in inversion-asymmetric systems we obtain results similar to Eqs. (5) and (6) with $`B_{}`$ replaced by $`k_{}`$. In particular, once again the spin vectors $`𝐒_{\alpha \pm }^{\mathrm{HH}}(𝐤_{})`$ of the spin-split HH states for a given in-plane wave vector $`𝐤_{}`$ are in general neither unit vectors nor they are antiparallel to each other.

We remark that previously it has proven useful to discuss spin splitting in quasi 2D HH systems by means of $`2\times 2`$ Hamiltonians acting only in the subspace of HH states Winkler (2003). On the other hand, in the present discussion it is crucial that we use the full $`4\times 4`$ Luttinger Hamiltonian to evaluate the spin orientation of quasi 2D systems. Here it should be kept in mind that the $`2\times 2`$ Hamiltonians are *effective* Hamiltonians that emerge from a perturbative unitary transformation of the full $`4\times 4`$ Luttinger Hamiltonian in order to decouple the motion of HH and LH states. To evaluate the spin orientation $`𝐒`$ we must either express $`𝐉`$ in the transformed basis or we must (re-)express the HH wave functions as four-component spinors as in Eqs. (1) and (2).

To confirm the qualitative results in Eq. (5) we show in Fig. 1 the self-consistently calculated spin orientation $`𝐒_{1\pm }^{\mathrm{HH}}`$ as a function of $`𝐤_{}`$ for a symmetric (100) GaAs-Al0.3Ga0.7As quantum well with hole density $`p=1\times 10^{11}`$ cm-2 and different values of an in-plane magnetic field $`B_{}`$ in $`x`$ direction. The left (right) column refers to the majority (minority) spin subband. The dashed lines show the Fermi contours. The numerical calculations follow Ref. Winkler (2003). They fully take into account the anisotropy of the crystal structure. We see in Fig. 1 that the magnitude and the orientation of $`𝐒_{1\pm }^{\mathrm{HH}}(𝐤_{})`$ depend sensitively on $`B_{}`$ and $`𝐤_{}`$, in sharp contrast to electron states. For small values of $`B_{}`$ and $`k_{}`$, the spin orientation vectors of the HH states in both spin subbands are essentially antiparallel to $`𝐁_{}`$, in agreement with Eq. (5a). For larger values of $`k_{}`$ we can have a finite angle between $`𝐒_{1\pm }^{\mathrm{HH}}(𝐤_{})`$ and $`𝐁_{}`$. Furthermore, the spin vectors $`𝐒_{1+}^{\mathrm{HH}}(𝐤_{})`$ and $`𝐒_1^{\mathrm{HH}}(𝐤_{})`$ for a given wave vector $`𝐤_{}`$ are in general neither parallel nor antiparallel. We remark that the $`𝐤_{}`$ dependence of the spin orientation can be studied analytically using the model underlying Eqs. (2) and (5). However, we do not reproduce here the lengthy formulas.

For the same system as in Fig. 1 we show in Fig. 2(a) the self-consistently calculated depopulation of the HH minority spin subband,

and Fig. 2(b) shows the $`x`$ component of the normalized spin polarization

$$𝝅=\frac{1}{p}\underset{\alpha ,\pm }{}𝑑k_{}^2\theta [E_{\alpha \pm }(𝐤_{})E_\mathrm{F}]𝐒_{\alpha \pm }(𝐤_{}).$$
(8)
Here $`E_{\alpha \pm }(𝐤_{})`$ is the spin-split subband dispersion and $`E_\mathrm{F}`$ is the Fermi energy. If $`𝐁_{}`$ is pointing in $`x`$ direction, the $`y`$ and $`z`$ components of $`𝝅`$ are zero. In Eq. (8) we have assumed that hole energies are negative and the temperature is zero. We evaluate Eq. (8) by means of analytic quadratic Brillouin zone integration Winkler (1993). It becomes clear from Fig. 2 that the depopulation of the HH minority spin subband due to $`B_{}`$ does not imply full spin polarization of the system. In Fig. 2 the minority spin subband is completely depopulated at $`B_D28.2`$ T while $`\pi (B_D)0.14`$. In agreement with Eq. (5a) we obtain a sign reversal of $`\pi `$ at $`B_{}21.4`$ T. The derivative of $`\pi (B_{})`$ is discontinuous at $`B_{}=B_D`$.

We have checked that the surprising results reported here are observed for a large variety of 2D HH systems, although the details depend on system parameters such as the density $`p`$, the width and the crystallographic growth direction of the quasi 2D system, and the material composition of the layers. In the calculations, we have neglected many-particle effects which will be important for a quantitative analysis of spin polarization in the density regime discussed here. However, to the best of our knowledge an appropriate theory for particles with $`j=3/2`$ is presently not available. We expect that the essential features in Figs. 1 and 2 will be confirmed by such a more quantitative theory.

The author appreciates stimulating discussions with E. Tutuc and M. Shayegan. This work was supported by BMBF.


text:
# Quantum-Implemented Selective Reconstruction of High-Resolution Images

## 1 Introduction

There is growing evidence that quantum-physical systems could be harnessed for information processing, including specifically image recognition, in two ways:

* by Turing-machine-based quantum computing using quantum logic gates ;
* by quantum processing similar to those in (oscillatory) associative neural nets  (cf. ).

This paper reports how it is possible to implement successful image recognition, as verified by our simulations, in a quantum holographic process . Since the natural fundamental quantum-wave dynamics is harnessed, it allows much easier and cheaper physical realization with much bigger sizes and resolutions of images than the mainstream quantum-computing approaches .

The main contribution of this paper is not to propose a generally-better image-recognition method, but to present its powerful alternative implementation into a quantum-wave medium (sec. 2), and to demonstrate its plausibility by computational experiments (sec. 3). Quantum-net’s capacities of connectivity, parallelism, storage, associativity, speed and miniaturization are enormous, even much greater than in classical holography .

In  it was shown how the Hopfield model with real-valued (thus not necessarily binary) activities of units / neurons, having linear (not sigmoid or signum) activation function, can be transformed into a quantum-holographic procedure  where the Hebbian memory-storage is replaced by multiple self-interferences of quantum plane-waves. This translation succeeded by the simplest variable exchange of the Hopfield’s real-valued variables with the complex-valued variables changing according to sinusoids (waves) (cf. ). Thereby, all input-to-output transformations are preserved. Thus, quantum-wave image recognition functions equivalently to Hopfield’s one, only the implementation is much miniaturized enabling almost infinitely large neural-like networks.

Since the opposite translation, i.e. digitalization of holography, was done in the sixties of the $`20^{th}`$ century to get the first computational associative memories, one might wonder what is new in the present proposal. The big experimental success of classical (optical, acoustic, microwave- etc., but also X-ray-, atom-, electron-) holography  is widely known, but not also the recent fast development of quantum optics  which gave birth to quantum holography  (good ”tutorial” in ). The latter promises to implement the well-known Hopfield model and its generalizations in a completely new framework where the former obstacles (memory-capacity limitations, problems with non-orthogonality of small-size inputs producing cross-talk) are very much reduced.

## 2 Web of quantum waves

Using neuro–quantum ”isomorphisms”, presented systematically in , and ”numbers-to-waves” translation, as in , we transform the Hopfield-like associative neural net into quantum formalism (details in ):

* Quantum wave-function $`\mathrm{\Psi }`$ acts as net’s state vector $`\stackrel{}{q}`$.
* Eigen-wave-functions $`\psi ^k`$ ($`k=1,\mathrm{},P`$) act as Hopfield’s pattern-bearing eigen-vectors (attractors) $`\stackrel{}{v}^k`$.
* The quantum Green-function propagator G replaces the Hebb memory matrix J.
* Thus, sum of self-interferences $`\psi ^k\psi ^k`$ of quantum waves $`\psi ^k`$ (that’s the ”hologram” G) implements the sum of auto-correlations of input-pattern configurations $`\stackrel{}{v}^k\stackrel{}{v}^k`$ (that’s the content-addressable associative memory J). ($``$ denotes tensor/outer product.)

The Hebb-equivalent expression for elements of G (i.e., the multiple cris-cross array $`_k\psi ^k\psi ^k`$ implementing matrix J) is:

$$G_{hj}=\underset{k=1}{\overset{P}{}}\psi _h^k(\psi _j^k)^{}$$
$`(1)`$
where $`h`$ and $`j`$ denote the unit / pixel / ”neuron” / quantum point at locations $`\stackrel{}{r}_1`$ and $`\stackrel{}{r}_2`$ at time $`t`$ ($`h,j=1,\mathrm{},N`$; $`N`$ can be almost infinite). The asterisk denotes complex conjugation (or, optically, phase conjugation).

After we have succeeded to encode patterns or images as eigen-states (attractors) $`\psi ^k`$ into the quantum system prescribed by eq. (1), we can reconstruct one (say, $`k_0^{th}`$) by presenting a new input similar to the $`k_0^{th}`$ stored one:

$$\mathrm{\Psi }_h^{output}=\underset{j=1}{\overset{N}{}}G_{hj}\mathrm{\Psi }_j^{input}=\underset{j=1}{\overset{N}{}}\left(\underset{k=1}{\overset{P}{}}\psi _h^k(\psi _j^k)^{}\right)\mathrm{\Psi }_j^{input}=$$

$$=\underset{k=1}{\overset{P}{}}\left(\underset{j=1}{\overset{N}{}}(\psi _j^k)^{}\mathrm{\Psi }_j^{input}\right)\psi _h^k\psi _h^{k_0}$$
$`(2)`$
describes the resulting selective retrieval (recognition) of image $`\stackrel{}{v}^{k_0}`$ encoded in $`\psi ^{k_0}`$. See detailed analysis in  or . Eq. (2) is in analogy with . In the quantum Dirac notation, eq. (2) is, using $`(\stackrel{}{a}\stackrel{}{b})\stackrel{}{c}=\stackrel{}{b},\stackrel{}{c}\stackrel{}{a}`$:

$$\mathrm{\Psi }^{output}=𝐆\mathrm{\Psi }^{input}=(\underset{k}{}\psi ^k\psi ^k)\mathrm{\Psi }^{input}=$$

$$=\underset{k}{}\psi ^k\mathrm{\Psi }^{input}\psi ^k\psi ^{k_0}$$
$`(3).`$
We assume that we can encode images $`\stackrel{}{v}^k`$ into quantum plane waves (i.e., propagating sinusoidally-changing probability-distribution for measuring a photon 111The $`k^{th}`$ mode of the photon has momentum $`\stackrel{}{p}^k`$ and energy $`E^k`$; $`\mathrm{}`$ is Planck’s constant; $`i=\sqrt{1}`$ at location $`\stackrel{}{r}`$ at time $`t`$):

$$\psi ^k(\stackrel{}{r},t)=A^k(\stackrel{}{r},t)e^{i\phi ^k(\stackrel{}{r},t)}=A^ke^{\frac{i}{\mathrm{}}(\stackrel{}{p}^k\stackrel{}{r}E^kt)}$$
$`(4).`$
We may choose the same constant amplitudes $`A`$, so that quantum phases (delays between wave-peaks) $`\phi `$ encode the whole information. Let us take $`A=1`$ (or $`A=1/\sqrt{N}`$ for convenient quantum normalization); so, $`A_j^k=1`$ (or another constant) for all $`k,j`$. (Various possibilities of amplitude and phase modulation see in .) The image-modulated laser-beam is thus: $`\psi ^k=(e^{i\phi _1^k},e^{i\phi _2^k},\mathrm{},e^{i\phi _N^k})`$ where the number of wave-front points (wave peaks) is $`N`$.

The ”isomorphism” of  allows us to exchange variables, $`\stackrel{}{v}^ke^{i\phi }`$, giving $`\psi _j^k=e^{i\phi _j^k}`$ instead of Hopfield-like $`\psi _j^k=v_j^k`$ (or $`\psi _j^k=A_j^k`$, respectively). With this exchange in equations (1) and (2), all the information-processing mathematics, verified by computer experiments of sec. 3 and , remains valid for sinusoid-encoded images also. Namely, because eq. (1) becomes

$$G_{hj}=\underset{k=1}{\overset{P}{}}e^{i\phi _h^k}e^{i\phi _j^k}=\underset{k=1}{\overset{P}{}}e^{i(\phi _h^k\phi _j^k)}$$
$`(5),`$
eq. (2) becomes

$$e^{i\phi _h^{output}}=\underset{j=1}{\overset{N}{}}\left(\underset{k=1}{\overset{P}{}}e^{i\phi _h^k}e^{i\phi _j^k}\right)e^{i\phi _j^{input}}=$$

$$=\underset{k=1}{\overset{P}{}}\left(\underset{j=1}{\overset{N}{}}e^{i\phi _j^k}e^{i\phi _j^{input}}\right)e^{i\phi _h^k}e^{i\phi _h^{k_0}}$$
$`(6).`$
If images are almost orthogonal, a wave carrying an image (those among many stored ones which is the most similar to the newly input one) is selectively reconstructed.

There is a non-local information-exchange involved in this holographic process, which in our quantum case exploits the quantum interference web (G) itself, not its static imprint onto a crystal plate as in classical holography .

Our information-processing result can be extracted from $`\psi ^{k_0}`$ using new quantum-optical (and computer-aided) techniques for ”ensemble”-measurement of observables or for quantum-holographic-(like) wavefront reconstruction :
$``$ quantum-phase estimation / engineering,
$``$ wave-packet / wave-function reconstruction / sculpting / engineering,
$``$ (coherent) quantum control / manipulation,
$``$ quantum tomography.

## 3 Computational experiments

The purpose of these experiments is just to verify the theory on those basic aspects of the real quantum processing which can be simulated. Real quantum-physical systems provide performance (much) beyond what has been simulated by us up to now, and beyond what is simulable at all.

All experiments were done on a Pentium-4 1.3-GHz PC using the following algorithm programmed in MatLab with Image Processing Toolbox:

* $`P`$ images with index $`k`$ were encoded into $`\stackrel{}{\stackrel{~}{v}}^k=(\stackrel{~}{v}_1^k,\mathrm{},\stackrel{~}{v}_N^k)`$ where a pixel’s greyness is described by $`\stackrel{~}{v}_j^k[0,255]`$ ($`j=1,\mathrm{},N`$).
* Images were preprocessed according to: $`v_j^k=\stackrel{~}{v}_j^k\frac{1}{N}_{j=1}^N\stackrel{~}{v}_j^k`$ for each $`k,j`$. The resulting vector $`\stackrel{}{v}^k`$ was then normalized to satisfy $`_{j=1}^N(v_j^k)^2=1`$. Such normalized $`\stackrel{}{v}^k`$ are assumed to be quantum-implemented into plane-wave/laser-beam $`\psi ^k`$.
* Memory matrix, Eq. (1), was calculated (”storage stage”).
* Later, in the ”selective reconstruction stage”, a new ”query / recall-key” input was inserted. The network reacted as described in Eq. (3), or equivalently in Eq. (6). The ”query-input” was completed (if partial initially) or corrected (if corrupted) based on memorized examples, and scaled back into -range.

Quality of reconstructed image $`\stackrel{}{v}`$ was measured with Peak Signal-to-Noise Ratio (in dB; for 255 grey-values):
$`PSNR=20\mathrm{log}_{10}\left(\frac{255}{RMSE}\right)`$ ;
$`RMSE=\sqrt{\frac{1}{N}_{j=1}^N(v_j^{original}v_j^{reconstructed})^2}`$.

We found that reconstruction-quality only slightly decreased with increasing number of images stored simultaneously, and that this behavior was similar regardless of the type of stored images and the type and rate of deviation of the query-image from the stored images. For two examples see Fig. 1. Compare these plots with Fig. 2 which demonstrates examples of ”image recovery” from occlusion or noise. Indeed, the capability of selective reconstruction using memory is almost the same for different rates of degradation (occlusion or corruption with noise) of the query-image or its deviation from the original stored image(s).

The performance is indeed holography-like — a small part of a hologram contains enough information about the whole pattern (stored in the hologram, our Eq. (5), in a parallel-distributed way) that the whole can be retrieved from the small part.

As evident from Fig. 2, the image which shared most pixels with the query-image, was selected from memory-matrix and reconstructed (”recognized”), being disrupted (merely) by cross-talk due to non-orthogonality of stored images. Such results, typical for associative nets and holography, were got also in the ”mixed-set experiment” (Fig. 3). Here, 3 very-different sets of 10 different-content images, i.e. with big inter-set differences and small intra-set differences, were simultaneously stored. Cross-talk backgrounds can be seen in Figs. 3(d-f), but the reconstructed images are not disrupted too much.

## 4 Conclusions

Our simulations confirm Hopfield-net’s capabilities. The novelty of our simulations is reconsideration of Hopfield-net’s characteristics in the age of powerful computers – early simulations of the eighties had a limited resolution of patterns rather than images. Moreover, our original proposal of quantum-wave implementation opens a possibility of nets having up to almost infinite size, and of processing of huge / high-resolution images. Therefore, Hopfield-net’s cross-talk and storage limitations do not manifest (too) much for our practical needs. The first problem, cross-talk, is reduced since images with a huge number of pixels are usually almost orthogonal. The second problem, memory-capacity of the Hopfield model is limited (to $`P0.14N`$), is much reduced with possibility of ”astronomically big” $`N`$. Since databases include limited-size images, we have not yet been able to demonstrate the benefits of (quantum) huge-image processing, but they are evident even from classical holography .

Instead of plane-waves, images could be encoded into Gabor wavelets  which are similar to quantum wave-packets. Other possible (great) improvements will be studied in the future.

Our proposal is enormously superior to other proposed quantum associative memories , based on the mainstream of the quantum computing science using quantum-implemented logic gates, in the sense of simplicity, miniaturization, natural physical realizability of associative processing, memory capacity and dimensionality of data (specifically, size and resolution of images). Models  are, however, more compatible with the mainstream attempts for an universal-purpose quantum computer, not merely for associative tasks which our model masters.

ACKNOWLEDGEMENTS: M.P. thanks for discussions to Prof. H.J. Caulfield, Prof. V. Bužek, Drs. C. Trugenberger and A. Vlasov. M.P. gratefully worked as EU Marie-Curie post-doc fellow (contract no. HPMF-CT-2002-01808).


text:
# Tracing the shock precursors in the L1448–mm/IRS3 outflows

## 1. Introduction

It is well–known that the SiO emission traces the material which has been processed by shocks. The SiO profiles present very broad linewidths that are centered at velocities different from that of the ambient cloud (Martín-Pintado, Bachiller, & Fuente 1992). As silicon is heavily depleted onto dust grains in the cold quiescent gas, the SiO abundances are extremely low in dark clouds ($``$10-12; Ziurys, Friberg, & Irvine 1989; Martín-Pintado et al. 1992). Nevertheless, Lefloch et al. (1998) reported the first detection of bright and narrow (with the same linewidth of the ambient gas) SiO lines at ambient velocities associated with the molecular outflows in the dark cloud NGC 1333 (t$`>`$104 years). Lefloch et al. (1998) suggested that the SiO emission arises from the postshock equilibrium gas after the interaction of powerful protostellar jets with local dense clumps. They concluded that the narrow SiO emission could be only detected in particular objects characterized by its high degree of clumpiness. Codella, Bachiller, & Reipurth (1999) confirmed the suggestion of Lefloch et al. (1998) that the ambient SiO emission arises from decelerated postshock material. In this scenario one would not expect to detect SiO ambient emission with narrow linewidths in very young molecular outflows.

The L1448 outflow is a very young outflow with a dynamical age of only $``$3500 years (Bachiller et al. 1990). Due its youth, this object provides an excellent laboratory to study the origin of the ambient SiO emission, because the time scale is too short for the postshock gas to reach the equilibrium (Chièze, Pineau des Forêts, & Flower 1998) and to decelerate to the radial velocities of those of the ambient gas.

The formyl radical (HCO) has not been detected so far in dark clouds (e.g. L183, where the derived HCO abundance is $`<`$2$`\times `$10-10; Schenewerk et al. 1988). However, gas–phase chemistry models predict even lower HCO abundances ($``$10-11) than the previously measured upper limits to this molecule (Leung, Herbst, & Huebner 1984).

In this Letter we present the first detection of very narrow SiO and HCO emissions at ambient velocities in the L1448–mm/IRS3 molecular outflows. The narrow lines and the radial velocities of the SiO emission do not seem to be related to decelerated postshock equilibrium gas but produced by the interaction of the magnetic and/or radiative shock–precursors with the ambient pre–shocked clumpy medium.

## 2. Observations and Results

The observations were carried out with the IRAM 30–m telescope at Pico Veleta (Spain). We observed the central position of L1448-mm (the driving source of the highly collimated outflow) and L1448–IRS3, and mapped the L1448–mm region observing towards the offsets (20,0), (–20,0), (0,20), and (0,–20). Table 1 summarizes the parameters of the observed transitions. The beam size of the telescope was $``$27” at $``$90 GHz and $``$18”at $``$140 GHz. The 3 mm and the 2 mm SIS receivers were tuned to single side band with image rejections $`>`$10 dB. We observed simultaneously the SiO J=2–1, H13CO+, HN13C, and HCO species using the frequency–switching mode with a frequency throw of 3 MHz. As spectrometers we used VESPA, which provided a velocity resolution of $``$0.14 km s-1. Typical system temperatures were 90–100 K. We observed the rest of the species (H13CN, CH3OH, N2H+, and the SiO J=3–2 transition) in position–switching mode. The system temperatures and the spectral resolutions were typically of 120–350 K and 0.13–0.20 km s-1, respectively. All the line intensities were calibrated in antenna temperature.

Fig. 1 shows all the spectra measured for the L1448–mm map and L1448–IRS3, and Table 1 shows the observed parameters for all the lines measured towards the central positions in L1448–mm and L1448–IRS3. For comparison, Fig. 1 shows, in the upper panels, the profiles of the blueshifted and redshifted shocked gas with broad, $``$50 km s-1, lines (see e.g., Martín–Pintado et al. 1992). The ambient profiles are shown in detail in the central panels of Fig. 1. The SiO and HCO ambient emissions show very narrow linewidths of $``$0.5 km s-1, smaller than those of the ambient cloud molecules like H13CO+, HN13C, and N2H+ which have widths of $``$1 km s-1. There are appreciable differences between the observed profiles of SiO, HCO, and the other molecules. While all ambient molecules towards L1448–mm show double peaked profiles with velocity peaks at 4.7 and 5.2 km s-1 (see vertical dashed lines in Fig. 1), the HCO and SiO lines are single gaussian narrow profiles with peaks at 4.7 and 5.2 km s-1 respectively (Table 1). VLA images of the ambient NH3 gas show that these two components correspond to distinct molecular clumps with different spatial distribution associated with the L1448–mm source (Curiel et al. 1999).

From our limited map, we notice that the detection of the SiO emission in the 5.2 km s-1 cloud is related to an increase of the ion abundance with respect to neutrals. Fig. 2 illustrates this effect by representing the SiO line intensities as a function of the ratio between the ion (H13CO+) and the neutral (HN13C) intensities for the 4.7 and the 5.2 km s-1 components. Although the data are scarce, the tendency of the SiO emission to increase with the ion to neutral intensity ratio is clear. Since the lines are optically thin, this shows that the larger SiO abundance in the ambient cloud is associated with an enhancement of the ion abundance relative to the neutrals. Furthermore, we find that the SiO profile towards L1448–mm (0,–20) presents a substantial broadening (redshifted with respect to the ambient gas) due to material which has already entered the shock. Although CH3OH is also a neutral molecule, its emission is shifted to the 5.2 km s-1 component. As discussed in section 4, these data are consistent with the fact that the SiO and CH3OH abundances are enhanced by grain mantle destruction produced by shocks (Millar, Herbst, & Charnley 1991).

Towards L1448–IRS3, we have measured the strongest narrow SiO emission in L1448. Again, the SiO and HCO lines are narrower (with linewidths $``$0.5 km s-1) than those of the ambient molecules. Toward this position all line profiles are single peaked, but the SiO emission shows an appreciable velocity shift relative to the emissions from ambient molecules ($``$0.2 km s-1; see vertical line in Fig. 1). Like in L1448–mm, this narrow emission is shifted in the same velocity direction as the shocked material which suggests that the ′′ambient′′ SiO is related to shocks.

## 3. Abundances

Assuming a kinetic temperature of 20 K (Curiel et al. 1999), Table 2 shows the derived H2 densities and the column densities of all molecules for the central region of L1448–mm and L1448–IRS3. The H2 densities derived from the SiO J=2–1 and J=3–2 lines are a few 105–106 cm-3. If we consider a HCO+ fractional abundance of $``$10-8 (Irvine et al. 1987), the derived SiO abundance towards L1448–mm (0,0), (–20,0), and (0,–20) is $``$10-11 for the 5.2 km s-1 component, an order of magnitude lower than what Lefloch et al. (1998) found in NGC 1333. The estimated upper limit to the SiO abundance for the 4.7 km s-1 cloud is $``$10-12 at L1448–mm (0,0), of the same order as that obtained in the cold quiescent clouds TMC1 and B1 (Ziurys et al. 1989; Martín–Pintado et al. 1992). Therefore, the SiO abundance in the 5.2 km s-1 component has been enhanced by more than one order of magnitude with respect to the quiescent gas of the 4.7 km s-1 cloud.

On the other hand, HCO has not been detected for the 5.2 km s-1 component (Fig. 1), with an upper limit of $`\chi `$(HCO)$``$4$`\times `$10-12 (Table 2). For the 4.7 km s-1 cloud, we derive a HCO abundance of few 10-11.

## 4. Discussion

HCO is predicted to be produced in dark clouds through ion–molecule reactions with abundances of 10-11 (Leung et al. 1984). Our detection of HCO and the abundances derived for this molecule ($``$10-11; Table 2) are in agreement with these models, suggesting that the 4.7 km s-1 material traces the quiescent gas. The SiO/HCO abundance ratio changes by more than two orders of magnitude from the 4.7 to the 5.2 km s-1 component of the ambient gas, indicating that the formation of both molecules are related to different mechanisms.

The spatial distribution and the kinematics of the narrow SiO emission suggest that this emission is associated with the molecular outflows towards L1448–mm and L1448–IRS3. However, the lines are very narrow and with the radial velocities of the ambient gas. We can rule out the idea that the ambient SiO emission is generated in the postshock equilibrium gas since the time scales for the shocked gas to decelerate to the ambient velocities are much larger ($``$104 years) than the dynamical age, and the SiO has not been detected in the ambient 4.7 km s-1 cloud.

In the following, we explore the possibility that the SiO ambient emission is a signature of the shock precursors associated with the jets driven by the young stars. There is a tendency for the ambient SiO emission to arise only from regions where the relative abundance of the ions to neutrals is larger than in the quiescent clumps (see Fig. 1). This indicates that the ions could have slipped to the 5.2 km s-1 component from the neutrals of the 4.7 km s-1 cloud, as expected if this narrow emission were produced by the magnetic precursor associated with C–type shocks. The magnetic precursor would force the plasma to stream through the neutral gas ahead the jump front, generating a velocity decoupling between the charged and neutral fluids (Draine 1980). As the charged particles were compressed in the first stages of the magnetic precursor, the ion density would initially increase in the 5.2 km s-1 component with respect to the quiescent 4.7 km s-1 cloud (see e.g. Flower et al. 1996).

Pilipp, Hartquist, & Havnes (1990) and Pilipp & Hartquist (1994) have considered grains as a separated fluid from the ion, neutral, and electron components. The behavior of the grains changes depending on several parameters like the $`H_2`$ density, mean charge, and grain size distribution. There are cases where the charged grains are coupled to the ions, leading to a small grain–neutral velocity drift just before the shock starts. The Chièze et al. (1998) time–dependent C–shock model shows that, even for time scales of $`<`$104 years, this velocity drift generates the sufficient friction to produce the ejection of small amounts of silicon from grains into the gas phase in the 5.2 km s-1 cloud through sputtering (Flower et al. 1996; Caselli, Hartquist, & Havnes 1997). Since these time scales are rather short, we should expect that silicon is released directly in the form of SiO, producing an abundance of 10-10 or 10-11 (Flower et al. 1996; Caselli et al. 1997). The lack of a velocity shift between the SiO and the ions for the 5.2 km s-1 component also indicates that we are observing the very early phases of the interaction of the shock precursor. This is supported by the data at the position (0,-20) where the broader SiO profile (see Fig. 1) suggests that in this case a large fraction of the neutrals has already entered the shock.

On the other hand, the detection of \[Si II\] 34.8 $`\mu `$m emission in L1448–IRS3 (Nisini et al. 2000) shows the existence of a photon–dominated region (PDR) which could be related to the radiative precursor of the J–type shocks (Shull & McKee 1979; Hollenbach & McKee 1989). In fact, the C–type shock models of Chièze et al. (1998) predict a jump discontinuity in the neutral flow in the early stages for time scales $``$104 years. The UV radiation field of the radiative precursor associated with such J–shock could affect the chemistry of the pre–shocked medium. It might be possible that the ambient SiO emission would be produced by the radiative precursor through the photodesorption of Si from the grains (Schilke et al. 2001). However, the SiO/HCO abundance ratio in the 5.2 km s-1 velocity component is very different from what is observed in PDRs (Schilke et al. 2001).

In summary, the linewidths, velocities, abundances, and spatial distribution of the SiO in L1448–mm/IRS3 could be explained by the radiative and/or magnetic shock precursors associated with the jets driven by protostars. Based on all these results, we speculate that we have probably detected the first evidence of the chemical effects on the ambient cloud produced by the precursors of interstellar shocks. Higher angular resolution observations of SiO, HCO, and tracers of the ambient gas are needed to clearly establish the origin of the ambient emission of SiO in L1448.

This work has been supported by the Spanish MCyT under projects number ESP2002–01627, AYA2002–10113–E and AYA2003–02785–E. We thank J.R. Goicoechea for useful comments and suggestions. We also thank an anonymous referee for helping us to improve the manuscript.