Brian Mcwilliams
email: brian.mcwilliams07@imperial.ac.uk

Giovanni Montana
email: g.montana@imperial.com

Predictive modeling with high-dimensional data streams: an on-line variable selection approach

In this paper we propose a computationally efficient algorithm for on-line variable selection in multivariate regression problems involving high dimensional data streams. The algorithm recursively extracts all the latent factors of a partial least squares solution and selects the most important variables for each factor. This is achieved by means of only one sparse singular value decomposition which can be efficiently updated on-line and in an adaptive fashion. Simulation results based on artificial data streams demonstrate that the algorithm is able to select important variables in dynamic settings where the correlation structure among the observed streams is governed by a few hidden components and the importance of each variable changes over time. We also report on an application of our algorithm to a multivariate version of the "enhanced index tracking" problem using financial data streams. The application consists of performing on-line asset allocation with the objective of overperforming two benchmark indices simultaneously.

I. INTRODUCTION

Streaming data arise in several application domains, including web analytics, healthcare monitoring and asset management, among others. In all such contexts, large quantities of data are continuously collected, monitored and analyzed over time. Often the main objective is to make real-time predictions by using the incoming streams as covariates in a regression model. In this work, we envisage a system that imports p input and q output data streams at discrete time points. The input data vector is denoted by x t ∈ R 1×p where the subscript refers to the time stamp and the dimension p may be very large. The output y t ∈ R 1×q may also be multivariate. A common task is to recursively estimate a linear regression function of form y t = f (x t ) which can be used to make future predictions, for instance at time t + 1. Our fundamental assumption is that, at any given time, only a few selected components of x t contain enough predictive power, and only those should be actively used to build the regression model. We embrace a penalized regression approach where the unimportant predictors are excluded from the model by forcing their coefficients to be exactly zero.

There are a number of statistical problems arising in this setting which we intend to tackle in this paper. Firstly, a decision has to be made on how to select the truly important predictive components on the input data streams that best explain the multivariate response in a computationally efficient manner. Secondly, since the components of x t may be highly correlated, variable selection arises in an ill-posed problem and special care is needed in order to deal with this difficulty. As will be clear later, we take a dimensionality reduction approach. Thirdly, the relationship between input and output streams is expected to change quite frequently over time, with the frequency of change depending on the specific application domain and nature of the data. This aspect requires the development of adaptive methods that are able to deal with possible non-stationarities and the notion of concept drift, that is the time-dependency of the underlying data generating process. To the best of our knowledge, little work has been done towards the development of a methodology that resolves all these three issues in a unified framework.

The literature concerning on-line variable selection methods for predictive modeling with data streams is still in its infancy. A search of the literature yielded two relatively recent works which address this issue within a penalized regression framework. The earlier method by [START_REF] Kim | Tracking of multivariate time-variant systems based on on-line variable selection[END_REF] proposes a modification to the least angle regression (LARS) algorithm of [START_REF] Efron | Least angle regression[END_REF] for L 1penalized regression, otherwise known as the Lasso, which allows it to be updated on-line. More recently, [START_REF] Anagnostopoulos | Online optimisation for variable selection on data streams[END_REF] developed an alternative approach to on-line L 1 -penalized regression based on recursive least squares. The Lasso is solved by using the shooting algorithm, a pathwise co-ordinate optimization algorithm [START_REF] Friedman | Pathwise coordinate optimization[END_REF]. The resulting procedure is related to adaptive recursive least squares algorithms which have been routinely applied, for instance, in the domain of adaptive filtering. Finally, we note how neither approach considers a multivariate response or the issue of multicollinearity among covariates.

In this work we aim to unify these problems into a single framework by proposing an efficient incremental and sparse partial least squares (PLS) algorithm for on-line variable selection and tracking of multivariate data streams. PLS regression is an extension of the multiple linear regression model and assumes the existence of a handful of latent factors explaining the variation observed in the data. It has the favorable properties in that it can be used to deal with situations where the data is multicollinear and in problems where the response is multivariate.

II. BRIDGE PARTIAL LEAST SQUARES REGRESSION

A. Partial least squares regression

Partial least squares (PLS) regression is a method of dimensionality reduction concerned with modeling the relationship between some input data X ∈ R n×p and the response or output Y ∈ R n×q [START_REF] Hoskuldsson | Pls regression methods[END_REF]. The assumption underlying PLS is that both X and Y are generated by a small number, R, of latent factors

X = R r=1 s (r) b (r) T + E, Y = R r=1 s (r) w (r) T + F (1)
where s (r) ∈ R n×1 are the latent factors and b (r) ∈ R p×1 and w (r) ∈ R q×1 are the factor loadings of X and Y , respectively. E and F are matrices of residuals with no assumed distribution. PLS finds the latent factors such that the covariance between input and output is maximized. In order to extract the full complement of latent factors, each one must be extracted sequentially. Once a factor has been extracted, a rank one deflation of the X and Y matrices is performed by subtracting the contribution of the current factor from the data, and a new iteration begins. The PLS literature is extensive and many methods exist for extracting the latent factors (see, for example [START_REF] Rosipal | Overview and recent advances in partial least squares[END_REF] for a recent review of PLS variants). The various algorithms usually differ beyond computation of the first latent factor by how the input and output data matrices are deflated.

In this work we focus on the commonly used PLS-2 algorithm [START_REF] Rosipal | Overview and recent advances in partial least squares[END_REF]. The algorithm iteratively finds R hidden factors of X such that S = XU where S = [s (1) , ...s (R) ]. U = [u (1) , ..., u (R) ] is a matrix of weights corresponding to the direction of maximal covariance between X and Y . Because it is assumed that X and Y are related through the hidden factors and the factors underlying X are a good predictor of Y , the response can be expressed in the form of a regression model

Ŷ = X β + F (2)
where β = Ũ Ŵ are the estimated coefficients. For all values of r, we define

M (r) = X (r) T Y (r)
that is the covariance matrix between input and output streams. The weight vector u (r) is found by solving:

ũ(r) = arg max u u T M (r) M (r) T u s.t. u = 1 (3)
which is the normalized eigenvector corresponding to the largest eigenvalue of M (r) M (r) T . Alternatively, this is the first left singular vector of the singular value decomposition (SVD) of M (r) . The loading vectors for both Y and X are found by performing univariate regressions

w (r) = s (r) T Y (r) s (r) T s (r) b (r) = s (r) T X (r) s (r) T s (r) (4)
After the extraction of the first factor, in order to extract subsequent factors X and Y must be deflated by subtracting the current latent factor to give X (r+1) = X (r) -s (r) b (r) T and Y (r+1) = Y (r) -s (r) w (r) T . The same procedure is then repeated until all factors are extracted. The deflation steps above are necessary because if rank(Y ) < rank(X), then the covariance matrix M M T will be rank deficient and so the number of PLS components which can be extracted without deflation will be limited to rank(Y ). For instance, in the case of univariate response, R separate SVD computations must be performed. This is the main limiting factor in developing an efficient on-line sparse PLS algorithm that we intent to remove. In order to circumvent this problem, we propose an approach that avoids the deflation steps altogether, thus requiring only one SVD computation for the extraction of all the latent factors.

B. Bridge PLS

Bridge PLS (BPLS) [START_REF] Gidskehaug | Bridge-PLS regression: two-block bilinear regression without deflation[END_REF] is a recent development which ensures that the full complement of PLS components may be extracted in one step by adding a ridge term to the eigenvalue problem. This ensures that the covariance matrix is full rank so we are no longer limited by the rank of Y in the number of components we are able to extract. This is a very important step as it opens the possibility for efficient on-line PLS implementations.

This goal is achieved by introducing a new covariance matrix

H = X T αI + (1 -α)Y Y T X = αX T X + (1 -α)M M T ( 5
)
where 0 ≤ α ≤ 1 is a ridge parameter which adds a small constant to the diagonal of Y Y T which prevents H from rank deficient. It can be noticed that H is a weighted sum between the covariance matrix of X and the covariance matrix of X and Y . When α = 0, this yields regular PLS and setting α = 1 yields a principal components regression.

All BPLS weights are then obtained in one step by solving the following modified PLS optimization problem

