text:
# A novel multigrid method for electronic structure calculations

## I Introduction

One of the goals of computational materials science is to calculate from first principles the various physical and chemical properties. This requires the solution of the electronic and ionic structures of the materials system in question. The density-functional theory (DFT) makes a huge step towards this goal by casting the untractable problem of many interacting electrons to that of noninteracting particles under the influence of an effective potential . The adiabatic approximation allows one to separate the ionic degrees of freedom from those of the electrons. However, in order to apply DFT in practice one has to resort to approximations for electron exchange and correlation such as the local-density approximation (LDA) or the generalized-gradient approximation (GGA). Moreover, in the case of systems consisting of hundreds or more atoms it is still a challenge to solve numerically efficiently for the ensuing Kohn-Sham equations.

The numerical solution of the Kohn-Sham equations is the concern of our present work. It deals with real-space (RS) methods, in which the values of the different functions are presented using three-dimensional point grids, and the partial differential equations are discretized using finite differences . The RS methods, as suggested by the name chosen, are contrasted with the popular plane-wave (PW) schemes . There are several aspects favouring the RS methods over the PW methods. Both of the methods are used in the context of pseudopotentials describing the electron-ion interactions, but only the RS can easily be used in all-electron calculations or with hard pseudopotentials of, i.e. first-row or transition metal atoms, because the RS grid can be refined in a natural way in the ion core regions (composite grids , adaptive coordinates ). Systems, such as surfaces, containing different length scales are more economically described in the RS than in the PW scheme because one needs not waste many grid points in the vacuum regions to describe the slowly varying tails of wave functions. In the RS methods periodic boundary conditions are not necessary. This leads to ease and accuracy in describing charged atomic clusters in contrast to PW methods requiring an artificial neutralizing background charge. Besides the above “physical arguments” there are also methodological and computational aspects favouring the RS methods. The RS methods allow a systematic convergence control by increasing the grid (or basis function) density. (The PW methods do also so by adjusting the cutoff energy of the plane wave expansions.) The so-called “order-N” methods , the computational cost of which scales linearly with the number of electrons, require localized real-space wave functions leading naturally to the employment of RS methods . The discretizations in the real-space grid can be made local and therefore parallelization can effectively use data decomposition in which different real-space regions are handled with different processing units and the communications between processing units will be mainly short-ranged .

More specifically, our choice for the numerical method is a multigrid scheme. Several approaches employing the multigrid idea within electronic structure calculations have appeared during recent years . The main idea of multigrid methods is that they avoid the critical slowing-down (CSD) phenomenon occuring when a partial differential equation discretized on a real space grid is solved with a simple relaxation method such as the Gauss-Seidel method. The discretization operators typically use information from a rather localized region of the grid at a time. Therefore the high frequency error of the length scale of the grid spacing is reduced very rapidly in the relaxation. However, once the high frequency error has effectively been removed, the very slow convergence of the low frequency components dominates the overall error reduction rate , i.e. CSD occurs. In multigrid methods one stops the relaxation on a given (fine) grid before CSD sets in and transfers the equation to a coarser grid (the so-called restriction operation) where the low-frequency components can be solved more efficiently. On the coarsest grid the problem is solved exactly or as accurately as possible, after which one interpolates (the so-called prolongation operation) the correction to finer grids, performing simultaneously relaxations in order to remove the high-frequency errors introduced in the interpolation.

The solution of the Poisson equation by multigrid methods is straightworward . This is because the error (or the correction needed) also obeys a Poisson equation and thus will be a smooth function to be presented and solved on the repeatedly coarser grids optimal to handle the lower frequencies. The solution of an eigenvalue problem, such as the Schrödinger equation, is a much more complicated task than that of the Poisson equation. The problem is no more linear because both the eigenfunction and the eigenvalue have to be solved simultaneously. Then the error no longer obeys the same equation as the solution. Also one has to solve for several eigenpairs (eigenvalues and corresponding eigenvectors). Moreover, the existence of both negative and positive eigenvalues makes the problem indefinite. This implies severe difficulties for many simple iterative methods which converge only in the case of a positive definite iteration matrix. In particular, it can easily be shown that using Gauss-Seidel relaxation for the Schrödinger equation the high frequency components typically converge as in the case of the Poisson equation, but the low frequency components may diverge, although the divergence may be slow . More complicated methods such as Kaczmarz relaxation are guaranteed to converge, but may have clearly inferior high frequency reduction rates, which are essential for the overall speed of multigrid methods. Other possible convergent methods include GMRES  which is considerably more complex than Gauss-Seidel relaxation.

A standard recipe for dealing with eigenproblems with multigrids is the full-approximation-storage (FAS) method originally described by Brandt . In FAS one solves for the entire problem on the coarse grids also and ends up in solving for a properly modified problem so that its solution can be used in correcting the fine grid solution. The FAS method may not be very straighworward to implement for the Schrödinger equation. It is also difficult to present some actual potential on the coarse levels accurately enough. However, some succesful applications of FAS have appeared in the context of electronic structure calculations by Beck et al.  and advanced strategies for FAS have been proposed .

Briggs et al.  employ a multigrid method in electronic structure calculations by linearizing the Schrödinger problem and presenting the potential contribution on the coarse levels by an error term (residual) only. Then on the coarse levels they solve effectively for the Poisson problem. Ancilotto et al.  modified the method by Briggs et al. by shifting to a full multigrid (FMG) scheme and by solving on the coarse grids a problem including a local potential term. The idea of FMG is to start the smoothing iterations from a coarse grid. Then the interpolation to a finer grid provides a good initial guess of the solution. The FMG scheme can accelerate the convergence remarkably with respect to the (above-described) V-cycle scheme in which one starts from the finest level. Fattebert  used a multigrid method with a block Galerkin inverse iteration (BGII) and GMRES in the relaxations. In the method, the current approximation is kept orthogonal against all the nearby states during the multigrid cycle. The inverse iteration converges for a given guess for the energy eigenvalue towards the nearest eigenvalue. In order to solve all the desired lowest eigenstates a good guess for the eigenvalue spectrum is needed in the beginning of iterations, but thereafter large computational savings may be expected because explicit orthogonalizations are not needed (at least between well-separated states).

A severe problem in the existing multigrid schemes for the Schrödinger equation is often that the coarse grids cannot well approximate the solutions of the coarse grid equations themselves. As a consequence the correction from coarse grids, no matter how accurately the equation is solved, may be ineffective in correcting the fine grid solution and as a result the overall process converges slowly. Therefore one has to restrict to the use of rather fine grids only and the convergence speed of the scheme is drastically lowered. In those multigrid methods, which use the potential also on the coarse grids the size of the coarsest grid has been typically of the order of 31x31x31 . Using the FAS method coarser grids are possible at least for systems with a small number of eigenstates solved . If a large number of eigenstates have to be solved problems may arise because the coarse grids may not be able to represent eigenstates with many nodes or the ordering of the states may change between the successive grids. To bypass these problems in FAS, rather complicated strategies are needed .

In order to avoid the coarse grid representation problems we utilize the so-called Rayleigh Quotient Multigrid (RQMG) method introduced by Mandel and McCormick . In this method the coarse grid relaxation passes are performed so that the Rayleigh quotient calculated on the fine grid will be minimized. In this way there is no requirement for the solution to be well represented on a coarse grid and the coarse grid representation problem is avoided. Mandel and McCormick introduced the method for the solution of the eigenpair corresponding to the lowest eigenvalue. We have generalized it to the simultaneous solution of a desired number of lowest eigenenergy states by developing a scheme which keeps the eigenstates separated by the use of a penalty functional, Gramm-Schmidt orthogonalization, and subspace rotations. Our generalization of the RQMG method is an attractive alternative for large-scale electronic structure calculations.

The Kohn-Sham equations have to be solved self-consistently, i.e. the wave functions solved from the single-particle equation determine via the density (solution of the Poisson equation and the calculation of the exchange-correlation potential) the effective potential for which they should again be solved. To approach this self-consistency requires an optimized strategy so that numerical accuracy of the wave functions and the potential increase in balance, enabling the most efficient convergence . In order to avoid the divergence of the self-consistency iterations, the mixing of the input and output solutions is needed. For this feedback procedure sophisticated schemes  and control strategies  have been presented.

The outline of the present paper is as follows. In Section II we represent shortly the most important ideas of the density functional theory. Section III is devoted for numerical methods, the most important of which is the Schrödinger equation solver developed. Also the strategy for the self-consistency iterations is discussed. In Section IV we demonstrate by the help of a couple of examples the performance of our scheme in calculating the electronic structures of small molecules and solid-state systems described by pseudopotentials. Section V summarises the work and gives outlines for the future developments.

## II The Kohn-Sham scheme

In the Kohn-Sham method for electronic structure calculations one solves for a set of equations self-consistently . In the following, we present the equations in the spin-compensated form. In practice, we have made the straightforward generalization using the spin-density functional theory. The set of equations reads as (atomic units with $`\mathrm{}=m_e=e=1`$ are used):

$$\left(\frac{1}{2}^2+V_{\mathrm{eff}}(𝐫)\right)\mathrm{\Psi }_i=ϵ_i\mathrm{\Psi }_i,$$
(1)
$$n(𝐫)=\underset{i}{\overset{N}{}}|\mathrm{\Psi }_i(𝐫)|^2,$$
(2)
$$V_{\mathrm{eff}}(𝐫)=V_{\mathrm{ion}}(𝐫)+V_\mathrm{H}(𝐫)+V_{\mathrm{XC}}(𝐫),$$
(3)
$$V_\mathrm{H}(𝐫)=\frac{n(𝐫^{})}{|𝐫𝐫^{}|}𝑑𝐫^{},$$
(4)
$$V_{\mathrm{XC}}(𝐫)=\frac{\delta E_{\mathrm{XC}}[n(𝐫)]}{\delta n(𝐫)}.$$
(5)
The first equation (1) is a Schrödinger equation for non-interacting particles in an effective potential $`V_{\mathrm{eff}}(𝐫)`$. For finite systems the wave functions are required to vanish at the boundaries of the computation volume. In the case of infinite periodic systems the complex wave functions have to obey the Bloch theorem at the cell boundaries. The electron density $`n(𝐫)`$ is obtained from a sum over the $`N`$ occupied states. The effective potential consists of an external potential $`V_{\mathrm{ion}}(𝐫)`$ due to ions (or nuclei in all-electron calculations), the Hartree potential $`V_\mathrm{H}(𝐫)`$ calculated from the electron density distribution, and the exchange-correlation potential $`V_{\mathrm{XC}}(𝐫)`$. In the examples of the present work we use the norm-conserving non-local pseudopotentials for the electron-ion interactions and the local-density approximation (LDA) for the exchange-correlation energy

$$E_{\mathrm{XC}}[n(𝐫)]=ϵ_{\mathrm{XC}}(n(𝐫))n(𝐫)𝑑𝐫,$$
(6)
and for the exchange-correlation potential

$$V_{\mathrm{XC}}(𝐫)=ϵ_{\mathrm{XC}}(n(𝐫))+n(𝐫)\frac{dϵ_{\mathrm{XC}}}{dn}_{|n=n(𝐫)}.$$
(7)
The Hartree potential is solved from the Poisson equation

$$^2V_\mathrm{H}(𝐫)=4\pi n(𝐫).$$
(8)
In practice, the electron density $`n(𝐫)`$ is substituted by the total charge density $`\rho (𝐫)`$, which includes the positive ionic (nuclear) charge neutralizing the system. In the case of finite systems, Dirichlet boundary conditions are used with the Coulomb potential values calculated using a multipole expansion. For periodic systems we fix the average Coulomb potential to zero and allow the periodic boundary conditions to result in the corresponding converged potential.

The self-consistent solution of the above Kohn-Sham equations leads to the ground state electronic structure minimizing the total energy