Ũ = arg max U U T HU s.t. U = 1 (6)
so that Ũ = [ũ (1) , ..., ũ(R) ] are the first R eigenvectors of H. The latent factors, S are then computed as X Ũ . The corresponding Y -loadings are Ŵ = S T S -1 S T Y . The final PLS regression coefficients are given by β = Ũ Ŵ . In our experiments we set α = 10 -5 so that H full rank, yet all PLS directions may be extracted accurately after computing the SVD of H only once; see [START_REF] Gidskehaug | Bridge-PLS regression: two-block bilinear regression without deflation[END_REF] for related discussions. The computational benefits gained by removing the necessity to perform R -1 additional SVD computations in the off-line case is a saving in computation time of O(Rnp 2 ). Reducing the PLS problem to a single SVD computation provides the key element for performing variable selection in an efficient way in both off-line and on-line scenarios.

III. NEW METHODS FOR SPARSE MODELLING

A. Sparse Bridge PLS

A regularized SVD method has recently been introduced by [START_REF] Shen | Sparse principal component analysis via regularized low rank matrix approximation[END_REF] as an efficient device to perform PCA with sparse loading vectors. The method relies on the best low rank approximation property of the SVD. Briefly, this is achieved by reformulating the PCA optimization problem as a regression between X and its best low rank approximation, which is solved by an SVD application. The loading vectors are then made sparse by applying a component-wise thresholding operation.

In this section we use the sparse SVD method of [START_REF] Shen | Sparse principal component analysis via regularized low rank matrix approximation[END_REF] in order to achieve an efficient variable selection algorithm within the Bridge PLS framework. We first calculate H as in Eq. ( 5) and define the SVD of H = U DV T . The bridge PLS criterion in Eq. 6 can be written as regression by whereby the criterion to be minimized is the residual sum of squares between H and its low rank approximation, as follows:

min ũ,ṽ H -ũṽ T 2 (7)
where ũ and ṽ ∈ R p×1 are restricted to be vectors with unit norm so that a unique solution may be obtained. It is known that the product of the first left and right singular vectors, u (1) v (1) is the best rank one approximation of H. Therefore Eq. 7 is solved by setting ũ = u (1) and ṽ = v (1) . We obtain sparse loadings by imposing a penalty on ũ and removing its scale constraint. Although we could use one of a number of penalty functions (see, for instance [START_REF] Friedman | Pathwise coordinate optimization[END_REF]), in this work we concentrate on the Lasso penalty, which places a restriction on the L 1 norm of ũ. This amounts to the following optimization problem:

min ũ,ṽ H -ũṽ T 2 + γ ũ s.t. ṽ = 1 ( 8
)
where γ is a parameter which controls the sparsity of the solution. If γ is large enough, it will force some variables to be exactly zero. The problem of Eq. ( 8) can be solved in an iterative fashion by first setting ũ = u (1) and ṽ = v (1) as before. Since ũ and ṽ are rank one vectors, the Lasso penalty can be applied as a component-wise soft thresholding operation on the elements of ũ (see, for instance, [START_REF] Friedman | Pathwise coordinate optimization[END_REF]). The sparse ũ are found by applying the threshold component-wise as follows:

ũ * = sgn H T ṽ H T ṽ -γ + ṽ * = H ũ * / H ũ * (9)
We then set ũ = ũ * and ṽ = ṽ * and iteratively apply Eq. ( 9) until ũ * -ũ < τ where τ is an arbitrarily small constant.

The procedure above allows all the PLS weight vectors to be extracted and made sparse at once without the need to recompute an SVD for each dimension.

B. Incremental Sparse Bridge PLS

In this section we develop the Sparse Bridge PLS algorithm to be used for variable selection in the streaming data setting. We call the resulting algorithm incremental Sparse Bridge PLS (iSB-PLS). In this case, we no longer assume we have access to the full data matrix X ∈ R n×p . Instead the data arrives sequentially at each time point, t, as x t ∈ R 1×p . Similarly, the response arrives is observable only at discrete time points as y t ∈ R 1×q .

The main challenge with applying the sparse Bridge PLS algorithm to streaming data is implementing an efficient method to calculate and update the SVD of H. Our solution to this problem consists in using the Adaptive SIM algorithm [START_REF] Erlich | Convergences of adaptive block simultaneous iteration method for eigenstructure decomposition[END_REF], a generalization of the power method which is able to adapt to changes in the data. When a new data point x t and its corresponding response y t arrives, we update the individual covariance matrices as follows

C t = λC t-1 + x T t x t M t = λM t-1 + x T t y t ( 10
)
where λ is a forgetting factor which exponentially downweights the contribution of past data points to the current covariance matrix. The Bridge PLS covariance matrix H t of Eq. ( 5) is constructed by summing the weighted PCA and PLS covariance matrices C t and M t M T t , which leads to

H t = αC t + (1 -α)M t M T t (11)
At each time point, the estimate of the eigenvectors of the covariance matrix, H are updated by performing one iteration of the SIM algorithm as follows:

Q = H t U t-1 U t = orth(Q) ( 12
)
where the function orth(Q) ensures that the columns of the matrix Q are mutually orthogonal. We use the Gram-Schmidt orthogonalization procedure which has a computational complexity of O(pR 2 ).

Once the weight vectors U t have been updated, they are made sparse using the a modified version of the iterative regularized SVD algorithm used for Sparse Bridge PLS in Section III-A. Since our algorithm is on-line and the solution is updated when a new data point arrives, we no longer iteratively apply the thresholding operation and instead apply it directly to the current estimate of the eigenvector. The simplified sparsification process for the r th weight vector is

u * = sgn u (r) |u (r) | -γ (r) + u * = u * u * (13)
The final steps of the Bridge PLS algorithm proceed as in the off-line case. The latent vectors S are computed as S = XU . However since the number of observations is effectively one, S will be an R-vector and the Y -loadings can be computed as W = Y T S/(S T S). The sparse PLS regression coefficients are β = U W so that the regression estimate at time t is ŷt = x t U W . Algorithm (1) details the resulting iSB-PLS procedure in full. In the initialization phase, we set U 0 = [u (1) 0 , ...u (R) 0 ] = I p×R to ensure that the initial estimates of the eigenvectors are mutually orthogonal. We also initialize m 0 = 0, γ = 0, and C 0 = 0. The forgetting factor, λ is chosen to be between between zero and one. When λ = 1, no data forgetting takes place, whereas λ = 0 has the effect of setting the sample size to the present data point only. Therefore, as the values of λ get close to zero, the algorithm becomes more adaptive and the selected variables may change more often.

Initialize U = I, m 0 = 0, γ = 0, C 0 = 0; Data: Input x t and output y t at time t Result: Sparse regression coefficients r) , q (1) ← a (1) ; u (r) ← q (r) / q (r) ; γ (r) ← findRoot(u (r) );