$`E_{\mathrm{tot}}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \mathrm{\Psi }_i^{}(𝐫)\left(\frac{1}{2}^2\right)\mathrm{\Psi }_i(𝐫)𝑑𝐫}+{\displaystyle \frac{1}{2}}{\displaystyle V_\mathrm{H}(𝐫)n(𝐫)𝑑𝐫}`$ (9)
$`+`$ $`{\displaystyle V_{\mathrm{ion}}(𝐫)n(𝐫)𝑑𝐫}+E_{\mathrm{XC}}+E_{\mathrm{ion}\mathrm{ion}},`$ (10)
where $`E_{\mathrm{ion}\mathrm{ion}}`$ is the repulsive interaction between the ions (nuclei) of the system. Instead of the self-consistency iterations the solution of the Kohn-Sham problem can be found by minimizing directly the total energy with respect to the wave function parameters, e.g. plane-wave coefficients . However, Kresse and Furthmüller  have found this scheme less efficient than the self-consistency iterations.

## III Numerical methods

### A Schrödinger equation solver

In our real space method we start from an initial guess for the effective potential and initial wave functions generated by random numbers in grid points. The wave functions and the Hartree potential are updated alternatingly towards self-consistency. The solution of the Poisson equation is a standard task for the multigrid scheme. If a reasonable guess for the Coulomb potential is not available, the FMG method will produce the solution starting from random numbers and requiring the work which scales linearly as a function of the size of the system ($`O(N)`$). During the Kohn-Sham iterations one can start from the present approximation of Coulomb potential and update it with respect to the new charge density by performing only a few V-cycles.

The solution of the wave functions is a much more complicated task than that of the Poisson equation because one has to solve an eigenvalue problem which in the state-of-the-art electronic structure calculations means the determination of several hundreds of eigenpairs. For this purpose we have developed a scheme based on RQMG method introduced by Mandel and McCormick for the solution of the eigenpair corresponding to the lowest eigenvalue. We begin by reviewing the basic principles of RQMG. This is most easily done in the framework of the so-called coordinate relaxation method. Thereafter we go through the modifications made in order to simultaneously solve for several eigenpairs.

Coordinate relaxation is a method of solving the discretized eigenproblem

$$Hu=\lambda Bu$$
(11)
by minimizing the Rayleigh quotient

$$\frac{u|H|u}{u|B|u}.$$
(12)
Above, $`H`$ and $`B`$ are matrix operators chosen so that the Schrödinger equation discretized on a real-space point grid with spacing $`h`$ is satisfied to a chosen order $`O(h^n)`$. In Eq. (12) $`u`$ is a vector containing the wave function values at the grid points. In the relaxation method, the current estimate $`u`$ is replaced by itself plus a multiple of some search vector $`d`$

$$u^{}=u+\alpha d,$$
(13)
and $`\alpha `$ is chosen to minimize the Rayleigh quotient. This leads to a simple quadratic equation for $`\alpha `$. (Find the minimum of the expression (15) below with respect to $`\alpha `$. In the case of a complex wave function one has to solve for the real and imaginary parts of $`\alpha `$ from a coupled pair of quadratic equations.) Moreover, the search vector $`d`$ is simply chosen to be unity in one grid point and to vanish in all other points. A complete coordinate relaxation pass is then obtained by performing the minimization at each point in turn and these passes can be repeated until the lowest state is found with desired accuracy.

Naturally, also the coordinate relaxation suffers from CSD because of the use of local information only in updating $`u`$ in a certain point. In order to avoid it one applies the multigrid idea. In the multigrid scheme by Mandel and McCormick the crucial point is that coarse grid coordinate relaxation passes are performed so that the Rayleigh quotient calculated on the fine grid will be minimized. In this way there is no requirement for the solution to be well represented on a coarse grid. In practice, a coarse grid search substitutes the fine grid solution by

$$u_f^{}=u_f+\alpha I_c^fd_c,$$
(14)
where the subscripts $`f`$ and $`c`$ stand for the fine and coarse grids, respectively, and $`I_c^f`$ a prolongation operator interpolating the coarse grid vector to the fine grid. The Rayleigh quotient to be minimized is then

$`{\displaystyle \frac{u_f+\alpha I_c^fd_c|H_f|u_f+\alpha I_c^fd_c}{u_f+\alpha I_c^fd_c|B_f|u_f+\alpha I_c^fd_c}}=`$ (15)
$`{\displaystyle \frac{u_f|H_fu_f+2\alpha I_f^cH_fu_f|d_c+\alpha ^2d_c|H_cd_c}{u_f|B_fu_f+2\alpha I_f^cB_fu_f|d_c+\alpha ^2d_c|B_cd_c}}.`$ (16)
The second form is obtained by relating the coarse grid operators, $`H_c`$ and $`B_c`$, with the fine grid ones, $`H_f`$ and $`B_f`$, by the Galerkin condition

$$\begin{array}{cc}H_c\hfill & =I_f^cH_fI_c^f\hfill \\ B_c\hfill & =I_f^cB_fI_c^f,\hfill \end{array}$$
(17)
and the restriction operator $`I_f^c`$ has to be the transpose of the prolongation operator

$$I_f^c=\left(I_c^f\right)^T.$$
(18)
The key point to note is that when $`H_fu_f`$ and $`B_fu_f`$ are provided from the fine grid to the coarse grid, the remaining integrals can be calculated on the coarse grid itself. Thus one really applies coordinate relaxation on the coarse grids to minimize the fine level Rayleigh quotient. This is a major departure from the earlier methods, which to some extent rely on the ability to represent the solution of some coarse grid equation on the coarse grid itself. Here, on the other hand, one can calculate the exact change in the Rayleigh quotient due to any coarse grid change, no matter how coarse the grid itself is. There is no equation whose solution would have to be representable.

Thus, in the Rayleigh quotient minimization multigrid (RQMG) algorithm the coordinate relaxation passes on each level keep track of the integrals in Eq. (15). Actually, on the finest level we use Gauss-Seidel relaxation, which very effectively smoothens the errors of the wavelength corresponding to the grid spacing. When calculating several eigenpairs Gauss-Seidel relaxation may also work as a residual minimization method. The idea is that the coarse grid-iterations with Gramm-Schmidt orthogonalization can provide the separation of the eigenstates so well that the subsequent finest level relaxations converge to the correct (nearest) eigenstates without orhogonalization. This requires that the effect of the coarse-level smoothings on the low-frequency components of the solutions overcomes the possible divergence tendency of these components caused by the Gauss-Seidel relaxation on the finest level.

Moreover, we discretize the original equation separately on each grid (discretization coarse grid approximation (DCA)) instead of using the Galerkin conditions of Eq. (17) This may in principle decrease the convergence rate and force a limit to the coarsest possible grid in order to avoid instability or divergence. However, we have observed this DCA implementation of RQMG to be quite stable and efficient. To avoid possible coarse-level instabilities occuring especially during the first few iteration cycles we may recalculate the Rayleigh quotient whenever coarse grid corrections are interpolated to a finer grid. Later when approaching the convergence the recalculation can be omitted.

For the matrix operators $`H`$ and $`B`$ we have used either high-order ($`O(h^4)`$ or higher) Mehrstellen or central difference stencils (CDS) . The use of high-order stencils reduces remarkably the density of grid points needed. The benefit of the Mehrstellen scheme is that more local information is used. The scheme leads to controlled accuracy and convergence properties and to more isotropic smoothing of the error in comparison with the use of CDS’s. The local nature enables also a more efficient parallel coding. As the prolongation operator $`I_c^f`$ we usually use trilinear interpolation and as the restriction operator $`I_f^c`$ its transpose, the so-called full-weighting operator, in which the coarse-grid values are chosen to be the averaged values of the surrounding fine grid points. The integrations are performed by the trapezoidal rule.

Next we consider the generalization of the RQMG method to the simultaneous solution of several ($`N`$) mutually orthogonal eigenpairs. The separation of the different states is divided into two or three subtasks. First, in order to make the coarse grid relaxations converge towards the desired state we apply a penalty functional scheme. Given the $`k`$ lowest eigenfunctions, the next lowest, $`(k+1)`$’th state is searched for by minimizing the functional

$$\frac{u_{k+1}|H|u_{k+1}}{u_{k+1}|B|u_{k+1}}+\underset{i=1}{\overset{k}{}}q_i\frac{u_i|u_{k+1}^2}{u_i|u_iu_{k+1}|u_{k+1}}.$$
(19)
The overlap integral in the penalty term is squared to make the penalty positive definite. The denominator is required to make the functional independent of the norms of $`u_i,i=1\mathrm{}k+1`$. The minimization of this functional is equivalent to imposing the orthonormality constraints against the lower $`k`$ states, when $`q_i\mathrm{}`$. By increasing the shifts $`q_i`$ any desired accuracy can be obtained, but in order to obtain a computationally efficient algorithm a reasonable finite value should be used, for example

$$q_i=(\lambda _{k+1}\lambda _i)+\mathrm{Q},$$
(20)
where $`Q`$ is a sufficiently large positive constant. In our test calculations $`Q`$ is of the order of $`Q=0.5\mathrm{}2`$ Ha.

We minimize the expression (19) simultaneously for all $`N`$ states. This simplifies the algorithm and enables a future parallelization over the eigenstates. Thus the current approximations are used for $`u_i,i=1\mathrm{}k`$. Moreover, changes in the $`u_i`$ during a given relaxation sweep are not used to update the penalty term in Eq. (19). This is sufficient, when the states are always ordered in the same way, in the order of increasing eigenvalue. In order to reduce computations, the $`B`$-innerproduct is actually used in calculating the penalty term integrals because the values of $`Bu`$ are readily available from the finer level. The substitution (14) is introduced in the functional (19) and the minimization with respect to $`\alpha `$ leads again to a quadratic equation. This time the coefficients contain terms due to the penalty part.

On the finest level, we do not apply the minimization of the penalty functional. The ideal situation would be if a residual minimization method, such as the Gauss-Seidel method, would keep the states calculated on the coarse levels separated. We found out in practical calculations that this is not true at least when the states are far from convergence. Therefore we have developed for the finest level a scheme which by employing Gramm-Schmidt orthogonalization and subspace rotation keeps the eigenstates orthogonal. The subspace rotation is a method to find the most optimally separated eigenvectors from the approximative ones. The major steps of the rotation are:

(i) Calculation of the Hamiltonian matrix elements between the current states:

$$\overline{H}_{i,j}=u_i|B^1H|u_j.$$
(21)
(ii) Calculation of the overlap matrix:

$$\overline{S}_{i,j}=u_i|u_j.$$
(22)
The use of matrix elements of Eqs. (21) and (22) leads to eigenvectors orthogonal in the desired Euclidian sense ($`I`$-orthogonal) and not in the sense of the $`B`$-innerproduct.

(iii) Diagonalization to find the optimal eigenvectors ($`u_k^{}=_j\overline{A}_{k,j}u_j`$) and corresponding eigenvalues ($`\lambda _k`$):

$$\underset{j}{}\overline{H}_{i,j}\overline{A}_{j,k}=\lambda _k\underset{j}{}\overline{S}_{i,j}\overline{A}_{j,k}.$$
(23)
In practice, we apply the approximation

$$u_i|B^1H|u_ju_i|u_j\frac{u_i|Hu_j}{u_i|Bu_j}.$$
(24)
The Gramm-Schmidt orthogonalization and the subspace rotations are organized so that the space of the eigenvectors is first divided to small clusters corresponding to close eigenvalues. The Gramm-Schmidt orthogonalization is then performed for each cluster at a time so that its eigenvectors become orthogonal against the eigenvectors of the clusters of lower eigenvalues. Then a subspace rotation is performed within the states belonging to the present cluster. The division to clusters reduces remarkably the cost of the subspace rotation. This is because the cost is proportional to $`O(N^3`$), where $`N`$ is the number of states rotated. Moreover, the subspace rotation requires the calculation of matrix elements which are more complicated than those for the simple Gramm-Schmidt orthogonalization.

According to our test calculations this subspace rotation scheme leads quite effectively to $`I`$-orthogonal eigenstates. This is seen as a convergence of the eigenvalue problem within the numerical accuracy, i.e. the residuals of different eigenstates vanish. In order to prove this one has, in practice, to introduce potential shifts which reduce the number of significant digits in the eigenvalues so that the numerical accuracy of the eigenvalue does not prevent to reach the numerical accuracy of the wavefunction, i.e. the vanishing residual. The error in the eigenvalue scales as the square of the residual. When applying the subspace rotation it is important to complete the highest eigenvalue cluster; otherwise the rotation may become inefficient.

The orthogonalization needed scales as $`O(N^3)`$. For small systems of several tens of eigenpairs this is not yet a problem. The algorithm is effective and the number of fine grid orthogonalizations remains quite plausible, for example, in comparison with the conjugate gradient search of eigenpairs employing only the finest grid . But for larger systems with hundreds of states it will be the bottleneck. One solution could be to rely on the finest level only on a residual minimization method when the initial stages of the iteration process have been performed and the solution is clearly on a stable track towards convergence.

### B Strategy for self-consistency iterations

The Kohn-Sham problem has to be solved self-consistently. This means that an optimal strategy is needed so that computing time is not wasted in the beginning of the self-consistency iterations to obtain unnecessarily accurate wave functions, because these will change during the later iterations due to the changes in the potential. Updating the potential, including the solution of the Poisson equation, is a much less time-consuming task than the update of all the wavefunctions. Therefore the potential update can be performed frequently .

The examples of this paper are small-molecule and bulk-solid systems described by pseudopotentials. The strategy used is schematically presented in Fig. 1. Similar strategies can certainly be applied in other kind of Kohn-Sham calculations, for example in those employing all-electron or jellium-type models. In the examples of this work the initial electron density is the superposition of the pseudoatom densities centered around given nuclear positions. From the superposition we calculate the initial effective potential, where the wavefunctions are solved accurately enough using the full-multigrid method. The FMG process is started from random numbers for the wavefunctions on the coarsest level. The accuracy of the wavefunctions is controlled by calculating the norms of the residuals of the eigenstates and it is finally improved by adding more V-cycles starting from the finest level. A certain accuracy is needed in order to initiate self-consistency iterations which converge without large density oscillations. Then the new electron density and the ensuing effective potential are calculated. The new potential is not directly fed into the next iteration but it is mixed in this place, as well as later between the self-consistency iterations, with the input potential of the iteration. We monitor the accuracy of the wave functions by calculating their residuals and require that the accuracy has improved from the previous iteration. Usually one V-cycle is sufficient for this, because the changes in the potential are small.

An important point is also to find a proper balance with respect to the pre- and postsmoothening sweeps on the different grid levels. Typically, on the finest level we made two pre- and postsmoothening sweeps whereas on the coarser grids their number is four. Actually, this means that on the finest level four successive sweeps are done if the potential is not updated. A potential update is always preceded and followed by two immediate smoothening sweeps.

## IV Tests

We test the performance of our scheme by calculating the self-consistent electronic structures of a CO<sub>2</sub> molecule as well as that of perfect bulk Si lattice with a supercell of 64 Si atoms. The former system is an example of the employment of Dirichlet boundary conditions and the use of “hard” pseudopotentials whereas the latter system represents the use of periodic boundary conditions and a supercell size typical in electronic structure calculations for point defects in solids.

The ions are described by pseudopotentials of the Kleinman-Bylander form ,

$`V_{ion}(𝐫)`$ $`=`$ $`{\displaystyle \underset{a}{}}V_{ion,loc}(|𝐫_a|)`$ (25)
$`+`$ $`{\displaystyle \underset{a,n,lm}{}}{\displaystyle \frac{1}{<\mathrm{\Delta }V_{lm}^a>}}|\mathrm{\Delta }V_{ion,l}(r_a)u_{lm}(𝐫_a)>`$ (26)
$`\times `$ $`<\mathrm{\Delta }V_{ion,l}(r_a^{})u_{lm}(𝐫_{}^{}{}_{a}{}^{})|,`$ (27)
where $`<\mathrm{\Delta }V_{lm}^a>`$ is a normalization factor,

$$<\mathrm{\Delta }V_{lm}^a>=u_{lm}(𝐫_a)\mathrm{\Delta }V_{ion,l}(r_a)u_{lm}(𝐫_a)d^3r,$$
(28)
and $`𝐫_a=𝐫𝐑_a`$, $`u_{lm}`$ are the atomic pseudopotential wave functions of angular and azimuthal momentum quantum numbers $`(l,m)`$, from which the $`l`$-dependent ionic pseudopotentials $`V_{ion,l}(r)`$ are generated using the Troullier-Martins scheme . The ion core is assumed to be spherically symmetric. $`\mathrm{\Delta }V_{ion,l}(r)=V_{ion,l}(r)V_{ion,loc}(r)`$ is the difference between the $`l`$-component of the ionic pseudopotential and the local ionic potential. We have chosen the s-component of the pseudopotential as the local component.

Because the functions $`u_{lm}(𝐫)`$ are short-ranged, operating on the wave-function by the nonlocal parts of the pseudopotential is in practice a multiplication by a sparse matrix. The numerical work required to compute this scales as the square of the number of atoms in the system, whereas in the conventional reciprocal-space formulation the work scales as the cube of the system size. The advantage of implementing the nonlocal pseudopotentials in real space has been noted also in the context of plane-wave methods .

In the previous multigrid implementations of the pseudopotential method , the nonlocal parts have only been employed on the finest grid. It is, however, straightforward to implement them also on the coarse levels, and we have found that this may increase the convergence rate and stability of the method.