β t at time t m t ←-λm t-1 + x T t y t ; C t ←-λC t-1 + x T t x t ; H t ←-αC t + (1 -α)m t m T t ; for r ← 1 to R do a (r) ← H t u (r) ; q (r) ← I p×p - r-1 k=1 u (k) u (k)T a (
u * ← sgn u (r) |u (r) | -γ (r) + ; u * ← u * u * ; u (r) t ← u * ; end s ← xU t ; w ← ys s T s ; β t ← U t s T ;

Algorithm 1: The iSB-PLS algorithm

The parameter γ controls the degree of sparsity. In some situations, such as in financial applications (e.g. Section IV-B), the user may wish to have direct control over the number of variables to be selected. In such a case, it is necessary to select a value of γ to induce the correct degree of sparsity in the solution. Our method of choice is a simple rootfinding algorithm: γ is assigned a value equal to the (p -θ) th largest component of |u| where p is the number of elements in u and θ is the desired number of selected variables. Applying the threshold operation with this value of γ will cause all but θ of the elements in u to become 0. The computational effort required for this sort operation is O(Rp log p) which makes it suitable for application in an on-line algorithm. In the on-line case the complexity introduced by the penalization function decreases from the off-line case as we operate only on a single data point at a time (i.e. n = 1). This makes the complexity of the penalization function at each time point O(Rp).

IV. EXPERIMENTAL RESULTS

A. Simulated data

In this section we report on a simulation experiment designed to demonstrate the performance of the iSB-PLS algorithm as an on-line variable selection method. The input is simulated by first introducing three hidden factors whose temporal evolution is governed by an autoregressive (AR) process of first order in the following way:

F t,j = δ j F t-1,j + ǫ t,j for t = 2 . . . , 400 (14)
where F t,j indicates the value of factor j at time t, starting with an arbitrary initial value at time t = 1, and independently for j = 1, 2, 3. The parameter j is the autoregressive coefficient for factor j, and we use δ 1 = 0.1, δ 2 = 0.4, δ 3 = 0.2. The error terms in each one of the three factors follow a normal distribution with variance set to 12.25 and means given by, respectively, 0, -1.5 and 1.5. Each input is generated as

x t,i = F t,j + η t η t ∼ N (0, 1) (15)
where x t,i indicates the values of data stream i at time t, for t = 1, . . . , 400 and i = 1, . . . , 60. The index j indicates that each stream depends only on a given time-varying hidden factor. Specifically, we create three groups of data streams by setting j = 1 for 1 ≤ i ≤ 20, j = 2 for 21 ≤ i ≤ 40 and j = 3 for 41 ≤ i ≤ 60.

In order to test the adaptive behavior of the iSB-PLS algorithm, we generate an univariate output by introducing time-dependent regression coefficients. Until time t = 100, all the variables associated with the first hidden factor strongly contribute to the output, and their regression coefficients are selected by sampling from a normal distribution centered at 10 and with low variance. Analogously, the variables associated with the second hidden factors have regression coefficients with mean 5 and with low variance. The variables associated with the third hidden factor are assigned zero coefficients. In order to introduce a non-stationary behavior, all the nonzero coefficients in the two groups of "active variables" are swapped at t = 101. At t = 301 until the end of the period, the first group of variables is assigned a zero coefficient and the group associated with the third hidden factor is assigned a coefficient sampled from a normal distribution centred around 10. In this way, the important predictors change over time and we expect these changes to be picked up in almost real-time by the algorithm. In this setting, we set R = 2 and the sparsity parameter γ is chosen automatically by the algorithm so that, at any given time, exactly 20 variables are selected. The forgetting factor λ is set to 0.98 to ensure a rapid adjustment when the coefficients switch while also keeping the switching frequency low to gain stability in the selected variables.

Figure 1 shows the results of a single run of this experiment. Clearly, the first PLS component is able to accurately select the most important group of variables. The second component always selects the second most important group of variables whilst mostly ignoring the group of variables selected by the first component. Neither component selects the inactive variables suggesting the algorithm is correctly able to distinguish important predictors from noise. As the coefficients switch, the algorithm only requires few data points before it detects the changes and adapts itself.

Figure 2 reports on the mean percentage of correctly selected variables in both components by the iSB-PLS algorithm in a Monte Carlo simulation consisting of 500 runs of this experiment. The solid line shows the mean percentage of correctly selected variables by the first and second PLS components. The shaded area shows the Monte Carlo error. It is clear that in the portions where the data is stationary, iSB-PLS will correctly select the important variables with very little error. In response to a change in the important factors, the percentage of correctly selected variables instantly decreases and quickly adapts to the new data. The algorithm eventually selects the correct variables after some settling time. However, during this time the variability of the result increases.

B. An application to index tracking

An example application of the iSB-PLS algorithm lies in the financial domain and is related to the index tracking problem. The objective of index tracking is to select a small portfolio of assets and determine weights, which represent a proportion of the total investment capital, so that the returns achieved by the portfolio track very closely those achieved by a benchmark index. Our application of a sparse algorithm to the portfolio selection and index tracking problem is supported by work in [START_REF] Brodie | Sparse and stable markowitz portfolios[END_REF] who propose sparse portfolios based on Lasso penalized regression. Furthermore the use of a latent factor model for index tracking is supported by evidence which suggests that the first principal component of index returns captures the market factor (see, for example [START_REF] Alexander | Sources of over-performance in equity markets: mean reversion, common trends and herding[END_REF]). Our framework unifies these two approaches by combining dimensionality reduction by projection onto latent factors with variable selection using a regularized regression. For this application, we use published data from the S&P and Nikkei indices as described in [START_REF] Beasley | An evolutionary heuristic for the index tracking problem[END_REF].

We have tested the iSB-PLS algorithm in a more involved setting where: (a) two indices (the S&P and the Nikkei) need to be simultaneously tracked, so the response is bivariate, and (b) both benchmark indices have been enhanced (see, for instance, [START_REF] Alexander | Equity indexing: Optimize your passive investments[END_REF]). "Enhanced" index tracking involves performing index tracking in the case where the target asset to be tracked are the index returns plus an additional 15% returns. The total number of available stocks is 323 and we set the portfolio size to 10. The forgetting factor is λ = 0.99 and we constrain the selected stock to be associated to the main latent factor only, so that R = 1, as in [START_REF] Alexander | Sources of over-performance in equity markets: mean reversion, common trends and herding[END_REF].

In order to assess whether our procedure selects and tracks the important variables over time, we compare its performance with the average returns obtained from a population of 1000 portfolios of the same size, with each portfolio being made of a randomly selected subset of assets. To make sure that the comparison is fair, the portfolio weights are also timevarying and are obtained by using a recursive least squares method with the same λ parameter. This comparison is made in order to determine whether the ability to update the portfolio composition in response to perceived changes in the market is really advantageous in an index tracking application.

Figure 3 shows the results of this test. It can be seen that iSB-PLS consistently overperforms both indices and selects a small portfolio achieving exactly the target annual returns of +15%. In comparison, the random portfolio underperforms the S&P index by 32.07% and the Nikkei by 8.42%. Our results suggest that the importance of certain stocks in the index is not constant over time so the ability to detect and adapt to these changes is certainly advantageous. Using a model that assumes a time-varying latent factor driving the asset returns is also advantageous in this setting, since its existence in real markets has been heavily documented in the financial literature. The bottom plot of Figure 3 is a heatmap illustrating how the make-up of the portfolio selected by iSB-PLS changes during the entire period. Specifically, it shows the existence of a few important stocks that are held for the majority of the period whereas other assets are picked and dropped more frequently throughout the period, further suggesting that it is advantageous to be able to adapt the constituents of a tracking portfolio.

V. CONCLUSIONS

In this work we have presented an on-line algorithm for variable selection in a multivariate regression context based on streaming data. As far as we are aware, this is the first such algorithm which combines dimensionality reduction and variable selection for data streams in a unified framework. From the simulation results we have shown that the algorithm is able to accurately select variables associated with the important factors underlying the data. In the case of nonstationary data where the important factors are changing, iSB-PLS is able to accurately track the changes.

iSB-PLS requires the specification of a number of parameters which are currently pre-specified by the user. The question of how to select, in an on-line and adaptive manner, the number of PLS components and the number of variables per component is an important one and we are currently working towards the development of self-tuning procedures.

There are several methods in the literature for automatically updating the individual model parameters. A mechanism for adapting the sparsity parameter, γ t at each time point was proposed by [START_REF] Anagnostopoulos | Online optimisation for variable selection on data streams[END_REF]. They achieve this by evaluating the Akaike information criterion (AIC) of the model with a value of γ t-1 , γ t-1 + c and γ t-1 -c, where c is some small constant. The value of γ t which is used at that time point is the one which minimizes the AIC. This method could be incorporated into iSB-PLS as a simple adaptive solution to the variable selection portion of the model selection problem.

A method to select the number of PLS projections online was proposed by [START_REF] Vijayakumar | Incremental online learning in high dimensions[END_REF] who use an approximation of leave-one-out cross validation. The algorithm initially sets the number of projections, R = 2 and recursively keeps track of a mean squared error term, e (r) t as a function of the number of components, using a forgetting factor in the following way e (r) t+1 = λe (r) t + (y t -ŷt ) 2 (16

)
where ŷt is the estimated response at time, t. If at time t + 1 adding a new PLS component causes a large enough reduction in error, the number of PLS components is increased. If adding the new component does not decrease the error enough, the number of PLS components is not changed. The forgetting factor λ has also been pre-selected, however a number of techniques exist for learning this parameter from the data in a streaming fashion. These techniques have been discussed in the literature concerning on-line learning of neural networks, as in [START_REF] Saad | On-Line Learning in Neural Networks[END_REF], and other time-varying processes, as in [START_REF] Niedźwiecki | Identification of time-varying processes[END_REF]. Furthermore, we are planning to apply these methods to related financial applications such as further extensions of index tracking for building market neutral portfolios and detecting market inefficiencies for algorithmic trading, as in [START_REF] Alexander | Equity indexing: Optimize your passive investments[END_REF] and [START_REF] Montana | Data stream mining for market-neutral algorithmic trading[END_REF], respectively. We are considering other applications in the field of text mining involving news feeds.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Results of test with simulated data. The top figure shows how, at any time, there are three blocks of data streams: active streams having larger (black) and smaller (gray) regression coefficients, and inactive streams (white) which only contributes to noise. Each block is related to a different hidden factor. The bottom figure shows the data streams selected on-line by each PLS component.

Fig. 3 .

 3 Fig. 3. Bivariate enhanced tracking (+15% annual returns) of the S&P and Nikkei indices using a dynamic portfolio of 10 stocks.


Adib Tavakoli

Kunal J Lulla

Tuomas Puurtinen

Ilari Maasilta

Eddy Collin

Laurent Saminadayar

Olivier Bourgeois

Specific Heat of Thin Phonon Cavities at Low Temperature: Very High Values Revealed by ZeptoJoule Calorimetry

Specific heat of phonon cavities is investigated in order to analyse the effect of phonon confinement on thermodynamic properties. The specific heat of free standing very thin SiN membranes in the low dimensional limit is measured down to very low temperatures (from 6 K to 50 mK). In the whole temperature range, we measured an excess of specific heat orders of magnitude bigger than the typical value observed in amorphous solids. Below 1 K, a cross-over in c p to a lower power law is seen, and the value of specific heat of thinner membranes becomes larger than that of thicker ones demonstrating a significant contribution coming from the surface.

We show that this high value of the specific heat cannot be explained by the sole contribution of 2D phonon modes (Lamb waves). The excess specific heat, being thickness dependent, could come from tunneling two level systems (TLS) that form in low density regions of amorphous solids located on the surfaces. We also show that the specific heat is strongly tuned by the internal stress of the membrane by orders of magnitude, giving unprecedentedly high values, making low stress SiN very efficient for energy storage at very low temperature.

The thermal properties of low-dimensional amorphous material nanowires, membranes etc... are at the heart of essential questions for many subjects such as bolometry, temperature sensing, or simply the understanding of phonon physics at very low temperature in dielectric nanosystems [1]. Even if the thermal transport in low-dimensional dielectric nanostructures has recently been the subject of numerous although difficult experiments at sub-kelvin temperatures [2][3][4][5], the specific heat, or the capacity of a body to rise in temperature following a supplied external energy, remains very poorly understood and experimentally almost unexplored [6].

At very low temperatures, two significant physical effects are expected to impact the thermal properties of thin membranes, driving away the specific heat far from the one of a conventional 3D solid material. The first is the effect of the confinement of the phonon modes in nanostructures thinner than the wavelength of phonons. In that kind of free standing membranes, it is expected that new phonon modes will appear having quadratic and not linear dispersion relations, modes whose contribution to the specific heat at low temperature is linear in T [6][7][8][9][10][11][12]. Such modes have been experimentally observed in the thermally relevant hypersonic frequency range [13]. Note that simple arguments would have naively predicted that the specific heat of a confined phonon gas in the low temperature limit should scale with T δ , δ being the effective dimension of the considered system (for a suspended membrane δ = 2 at low enough temperature).

Secondly, it is well known that the thermal behaviour of amorphous (glassy) materials differs significantly from their crystalline counterparts. For more than forty years, thermal properties of such glassy materials have been explained by the possible presence of dynamic defects acting like tunneling two-level-systems (TLS) [14][15][16]. Below a certain temperature, usually around 1 K, these defects add a sizeable contribution to the specific heat through additional degrees of freedom providing also a linear variation of c p in temperature. These dynamic defects are expected to scatter phonons, thus acting on the thermal transport of such low dimensional structures [4,5].

Consequently, very particular behaviour of the phonon specific heat is expected at very low temperature in quasi-2D self-suspended amorphous materials. Such systems should thus allow to address the question of the main mechanism governing the thermal properties of confined glassy material when it moves from a standard system described by nuumnuum elastic theory to a system dominated by TLS, or the question of how the internal characteristics of an amorphous material such as stress can influence thermal properties [17][18][19]. Apart from the fundamental aspect of differentiating between both theoretical models, since specific heat governs the characteristic time scales of low temperature thermal detectors which often employ suspended ultra-thin membranes, it is essential to solve this controversy by doing direct experiments on them.

Here, we report highly sensitive specific heat measurements of thin membranes in the quasi-2D limit. They were performed on thin silicon nitride suspended structures at sub-Kelvin temperatures using zeptoJoule calorimetry. By comparing results obtained for two different thicknesses (100 nm and 300 nm), it is shown that above 1 K the specific heat is dominated by 3D bulk contributions whereas for T < 1 K surface contributions dominate.

Moreover, even if we observe a close to linear behaviour of the specific heat as a function of the temperature below 1 K, we show that the 'standard' elastic continuum model (Lamb modes) cannot explain even the order of magnitude of the measured specific heat. The presence of TLS in the amorphous dielectric is the most probable scenario to explain this anomalously high specific heat. The absolute value of the specific heat measured in our experiment is up to five orders of magnitude larger than that measured in a crystalline bulk material, making the used SiN thin membranes the system with the largest specific heat in this range of temperatures. Specific heat of suspended thin membranes of silicon nitride has been measured down to very low temperatures using a highly sensitive technique. Various types of silicon nitride films have been used to fabricate the sensors: stoichiometric Si 3 N 4 having a high internal tensile stress of 0.95 GPa, and two non-stoichiometric silicon nitrides, one low stress (0.2 GPa) and one super low stress (< 0.1 GPa). We used 100 nm and 300 nm thick SiN films commercially available deposited by low pressure chemical vapor process on top of a Si substrate. The thicknesses have been chosen to be of the same order of magnitude as the dominant phonon wavelength to ensure that only the 2D phonon modes are occupied in the membrane. Here we define the dominant phonon wavelength λ dom to correspond to the wavelength at which the Planck distribution (spectral energy density per volume) for bosons is maximum [20]. For 3D phonons, the dominant phonon wavelength is given by:

λ dom = hv s,i 2.82k B T (1)
with v s,i the velocity of sound for mode i, h is the Planck constant and k B the Boltzmann constant; for SiN v s = 10000 m s -1 , λ dom = 100 nm around 1 K. By fabricating 100 nm thick membranes, the 2D phonon regime for thermal properties should be accessible at the lowest temperature of the experiment (below 0.1 K).

The thermal sensor is made of a self-supported SiN membrane of size 40 µm by 40 µm defining the dimensions of the quasi-2D phonon cavity suspended by eight arms (6 µm wide and 100 µm long). The long suspending arms lead to a very good thermal insulation of the membranes from the Si frames [5]. These arms serve also as mechanical support for the electrical leads of the transducers (see Fig. 1). A copper heater and a niobium nitride thermometer (NbN) are structured by laser lithography on the SiN membrane as shown in Fig. 1 a, the electrical leads are made with superconducting materials (NbTi). The temperature coefficient of resistance of the NbN calculated through α = 1 R dR dT can be higher than -1 K -1 , a particularly high value for thin film technology [21].

An example of a measurement of heat capacity is presented in Fig. 2. It is measured using an ac-calorimetry technique, a dynamic method based on modulating the temperature being efficient and highly sensitive especially at low temperatures. In recent papers [START_REF] Bourgeois | [END_REF][24][25][26],

we have demonstrated on these sensors a sensitivity for specific heat measurement of 2 ×

10 -16 J K -1 √ Hz -1
, for an oscillation of temperature of 3 mK, giving an unprecedented sensitivity in energy of ±10 -19 J √ Hz -1 , i.e. ± 100 zeptoJoule per √ Hz -1 at 100 mK. The sample holder is installed in a vacuum chamber (< 2 × 10 -6 mbar before cryopumping) on a stage regulated with a precision better than 20 µK at a temperature of 100 mK measured using a calibrated Speer carbon thermometer [21]. This set-up allows the measurement of the specific heat of 2D phonon cavities from below 100 mK to above 7 K as it can be seen in Fig. 2.

Specific heat measurements done on two sets of membranes of high-stress SiN of thicknesses 100 nm and 300 nm are presented in Fig. 3 in different units: 1-in units normalized by the surface area J m -2 K -1 highlighting the surface effects in panel (a), and 2-in the regular units J g -1 K -1 for volumetric specific heat highlighting bulk effects in panel (b). As expected, the specific heat (surface or volumetric) decreases when lowering the temperature as less and less degrees of freedom are excited. The total heat capacity signal is dominated by the SiN membrane contribution, as the heat capacities of the copper heater and the NbN thermometer are at least two orders of magnitude smaller over the whole temperature range of the experiments and will be neglected in the following discussion.

If we focus first on the low temperature part, one clearly sees that below 1 K, both thin and thick membranes have approximately the same surface specific heat (see Fig. 3a), whereas in the high temperature limit (above 1 K), it's the specific heat per unit volume which is almost identical for both types of samples (see Fig. 3b). This indicates that in the low temperature range, the specific heat of the samples seem to be governed by the surfaces, whereas at higher temperatures, the specific heat is dominated by volumetric effects. One basis for interpretation resides in comparing the thickness of the samples with the dominant wavelength of the phonons. Indeed, a possible change between a 3D volume effect and a 2D surface effect arises when the phonon dominant wavelength becomes bigger than the membrane thickness. As mentioned earlier, by using Eq. 1, we obtained that λ dom is on the order of few hundreds of nanometer around 1 K. One can conclude that at temperatures lower than that the phonon specific heat is dominated by the lowest 2D phonon modes, as their wavelength is much larger than the thickness of the sample. As the temperature increases, λ dom decreases and one recovers the standard three dimensional behaviour for the specific heat. Considering these two limits, we can infer that we expect a 2D-3D dimensional cross-over in the specific heat in this system, the cross-over temperature being below 1 K.

Then, the following questions need to be addressed: can we attribute the change of specific heat with temperature to a phonon dimensionality effect? does it appear with the correct order of magnitude? In the following, we will estimate the Debye 3D, Debye 2D and elastic continuum model (ECM) of specific heat and compare the different orders of magnitude obtained to sort out the various physical mechanisms at play.

The Debye 3D model of specific heat is characterized by its temperature variation following a cubic power law C ∝ T 3 as given by the formula [START_REF] Cleland | Foundations of Nanomechanics[END_REF]:

c Debye3D = 2π 2 k 4 B ρ15 3
2 v 3 T + 1 v 3 L T 3 (2)
where ρ is the mass density in g m -3 , is the reduced Planck constant, and k B is the Boltzmann constant, v T and v L the transversal and longitudinal speeds of sound respectively 6200 m s -1 and 10300 m s -1 for silicon nitride [7]; the specific heat c Debye3D is in J g

-1 K -1
. This model has demonstrated over the years its great ability to fit c p experimental values of crystalline electrically insulating materials [START_REF] Rosenberg | Low temperature Solid State Physics[END_REF][START_REF]Specific Heats at Low Temperatures[END_REF]. It can be adapted to 2D systems by integrating over the wavevectors parallel to the membrane surface in a planar geometry (keeping the regular longitudinal and transverse 3D bulk phonon modes):

c Debye2D = 3ζ(3)k 3 B π 2 2 v 2 T + 1 v 2 L T 2 (3)
where ζ is the Riemann zeta function. This area dependent specific heat in J m -2 K -1

does not depend on the thickness of the 2D membrane and has a lower temperature exponent as compared to the Debye 3D model, the power law being given by the dimension of the system. This Debye 2D model is compared to the experimental surface specific heat as measured on SiN membranes in Fig. 3a. Even if the Debye 2D specific heat is larger than the Debye 3D one, significant discrepancies still exist. Indeed, the measured absolute value of c p is in strong disagreement with the calculated ones, being few orders of magnitude bigger.

A refinement of the model seems necessary. This is done by taking into account the modified phonon modes present in thin membranes. In the following, the specific heat of thin free-standing membranes at low temperature is obtained through the calculation of phonon modes of a phonon cavity using the elastic continuum model [7,[9][10][11][12]. In this model, some bulk polarizations will couple together at the surfaces generating new set of vibrational modes: dilatational waves (or symmetric Lamb waves), flexural waves (or asymmetric Lamb waves) and horizontal shear waves [10,11,[START_REF] Graff | Wave motion in elastic solids[END_REF]. The modified 2D dispersion relations of these modes have important consequences on the low temperature specific heat. Again, by integrating over the wavevector parallel to the membrane surface q || , the total surface specific heat can be written as:

c ECM = A k B T 2 2π σ ∞ m=0 ∞ 0 dq || q || ( ω m,σ ) 2 exp(β ω m,σ ) [exp(β ω m,σ ) -1] 2 (4)
where β = (k B T ) -1 , σ represents the different modes and m is the summation over all the branches of a mode.

According to the calculations at low enough temperature, for membranes much thinner than λ dom , the contribution of the flexural waves becomes dominant. In contrast to the Debye 2D model, the elastic continuum theory then predicts a linear variation of specific heat (not quadratic) with temperature at very low temperature [7,9,11]. This linear dependence comes from the quadratic dispersion relation of the dominant flexural modes.

The temperature at which this 3D-2D cross-over appears depends on the thickness of the phonon cavity via the dependence on temperature of the dominant phonon wavelength. The surface specific heat, as evaluated by ECM, is shown in Fig. 3a All we have considered so far are in principle valid only for crystalline systems. As SiN thin films are amorphous and disordered, some specific thermal properties may be expected.

One of the most important properties of non-crystalline materials is that they contain a large amount of two-level systems. Two level systems are atoms or groups of atoms that have two metastable configurations in their energy landscape. At low enough temperature these atoms or groups of atoms can switch between these two states only by quantum tunneling.

An important ingredient when considering the physics of these TLS is that a uniform and broad spectral distribution for time constants or equivalently barrier heights is assumed.

The exact origin of these two level systems is still heavily debated, though they can be seen as extra degrees of freedom in the solid [19,[START_REF] Queen | [END_REF][32][33][34][35][36][37]. At low temperature, phonons can interact strongly with the TLSs, having not only consequences on the specific heat but also on the phonon thermal transport [4,5,16]. The contribution of the TLSs to the specific heat is given by [15,19,32]:

c T LS = π 2 6 n 0 ρ k 2 B T = αT (5)
with ρ the mass density in g m -3 of the material and n 0 the density of TLS; α is in J g -1 K -2 . c T LS is the specific heat linked to low energy excitation specific to glassy materials, it will dominate the c p signal at low temperature. Another contribution to specific heat in amorphous solids can come from extra non-propagating modes; this is usually modeled by an excess specific heat having a cubic power law [32]:

c ex = γ T 3 (6)
γ being a factor determined experimentally. For most disordered materials, γ lies in the range 1 to 10×10 -7 J g -1 K -4 . This excess specific heat has been observed in most disordered materials [16,32]. Non-propagating modes and TLS are both associated with low density regions specific to disordered materials where voids or floppy bonds are present.

The total specific heat of amorphous materials will be given by the relation

c tot = c T LS +c ex = αT + γT 3 .
We can now compare the contribution of these two different mechanisms in disordered solids to the measured specific heat. In Fig. 3b, we have reported our experimental data for the volumetric specific heat and plotted the theoretical contributions to c p , the excess specific heat and the specific heat due to two level systems, as obtained from Eq. 5 and Eq. 6;

the prefactors γ in Eq. 6 and the density of two level system n 0 in Eq. 5 being the adjustable parameters. We have also plotted data taken by other groups on 50 nm and 200 nm thick SiN membranes [16,[START_REF] Queen | [END_REF]. By comparing the different experimental data, we can first notice that the specific heat we observe is larger by more than two orders of magnitude than the one measured on a-SiN membranes .

In order to reproduce the experimental data using the TLS model, we have to consider the specific heat due to localized soft modes (c ex , Eq. 6) as well as the specific heat due to two level systems (c T LS , Eq. 5), taking for the adjustable parameters γ ≈ 10 -4 J g -1 K -4

and n 0 ≈ 1.5 × 10 49 J -1 m -3 . Using these parameters, we can reproduce quantitatively the experimental data as shown in Fig. 3b with the red line fit. The membrane specific heat is dominated over the whole temperature range of the experiment by the amorphous nature of the SiN: above 1 K, the specific heat is interpreted to be dominated by localized soft modes, whereas in the low temperature range, the specific heat is almost entirely due to the presence of two level systems. The density of TLS in the 100 nm thick membrane is n 0 ≈ 1.5 × 10 49

J -1 m -3 a relatively high value as for 300 nm n 0 ≈ 1.5 × 10 48 J -1 m -3 , lying slightly above the measured range in amorphous silicon n 0 ≈ 10 44 J -1 m -3 to n 0 ≈ 10 48 J -1 m -3 [32].

To confirm this scenario of TLS dominated c p , it is interesting to try to control the density of the two level systems in the material; this can be done by changing the mass density of the material by the following mechanism: when the composition of SiN is changed this affects the internal stress along with the mass density. Since the TLS are thought to result from dangling bonds arising from voids or surfaces, varying the mass density and stress should affect the TLS density and with it the specific heat. We consequently expect a higher TLS density in lower mass density materials having lower mechanical stress. Therefore, we decided to reproduce the same c p experiments on high-stress (HS), low-stress (LS) and super low-stress (SLS) amorphous silicon nitride membranes to probe the effect of various TLS density on the specific heat.

The temperature dependence of the specific heat of the different 100 nm thick SiN membranes is depicted in Fig. 4. Using Eq. 5 we can extract from our measurements an order of magnitude for the concentration of TLS in our membranes, summarized in the Table I.

This table displays the composition, the mass density, the TLS density and the stress for each composition as given by the supplier [START_REF]The SiN/Si wafers have been pruchased from Si-Mat company[END_REF]. The composition and the mass densities were determined by Rutherford Backscattering Spectrometry (RBS).

The main observation from Fig. 4 is that by decreasing the stress and the mass density between HS and SLS SiN, the specific heat is increased by two orders of magnitude.

These values are surprisingly high in comparison to what has been measured until now on amorphous dielectric materials [19,32,[START_REF] Pohl | Low Temperature Specific Heat of Glasses in Amorphous Solids[END_REF]. We suppose that such high value of the specific heat can be attributed to the presence of two level systems and localized modes for each internal stress. The unit of α is J g -1 K -2 and for γ J g is the mass density of material measured by Rutherford Backscattering Spectrometry, n 0 is the TLS density of states as extracted from the measurements using Eq. 5, and the stress for each composition, as given by the supplier [START_REF]The SiN/Si wafers have been pruchased from Si-Mat company[END_REF].

in amorphous SiN. One remarkable feature appears when inspecting Table I: it seems that there is a correlation between the density of TLS and the mass density. Indeed, there is a factor of twenty between the density of TLS in HS and SLS a-SiN x . In between these two samples, the mass density varies by more than 30%. We guess that this variation of mass density is due to the appearance of low density regions in the amorphous SiN located on the surface; it is precisely these low density regions which are favourable for the formation of TLS; this could explain the difference by orders of magnitudes of specific heat between SLS and HS SiN through a much higher TLS density. Such a behaviour has never been reported before; more investigations will be required to elucidate its true origin.

In conclusion, we have measured specific heat of thin SiN membranes from 6 K down to very low temperature (50 mK) to investigate the heat capacity of 2D phonon cavities.

In the whole temperature range, the specific heat does not exhibit a single power law as a function of temperature but a cross-over between two different regimes around 1 K. As it has been suggested first by Gusso et al. [9], this cross-over could have been correlated to a 3D-2D transition when the phonon wavelength becomes of the order of the thickness of the membrane. However, first, the position of the cross-over in temperature is not in agreement with the ECM model, it is at too high temperature and second the absolute value of the measured specific heat is much larger than that predicted from our ECM analysis in terms of surface specific heat. This is ruling out that, in SiN membranes, the contribution of the Lamb modes to c p can explain such very high values of specific heat. Rather, we have

FIG. 1 .

 1 FIG. 1. Suspended membrane-based zeptoJoule nanocalorimeter. (a) A top view of the 2D nanocalorimeter made of amorphous SiN that is functionalized with two transducers: a heater and a thermometer. (a), (b), and (c) are the angled views of the nanocalorimeter that is suspended by eight suspension beams.

FIG. 2 .

 2 FIG. 2. Heat capacity of a 100 nm thin membrane made out of high-stress Si 3 N 4 . The measurement of the heat capacity is shown over a temperature range of 0.07 K to 7 K. The noise has been measured at 100 mK as shown in inset. The full width at half maximum shows that the measurement is done with a sensitivity in energy of ± 100 zeptoJoule.

FIG. 3 .

 3 FIG. 3. The specific heat of thin membranes made of high-stress Si 3 N 4 . (a) The heat capacity normalized by the surfaces area of the membranes, in comparison to the predicted specific heat for 2D phonon cavity for two sets of membranes of thicknesses 100 nm and 300 nm. The green lines are the prediction of elastic continuum theory, and the red line is the 2D Debye model. The two black straight lines are guides to the eyes for the temperature power laws T and T 2 (b) The heat capacity normalized by volume, same data as in (a). The blue line is the excess of specific heat C exfrom Eq. 6, the green dashed line is the specific heat of TLSs C T LS ∝ T , and the red line is the total estimated specific heat from the sum of Eq. 5 and Eq. 6 for an amorphous solid from which the TLS density is extracted. The black line is the predicted specific heat from the de Debye 3D model. The green and orange dots are the specific heat of amorphous SiN of thickness 50 nm and 200 nm respectively, measured by Queen et al[START_REF] Queen | [END_REF] plotted for comparison.

  only above 0.5 K due to its very low value. The measured c p (T ) and the calculated one by ECM have very dissimilar absolute values showing that the elastic continuum model do not give a good account of what has been measured by few orders of magnitude.

FIG. 4 .

 4 FIG. 4. The temperature dependence of the specific heat for different internal stress of SiN x . The samples have all the same geometry with a thickness of 100 nm. The specific heat measurements have been done for two sets of membranes of each kind of SiN: blue, red and green dots are respectively the specific heat of high-stress (HS), low-stress (LS) and super low-stress (SLS) amorphous silicon nitride. In (a) the specific heat vs T is plotted in a log-log scale and in (b) c p /T versus T 2 to highlight the α (intercept with y axis) and γ (the slope of the linear fit) values

TABLE I .

 I Si 43 N 57 Si 46 N 54 Si 48 N 52 Parameters for high-stress (HS), low-stress (LS), and super low-stress (SLS) SiN x . ρ

	SiN x	HS	LS	SLS
	composition ρ (g cm -3 )	3.2	2.9	2.4±0.2
	n 0 (J -1 m -3 ) 1.5×10 49 3.6 ×10 49 3.3 ×10 50
	stress (GPa)	0.85	0.2	<0.1

-1 K -4 .

demonstrated that the exceptionally high magnitude of the specific heat likely results from the presence of a large density of two level systems, as expected for the amorphous disordered systems as the ones used in this work (SiN membranes). Our modeling indicates that above 1 K the specific heat is dominated by an excess bulk-like contribution, possibly due to nonpropagating modes, and by two-level systems at lower temperatures. Using different SiN compositions, we have also shown that by changing the density of two-level systems through the mass density of the SiN material, unprecedentedly high values of specific heat can be reached [19,[START_REF] Pohl | Low Temperature Specific Heat of Glasses in Amorphous Solids[END_REF][START_REF] Perez-Castaneda | [END_REF]. The formation of TLS is enhanced by to the presence of low-density regions in the SiN, and the presence of large surfaces due to the membrane geometry of the samples. These results will have a significant impact on the performances of low temperature sensors.

The authors thank the technical support provided by the Nanofab, Pôle Capteur, Pôle électronique, Pôle de cryogénie facilities in Institut Néel. We thank Cyril Bachelet and Claire Marrache-Kikuchi from Irene Joliot-Curie Lab-CNRS for their help in the RBS


Stefan Lalkovski

P Van Isacker

IBM-1 calculations towards the neutron-rich nucleus 106 Zr

Keywords:

The neutron-rich N = 66 isotonic and A = 106 isobaric chains, covering regions with varying types of collectivity, are interpreted in the framework of the interacting boson model. Level energies and electric quadrupole transition probabilities are compared with available experimental information. The calculations for the known nuclei in the two chains are extrapolated towards the neutron-rich nucleus 106 Zr.

I. INTRODUCTION

In the last decade the neutron-rich nuclei in the 40 ≤ Z ≤ 50 region have attracted both theoretical and experimental attention. They were extensively studied via spontaneous or induced fission reactions. Nuclei from this region of Segré chart exhibit vibrational, transitional, and rotational types of collectivity. The neutron-rich palladium and ruthenium nuclei, for example, show a typical transitional behavior while molybdenum isotopes exhibit a vibrational to near-rotational evolution [1,2]. Such changes in the degree of collectivity are even stronger in the zirconium isotopic chain [1,2] where the structure evolves from near-vibrational in 98 Zr to rotational-like in 104 Zr, the latter being the most deformed nucleus of the neutron-rich Z = 40 isotopes. As the next zirconium isotope 106 Zr is an N = 66 mid-shell nucleus, an even larger deformation can be expected.

A detailed analysis of some spectroscopic observables, such as the ratio R 4/2 of excitation energies of the first 2 + and 4 + levels [1,2] or the amplitude of even-odd staggering in the γ band [2], shows that mid-shell effects in palladium isotopes occur in 114,116 Pd, i.e. two and four neutrons beyond N = 66. In the ruthenium isotopes these effects arise precisely at N = 66 while in the molybdenum chain they occur in 106 Mo, i.e. two neutrons before mid-shell. Furthermore, it was shown that the degree of collectivity decreases in the heavier N = 66, 68 molybdenum isotopes 108,110 Mo. In spite of the fact that the experimental information on the heavy zirconium nuclei is poor, it is clear that the degree of collectivity increases towards 104 Zr [1] which has 64 neutrons and is the heaviest zirconium isotope known to date. The shift of the mid-shell in the various isotopic chains was explained in Ref. [1] as resulting from the filling of the neutron h 11/2 and proton g 9/2 intruder orbitals.

Self-consistent calculations show that in the region of neutron-rich nuclei approaching the neutron-drip line, the single-particle shell structure may significantly change due to the diffuseness of the neutron density [3]. Also, as widely discussed in the literature, due to a predicted weakening of the spin-orbit force, new magic numbers may be expected in neutron-rich nuclei. For in-stance, in a study of N = 82 isotones, evidence for shell quenching was presented in Ref. [4]. On the other hand, Jungclaus et al. [5] argued that from the isomeric decay study in 130 Cd there is no direct evidence for the shell quenching in the region of the heavy cadmium nuclei.

The lower-Z N = 82 nuclei, however, have not yet been studied experimentally and the suppression of the shell effects in this region is still an open question. Properties of these nuclei were predicted within a relativistic Hartree-Fock Bogoliubov (HFB) approach [6]. The calculations, performed with different parametrizations, show that the N = 82 shell gap persists in the heavy palladium and ruthenium nuclei but that a weakening of this gap is expected in the zirconium isotopes. For example, the two-neutron separation energies from 118 Zr to 126 Zr, obtained via HFB calculations with the SkP interaction, show no discontinuity at N = 82 but a rather smooth behavior. As discussed above, nuclear collective behavior is expected to enhance towards the mid-shell as a result of the increase in valence particle number. If, however, the HFB+SkP scenario for the zirconium isotopes turns out to be valid and N = 82 vanishes or is weakened in the neutron-rich region, then the determination of valenceparticle number is fraught with ambiguity, leading to a different behavior of nuclear collectivity in this isotopic chain. The aim of the present paper, therefore, is to predict the spectroscopic properties concerning the N = 66 mid-shell zirconium nucleus 106 Zr, based on the assumption that N = 82 remains a magic number in the heavy zirconium isotopes. Comparison with results of future experiments on this nucleus will then reveal whether this hypothesis is borne out or not.

II. MODEL

Often the only available information about neutronrich nuclei comes from prompt γ-ray spectroscopy. The observed levels are grouped into ∆J = 2 sequences, corresponding to the ground-state band, and/or into ∆J = 1 sequences, based on a 2 + 2 level which usually is interpreted as a quasi-γ band [7]. Such structures naturally appear in the framework of the interacting boson model (IBM) [START_REF] Iachello | The Interacting Boson Model[END_REF] which has been shown to be successful in the description of nuclear collective properties.

The IBM in its first version, known as IBM-1, is based on the assumption that nuclear collectivity can be expressed in terms of s and d bosons [START_REF] Iachello | The Interacting Boson Model[END_REF]. The model Hamiltonian is constructed from a set of 36 operators, bilinear in the boson creation and annihilation operators and generating the U(6) Lie algebra. Dynamical symmetries occur if the Hamiltonian can be written as a combination of invariant (or Casimir) operators of specific subalgebras of U(6) [START_REF] Iachello | Lie Algebras and Applications[END_REF] and three such cases occur, namely the spherical vibrational limit U(5), the deformed limit SU(3), and γ-soft limit SO (6). These dynamical symmetries generate energy spectra with states that are labeled by the irreducible representations of the algebras in the respective chains that reduce the dynamical algebra U(6) into the symmetry algebra SO(3) of rotations in threedimensional space. The different limits thus correspond to nuclei with distinct collective properties. For example, in the SO( 6) limit the ground-state and γ bands lie in the same representation while in the SU(3) limit they are in different ones, leading to forbidden or weaker inter-band transitions in the latter case.

In the A ≈ 110 nuclei of interest here strong transitions between the two bands are observed indicating that none of these nuclei can be interpreted in the exact SU(3) limit, but rather a transitional behavior should be expected. This can be achieved in the IBM-1 by the use of the full Hamiltonian which reads [START_REF] Iachello | The Interacting Boson Model[END_REF]

Ĥ = ǫ d nd + κ Qχ • Qχ + κ ′ L • L + c 3 T3 • T3 + c 4 T4 • T4 , (1)
where nd

≡ d † • d and Lµ ≡ √ 10[d † × d] (1)
µ are the d-boson number and the angular momentum operators, respectively. Furthermore, the quadrupole operator is defined as Qχ

µ ≡ [d † × s] (2)
µ + χ[d † × d] (2)
µ , while the last two terms in the Hamiltonian involve the operators T3,µ

≡ [d † × d] (3) µ and T4,µ ≡ [d † × d] (4)
µ . In the present work there is no need for the rotational term L • L and we take κ ′ = 0 throughout. The total number of bosons N = n s + n d is taken as half the number of valence particles or holes, counted from the nearest closed-shell configuration, following the prescription of Ref. [START_REF] Casten | Nuclear Structure from a Simple Perspective[END_REF]. As such, IBM-1 calculations are indirectly related to the underlying shell structure. Finally, in the consistent-Q formalism [START_REF] Warner | [END_REF] the operator for electric quadrupole transitions is introduced as Tµ

(E2) = e b Qχ µ
, where e b is the boson effective charge and Q χ µ is the quadrupole operator with the same parameter χ as in the Hamiltonian (1).

III. RESULTS AND DISCUSSION

Since the addition or subtraction of a few nucleons may change significantly a nuclear spectrum, two nuclei with different numbers of valence neutrons and protons but with the same total number of valence nucleons- and hence the same number of bosons-may display different types of collectivity. In adjusting the parameters of the Hamiltonian (1) to spectra observed in a given region of the nuclear chart, it is therefore not sufficient to assume them to be constant for all nuclei since that would lead to identical spectra for nuclei with the same N . Instead, some dependence of the parameters on the separate neutron and proton boson numbers N ν and N π should be imposed. To avoid the complexity of determining the correct functional dependence on N ν and N π for each of the Hamiltonian parameters, we have followed the simpler procedure of separately considering chains of isotopes, isotones, or isobars that end with the nucleus of interest, 106 Zr, and fitting each chain with an independent set of parameters. This method avoids the problem of having identical spectra for constant boson number N (such nuclei belong to different chains) and, in addition, since the extrapolation to 106 Zr obtained from the different chains may vary, it gives an idea on the possible error of the prediction.

In the present application to 106 Zr, the N = 66 isotonic chain includes γ-unstable 112 Pd, more triaxial 110 Ru, and near-rotational 108 Mo. The second set of nuclei belongs to the A = 106 isobaric chain [16] extending from vibrational-like 106 Pd to near-rotational 106 Mo. Unfortunately, the zirconium isotopic chain does not allow the determination of a unique set of model parameters since only yrast data are available for the neutron-rich isotopes 100,102,104 Zr. Consequently, only the isotonic and isobaric chains are considered in the fit. Nuclei in a given chain are distinguished by the number of bosons (i.e., particle or hole pairs), counted from the nearest closed-shell configurations which are Z = 28 or Z = 50 for the protons, and N = 50 or N = 82 for the neutrons. Figure 1 summarizes known nuclei in the neighbourhood of 106 Zr. Nuclei with known excited levels are distinguished (full lines) from those of which only ground-state properties are known (dashed lines). Furthermore, the nuclei used in the fit are high-lighted with thick lines.

The model parameters ǫ d , κ, c 3 , and c 4 , listed in Table I, are obtained after a least-squares fit to the experimental level energies. Since its effect on the excitation energies is weak, the parameter χ is varied independently to give a fair description of E2 transition rates in 106 Pd within the consistent-Q formalism, and other nuclei are calculated with the same value. Although within a given chain parameters may vary with boson number, any allowed variation does not markedly improve the quality of the fit. Consequently, to keep the number of parameters to a minimum, they are taken constant for a given chain and any structural evolution within that chain is due to the changing boson number N .

The quality of the fit is summarized by the root-meansquare (rms) deviation between the experimental and calculated level energies. The rms values obtained in the two independent fits are given in Table I. The experimental [12,13,14,15,16] and theoretical energies are listed in Tables II and III for the N = 66 and A = 106 chains, respectively. The nucleus 106 Pd is known in greater detail and is shown separately in Table IV. In both the isotopic and isobaric chains, the observed energy of the 2 + 1 state decreases when approaching the nucleus 106 Zr. The same trend is observed for the calculated energies as a function of boson number. The experimental energy of the γ-band head also decreases with increasing boson number and a similar evolution is observed for the 2 + 2 theoretical energies. According to the systematics of the less neutron-rich nuclei, the 2 + 1 level in 106 Zr should thus lie around 140 keV while the 2 + 2 level should be in the region 600-650 keV. Furthermore, the two different fits give approximately the same γ-band energies. Even though the two fitted parameter sets are close (see Table I), the N = 66 set gives rise to a slightly faster increase of the ground-state band energies than is the case for the A = 106 set.

The displacement of odd-spin with respect to evenspin γ-band levels can be quantified with the following function:

S(J) = E(J) - (J + 1)E(J -1) + JE(J + 1) 2J + 1 , (2)
where E(J) is the energy of the level with angular momentum J. The experimental and theoretical staggering in the N = 66 isotones is plotted in Fig. 2; to gauge any systematic behavior, plots are drawn to the same scale.

In general, the staggering amplitude is high for spherical nuclei and decreases towards the mid-shell region. The N = 66 isotones do not constitute an exception to the systematic trends established in Ref. [2]. The γ-unstable nucleus 112 Pd (10 bosons) shows the highest amplitude in S(J) among the N = 66 isotones. The γ-band staggering is less pronounced in 110 Ru (11 bosons) while in 108 Mo (12 bosons) it is highly suppressed. From these systematics trends a γ band with a very low staggering amplitude is expected in 106 Zr, namely S(J) ≈ 10 keV, which would make 106 Zr the nucleus with the lowest staggering amplitude observed in the 40 ≤ Z ≤ 50 region. Also electric quadrupole transition probabilities have been considered for the neutron-rich nuclei belonging to the N = 66 and A = 106 chains. The calculated values are listed in Tables V and VI along with the few known experimental values [16]. The boson effective charge e b = 0.11 eb is obtained after a least-squares fit to the B(E2) values measured in 106 Pd and this value is used for all other nuclei. The results show an increase of the B(E2; 2 + 1 → 0 + 1 ) value from the near-vibrational 106 Pd towards the near-rotational 108 Mo. This B(E2) value is predicted highest for 106 Zr and the calculations for the two different chains give a consistent value of B(E2; 2 + 1 → 0 + 1 ) ≈ 0.50-0.55 e 2 b 2 . Furthermore, we note from Table VI the behavior predicted for the B(E2; 2 + 2 → 0 + 1 ) value which is small in 106 Pd and 106 Zr but becomes large for the isotopes in between, a characteristic feature of transitional nuclei.

IV. CONCLUSIONS

The aim of this paper was to illustrate the usefulness of the interacting boson model in its simplest version, the IBM-1, for predicting properties of exotic nuclei. The usual difficulty that arises with such attempts is that the variations of the model's parameters with the valence neutron and proton numbers are not known, pre- This method was applied to the neutron-rich members of the N = 66 isotonic and A = 106 isobaric chains of which level energies and electric quadrupole transition probabilities were fitted with IBM-1. The two chains in-tersect at 106 Zr which allowed the prediction of this nucleus' excitation energies and electric quadrupole transition properties.

As a final remark we emphasize that the IBM-1 is a valence-nucleon model and that the extrapolations as described here crucially depend on the definition of neighboring closed-shell configurations. Thus, the results on 106 Zr indirectly involve the assumption of the persistence of magic numbers in this region of neutron-rich nuclei.

V. ACKNOWLEDGMENTS S.L. acknowledges the Bulgarian National Science Fund for a Research Fellowship Grant. This work

FIG. 1 :

 1 FIG.1: Even-even nuclei in the region close to 106 Zr. All nuclei in boxes with full lines have known excited levels while those in boxes with broken lines are only known in their ground state. The nuclei used in the fit are in boxes with thick lines.

106 FIG. 2 :

 1062 FIG. 2: Experimental (full symbols connected with full lines) and theoretical (open symbols connected with dashed lines) γ-band staggering S(J) in the N = 66 isotones. In 106 Zr both near-identical predictions are shown based on the extrapolation of the N = 66 isotones (open squares) and of the A = 106 isobars (open diamonds).

TABLE I :

 I Parameters and rms deviation (in keV) for the N = 66 isotones and the A = 106 isobars.

		ǫ d	κ	c3	c4	χ a rms
	N = 66	1171 -24.3 -70.8 -134.9 -0.30	58
	A = 106 1053 -25.6 -70.6 -118.3 -0.30	143

a Dimensionless.

TABLE II :

 II Experimental and theoretical excitation energies (in keV) of levels in neutron-rich N = 66 isotones.

	nucleus	J π gsb	E ex gsb	E th gsb	J π γ	E ex γ	E th γ

TABLE IV :

 IV Experimental and theoretical excitation energies (in keV) of levels in 106 Pd.

	nucleus	J π	E1	E2	E3	E4	E5
	106 Pd	0 + ex	0	1134	1706	2001	2278
	(N = 7)	0 + th 2 + ex	0 512	1017 1128	1384 1562	2080 1909	2297 2242
		2 + th 3 + ex 3 + th 4 + ex	519 1558 1483 1229	982 1932	1575 2077	1859 2283	2004
		4 + th 5 + ex	1153 2366	1617	2010	2213
		5 + th 6 + ex	2141 2077
		6 + th	1896

TABLE V :

 V Theoretical B(E2) values (in units e 2 b 2 ) for transitions in neutron-rich N = 66 isotones.

	transition	112 Pd	110 Ru	108 Mo	106 Zr
	2 + 1 → 0 + 1 4 + 1 → 2 + 1 0 + 2 → 2 + 1 2 + 2 → 0 + 1	0.27 0.43 0.22 0.015	0.35 0.54 0.17 0.024	0.45 0.65 0.057 0.025	0.53 0.75 0.003 0.020

TABLE VI :

 VI Experimental and theoretical B(E2) values (in units e 2 b 2 ) for transitions in neutron-rich A = 106 isobars.

	transition	106 Pd		106 Ru	106 Mo 106 Zr
		Expt	Theo	Theo	Theo	Theo
	2 + 1 → 0 + 1 4 + 1 → 2 + 1 0 + 2 → 2 + 1 2 + 2 → 0 + 1	0.134 (4) 0.211 (17) 0.22 0.14 0.105 (20) 0.13 0.0035 (3) 0.0019 0.0097 0.023 0.23 0.37 0.36 0.55 0.16 0.092	0.54 0.77 0.0044 0.021
	was partially supported by the Agence National de
	Recherche, France, under contract nr ANR-07-BLAN-
	0256-03.

cluding a reliable extrapolation to unknown regions of the nuclear chart. Our method proposes to circumvent this problem by studying the structural evolution in three different chains of nuclei, namely isotopic, isotonic, and isobaric ones, which cross at the nucleus of interest at the outskirts of the region of stable nuclei. Predictions are obtained by extrapolating the different chains to the exotic nucleus in question, and, in addition, a comparison of these extrapolations gives an idea of the errors involved.