The CO<sub>2</sub> molecule is placed diagonically in the center of a cubic computation volume of the size of (12.6 a$`{}_{0}{}^{})^3`$. Experimental bond lengths are used. Dirichlet boundary conditions are used so that the potential values outside the cube are obtained from a multipole expansion of the charge density. The point mesh used is 63<sup>3</sup>, giving the grid spacing $`h`$ = 0.20 a<sub>0</sub>. The Mehrstellen discretization by Briggs et al.  is used.

In this calculation we used a mixing scheme, where the new effective potential $`V_{in}^{i+1}`$ is obtained from the input and output potentials according to

$$V_{in}^{i+1}=(1\kappa )V_{in}^i+\kappa V_{out}^i.$$
(29)
The convergence of the self-consistency iterations employing the strategy described above (Fig. 1) is shown in Fig. 2. The deviation of the total energy from the converged value is given as a function of self-consistency iteration steps performed. The zeroth iteration is a full-multigrid solution for the wavefunctions in the initial potential. Two V-cycles are performed on the finest level at this point. The effective potential obtained from the output electron density was mixed with the initial potential using the feedback $`\kappa =0.4`$. Next, at iteration one, the wave-functions are relaxed in this new potential using one V-cycle. From this point on, the four curves in the figure give the convergence with different values of the feedback parameter $`\kappa `$. One V-cycle per one self-consistency iteration step is done. A wide range of values for $`\kappa `$ gives satisfactory convergence indicating a robust behaviour for the scheme. The accuracy of 1 meV, which is sufficient in practical calculations, is reached after three or four V-cycles. The implementation of the non-local parts of the pseudopotential on the coarse levels is found to speed up the convergence especially in this region. From Fig. 2 we obtain an average convergence rate of approximately one decade per self-consistency iteration. This is of the same order as those reported by Wang and Beck  in their FAS scheme or by Kresse and Furthmüller  in their pseudopotential scheme employing self-consistency iterations. The convergence rate of one decade per self-consistency iteration is better than that obtained by Ancilotto et al.  in the FMG scheme and much better than the rate reached in the linearized multigrid scheme by Briggs et al. .

We have solved for the electronic structure of perfect Si lattice described by a supercell of 64 Si ions. The lattice constant of 20.38 a<sub>0</sub> used is the equilibrium value obtained in a plane-wave calculation, with which we have compared our real-space results. The first Brillouin zone is sampled in this test using the $`\mathrm{\Gamma }`$-point only. The point mesh used for the wave-functions is 64<sup>3</sup>, giving the grid spacing $`h`$ = 0.32 a<sub>0</sub>. For the densities and potentials we use a finer grid of 128<sup>3</sup> points. The other numerical parameters and the iteration strategy are the same as in the CO<sub>2</sub> test. The resulting valence electron density on the (110)-plane is given in Fig. 3. The area of the figure corresponds to the extent of the supercell. One notes that exactly the same features are reproduced at the equivalent points in different regions of the supercell. This means that a fully converged result has been found. We have compared the results of our real-space code to those obtained using the plane-wave method. The energy cutoff, 18 Ry, of the plane-wave expansion was chosen so that it results in a real-space point mesh of 64<sup>3</sup>, i.e. it is the same as in our real-space calculation. The widths of the valence band and band gaps obtained by the two methods agree with an accuracy of 3 meV. In the case of degenerate eigenstates the real-space code results in degenerate eigenenergies with an accuracy better than 1 meV. The convergence towards to the self-consistent solution occurs similarly as for the CO<sub>2</sub> molecule in Fig. 2. Thus, the convergence process seems to be independent of the size of the system.

## V Summary and outlook

In this work we have generalized the RQMG method introduced by Mandel and McCormick  for the simultaneous solution of a desired number of lowest eigenenergy states. This approach can be viewed as belonging to a third group of multigrid methods, in addition to FAS and the techniques where the eigenproblem is linearized (e.g. inverse iteration. In principle, one can use arbitrarily coarse grids in RQMG, whereas in the other multigrid methods one has to be able to represent all the states on the coarsest grid.

We have demonstrated the feasibility of the method by electronic structure calculations for the CO<sub>2</sub> molecule and bulk Si described by pseudopotentials. Our strategy for the self-consistent solution consists of a full-multigrid solution for the wave-functions in the initial potential, and subsequent self-consistency iterations. Less than five V-cycles are generally sufficient for practically sufficient accuracy. The cpu-times required for the FMG and SCF steps are roughly equal.

We have applied the method also in two-dimensional problems for quantum dots employing the current-spin-density functional theory, in three-dimensional cylindrically symmetric systems, and also for calculation of positron states in solids.

We believe that our method will eventually compete with the standard plane-wave methods for electronic structure calculations. However, some straightforward programming is still required. For calculations, where the optimization of the ionic structure is necessary, the Hellmann-Feynman forces will be implemented. In order to remove the spurious dependence of the total energy on the position of the atoms with respect to the grid points, Fourier-filtering of the pseudopotentials is required. Complex wave functions for any k-point are easily implemented, and are already in use in two-dimensional geometries.

Parallelization over k-points can be done easily. One only needs to communicate the electron density and effective potential at the end of each V-cycle. During the RQMG V-cycle, the states are all relaxed simultaneously and independently of each other. Therefore parallelization over states is natural and easy to implement. However, for larger systems the Gramm-Schmidt orthogonalization becomes very inefficient in a state-parellel code. The most efficient and yet straightforward choice is real-space data decomposition, where each processor is mapped to a specific region of space.

###### Acknowledgements.

M.H. acknowledges stimulating discussions with Prof. S. F. McCormick and Prof. T. L. Beck. T.T. acknowledges financial support by the Vilho, Yrjö and Kalle Väisälä foundation. This research has been supported by the Academy of Finland through its Centre of Excellence Programme (2000 - 2005).


format: text

text:
# Tensor algebras, displacement structure, and the Schur algorithm

## 1. Introduction

It was recently showed in  that the tensor algebras have displacement structure. The goal of this paper is to further develop this connection. The displacement structure theory, as it was initiated in the paper , consists in recursive factorization of matrices with additional structure encoded by a so-called displacement equation. Many applications of this theory were found in applied fields, like in the study of wave propagation in layered media, or filtering and modeling of nonstationary processes (see  for a survey). On the mathematical side, there are applications to bounded interpolation in upper triangular algebras (including the classical setting of bounded analytic functions on the unit disk) and to triangular factorizations of operators.

On the other hand, the tensor algebra is a classical concept and in  it is showed that there is a connection expressed by the fact that any element $`T`$ of the unit ball of the algebra (with respect to an appropriate operator norm) satisfies the equation

$$A\underset{k=1}{\overset{N}{}}F_kAF_k^{}=GJG^{},$$

where $`A=IT^{}T`$ and the so-called generators $`F_k`$, $`G`$ and $`J`$ are described in the next section. This equation is of displacement type so we can develop a theory for the tensor algebra that parallels the one dimensional case. It is the main goal of this paper to describe several such results.

The paper is organized as follows. In Section 2 we review the basic elements relevant to the tensor algebra and to the so-called noncommutative analytic Toeplitz algebra which is a natural extension of the tensor algebra introduced in , and studied in other papers such as , , , . We describe an isometric isomorphic representation of this later algebra by a certain algebra $`𝒰_𝒯()`$ of upper triangular bounded operators and this represention allows us to rely on the results and ideas in the displacement structure theory. In this section we also show that a certain class of positive definite kernels on the free semigroup with $`N`$ generators, introduced in , also has displacement structure of the type mentioned above and in Corollary 2.5 we explain the connection between this class of positive definite kernels and the algebra $`𝒰_𝒯()`$ which is well-known for $`N=1`$. In Section 3 we obtain results corresponding to so-called scattering experiments. Theorem 3.1 is the main result in this direction and as a by-product we can obtain Schur type parametrizations for $`𝒰_𝒯()`$ and the realization of the elements of the algebra as transfer maps of a certain class of nonstationary linear systems. Section 4 develops a Schur type algorithm for $`𝒰_𝒯()`$. We plan to explore approximation properties (Szegö type theory) related to this algorithm in a future paper. This paper is part of .

## 2. Preliminaries

We briefly review the construction of the tensor algebras and we explain their displacement structure.

### 2.1. Tensor Algebras

We introduce some notation relevant to tensor algebras and we describe their representation as algebras of upper triangular operators (see, for instance, ). The associative tensor algebra $`𝒯()`$ generated by the complex vector space $``$ is defined by the algebraic direct sum

$$𝒯()=_{k0}^k,$$

where

$$^k=\underset{kfactors}{\underset{}{\mathrm{}}}$$

is the $`k`$-fold algebraic tensor product of $``$ with itself. The elements of $`𝒯()`$ are terminating sequences $`x=(x_0,x_1,\mathrm{},x_k,0,\mathrm{})`$, where $`x_p^p`$ is called the $`p`$th homogeneous component of $`x`$ (note that $`^0=`$, $`^1=`$). The addition and multiplication of elements in $`𝒯()`$ are defined componentwise:

$$(x+y)_n=x_n+y_n$$

and

$$(xy)_n=\underset{k+l=n}{}x_ky_l,(x_0y_n=y_nx_0=x_0y_n).$$

In order to deal with displacement structure matters it is convenient to use the fact that $`𝒯()`$ can be represented by upper triangular matrices of a special type. For a certain simplicity, from now on we take $`=^N`$, $`N1`$. Then we denote by $`()`$ the full Fock space associated to $``$, that is, the Hilbert space

$$()=_{k0}^k$$

obtained by taking the Hilbert direct sum of the spaces $`^k`$, $`k0`$, on which we consider the tensor Hilbert space structure induced by the Euclidean norm on $`=^N`$. If $`\{e_1,e_2,\mathrm{},e_N\}`$ is the standard basis of $`^N`$, then $`\{e_{i_1}\mathrm{}e_{i_k}i_1,\mathrm{},i_k\{1,2,\mathrm{},N\}\}`$ is an orthonormal basis for $`^k`$. It is more convenient to write $`e_\sigma `$ instead of $`e_{i_1}\mathrm{}e_{i_k}`$, where $`\sigma =i_1\mathrm{}i_k`$ is viewed as a word in the unital free semigroup $`𝔽_N^+`$ with $`N`$ generators $`1`$, …,$`N`$. The empty word is the identity element of $`𝔽_N^+`$. The length of the word $`\sigma `$ is denoted by $`|\sigma |`$. We also use to write the elements of $`𝔽_N^+`$ in lexicographic order. The space $`^{k+1}`$ can be identified with the direct sum of $`N`$ copies of $`^k`$. More precisely, we notice that $`\{(\delta _{jl}e_\sigma )_{l=1}^Nj\{1,\mathrm{},N\},|\sigma |=k\}`$ is an orthonormal basis for $`\underset{Nterms}{\underset{}{^k\mathrm{}^k}}=(^k)^N`$ ($`\delta _{jl}`$ is the Kronecker symbol). We deduce that the mapping

$$\varphi _k(((\delta _{jl}e_\sigma )_{l=1}^N))=e_{j\sigma },j=1,\mathrm{},N,$$

extends to a unitary operator from $`(^k)^N`$ onto $`^{k+1}`$ and the formula

$$\psi _k=\varphi _{k1}\varphi _{k2}^N\mathrm{}\varphi _1^{N^{k2}},k1,$$

gives a unitary operator from $`_k=^{N^{k1}}`$ onto $`^k`$ (the notation $`T^l`$ means the direct sum of $`l`$ copies of the operator $`T`$). We finally obtain a unitary operator $`\psi :_{k=0}^{\mathrm{}}_k()`$ by the formula $`\psi =_{k=0}^{\mathrm{}}\psi _k`$, where $`\psi _0`$ is the identity on $`=_0`$. Each $`_k`$, $`k2`$, is also identified with the direct sum of $`N`$ copies of $`_{k1}`$,

$$_k=\underset{Nterms}{\underset{}{^{N^{k2}}\mathrm{}^{N^{k2}}}},$$

but the explicit mention of the identification can be omitted in this case.

An operator $`T(_{k=0}^{\mathrm{}}_k)`$ has a matrix representation $`T=[T_{ij}]_{i,j=0}^{\mathrm{}}`$ with $`T_{ij}(_j,_i)`$. We say that $`T`$ is upper triangular provided that $`T_{ij}=0`$ for $`i>j`$. We denote by $`𝒰_𝒯^0()`$ the set of upper triangular operators $`T=[T_{ij}]_{i,j=0}^{\mathrm{}}`$ in $`(_{k=0}^{\mathrm{}}_k)`$ with the property that for $`ij`$ and $`i,j1`$,

(2.1)
$$T_{ij}=T_{i1,j1}^N$$

and $`T_{0j}=0`$ for $`j>k`$ and some $`k`$. We notice that the entries $`T_{0j}`$, $`j0`$, determine the matrix $`T`$. Also, a simple calculation shows that $`𝒰_𝒯^0()`$ is an associative algebra.

We define an algebra isomorphism

(2.2)
$$\mathrm{\Phi }:𝒯()𝒰_𝒯^0()$$

as follows: for $`x𝒯()`$, $`x=(x_0,x_1,\mathrm{})`$, there are complex numbers $`c_\sigma `$, $`\sigma 𝔽_N^+`$, such that

(2.3)
$$x_j=\underset{|\sigma |=j}{}c_\sigma e_\sigma .$$

For $`j0`$, $`T_{0j}`$ is the row matrix $`[c_\sigma ]_{|\sigma |=j}`$. Then $`T_{0j}=0`$ for sufficiently large $`j`$ and we can define $`T(_{k=0}^{\mathrm{}}_k)`$ by using (2.1). Set $`\mathrm{\Phi }(x)=T`$. Clearly $`\mathrm{\Phi }`$ is a bijection. In order to see that $`\mathrm{\Phi }`$ is also an algebra morphism it is convenient to introduce the following operators. First, the operator $`C_j^{}`$, $`j=1,\mathrm{},N`$, is an element of $`𝒰_𝒯^0()`$ defined by the formula: $`C_j^{}=\mathrm{\Phi }(0,e_j,0,\mathrm{})=[T_{lk}^j]_{l,k=0}^{\mathrm{}}`$, where $`T_{0k}^j=0`$ for $`k1`$ and $`T_{01}^j=e_j`$. Then $`C_j^+=(C_j^{})^{}`$, $`j=1,\mathrm{},N`$. If we use the notation $`C_\sigma ^{}`$ to denote the operator $`C_{i_1}^{}\mathrm{}C_{i_k}^{}`$, where $`\sigma =i_1\mathrm{}i_k`$, then we deduce that each $`T𝒰_𝒯^0()`$ has a representation

(2.4)
$$T=\underset{\sigma 𝔽_N^+}{}c_\sigma C_\sigma ^{}.$$

Now the required properties of $`\mathrm{\Phi }`$ follow from the easily verifiable fact that

$$\mathrm{\Phi }((0,e_j,0,\mathrm{})(0,e_k,0,\mathrm{}))=C_j^{}C_k^{}.$$

It appears to be quite natural to consider the algebra $`𝒰_𝒯()`$ of all those operators in $`(_{k=0}^{\mathrm{}}_k)`$ that satisfy the condition (2.1). It turns out that $`𝒰_𝒯()`$ is isometrically isomorphic to the so-called non-commutative analytic Toeplitz algebra, recently studied in papers like , , . We also denote by $`𝒮()`$ the Schur class of all contractions in $`𝒰_𝒯()`$. We can extend this setting as folows. Let $`_1`$ and $`_2`$ be Hilbert spaces, then $`𝒰_𝒯(,_1,_2)`$ denotes the set of the operators $`T=[T_{ij}]_{i,j=0}^{\mathrm{}}`$ in $`(_{k=0}^{\mathrm{}}(_k_1),_{k=0}^{\mathrm{}}(_k_2))`$ obeying (2.1). Instead of $`𝒰_𝒯(,,)`$ we write $`𝒰_𝒯(,)`$. Also $`𝒮(,_1,_2)`$ denotes the set of the contractions in $`𝒰_𝒯(,_1,_2)`$. We notice that $`𝒰_𝒯()=𝒰_𝒯(,,)`$.

### 2.2. Displacement Structure

The systematic study of the displacement structure was initiated in , and since then manifold applications have been found in the fields of mathematics and electrical engineering (for a recent surveys, see ). Two main themes of this theory are about scattering experiments and recursive factorizations. A matrix $`R`$ has displacement structure with respect to the generators $`F`$ and $`G`$ provided that

(2.5)
$$AFAF^{}=GJ_{pq}G^{},$$

where $`J_{pq}=I_pI_q`$ is a signature matrix that specifies the displacement inertia of $`A`$. In many applications, $`r=p+q`$ is much smaller than the size of $`A`$. From now on we will write $`J`$ instead of $`J_{11}`$. There is also a version of this theory allowing $`A`$ to depend on a (usually integer) paramater, see  or . A major result concerning matrices with displacement structure is that, under suitable assumptions on the generators, the succesive Schur complements of $`A`$ inherit a similar structure. This allows to write the Gaussian elimination for $`A`$ only in terms of generators (see , ).

We now describe a connection between displacement structure and $`𝒰_𝒯()`$. Let $`C_k=[T_{ij}^k]_{i,j=0}^{\mathrm{}}`$, $`k=1,\mathrm{},N`$, be isometries defined by the formulae: $`T_{ij}=0`$ for $`ij+1`$ and $`T_{i+1,i}`$ is a block-column matrix consisting of $`N`$ blocks of dimension $`dim_i`$, all of them zero except for the $`k`$th block which is $`I__i`$. These operators are closely related to the operators $`C_k^+`$, in the sense that there is a unitary operator $`u`$ such that $`C_k=uC_k^+u^{}`$, $`k=1,\mathrm{},N`$. The following result was noticed in .

###### Theorem 2.1.

Let $`T=[T_{ij}]_{i,j=0}^{\mathrm{}}𝒮()`$ and $`A=IT^{}T`$. Then

(2.6)
$$A\underset{k=1}{\overset{N}{}}C_kAC_k^{}=GJG^{},$$

where the generator $`G`$ is given by the formula:

$$G=\left[\begin{array}{cc}1& T_{00}^{}\\ 0& T_{01}^{}\\ \mathrm{}& \mathrm{}\end{array}\right].$$

The equation (2.6) can be rewritten in the form of a time variant displacement equation by setting $`F=[C_1,\mathrm{}C_N]`$. We now discuss some additional examples.

#### 2.2.1. Positive definite kernels on $`𝔽_N^+`$

A mapping $`K:𝔽_N^+\times 𝔽_N^+`$ is called a positive definite kernel on $`𝔽_N^+`$ provided that

$$\underset{i,j=1}{\overset{k}{}}K(\sigma _i,\sigma _j)\lambda _j\overline{\lambda }_i0$$

for each positive integer $`k`$ and each choice of words $`\sigma _1,\mathrm{},\sigma _k`$ in $`𝔽_N^+`$ and complex numbers $`\lambda _1,\mathrm{},\lambda _k`$. Without loss of generality, we assume $`K(\mathrm{},\mathrm{})=1`$. The semigroup $`𝔽_N^+`$ acts on itself by juxtaposition and we assume that the positive definite kernel $`K`$ is invariant under this action, that is

$$K(\tau \sigma ,\tau \sigma ^{})=K(\sigma ,\sigma ^{}),\tau ,\sigma ,\sigma ^{}𝔽_N^+.$$

By the invariant Kolmogorov decomposition theorem, see e.g., , Ch. II, there exists an isometric representation $`u`$ of $`𝔽_N^+`$ on a Hilbert space $`𝒦`$ and a mapping $`v:𝔽_N^+𝒦`$ such that

$$K(\sigma ,\tau )=v(\tau ),v(\sigma )_𝒦,$$

$$u(\tau )v(\sigma )=v(\tau \sigma )$$

for all $`\sigma ,\tau 𝔽_N^+`$, and the set $`\{v(\sigma )\sigma 𝔽_N^+\}`$ is total in $`𝒦`$. The mapping $`v`$ is unique up to unitary equivalence and $`u`$ is uniquely determined by $`v`$. Also, the previous representation of $`K`$ implies, as noticed in , that $`u(1)`$, $`\mathrm{}`$, $`u(N)`$ are isometries with orthogonal ranges if and only if $`K(\sigma ,\tau )=0`$ for all pairs $`(\sigma ,\tau )`$ with the property that there is no $`\alpha 𝔽_N^+`$ such that $`\sigma =\alpha \tau `$ or $`\tau =\alpha \sigma `$. In this case, $`K`$ is called in  a multi-Toeplitz kernel and we will use the same terminology. We can define

(2.7)
$$S_{ij}=[K(\sigma ,\tau )]_{|\sigma |=i,|\tau |=j}(_j,_i).$$

By induction on the length of words and using properties of the lexicographic order on words we check that if $`K`$ is a positive definite multi-Toeplitz kernel, then for all $`i,j1,`$

(2.8)
$$S_{ij}=S_{i1,j1}^N.$$

This relation is an analogue of (2.1), therefore we expect that positive definite multi-Toeplitz kernels have displacement structure.

###### Theorem 2.2.

Let $`S=[S_{ij}]_{i,j=0}^{\mathrm{}}`$ be the array of operators associated to a positive definite multi-Toeplitz kernel on $`𝔽_N^+`$. Then

(2.9)
$$S\underset{k=1}{\overset{N}{}}C_kSC_k^{}=GJG^{},$$

where the generator $`G`$ is given by the formula:

$$G=\left[\begin{array}{cc}1& 0\\ S_{01}^{}& S_{01}^{}\\ \mathrm{}& \mathrm{}\end{array}\right].$$

The equality above can be interpreted in the sense of an equality of unbounded operators with dense domain made of the terminating sequences in $`_{l=0}^{\mathrm{}}_l`$.

###### Proof.

We write $`S=L+Q`$, where $`L=[L_{ij}]_{i,j=0}^{\mathrm{}}`$ and $`Q=[Q_{ij}]_{i,j=0}^{\mathrm{}}`$ are the lower and, respectively, the upper parts of $`S`$, such that $`L_{ij}=S_{ij}`$ for $`i>j`$ and $`L_{ij}=0`$ for $`ij`$, while $`Q_{ij}=S_{ij}`$ for $`ij`$ and $`Q_{ij}=0`$ for $`i>j`$. Since, formally (that is, computing the entries of the arrays that are involved according to the usual rule of multiplication of operators), $`LC_k=C_kL`$ and $`QC_k^{}=C_k^{}Q`$ for all $`k=1,\mathrm{},N`$, we deduce that

$$\begin{array}{ccc}S_{k=1}^NC_kSC_k^{}\hfill & =\hfill & L+Q_{k=1}^NC_k(L+Q)C_k^{}\hfill \\ & & \\ & =\hfill & L+Q_{k=1}^NC_kLC_k^{}_{k=1}^NC_kQC_k^{}\hfill \\ & & \\ & =\hfill & L+Q_{k=1}^NLC_kC_k^{}_{k=1}^NC_kC_k^{}Q.\hfill \end{array}$$

But $`I=_{k=1}^NC_kC_k^{}=P__0`$, the orthogonal projection on $`_0`$, so we can deduce that

$$\begin{array}{ccc}S_{k=1}^NC_kSC_k^{}\hfill & =\hfill & L+QL(IP__0)(IP__0)Q\hfill \\ & & \\ & =\hfill & LP__0+P__0Q=GJG^{}.\hfill \end{array}$$

We now consider truncations of a positive definite multi-Toeplitz kernel on $`𝔽_N^+`$. If $`K`$ is such a kernel and $`S_{ij}`$ are associated to $`K`$ by the formula (2.7), then we define $`K_n=[S_{ij}]_{i,j=0}^n`$. These are positive operators on $`_{l=0}^n_l`$ and their entries satisfy (2.8). We also introduce the operators $`C_{k,n}=P_{_{l=0}^n_l}C_k_{l=0}^n_l`$. We obtain another example of a displacement equation.

###### Corollary 2.3.

The positive operator $`K_n`$ satisfies the displacement equation

(2.10)
$$K_n\underset{k=1}{\overset{N}{}}C_{k,n}K_nC_{k,n}^{}=GJG^{},$$

where the generator $`G`$ is given by the formula:

$$G=\left[\begin{array}{cc}1& 0\\ S_{01}^{}& S_{01}^{}\\ \mathrm{}& \mathrm{}\\ S_{0n}^{}& S_{0n}^{}\end{array}\right].$$

###### Proof.

This is an immediate consequence of Theorem 2.2. ∎

It turns out that the operators satisfying (2.10) have a certain explicit description. For $`x_{l=0}^n_l`$ we denote by $`U(x)`$ the upper triangular operator on $`_{l=0}^n_l`$ satisfying (2.8) and whose first row is $`x`$.

###### Theorem 2.4.

A positive operator $`A`$ on $`_{l=0}^n_l`$ satisfies (2.10) if and only if $`A=U(x)^{}U(x)U(y)^{}U(y)`$ for two vectors $`x,y_{l=0}^n_l`$ with the property that there is another vector $`z_{l=0}^n_l`$ such that $`U(y)=U(z)U(x)`$ and $`U(z)`$ is a contraction.

###### Proof.

If $`A=U(x)^{}U(x)U(y)^{}U(y)`$ then we use the fact that for any $`x_{l=0}^n_l`$ and $`k=1,\mathrm{},N`$, $`C_{n,k}U(x)^{}=U(x)^{}C_{n,k}`$, and deduce (2.10) as in the proof of Theorem 2.2. Conversely, if $`A`$ satisfies (2.10), then we use the fact that $`C_{\sigma ,n}=0`$ for $`|\sigma |>n`$ in order to deduce that

$$A=GJG^{}+\underset{|\sigma |=1}{\overset{n}{}}C_{\sigma ,n}GJG^{}C_{\sigma ,n}^{}.$$

Hence, if $`G=[x^{},y^{}]`$ for some $`x,y_{l=0}^n_l`$, then $`A=U(x)^{}U(x)U(y)^{}U(y)`$. Since $`A`$ is positive, it follows that there exists a contraction $`T`$ such that $`U(y)=TU(x)`$. We notice that we can assume, without loss of generality, that $`x_00`$ and then we deduce that $`C_{n,k}T^{}=T^{}C_{n,k}`$ for $`k=1,\mathrm{},N`$, which implies that $`T=U(z)`$ for some $`z_{l=0}^n_l`$. ∎

We mention a simple but useful consequence of this result. Let $`K`$ be a positive definite multi-Toeplitz kernel on $`𝔽_N^+`$ with $`K(\mathrm{},\mathrm{})=1`$ and let $`S_{ij}`$, $`i,j0`$, be the operators associated to $`K`$ by (2.7). Let

$$x_n=\left[\begin{array}{c}1\\ S_{01}^{}\\ \mathrm{}\\ S_{0n}^{}\end{array}\right],y_n=\left[\begin{array}{c}0\\ S_{01}^{}\\ \mathrm{}\\ S_{0n}^{}\end{array}\right],$$

then $`K_n=U(x_n)^{}U(x_n)U(y_n)^{}U(y_n)`$ and by the previous result there exists $`z_n_{l=0}^n_l`$ such that $`U(y_n)=U(z_n)U(x_n).`$ In fact, $`z_n`$ is uniquely determined and also $`P_{_{l=0}^n_l}U(z_{n+1})_{l=0}^n_l=U(z_n)`$ for all $`n1`$. It follows that there exists a unique $`T(K)𝒮()`$ such that $`P_{_{l=0}^n_l}T_{l=0}^n_l=U(z_n)`$ and $`T(K)_{00}=0`$. Denote by $`𝒮_0()`$ the set of those elements $`T𝒮()`$ with $`T_{00}=0`$.

###### Corollary 2.5.

The mapping $`KT(K)`$ gives a one-to-one correspondence between the set of positive definite multi-Toeplitz kernels on $`𝔽_N^+`$ with $`K(\mathrm{},\mathrm{})=1`$ and $`𝒮_0()`$.

#### 2.2.2. Pick kernels.

This type of kernels were recently introduced in , , and they were shown to have a certain universality property with respect to a Nevanlinna-Pick problem, .

Let $`^N`$ be the open unit ball in $`^N`$. For distinct points $`\lambda _1`$, $`\mathrm{}`$, $`\lambda _L`$ in $`^N`$ and complex numbers $`b_1`$, $`\mathrm{}`$, $`b_L`$, we define the Pick kernel by the formula

$$R=\left[\frac{1b_j\overline{b_l}}{1\lambda _j,\lambda _l}\right]_{j,l=1}^L,$$

where $`\lambda _j,\lambda _l`$ denotes the Euclidean inner product in $`^N`$.

We can define

$$F_j=\text{diag}(\lambda _l^{(j)})_{l=1}^L,j=1,\mathrm{},N,$$

the diagonal matrix with the diagonal made of the $`j`$th components of the points $`\lambda _1`$, $`\mathrm{}`$, $`\lambda _L`$. Also, we set

$$G=\left[\begin{array}{cc}1& b_1\\ 1& b_2\\ \mathrm{}& \mathrm{}\\ 1& b_L\end{array}\right].$$

We can check by direct computation that

(2.11)
$$R\underset{k=1}{\overset{N}{}}F_kRF_k^{}=GJG^{}.$$

It was shown in  how to use this equation in order to solve some multidimensional Nevanlinna-Pick type problems as those in , .

## 3. Scattering Experiments

One of the main results in displacement structure theory refers to the possibility of associating a scattering operator to the data given by the generators of the equation (2.5). In this section we prove a similar result for displacement equations of the type of (2.11) and we explain the connection with a special class of time variant linear systems.

Let $`_1`$, $`_2`$, and $`𝒢`$ be Hilbert spaces. We consider a displacement equation of the form

(3.1)
$$A\underset{k=1}{\overset{N}{}}F_kAF_k^{}=GJG^{},$$

where $`F_k(𝒢)`$, $`k=1,\mathrm{},N`$, are given contractions on the Hilbert space $`𝒢`$. Also $`G=\left[\begin{array}{cc}U& V\end{array}\right](_1_2,𝒢)`$ and $`J=I__1I__2`$. The wave operators associated to (3.1) are introduced by the formulae: $`U_{\mathrm{}}^{}=\left[U_k\right]_{k=0}^{\mathrm{}}`$, $`V_{\mathrm{}}^{}=\left[V_k\right]_{k=0}^{\mathrm{}}`$, where $`U_k=\left[F_\sigma U\right]_{|\sigma |=k}:_k_1𝒢`$, $`V_k=\left[F_\sigma V\right]_{|\sigma |=k}:_k_2𝒢`$. We will assume that both $`U_{\mathrm{}}`$ and $`V_{\mathrm{}}`$ are bounded and also that $`lim_k\mathrm{}_{|\sigma |=k}F_\sigma ^{}g=0`$ for all $`g𝒢`$. Under these assumptions we deduce that (3.1) has a unique solution given by

(3.2)
$$A=U_{\mathrm{}}^{}U_{\mathrm{}}V_{\mathrm{}}^{}V_{\mathrm{}}.$$

The next result was already noticed in , but here we present a different approach based on system theoretic ideas.

###### Theorem 3.1.

The solution (3.2) of the displacement equation (3.1) is positive if and only if there exists $`T𝒮(,_1,_2)`$ such that $`V_{\mathrm{}}=TU_{\mathrm{}}`$.

###### Proof.

Assume $`A=U_{\mathrm{}}^{}U_{\mathrm{}}V_{\mathrm{}}^{}V_{\mathrm{}}0`$ and let $`A=LL^{}`$ be a factorization of $`A`$ with $`L(,𝒢)`$ for some Hilbert space $``$. From (3.1) we deduce that

$$LL^{}+VV^{}=\underset{k=1}{\overset{N}{}}F_kLL^{}F_k^{}+UU^{}.$$

In matrix form,

(3.3)
$$\left[\begin{array}{cc}L& V\end{array}\right]\left[\begin{array}{c}L^{}\\ V^{}\end{array}\right]=\left[\begin{array}{cccc}F_1L& \mathrm{}& F_NL& U\end{array}\right]\left[\begin{array}{c}L^{}F_1^{}\\ \mathrm{}\\ L^{}F_N^{}\\ U^{}\end{array}\right].$$

Define $`A^{}=\left[\begin{array}{cc}L& V\end{array}\right]`$ and $`B^{}=\left[\begin{array}{cccc}F_1L& \mathrm{}& F_NL& U\end{array}\right]`$, then we deduce from (3.3) that there exists an unitary operator $`\theta _0(\overline{(B)},\overline{(A)})`$ such that $`A=\theta _0B.`$ It follows that there exist Hilbert spaces $`_1`$, $`_2`$, and an unitary extension $`\theta (^N_1_1,_2_2)`$ of $`\theta _0`$, hence this extension satisfies the relation

(3.4)
$$\left[\begin{array}{c}A\\ 0__2\end{array}\right]=\theta \left[\begin{array}{c}B\\ 0__1\end{array}\right].$$

Let $`\theta _{ij}`$, $`i\{1,2,3\}`$, $`j\{1,2,\mathrm{},N+2\}`$, be the matrix coefficients of $`\theta `$. It is convenient to rename some of these coefficients. Thus, we set

$$\begin{array}{cccc}X_k=\theta _{1k},& k=1,\mathrm{}N,& & Z=\theta _{1,N+1},\\ & & & \\ Y_k=\theta _{2k},& k=1,\mathrm{}N,& & W=\theta _{2,N+1}.\end{array}$$

From (3.4) we deduce that

$$L^{}=\underset{k=1}{\overset{N}{}}X_kL^{}F_k^{}+ZU^{}$$

and

$$V^{}=\underset{k=1}{\overset{N}{}}Y_kL^{}F_k^{}+WV^{}.$$

By induction we deduce that

(3.5)
$$V^{}=WU^{}+\underset{k=1}{\overset{N}{}}\underset{|\sigma |=0}{\overset{n}{}}Y_kX_\sigma ZU^{}F_{k\sigma }^{}+\underset{|\tau |=n+1}{}Q_\tau L^{}F_\tau ^{},$$

where $`Q_\tau `$ are monomials of length $`|\tau |`$ in the variables $`X_1`$, $`\mathrm{}`$,$`X_N`$, $`Y_1`$, $`\mathrm{}`$, $`Y_N`$. Since $`\theta `$ is unitary it follows that all $`Q_\tau `$ are contractions.

We define $`T_{00}=W`$ and for $`j>0`$,

$$T_{0j}=[Y_kX_\sigma Z]_{|\sigma |=j1;k=1,\mathrm{},N}.$$

Then we define $`T_{ij}`$, $`i>0`$, $`ji`$, by the formula (2.1) and $`T_{ij}=0`$ for $`i>j`$. We show that $`T=[T_{ij}]_{i,j=0}^{\mathrm{}}`$ belongs to $`𝒮(,_1,_2)`$. To that end, we introduce the notation: for $`n0`$,

$$\begin{array}{ccc}X_1(n)=\left[\begin{array}{ccc}X_1& \mathrm{}& X_N\end{array}\right]^{N^n},\hfill & & Y_1(n)=\left[\begin{array}{ccc}Y_1& \mathrm{}& Y_N\end{array}\right]^{N^n},\hfill \\ & & \\ X_2(n)=Z^{N^n},\hfill & & Y_2(n)=W^{N^n}.\hfill \end{array}$$

We check by induction that $`T_{00}=Y_2(0),`$ $`T_{01}=Y_1(0)X_2(1)`$, and for $`j2`$,

$$T_{0j}=Y_1(0)X_1(1)\mathrm{}X_1(j1)X_2(j).$$

Consequently, $`T^{}`$ is the transfer map of the linear time variant system

$$\{\begin{array}{ccc}\hfill x(n+1)& =& X_1(n)^{}x(n)+Y_1(n)^{}u(n)\hfill \\ & & \\ \hfill y(n)& =& X_2(n)^{}x(n)+Y_2(n)^{}u(n),n0.\hfill \end{array}$$

Since each matrix $`\left[\begin{array}{cc}X_1(n)& X_2(n)\\ Y_1(n)& Y_2(n)\end{array}\right]`$, $`n0`$, is a contraction, it follows that $`T^{}`$ is a contraction, hence $`T`$ is a contraction and $`T𝒮(,_1,_2)`$. Also, since $`lim_k\mathrm{}_{|\sigma |=k}F_\sigma ^{}g=0`$ for all $`g𝒢`$, we deduce from (3.5) that $`V_{\mathrm{}}=TU_{\mathrm{}}`$. ∎

This result is a generalization of Theorem 3.1. It can be used to solve interpolation and moment problems in several non commuting variables as suggested in . The proof given here emphasizes the fact that this type of problems appears as an interesting particular case of similar questions in time variant setting (as described, for instance, in ). Next we develop this idea by giving a parametrization of the elements in $`𝒮(,_1,_2)`$.

### 3.1. Schur Parameters

We introduce a Schur type parametrization for $`𝒮(,_1,_2)`$ and then give a number of applications. A first instance of such a result might be considered Euler’s description of $`SO(3)`$. Schur’s classical result in  gives a parametrization of the contractive holomorphic functions on the unit disk and Szegö’s theory of orthogonal polynomials provides an alternative to this result for probability measures on the unit circle.

Here we provide a refinement of a parametrization of $`𝒮(,_1,_2)`$ given in . We can deal with a slightly more general situation than $`𝒮()`$. For a contraction $`C`$ we denote by $`𝒟_C`$ its defect space, that is the closure of the range of the defect operator $`D_C=(IC^{}C)^{1/2}`$. Let $`\mathrm{\Pi }(_1,_2)`$ be the family of sequences $`\{\gamma _\sigma \}_{\sigma 𝔽_N^+}`$ with the properties: $`\gamma _{\mathrm{}}(_1,_2)`$ and for $`|\sigma |>0`$,

(3.6)
$$\gamma _{k\sigma }(𝒟_{\gamma _\sigma },𝒟_{\gamma _{k\sigma 1}^{}}),$$

where $`k=1,\mathrm{}N`$ and $`k\sigma 1`$ denotes the predecessor of $`k\sigma `$ with respect to the lexicographic order.

###### Theorem 3.2.

There exists a one-to-one correspondence between $`𝒮(,_1,_2)`$ and $`\mathrm{\Pi }(_1,_2)`$.

###### Proof.

A simple proof can be obtained as an application of the Schur parametrization of upper triangular contractions in . Let $`T𝒮(,_1,_2)`$ and for each $`j0`$, $`T_{0j}`$ is the row matrix $`[c_\sigma ]_{|\sigma |=j}`$. We use Theorem 2.2.1 in  in order to associate to $`T`$ a family $`\mathrm{\Gamma }`$ of Schur parameters. The Schur parameters of the contraction

$$\left[\begin{array}{cc}T_{00}& T_{01}\\ 0& T_{11}\end{array}\right]=\left[\begin{array}{ccccc}c_{\mathrm{}}& c_1& \mathrm{}& & c_N\\ 0& c_{\mathrm{}}& 0& \mathrm{}& 0\\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ 0& \mathrm{}& & 0& c_{\mathrm{}}\end{array}\right]$$

are given by $`\mathrm{\Gamma }_{00}=\gamma _{\mathrm{}}=c_{\mathrm{}}`$, $`\mathrm{\Gamma }_{11}=\mathrm{\Gamma }_{00}^N`$ and $`\mathrm{\Gamma }_{01}(𝒟_\gamma _{\mathrm{}}^N,𝒟_\gamma _{\mathrm{}})`$.

For the next step we associate Schur parameters to the contraction

$$\left[\begin{array}{ccc}T_{00}& T_{01}& T_{02}\\ 0& T_{11}& T_{12}\\ 0& 0& T_{22}\end{array}\right].$$

We notice that

$$\left[\begin{array}{cc}T_{11}& T_{12}\\ 0& T_{22}\end{array}\right]=\left[\begin{array}{cc}T_{00}^N& T_{01}^N\\ 0& T_{11}^N\end{array}\right].$$

By the formula (2.2.3) in , the entries of the matrix $`\left[\begin{array}{cc}T_{11}& T_{12}\\ 0& T_{22}\end{array}\right]`$ are polynomials in the Schur parameters, their adjoints and their defect operators. Consequently, the contraction associated to the Schur parameters $`\mathrm{\Gamma }_{00}^N`$, $`\mathrm{\Gamma }_{01}^N`$ and $`\mathrm{\Gamma }_{11}^N`$ is precisely $`\left[\begin{array}{cc}T_{00}^N& T_{01}^N\\ 0& T_{11}^N\end{array}\right].`$ From the uniqueness of the Schur parameters, we deduce that the Schur parameters of $`\left[\begin{array}{cc}T_{11}& T_{12}\\ 0& T_{22}\end{array}\right]`$ are precisely $`\mathrm{\Gamma }_{00}^N`$, $`\mathrm{\Gamma }_{01}^N`$ and $`\mathrm{\Gamma }_{11}^N`$. By induction, we deduce that we can associate to $`T`$ the family of Schur parameters $`\mathrm{\Gamma }=\{\mathrm{\Gamma }_{ij}\}_{ij}`$ such that each $`\mathrm{\Gamma }_{ij}`$ satisfies (2.1). We notice that $`\mathrm{\Gamma }_{0j}:𝒟_{\mathrm{\Gamma }_{0,j1}}^N𝒟_{\mathrm{\Gamma }_{0,j1}^{}}`$ for $`j1`$. Also, each $`\mathrm{\Gamma }_{0j}`$ is a row contraction, $`\mathrm{\Gamma }_{0j}=[\gamma _\sigma ^0]_{|\sigma |=j}`$. By Proposition 1.4.2 in ,

$$\mathrm{\Gamma }_{01}=\left[\begin{array}{cccc}\gamma _1& D_{\gamma _1^{}}\gamma _2& \mathrm{}& D_{\gamma _1^{}}\mathrm{}D_{\gamma _{N1}^{}}\gamma _N\end{array}\right],$$

where $`\gamma _k(𝒟_\gamma _{\mathrm{}},𝒟_{\gamma _{k1}^{}})`$, $`k=1,\mathrm{},N`$. Also, there exist unitary operators $`u_1:𝒟_{\mathrm{\Gamma }_{01}}_{k=1}^N𝒟_{\gamma _k}`$ and $`v_1:𝒟_{\mathrm{\Gamma }_{01}^{}}𝒟_{\gamma _N^{}}`$. Next, we have that $`\mathrm{\Gamma }_{02}:𝒟_{\mathrm{\Gamma }_{01}}^N𝒟_{\mathrm{\Gamma }_{01}^{}}`$ and we can define the row contraction $`\stackrel{~}{\mathrm{\Gamma }}_{02}=[\gamma _\sigma ^0]_{|\sigma |=2}`$ from $`\left(_{k=1}^N𝒟_{\gamma _k}\right)^N`$ into $`𝒟_{\gamma _N^{}}`$ by the formula $`\stackrel{~}{\mathrm{\Gamma }}_{02}=v_1\mathrm{\Gamma }_{02}(u_1^N)^{}`$. As above,

$$\stackrel{~}{\mathrm{\Gamma }}_{02}=\left[\begin{array}{cccc}\gamma _{11}& D_{\gamma _{11}^{}}\gamma _{12}& \mathrm{}& D_{\gamma _{11}^{}}\mathrm{}D_{\gamma _{N,N1}^{}}\gamma _{NN}\end{array}\right],$$

where $`\gamma _{k\tau }(𝒟_{\gamma _\tau },𝒟_{\gamma _{k\tau 1}^{}})`$, $`k=1,\mathrm{},N`$, $`|\tau |=1`$. We proceed by induction in order to construct the whole family of Schur parameters $`\{\gamma _\sigma \}_{\sigma 𝔽_N^+}`$ with $`\gamma _{\mathrm{}}(_1,_2)`$ and obeying (3.6). ∎

The parameters $`\{\mathrm{\Gamma }_{0j}\}_{j0}`$ are the Schur parameters associated to $`T`$ in . The above proof explain their construction within the framework of . The new parameters $`\gamma _\sigma `$ are more like ”true” Schur parameters in the sense that they are contractions subject to the only restrictions given by the compatibility relations (3.6).

A similar parametrization can be obtained for the set of positive multi-Toeplitz kernels. Denote by $`\mathrm{\Pi }_0(_1,_2)`$ the set of families $`\{\gamma _\sigma \}_{\sigma 𝔽_N^+}`$ in $`\mathrm{\Pi }(_1,_2)`$ with $`\gamma _{\mathrm{}}=0`$.

###### Theorem 3.3.

There exists a one-to-one correspondence between the set of $``$-valued positive multi-Toeplitz kernels with $`K(\mathrm{},\mathrm{})=1`$ and $`\mathrm{\Pi }_0()`$.

###### Proof.

This is a consequence of Theorem 3.2 and Corollary 2.5. A direct proof can be obtained as an application of Theorem 1.5.3 in . The parametrization given here is a refinement of one given in . ∎

An application gives another representation of the class $`𝒮(,_1,_2)`$. Motivated by the proof of Theorem 3.1 we consider the class $`𝒮_N(_1,_2,)`$ of systems

$$\mathrm{\Sigma }\{\begin{array}{ccc}\hfill x(n+1)& =& A(n)x(n)+B(n)u(n)\hfill \\ & & \\ \hfill y(n)& =& C(n)x(n)+D(n)u(n),n0,\hfill \end{array}$$

where $`_1`$, $`_2`$, $``$ are Hilbert spaces, $`A(0)(,^N)`$, $`B(0)(_2,^N)`$, $`C(0)(,_1)`$, $`D(_2,_1)`$, and $`A(n)=A(n1)^N`$, $`B(n)=B(n1)^N`$, $`C(n)=C(n1)^N`$, $`D(n)=D(n1)^N`$ (these relations are all of the type of (2.1)). Also we assume that for all $`n0`$, the operator $`\left[\begin{array}{cc}A(n)& B(n)\\ C(n)& D(n)\end{array}\right]`$ is a coisometry. The evolution of the system $`\mathrm{\Sigma }`$ can be encoded by the transfer map $`T_\mathrm{\Sigma }`$ that shows how the infinite input sequence $`(u(0),u(1),\mathrm{})`$ is transformed by the system into an output sequence $`(y(0),y(1),\mathrm{})`$ provided that $`x(0)=0`$. It is easy to see that the transfer map $`T_\mathrm{\Sigma }`$ is a lower triangular contraction with the matrix coefficients $`T_\mathrm{\Sigma }(n,n)=D(n)`$, $`T_\mathrm{\Sigma }(n+1,n)=C(n+1)B(n)`$ and for $`j>1`$, $`T_\mathrm{\Sigma }(n+j,n)=C(n+j)A(n+j1)\mathrm{}A(n+1)B(n)`$.

###### Theorem 3.4.

$`T`$ belongs to $`𝒮(,_1,_2)`$ if and only if $`T^{}`$ is the transfer map of a system in $`𝒮_N(_1,_2,)`$ for some Hilbert space $``$.

###### Proof.

We use the construction from the proof of Theorem 2.3.3 in . Thus, let $`\{\mathrm{\Gamma }_{ij}\}_{ij}`$ be the Schur parameters associated to $`T`$ in the proof of Theorem 3.2. For each $`i0`$ we consider the row contraction

$$L_i=\left[\begin{array}{cccc}\mathrm{\Gamma }_{ii}& D_{\mathrm{\Gamma }_{ii}^{}}\mathrm{\Gamma }_{i,i+1}& D_{\mathrm{\Gamma }_{ii}^{}}D_{\mathrm{\Gamma }_{i,i+1}^{}}\mathrm{\Gamma }_{i,i+2}& \mathrm{}\end{array}\right].$$

Let $`V_i`$ be the isometry associated to $`L_i`$ by the formula 1.6.3 in , that is,

(3.7)
$$V_i=\left[\begin{array}{cccc}\mathrm{\Gamma }_{ii}& D_{\mathrm{\Gamma }_{ii}^{}}\mathrm{\Gamma }_{i,i+1}& D_{\mathrm{\Gamma }_{ii}^{}}D_{\mathrm{\Gamma }_{i,i+1}^{}}\mathrm{\Gamma }_{i,i+2}& \mathrm{}\\ D_{\mathrm{\Gamma }_{ii}}& \mathrm{\Gamma }_{ii}^{}\mathrm{\Gamma }_{i,i+1}& \mathrm{\Gamma }_{ii}^{}D_{\mathrm{\Gamma }_{i,i+1}^{}}\mathrm{\Gamma }_{i,i+2}& \mathrm{}\\ 0& D_{\mathrm{\Gamma }_{i,i+1}}& \mathrm{\Gamma }_{i,i+1}^{}\mathrm{\Gamma }_{i,i+2}& \mathrm{}\\ 0& 0& D_{\mathrm{\Gamma }_{i,i+2}}& \\ \mathrm{}& \mathrm{}& & \mathrm{}\end{array}\right]=\left[\begin{array}{cc}\mathrm{\Gamma }_{ii}& B(i)^{}\\ C(i)^{}& A(i)^{}\end{array}\right].$$

$`V_0`$ is defined on the space $`_{k=0}^{\mathrm{}}𝒟_{\mathrm{\Gamma }_{0k}}^N`$ with values in in $`_{k=0}^{\mathrm{}}𝒟_{\mathrm{\Gamma }_{0k}}`$. We set $`=_{k=0}^{\mathrm{}}𝒟_{\mathrm{\Gamma }_{0k}}`$ and let $`\varphi :(_{k=0}^{\mathrm{}}𝒟_{\mathrm{\Gamma }_{0k}})^N_{k=0}^{\mathrm{}}𝒟_{\mathrm{\Gamma }_{0k}}^N`$ be the unitary operator defined by the formula $`\varphi (_{l=1}^N_{k=0}^{\mathrm{}}u_{kl})=_{k=0}^{\mathrm{}}_{l=1}^Nu_{kl}`$. Then we define

$$\left[\begin{array}{cc}\overline{D}^{}(i)& \overline{B}^{}(i)\\ \overline{C}^{}(i)& \overline{A}^{}(i)\end{array}\right]=\left[\begin{array}{cc}\mathrm{\Gamma }_{ii}& B^{}(i)\varphi ^{N^i}\\ (\varphi ^{N^{i1}})^{}C^{}(i)& (\varphi ^{N^{i1}})^{}A^{}(i)\varphi ^{N^i}\end{array}\right].$$

By Theorem 2.3.3 in , $`T_{i,i+1}^{}=C(i+1)B(i)=\overline{C}(i+1)\overline{B}(i)`$ and for $`j>i`$,

(3.8)
$$T_{ij}^{}=C(j)A(j1)\mathrm{}A(i+1)B(i)=\overline{C}(j)\overline{A}(j1)\mathrm{}\overline{A}(i+1)\overline{B}(i).$$

Since the elements of the family $`\mathrm{\Gamma }`$ satisfy (2.1) and the matrix coefficients of $`V_i`$ are monomials in $`\mathrm{\Gamma }_{ij}`$, their adjoints and defect operators, it follows that the system associated to the families $`\{\overline{A}(n)\}_{n0}`$, $`\{\overline{B}(n)\}_{n0}`$, $`\{\overline{C}(n)\}_{n0}`$ and $`\{\overline{D}(n)\}_{n0}`$ belongs to $`𝒮(,)`$ and $`T^{}`$ is the transfer map of this system. ∎

## 4. Recursive Constructions

As noticed in the proof of Theorem 3.4, the Schur parameters introduced in Theorem 3.2 are basically given by the parametrization of the unitary group. A convenient way to encode the dependence of the matrix coefficients of $`T`$ on the Schur paramaters is given by a continued fraction algorithm introduce for $`N=1`$ by Schur in . Here we describe such an algorithm for $`N>1`$. Of course, the main dificulty consists in keeping all the contructions within the tensor algebra. Due to the fact that there is a standard filtration on the tensor algebra, the contructions do not go through as in the one dimensional case, still we can obtain a sort of graded Schur type algorithm. For the algebra of upper triangular matrices, a Schur type algorithm was known (as a consequence of a similar algorithm for the spectral factorization of a family of matrices with displacement structure, see  or  for a detailed discussion). If we use it directly to our situation, it will produce the paramaters $`\{\gamma _\sigma \}_{\sigma 𝔽_N^+}`$ of Theorem 3.2, but the objects produced at each step are no longer in the tensor algebra. In this section we show how to remedy this situation.

We first notice another useful representation of a transfer map (we will assume in this section that $`_1=_1=`$). It is convenient to extend the symbol map (2.2). Thus, for $`T(_{k=0}^{\mathrm{}}_k_1,_{k=0}^{\mathrm{}}_k_2)`$ an upper triangular operator satisfying (2.1), we write

$$\mathrm{\Phi }(\left[T_{0j}\right]_{j=0}^{\mathrm{}})=T.$$

Now, if $`T_\mathrm{\Sigma }`$ is the transfer map of a system $`\mathrm{\Sigma }`$ in $`𝒮_N()=𝒮_N(,,)`$, then

$$\begin{array}{ccc}\hfill T_\mathrm{\Sigma }^{}& =& \mathrm{\Phi }(D(0)^{},0,\mathrm{})+\mathrm{\Phi }(0,B(0)^{},0,\mathrm{})\hfill \\ & & \\ & & \times (I\mathrm{\Phi }(0,A(0)^{},0,\mathrm{}))^1\mathrm{\Phi }(C(0)^{},0,\mathrm{}),\hfill \end{array}$$

provided that $`\mathrm{\Phi }(0,A(0)^{},0,\mathrm{})<1`$. We notice that this is the case if and only if $`A(0)<1`$. Let $`A(0)^{}=\left[A_j(0)^{}\right]_{j=1}^N`$, $`B(0)^{}=\left[B_j(0)^{}\right]_{j=1}^N`$. Then

$$\mathrm{\Phi }(0,A(0)^{},0,\mathrm{})=\underset{j=1}{\overset{N}{}}\mathrm{\Phi }(A_j(0)^{},0,\mathrm{})C_j^{}$$

and

$$\mathrm{\Phi }(0,B(0)^{},0,\mathrm{})=\underset{j=1}{\overset{N}{}}\mathrm{\Phi }(B_j(0)^{},0,\mathrm{})C_j^{},$$

so that

(4.1)
$$\begin{array}{ccc}\hfill T_\mathrm{\Sigma }^{}& =& \mathrm{\Phi }(D(0)^{},0,\mathrm{})+(_{j=1}^N\mathrm{\Phi }(B_j(0)^{},0,\mathrm{})C_j^{})\hfill \\ & & \\ & & \times (I_{j=1}^N\mathrm{\Phi }(A_j(0)^{},0,\mathrm{})C_j^{})^1\mathrm{\Phi }(C(0)^{},0,\mathrm{}),\hfill \end{array}$$

provided that $`A(0)<1`$. Next, we prove a refinement of (4.1) which is the key step of a Schur type algorithm. We introduce the translation operator given by the matrix

$$S=\left[\begin{array}{cccc}0& 0& & \mathrm{}\\ 1& 0& & \mathrm{}\\ 0& 1& 0& \mathrm{}\\ \mathrm{}& & \mathrm{}& \end{array}\right].$$

###### Theorem 4.1.

If $`T`$ belongs to $`𝒮()`$ then there exists an unique $`T_1𝒮(,𝒟_{T_{00}}^N,𝒟_{T_{00}^{}})`$ such that

$$\begin{array}{ccc}\hfill T& =& \mathrm{\Phi }(T_{00},0,\mathrm{})+\mathrm{\Phi }(D_{T_{00}^{}},0,\mathrm{})T_1S\hfill \\ & & \\ & & \times (I\mathrm{\Phi }(T_{00}^{},0,\mathrm{})T_1S\mathrm{\Phi }(D_{T_{00}},0,\mathrm{}).\hfill \end{array}$$

###### Proof.

Let $`\{\mathrm{\Gamma }_{ij}\}_{ij}`$ be the Schur parameters associated to $`T`$ as in the proof of Theorem 3.2. Then $`|\mathrm{\Gamma }_{00}|1`$. If $`|\mathrm{\Gamma }_{00}|=1`$, then we must have $`\mathrm{\Gamma }_{ij}=0`$ for $`j>i`$, so that $`T=\mathrm{\Phi }(\mathrm{\Gamma }_{00},0,\mathrm{})`$ and $`T_1=\mathrm{}=T_N=0`$. Assume $`|\mathrm{\Gamma }_{00}|<1`$. Then $`T^{}`$ is the transfer map of the system associated to the operators $`V_i`$ defined by (3.7). Thus, we have

$$\begin{array}{ccc}\hfill V_0& =& \left[\begin{array}{ccc}\mathrm{\Gamma }_{00}& B(0)^{}& \\ C(0)^{}& A(0)^{}& \end{array}\right]=\left[\begin{array}{cccc}\mathrm{\Gamma }_{00}& D_{\mathrm{\Gamma }_{00}^{}}\mathrm{\Gamma }_{01}& D_{\mathrm{\Gamma }_{00}^{}}D_{\mathrm{\Gamma }_{01}^{}}\mathrm{\Gamma }_{02}& \mathrm{}\\ D_{\mathrm{\Gamma }_{00}}& \mathrm{\Gamma }_{00}^{}\mathrm{\Gamma }_{01}& \mathrm{\Gamma }_{00}^{}D_{\mathrm{\Gamma }_{01}^{}}\mathrm{\Gamma }_{02}& \mathrm{}\\ 0& D_{\mathrm{\Gamma }_{01}}& \mathrm{\Gamma }_{01}^{}\mathrm{\Gamma }_{02}& \mathrm{}\\ 0& 0& D_{\mathrm{\Gamma }_{02}}& \\ \mathrm{}& \mathrm{}& & \mathrm{}\end{array}\right]\hfill \\ & & \\ & =& \left[\begin{array}{cccc}\mathrm{\Gamma }_{00}& D_{\mathrm{\Gamma }_{00}^{}}& 0& \mathrm{}\\ D_{\mathrm{\Gamma }_{00}}& \mathrm{\Gamma }_{00}^{}& 0& \mathrm{}\\ 0& 0& I& \\ \mathrm{}& \mathrm{}& & \mathrm{}\end{array}\right]\left[\begin{array}{cccc}I& 0& 0& \mathrm{}\\ 0& \mathrm{\Gamma }_{01}& D_{\mathrm{\Gamma }_{01}^{}}\mathrm{\Gamma }_{02}& \mathrm{}\\ 0& D_{\mathrm{\Gamma }_{01}}& \mathrm{\Gamma }_{01}^{}\mathrm{\Gamma }_{02}& \\ \mathrm{}& \mathrm{}& & \mathrm{}\end{array}\right]\hfill \\ & & \\ & =& \left[\begin{array}{ccc}\mathrm{\Gamma }_{00}& D_{\mathrm{\Gamma }_{00}^{}}& 0\\ D_{\mathrm{\Gamma }_{00}}& \mathrm{\Gamma }_{00}^{}& 0\\ 0& 0& I\end{array}\right]\left[\begin{array}{ccc}I& 0& 0\\ 0& \mathrm{\Gamma }_{01}& B^{}(0)^{}\\ 0& C^{}(0)^{}& A^{}(0)^{}\end{array}\right]\hfill \\ & =& \left[\begin{array}{ccc}\mathrm{\Gamma }_{00}& D_{\mathrm{\Gamma }_{00}^{}}\mathrm{\Gamma }_{01}& D_{\mathrm{\Gamma }_{00}^{}}B^{}(0)^{}\\ D_{\mathrm{\Gamma }_{00}}& \mathrm{\Gamma }_{00}^{}\mathrm{\Gamma }_{01}& \mathrm{\Gamma }_{00}^{}B^{}(0)^{}\\ 0& C^{}(0)^{}& A^{}(0)^{}\end{array}\right].\hfill \end{array}$$

We assume first that $`A(0)<1`$. Then formula (4.1) gives

$$\begin{array}{ccc}\hfill T& =& \mathrm{\Phi }(\mathrm{\Gamma }_{00},0,\mathrm{})+\mathrm{\Phi }(0,B(0)^{},0,\mathrm{})\hfill \\ & & \\ & & \times \left(I\mathrm{\Phi }(0,A(0)^{},0,\mathrm{})\right)^1\mathrm{\Phi }(C(0)^{},0,\mathrm{})\hfill \\ & & \\ & =& \mathrm{\Phi }(\mathrm{\Gamma }_{00},0,\mathrm{})+\left[\begin{array}{cc}\mathrm{\Phi }(0,D_{\mathrm{\Gamma }_{00}^{}}\mathrm{\Gamma }_{01},0,\mathrm{})& \mathrm{\Phi }(0,D_{\mathrm{\Gamma }_{00}^{}}B^{}(0)^{},0,\mathrm{})\end{array}\right]\hfill \\ & & \\ & & \times \left[\begin{array}{cc}I\mathrm{\Phi }(0,\mathrm{\Gamma }_{00}^{}\mathrm{\Gamma }_{01},0,\mathrm{})& \mathrm{\Phi }(0,\mathrm{\Gamma }_{00}^{}B^{}(0)^{},0,\mathrm{})\\ \mathrm{\Phi }(0,C^{}(0)^{},0,\mathrm{})& I\mathrm{\Phi }(0,A^{}(0)^{},0,\mathrm{})\end{array}\right]^1\left[\begin{array}{c}\mathrm{\Phi }(D_{\mathrm{\Gamma }_{00}},0,\mathrm{})\\ 0\end{array}\right].\hfill \end{array}$$

Since $`A(0)<1`$, we also have $`A^{}(0)<1`$ so that we can use a well-known formula for the inversion of a $`2\times 2`$ matrix in order to deduce that

$$\begin{array}{ccc}\hfill T& =& \mathrm{\Phi }(\mathrm{\Gamma }_{00},0,\mathrm{})+\mathrm{\Phi }(D_{\mathrm{\Gamma }_{00}^{}},0,\mathrm{})(\mathrm{\Phi }(0,\mathrm{\Gamma }_{00},0,\mathrm{})\mathrm{\Delta }^1\hfill \\ & & \\ & & +\mathrm{\Phi }(0,B^{}(0)^{},0,\mathrm{})(I\mathrm{\Phi }(0,A^{}(0)^{},0,\mathrm{}))^1\hfill \\ & & \\ & & \times \mathrm{\Phi }(0,C^{}(0)^{},0,\mathrm{})\mathrm{\Delta }^1)\mathrm{\Phi }(D_{\mathrm{\Gamma }_{00}},0,\mathrm{})\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill \mathrm{\Delta }& =& I\mathrm{\Phi }(0,\mathrm{\Gamma }_{00}^{}\mathrm{\Gamma }_{01},0,\mathrm{})\mathrm{\Phi }(0,\mathrm{\Gamma }_{00}^{}B^{}(0)^{},0,\mathrm{})\hfill \\ & & \\ & & \times (I\mathrm{\Phi }(0,A^{}(0)^{},0,\mathrm{}))^1\mathrm{\Phi }(0,C^{}(0)^{},0,\mathrm{})\hfill \\ & & \\ & =& I\mathrm{\Phi }(0,\mathrm{\Gamma }_{00}^{},0,\mathrm{})(\mathrm{\Phi }(0,\mathrm{\Gamma }_{00},0,\mathrm{})\hfill \\ & & \\ & & +\mathrm{\Phi }(0,B^{}(0)^{},0,\mathrm{})(I\mathrm{\Phi }(0,A^{}(0)^{},0,\mathrm{}))^1\mathrm{\Phi }(0,C^{}(0)^{},0,\mathrm{}))\hfill \end{array}$$

Since $`\mathrm{\Phi }(0,\mathrm{\Gamma }_{00}^{},0,\mathrm{})<1`$, $`\mathrm{\Delta }`$ is invertible and the previous calculation makes sense. Also, we define

$$R=\mathrm{\Phi }(0,\mathrm{\Gamma }_{01},0,\mathrm{})+\mathrm{\Phi }(0,B^{}(0)^{},0,\mathrm{})(I\mathrm{\Phi }(0,A^{}(0)^{},0,\mathrm{}))^1\mathrm{\Phi }(0,C^{}(0)^{},0,\mathrm{}).$$

Since $`\left[\begin{array}{cc}\mathrm{\Gamma }_{01}& B^{}(0)^{}\\ C^{}(0)^{}& A^{}(0)^{}\end{array}\right]`$ is a contraction (in fact, it is an isometry), we can easily check that $`R𝒮_0()`$. Actually, we can notice that $`R`$ is obtained from the Schur parameters $`\{\stackrel{~}{\mathrm{\Gamma }}_{ij}\}_{ij}`$, where $`\stackrel{~}{\mathrm{\Gamma }}_{ii}=0`$ and $`\stackrel{~}{\mathrm{\Gamma }}_{ij}=\mathrm{\Gamma }_{ij}`$ for $`j>i`$. We obtain that

(4.2)
$$T=\mathrm{\Phi }(\mathrm{\Gamma }_{00},0,\mathrm{})+\mathrm{\Phi }(D_{\mathrm{\Gamma }_{00}^{}},0,\mathrm{})R(I\mathrm{\Phi }(\mathrm{\Gamma }_{00}^{},0,\mathrm{})R)^1\mathrm{\Phi }(D_{\mathrm{\Gamma }_{00}},0,\mathrm{}).$$

Since $`\mathrm{\Phi }(\mathrm{\Gamma }_{00}^{},0,\mathrm{})<1`$, the previous formula makes sense without our assumption that $`A(0)<1`$. In fact, an approximation argument shows that (4.2) holds without that assumption. Thus, let $`\rho (0,1)`$ and consider $`T_\rho ^{}`$ the transfer map of the system $`\{\rho A(n),B(n),C(n),D(n)\}`$ (not necessarely in $`𝒮_N()`$). All the previous calculations go through and we obtain a contraction $`R_\rho `$ such that (4.2) holds for $`T_\rho `$. Letting $`\rho 1`$ we obtain (4.2). ∎

This result leads to the following Schur type algorithm for a $`T𝒮()`$.

$$\{\begin{array}{c}T_0=T,\gamma _n=T_{n,00},n0,\hfill \\ \\ T_{n+1}=\text{the unique solution of the equation}\hfill \\ T_n=\mathrm{\Phi }(\gamma _n,0,\mathrm{})\hfill \\ +\mathrm{\Phi }(D_{\gamma _n^{}},0,\mathrm{})T_{n+1}S^{N^n}(I\mathrm{\Phi }(\gamma _n^{},0,\mathrm{})T_{n+1}S^{N^n})^1\mathrm{\Phi }(D_{\gamma _n},0,\mathrm{}).\hfill \end{array}$$

The algorithm generates a sequence of elements $`T_n𝒮(,𝒟_{T_{n1}^N},𝒟_{T_{n1}^{}})`$ and it is easily seen that $`\gamma _0=\mathrm{\Gamma }_{0n}`$, $`n0`$. Therefore $`\{\gamma _n\}_{n0}`$ uniquely determine $`T`$ and the algorithm encodes the dependence of $`T`$ on these parameters.