text:
# Characterization of the channel diffusion in a sodium tetrasilicate glass via molecular-dynamics simulations

## Abstract

We study the structural and dynamical characteristics of the sodium atoms inside and outside the “diffusion channels” in glassy Na2O-4SiO2 (NS4) using classical molecular dynamics. We show that on average neither energetic arguments nor local environment considerations can explain the increased density of sodium atoms inside the subspace made of the channels. Nevertheless we show that at low temperature the mean square displacement of the sodium atoms inside this subspace is significantly larger than the one of the atoms outside the channels.

The diffusion of alkali atoms in silicate glasses is an important matter under investigation since several years. In particular the properties of sodium atoms inside the amorphous tetrahedral network of silica have been the topic of both experimental work  and molecular-dynamics simulations (MD)  since a simple glass like Na2O-4SiO2 (NS4) can be used as a prototype for more complicated glasses. The question of how the alkali atoms diffuse inside the glassy network is still a matter of debate and a generally accepted theory of ion transport in glasses is still missing . In a previous MD study of NS4 we have shown that the ions follow preferential pathways (“channels”) inside the glassy matrix . Nevertheless contrarily to the popular idea proposed by Greaves  and developed in further studies , these channels are neither static nor due to a microsegregation of the sodium atoms but have to be seen dynamically in the sense that the channels are those regions of space in which a great number of sodiums have passed during a given simulation time. The existence of these channels gives rise to a pre-peak in the structure factor at around $`q=0.95`$ Å-1 seen both in experiments  and classical MD simulations . In this last paper Horbach et al. show that the slow dynamics of the sodium atoms is closely related to the one of the underlying silica matrix which is coherent with the fact that there exists a strong correlation between the channels and the location of the non-bridging oxygen atoms, as shown in previous MD studies .
Once the existence of the channels is established, the next step is naturally to analyze their characteristics and to determine why the sodium atoms take these preferential pathways. This is the aim of the present study. In that direction we analyze the potential energy and the local structure of the sodium atoms whether they are INside the channels (NaIn) or OUTside the channels (NaOut). It is indeed generally believed that the diffusion of the ions occurs via “hopping motions between well-defined potential minima”  which we should be able to detect in our simulations. We analyze also the time evolution of the sodium densities inside and outside the channels which permits us to detect the dynamics of the channels as a function of temperature and to quantify the differences between the sodium densities, differences that so far have only been suggested . To elucidate the diffusion mechanisms of the Na atoms we determine a characteristic “residence” time inside the channels and we study the mean square displacements versus temperature of the NaIn and NaOut atoms.

Our system is made of 648 particles, namely $`N_{\mathrm{Na}}=`$ 86 sodium, 173 silicon and 389 oxygen atoms, confined in a cubic simulation box of edge length $`L=20.88`$ Å to which we apply periodic boundary conditions. The density is thus the experimental one, i.e. 2.38 g cm-3 . The interaction potential used in our simulations is a modified version of the one proposed by Kramer et al  and its complete description can be found in . It has been shown in previous studies that this potential is able to describe reliably many structural and dynamical properties of different sodium silicate melts  and especially NS4 . We have used as initial structure a $`\beta `$-cristobalite crystal in which we have randomly substituted the appropriate number of SiO4 tetrahedra by Na2O3 “molecules”. Subsequently the system is melted and equilibrated in the liquid phase at high temperature, i.e. 4000 K, for 50000 time steps ($`35`$ ps). Then it is cooled down very rapidly with a linear cooling schedule at a quench rate of $`2.3\times 10^{14}`$ K s-1. During the quenching process, the configurations of the system (positions and velocities of all the atoms) are saved at different temperatures ($`T`$ 4000, 3100, 2500, 2300, 1900 and 1700 K). These configurations are subsequently used as starting points of production runs of 2 million steps ($`2.8`$ ns) performed in the micro-canonical ensemble \[($`N`$$`V`$$`E`$$`=`$ const\]. During these production runs we save 2000 configurations equally spaced in time. The glass transition temperature of the system can be roughly estimated from the location of the bend in a plot of the potential energy vs $`T`$ to be around $`T_g`$ 2400 K. At each temperature the results have been averaged over three independent samples.

As mentioned above, we have shown in previous studies  that at low temperature the sodium atoms diffuse through the (quasi-)frozen silica matrix within a subset of the total available space. This subset is made of pockets (typical size 3-6 Å) connected together via small pathways (similar in a sense to a neuronal network) and this network is what we call channels. Our aim in this work is to characterize the channel diffusion of the sodium atoms, i.e. to study the possible structural and dynamical differences between the NaIn and the NaOut atoms. In order to determine the channels, we divide the simulation box in $`N_{tot}=20^3`$ small distinct cubes of volume $``$ 1 Å3 and determine the number density of the Na atoms in each of these cubes during the 2.8 ns of the simulation. Then we consider only the upper 10 % of all the visited cubes. These cubes form “the core” of the channels (the space in which the sodium atoms have been the more often) or shortly “the channels”. More details can be found in Ref. . Once the channels have been defined via this time integration, we can go back and determine for every time $`t`$, NaIn and NaOut and start our investigations.
Why do the Na atoms visit a specific fraction of the total available space more frequently at low temperature? The first explanation that comes to mind is that the visited sites are more favorable energetically. To check this idea, we have calculated the individual potential energy distributions of the sodium atoms at 1700 K whether they are inside or outside the channels. The two distributions have been averaged over the 2000 configurations saved at 1700 K and over the 3 independent samples and are represented in Fig. 1. It is obvious in Fig. 1 that there is no significant difference between the distribution of NaIn and NaOut: both are gaussians centered approximately at $``$3.8 eV with a Half Width at Half Maximum of $``$ 0.6 eV. It appears therefore that the energy is not the driving force explaining why the Na atoms are inside the channels. One may argue that the time average can smear out the effect but in our previous study we have shown that the sodium atoms visit sites previously occupied by other sodium atoms  and therefore if these sites were favorable energetically it should appear in Fig. 1. We have also calculated the energy distributions of Si and O atoms (not shown in Fig. 1 ) and found that they have an average energy of $``$35 eV and $``$12 eV respectively.
In fact, the individual potential energy distributions reflect the whole environment of the Na atoms. Therefore even if there exists a structural difference at short range, it will not be visible in Fig. 1 because of all the other (long-range) interactions. It is hence justified to look more precisely at the environment of the sodium atoms via the radial pair distribution function and the integrated number of neighbors. For a given $`\alpha \beta `$ pair they are defined by:

$`g(r)_{\alpha \beta }={\displaystyle \frac{V}{4\pi r^2N_\alpha dr}}dn_\beta `$ (1)
and,

$`N(r)_{\alpha \beta }={\displaystyle \frac{N_\alpha }{V}}{\displaystyle _0^r}4\pi r^2g(r)_{\alpha \beta }𝑑r`$ (2)
Note that the number of first nearest $`\beta `$ neighbors around species $`\alpha `$, is given by the value of $`N(r)_{\alpha \beta }`$ at the first minimum of $`g_{\alpha \beta }(r)`$. In Fig. 2 we have represented $`g(r)`$ and $`N(r)`$ at $`T`$ = 1700 K for (a) Na$`{}_{\mathrm{Out}}{}^{}`$Na, Na$`{}_{\mathrm{In}}{}^{}`$Na and Na$``$Na pairs, (b) Na$`{}_{\mathrm{Out}}{}^{}`$O, Na$`{}_{\mathrm{In}}{}^{}`$O and Na$``$O pairs and (c) Na$`{}_{\mathrm{Out}}{}^{}`$Si, Na$`{}_{\mathrm{In}}{}^{}`$Si and Na$``$Si pairs. The major differences can be seen in (a) where it appears that a NaIn has on average 1 supplemental Na neighbor compared to a NaOut atom. This goes together with a reduction of the number of nearest Si neighbors (c). On the contrary the local oxygen environment (b) seems to be the same for the sodiums inside and outside the channels. Finally concerning the nearest neighbor distances with the other species there is no significant change between NaIn and NaOut except a slight decrease (0.1 Å) of the Na$`{}_{\mathrm{In}}{}^{}`$Na distance compared to the distance between NaOut and the other sodium atoms. It appears therefore that the local structure of the sodium atoms inside and outside the channels is not significantly different in order to explain the existence of the channels. Nevertheless it appears from the study of the radial pair distribution functions that on average the NaIn atoms have more Na neighbors than the NaOut atoms. Therefore it is justified to calculate at any instant $`t`$ the number of sodiums inside the channels ($`N_{\mathrm{Na}}^{\mathrm{In}}(t)`$) and the number of sodiums outside the channels ($`N_{\mathrm{Na}}^{\mathrm{Out}}(t)`$). These quantities are shown in Fig. 3 at 1700 K (a), 2300 K (b) and 3100 K (c). A first information can be obtained from the shape of the curves at the different temperatures. At 1700 K (a), the fact that the curves become flat after $`0.5`$ ns indicates that the channels are well defined through the frozen silica matrix and are made almost of the same small cubes during the rest of the simulation. At 2300 K (b), the situation is rather different since $`N_{\mathrm{Na}}^{\mathrm{In}}`$ goes through a maximum (and consequently $`N_{\mathrm{Na}}^{\mathrm{Out}}`$ goes through a minimum) after $`1.4`$ ns. This is a signature of the residual motion of the silica matrix which is only quasi-frozen and hence a manifestation of the motion of the channels. In fact, the channels that we have defined in our calculation at that temperature are the common part of the “moving” channels we could have defined over consecutive fractions of the whole simulation. It is therefore reasonable to assume that the real channels at $`1.4`$ ns (half the simulation length) are closest to the channels that we have defined over the whole simulation. At 3100 K (c), the fact that the curves are flat indicates that the channels do not really exist anymore. What we have defined as channels at that temperature are only randomly distributed clusters of density fluctuations. Actually, the matrix is melted and the sodium atoms may visit the entire accessible space (the system is ergodic) and hence the location of the clusters is completely random and depends crucially on the length of the run. Thus, the definition of the channels is quite artificial at high temperature. Nevertheless when we compare the values of $`N_{\mathrm{Na}}^{\mathrm{In}}`$ after $`1.4`$ ns for the different temperatures we do not see an increase of $`N_{\mathrm{Na}}^{\mathrm{In}}`$ with increasing temperature contrarily to what was suggested by Meyer et al.  to explain an increase of the quasielastic amplitude in the dynamic structure factor of sodium disilicate \[NS2\] melts. Finally we see in Fig. 3 that for all the temperatures the major part of the sodium atoms is outside the channels. Nevertheless as shown in  the total volume of the channels is small ($``$ 10 % at 2000 K) therefore the sodium densities can be very different inside and outside the channels. In Fig. 4 we have represented the sodium density $`\rho _{\mathrm{Na}}`$ inside and outside the channels (which is simply $`N_{\mathrm{Na}}^{\mathrm{In}}/V_{chan}`$ and $`N_{\mathrm{Na}}^{\mathrm{Out}}/(L^3V_{chan})`$ where $`V_{chan}`$ is the volume of the channels) as a function of $`T`$. It is immediately noticeable that the sodium density inside the channels is higher than the one outside and in particular at low temperature (below $`T_g`$). So far we knew that the channels are a subspace highly visited by the sodium atoms (they have been defined this way) but now we can conclude in addition that they are also a subspace of high sodium density. We emphasize that a higher sodium concentration does not imply necessarily a clustering of alkalis . Actually we have shown that no clustering can be observed on a single snapshot of the simulation box . It is also noticeable that the two curves converge towards the limit of uniform density at high temperature. This is another evidence that at high temperature the channels do not exist and that what we have defined as channels are randomly distributed clusters.
Once we know the position of the sodium atoms with respect to the channels it is naturally of interest to calculate a mean “residence” time of the Na atoms within the channels. This can be done by determining the probability $`P(0,t)`$ of a sodium atom inside the channels at $`t=0`$ to be inside the channels after a time $`t`$. $`P(0,t)`$ is given by:

$`P(0,t)={\displaystyle \frac{1}{N_{\mathrm{Na}}^{\mathrm{In}}(0)}}{\displaystyle \underset{i=1}{\overset{N_{\mathrm{Na}}^{\mathrm{In}}(0)}{}}}n_i(t)`$ (3)
with:

$`\{\begin{array}{c}n_i=0\text{ outside the channels}\hfill \\ n_i=1\text{ inside the channels}\hfill \end{array}`$ (6)
This probability is represented in Fig. 5 at different temperatures. Since we save the configurations every 1.4 ps (1000 MD steps), the fast decay from 1 at shorter times can not be seen in Fig. 5 but since the channels have been determined with this “time step” we have decided to calculate $`P(0,t)`$ in the same conditions. This means therefore that our probability is slightly overestimated since we do not take into account the vibrational effects.
In Fig. 5 we see that $`P(0,t)`$ decreases and converges towards the long time limit $`N_{\mathrm{Na}}^{\mathrm{In}}/N_{\mathrm{Na}}`$ very rapidly at high temperature while at low temperature ($`T2300`$ K) this limit is not reached after 1.5 ns. After subtracting the long time limit, these curves can be fitted reasonably well by an exponential function with a time constant $`\tau `$ which is represented versus 1/T in a linear-log plot in the inset of Fig. 5. From this inset it is clear that $`\tau (T)`$, the typical residence time, shows an Arrhenius behavior with an activation energy around 1.45 eV. This activation energy is larger than the activation energy of the Na diffusion process and other more specific characteristics of the Na displacement (hops between preferential sites, forward-backward jumps, decorrelation mechanism) which is around 1.3 eV . This seems to indicate that the mechanism governing the expulsion of the Na ions outside the channels is energetically more costly than the other mechanisms controlling the mean displacement of the ions.
Finally we have decided to investigate the influence of the location of the sodium atoms on their diffusion by calculating “conditional” mean squared displacements (MSD) ($`R^2(t)=\left|\stackrel{}{r}(t)\stackrel{}{r}(0)\right|^2`$). On the one hand we have calculated $`d^2(t)_{\mathrm{In}/\mathrm{Out}}`$ which is the “MSD” of the atoms inside/outside the channels at time $`t=0`$ and $`t0`$ (this is a reasonable calculation since we know that $`P(0,t)`$ remains non zero over the length of the simulation). These quantities are represented in Fig. 6 for (a) $`1700`$ K, (b) $`1900`$ K and (c) $`2300`$ K. On the other hand we have calculated $`R^2(t)_{\mathrm{In}/\mathrm{Out}}`$ strictly speaking: in that case the atoms are continuously inside/outside the channels between $`0`$ and $`t`$. In this latter case the number of steps an Na atom is consecutively inside/outside the channels is small especially at high temperature therefore this quantity is only represented for $`1700`$ K in the inset of Fig. 6 (a). From the results represented in Fig. 6 we see that at low temperature $`(T=`$1700, 1900 K) $`d_{\mathrm{In}}^2`$ is always larger than $`d_{\mathrm{Out}}^2`$ (at long times the two curves collapse since the history of the In and Out atoms becomes the same on average). This difference is maximum after $`0.1`$ ns when $`d_{\mathrm{In}}^2`$ is twice as important than $`d_{\mathrm{Out}}^2`$. This indicates clearly that the Na atoms are more mobile inside the channels. This result is coherent with the findings of Horbach et al.  who state that the “sodium atoms move quickly between preferential sites” but has never been quantified directly so far. It is confirmed in the inset of Fig. 6 (a) by the difference between the MSD of the Na atoms that are moving only inside the channels and the MSD of the sodium atoms moving only outside the channels. At $`T=2300`$ K (Fig. 6 (c)), we see that the difference between $`d_{\mathrm{In}}^2`$ and $`d_{\mathrm{Out}}^2`$ vanishes even though the channels still exist. Nevertheless as shown in Fig. 3 at that temperature the channels are not static anymore and therefore the distinction between In and Out has no clear meaning. At even higher temperatures (not shown), this is of course also true, and the two curves are superimposed. We can hence conclude that when the channels are formed and “static” (lowest temperatures), the sodium diffusivity is higher when the Na atoms are inside the channels.

In summary, with the use of classical molecular dynamics simulations on Na2O-4SiO2 systems we have analyzed in more detail the properties of the sodium trajectories inside the glassy matrix and especially the differences between the Na atoms inside the diffusion channels and those outside of these channels (the existence of these channels defined as the fraction of space mostly visited by the Na atoms has been shown in a previous study ).
Firstly we have shown that the average potential energy of the Na atoms is the same whether they are In or Out. This is also true concerning their local environment except a slightly higher Na coordination number for the sodium atoms inside the channels. This is coherent with the fact that, in comparison with the rest of the system, the concentration of sodium atoms at low temperature (T $``$ 3000 K) is higher inside the channels. At high temperature the matrix can not be considered as frozen and therefore the channels do not exist anymore. We have defined a characteristic residence time inside the channels $`\tau (T)`$ that shows an Arrhenius dependence with an activation energy around 1.45 eV. This value is slightly higher than the characteristic energy of the other activated processes present in this system. Finally, we have shown that, at low temperature when the channels can be considered as frozen, the mobility of the sodiums inside the channels is larger than the one of the Na atoms in the rest of the system: in that sense one can really speak about channel diffusion.
In order to find the origin of this effect and in particular the role of the SiO2 matrix, one might change the vibrational characteristics of the matrix or those of the cation and see how these changes affect the sodium dynamics. These are investigations currently under way.

We thank W. Kob for interesting discussions. Part of the numerical calculations were performed at “Centre Informatique National de l’Enseignement Supérieur” in Montpellier, France.


format: text

text:
# Effects of the magnetic fields on the protoneutron star structure.

The equation of state of a protoneutron star matter including the effects of magnetic field, temperature, nuclear asymmetry and trapped neutrinos are studied. Using the Oppenheimer-Volkoff-Tolman equation global parameters of the protoneutron star are obtained. Presence of the nuclear asymmetry and magnetic field make the protoneutron star more compact.

: J. Phys. G: Nucl. Phys.

Properties of dense matter in strong magnetic field has been the subject of investigations in astrophysics of white dwarfs, protoneutron and neutron stars. It is motivated by the fact that magnetic fields of the order of $`10^8`$ $`G`$ \- $`10^{18}`$ $`G`$ are known to exist in many cases of these objects. Properties of protoneutron stars, under various assumptions concerning composition and equation of state of hot lepton rich matter were studied by many authors (Strobel et al. 1999, Bombaci et al. 1995, Takatsuka 1995, Mańka et al. 2001) . Matter inside a protoneutron star is highly degenerate and chemical potential of its constituents are a few hundreds of $`MeV`$. The strength of magnetic field of a protoneutron star changes from $`10^8`$ $`G`$ at the surface up to $`10^{18}`$ $`G`$ in the center of the star. The interior magnetic fields are a few orders of magnitude larger than the surface ones .
The aim of this letter is to study the basic properties of the protoneutron stars, such as masses and radii. The construction of protoneutron star model is based on various realistic equation of states and results in the general picture of protoneutron star interiors. Thus the proper form of the equation of state is essential in determining protoneutron star properties, such as the mass-range and the mass-radius relation. The complete and more realistic description of a protoneutron star requires taking into consideration effects of finite temperature and nonzero magnetic fields. Using the relativistic mean-field theory approach the adequate form of the equation of state enlarged by contributions coming from magnetic field and temperature is constructed and serves as an input to the Oppenheimer-Volkoff-Tolman equations. The theory considered here comprises electrons, neutrinos, scalar, vector-scalar and vector-isovector mesons. The RMF theory implies that the nucleon interactions appear through the exchange of meson fields . This theory is very useful in describing nuclear matter and finite nuclei. Its extrapolation to large charge asymmetry is of considerable interest in nuclear astrophysics and particulary in constructing protoneutron and neutron star models where extreme conditions of isospin are realized. The Lagrangian of the theory contains baryon, meson and lepton degrees of freedom and can be represented as the sum

$$=_B+_L+_{}+_G+_{QED},$$
(1)
where $`_B`$, $`_L`$, $`_{}`$, $`_G`$ describe the baryonic, leptonic, mesonic and gravitational terms, respectively. The $`_{QED}`$ is the Lagrangian density function of the QED theory. The fermion fields are composed of neutrons, protons, electrons and neutrinos. In strong magnetic field, contributions from the anomalous magnetic moments of protons and neutrons ($`\mu _b`$) must be also considered . The baryonic term of the Lagrangian density function is given by

$`_B=i\overline{\psi }\gamma ^\mu D_\mu \psi \overline{\psi }(Mg_\sigma \phi )\psi +\mu _b\overline{\psi }\sigma _{\mu \nu }F^{\mu \nu }\psi ,`$ (2)
where the covariant derivative $`D_\mu `$ is defined as

$`D_\mu =_\mu +{\displaystyle \frac{1}{2}}ig_\rho \rho _\mu ^a\sigma ^a+ig_\omega \omega _\mu +iQA_\mu .`$ (3)
The anomalous magnetic moments introduced via the minimal coupling of the nucleons to the electromagnetic field tensor for each baryon give contributions to the Lagrangian density function which are contained in the $`\mu _b\overline{\psi }\sigma _{\mu \nu }F^{\mu \nu }\psi `$ term , where $`\sigma _{\mu \nu }=\frac{i}{2}\{\gamma _\mu ,\gamma _\nu \}`$.

The leptonic part of the total Lagrangian has the form

$`_L=i{\displaystyle \underset{f}{}}\overline{L}_f\gamma ^\mu \stackrel{~}{D}_\mu L_f{\displaystyle \underset{f}{}}g_f(\overline{L}_fHe_{Rf}+h.c),`$ (4)
where $`\stackrel{~}{D}_\mu =_\mu ieA_\mu `$ and $`e`$ is the electron charge. The vector potential is given by $`A_\mu =\{A_0=0,A_i\}`$ where

$`A_i={\displaystyle \frac{1}{2}}\epsilon _{ilm}x^lB_0^m.`$
The gauge in which uniform magnetic field $`B`$ lies along the z-axis was chosen $`B_0^m=(0,0,B_z)`$. The mesonic part of the Lagrangian density function consisting of isoscalar (scalar $`\sigma `$, vector $`\omega `$) and isovector (vector $`\rho `$) mesons is given by

$`_{}={\displaystyle \frac{1}{2}}_\mu \phi ^\mu \phi U(\phi ){\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }{\displaystyle \frac{1}{4}}W_{\mu \nu }W^{\mu \nu }{\displaystyle \frac{1}{2}}M_\omega ^2\omega _\mu \omega ^\mu `$
$`{\displaystyle \frac{1}{4}}c_3(\omega _\mu \omega ^\mu )^2{\displaystyle \frac{1}{4}}R_{\mu \nu }^aR^{a\mu \nu }{\displaystyle \frac{1}{2}}M_\rho ^2\rho _\mu ^a\rho ^{a\mu },`$
where $`F_{\mu \nu }`$ is the electromagnetic stress tensor, $`W_{\mu \nu }`$ and $`R_{\mu \nu }^a`$ are vector meson fields strength. The potential function $`U(\phi )`$ has a very well known form introduced by Boguta and Bodmer . The parameters used in this model are collected in Table 1.

In the mean-field approximation baryon currents in field equations are replaced by their ground states expectation values. Mesons fields are replaced by their mean values $`\sigma =<\phi >`$, $`w_\mu =<\omega _\mu >=\delta _{\mu ,0}w`$ and $`r_\mu ^a=<\rho _\mu ^a>=\delta ^{a,3}\delta _{\mu ,0}r`$. The Dirac equation for the nucleon quasiparticle has the form

$$(i\gamma ^\mu \overline{D}_\mu m_{eff})\psi =0$$
(5)
with the covariant derivative given by

$$\overline{D}_\mu =_\mu +\frac{1}{2}ig_\rho r_\mu ^a\sigma ^a+ig_\omega w_\mu +iQA_\mu $$
(6)
The effective nucleon mass $`m_{eff}`$ equals

$$m_{eff}=M\delta =Mg_s\sigma .$$
(7)
This redefine the proton and neutron chemical potentials

$`\mu _p=\epsilon _p+g_\omega w+{\displaystyle \frac{1}{2}}g_\rho r`$ (8)
$`\mu _n=\epsilon _n+g_\omega w{\displaystyle \frac{1}{2}}g_\rho r`$ (9)
where $`\epsilon _f=\sqrt{k^2+m_{eff}^2}`$ ($`f=\{p,n\}`$). Having obtained the effective nucleon chemical potentials one can describe the chemical equilibrium of the system which is imposed through relations between chemical potentials. Protoneutron star matter is assumed to be in beta equilibrium and electrically neutral

$`p+en+\nu _e.`$ (10)
This can be expressed as a relation between the chemical potentials of the protoneutron star constituents

$`\mu _{\nu _e}=\mu _e+\mu _p\mu _n.`$ (11)
Using relations (8) and (9) the neutrino chemical potential protoneutron star matter is obtained

$`\mu _{\nu _e}=\mu _e+\epsilon _p\epsilon _n+g_\rho r.`$
The charge neutrality means that $`n_e=n_p`$. The assumption that only $`\nu _e`$ are captured inside the star core was made. The proton fraction is defined by the physical conditions in the star and is determined by beta equilibrium and by nuclear asymmetry. Effects due to muon neutrino and antineutrino were completely ignored because muons cannot be produced in the low-density medium. In this letter the effects of moderate magnetic field on the equation of state of a relativistic, degenerate gas is considered. The well known dispersion relation for baryons takes the form

$`E_{n,p_z,s}^b=\sqrt{p_z^2+(\sqrt{m_{eff}^2+2nQB_z}+s\mu _bB_z)^2}`$ (12)
and for small values of magnetic fields strength the equation (12) is given by

$`E_{n,p_z}^b\sqrt{p_z^2+m_{eff}^2+2nQB_z}.`$
Along the field the particle motion is free and quasi-one-dimensional with the modified density of states. The fermion density of states in the absence of magnetic field is replaced by the sum

$`2{\displaystyle \frac{d^3p}{(2\pi )^3}}{\displaystyle \underset{s}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}[2\delta _{n0}]{\displaystyle \frac{eB_z}{(2\pi )^2}𝑑p_z},`$
where the symbol $`\delta _{n0}`$ denotes the Kronecker delta  and thus the spin degeneracy equals 1 for the ground ($`n=0`$) Landau level and 2 for $`n1`$. The total pressure of the system can be described as the sum of pressure coming from fermions, mesons plus corrections coming from electromagnetic field

$`P=P_f+P_{}+P_{QED}.`$
The contributions from the same constituents form the energy density

$`\epsilon =\epsilon _f+\epsilon _{}+\epsilon _{QED},`$
where $`\epsilon _{QED}=P_{QED}=\frac{1}{2}B^2`$. The fermion pressure is defined as

$`P_f={\displaystyle \frac{\gamma _im_i^4}{4\pi ^2}}{\displaystyle \underset{s}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}[2\delta _{n0}](I_{1,2,+}(z_i/t_i,(\sqrt{\delta _i^2+2\gamma _in}+s{\displaystyle \frac{\mu _b}{Q_i}}\gamma _i)^2)`$
$`+I_{1,2,}(z_i/t_i,(\sqrt{\delta _i^2+2\gamma _in}+s{\displaystyle \frac{\mu _b}{Q_i}}\gamma _i)^2)`$
where the Fermi integral

$$I_{\lambda ,\eta \pm }(u,\alpha )=\frac{(\alpha +x^2)^{\lambda /2}x^\eta dx}{e^{\left(\sqrt{\alpha +x^2}u\right)}+1}$$

was used, $`z_i=\mu _0/m_i`$, $`t_i=k_BT_0/m_i`$, $`u_i=z_i/t_i`$ and $`\gamma _i=B_z/B_{c,i}=Q_iB/m_i^2`$. The critical magnetic field strength for electron equals $`B_c=m_e^2/e`$= $`4.414\times 10^{13}`$ $`G`$ and for proton $`B_c10^{20}`$ $`G`$. The index $`i`$ denotes fermions.

The fermion energy density is defined with the use of the Fermi integral

$`\epsilon _f={\displaystyle \frac{\gamma _im_i^4}{4\pi ^2}}{\displaystyle \underset{s}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}[2\delta _{n0}](I_{1,0,+}(z_i/t_i,(\sqrt{\delta _i^2+2\gamma _in}+s{\displaystyle \frac{\mu _b}{Q_i}}\gamma _i)^2)`$
$`+I_{1,0,}(z_i/t_i,(\sqrt{\delta _i^2+2\gamma _in}+s{\displaystyle \frac{\mu _b}{Q_i}}\gamma _i)^2).`$
The mass-radius relation for the protoneutron star with and without magnetic field is presented in figure 1. The temperature is fixed and equals 30 $`MeV`$, the lepton fraction is assumed $`Y_L=0.4`$ and the proton fraction $`Y_P0.3`$. These conditions are adequate for the trapped neutrino case. In this figure the mass-radius relations for protoneutron star with and without field $`\rho `$ are depicted. The magnetic field causes that these objects are more compact and similar to a ordinary neutron star. Comparing protoneutron star parameters namely masses and radii one can go to the conclusion that protoneutron star model calculated in the presence of magnetic field results with smaller radii and bigger masses than the one obtained in the absence of magnetic field. The neutrino chemical potential as a function of the central density is shown in figure 2. One can compare the influence of the nuclear asymmetry which is driven by the meson $`\rho `$. Its presence in the protoneutron star matter diminishes the value of neutrino chemical potential and decreases the protoneutron star mass. This feature strongly depends on the nuclear interactions inside a protoneutron star. Figure 3 depicts the comparision of neutrino chemical potentials calculated with and without the influence of asymmetry with the electron one. This picture shows the profiles of the mentioned above chemical potentials inside the star. The nuclear asymmetry diminishes the value of neutrino chemical potential inside the star and thus lowers its density. The influence of asymmetry is the most distinctive in the high density region inside a protoneutron star. It was examined for moderate strength of magnetic field. In conclusion, the magnetic field and nuclear asymmetry make the protoneutron star more compact astrophisics object. ————————————————–


text:
# 1 Introduction

## 1 Introduction

The notion of strict positive realness (SPR) of transfer functions plays an important role in absolute stability theory, adaptive control and system identification\[1-23\]. Motivated by Kharitonov’s seminal theorem on the robust stability for a family of polynomials, a number of recent papers has concentrated on the strict positive realness for a family of transfer functions. In the spirit of Kharitonov, the robust SPR analysis and design problems were first formulated by Dasgupta and Bhagwat. They showed that every transfer function in an interval transfer function family is strictly positive real if and only if sixteen prescibed vertex transfer functions in this family are strictly positive real. The sixteen critical vertex transfer functions can be constructed explicitly using Kharitonov’s four vertex polynomials. This result was subsequently improved by Chapellat and Bhattacharyya, Wang and Huang, where only eight out of the sixteen critical vertex transfer functions need to be checked\[11-12\]. For a family of transfer functions with affine linearly correlated perturbations, or more generally, multilinearly correlated perturbations, Dasgupta, Anderson et al. showed that it suffices to check all vertices in order to ensure the strict positive realness of the entire family. By resort to the concept of positive polynomial pairs and root interlacing properties, Hollot and Huang solved the robust SPR design problem for low order and structured families\[9-10\]. Anderson et al. considered the general robust SPR design problem, and by using the Hilbert transform, provided a constructive method. Betser and Zeheb made some further improvements.

This paper characterizes SPR regions in coefficient space and presents a robust design method for SPR transfer functions. We first introduce the concepts of SPR regions and weak SPR regions and give a complete characterization of them. We show that the SPR region associated with a fixed polynomial is unbounded, whereas the weak SPR region is bounded. We then prove that the intersection of several weak SPR regions associated with different polynomials can not be a single point. Furthermore, we show how to construct a point in the SPR region from a point in the weak SPR region. Based on these theoretical development, we propose an algorithm for robust design of SPR transfer functions. This algorithm works well for both low order and high order polynomial families. Illustrative examples are provided to show the effectiveness of this algorithm.

## 2 Preliminaries

Denote $`P^n`$ as the $`n`$-th order real polynomial family, $`R^n`$ as the $`n`$ dimensional real field, and $`H^nP^n`$ as the set of all $`n`$-th order Hurwitz stable polynomials.

In the following definitions,$`b()P^m,a()P^n,`$ and $`p(s)=b(s)/a(s)`$ is a rational function.

Definition 1 $`p(s)`$ is said to be strictly positive real(SPR),denoted as $`p(s)`$SPR,if $`b(s)P^n,a(s)H^n,`$ and Re$`[p(j\omega )]>0,`$ $`\omega R.`$

Definition 2 $`p(s)`$ is said to be weak SPR (WSPR),denoted as $`p(s)`$WSPR,if $`b(s)P^{n1},a(s)H^n,`$ and Re$`[p(j\omega )]>0,`$ $`\omega R.`$

Definition 3 Given $`a(s)H^n,`$ the set of the coefficients (in $`R^{n+1}`$) of all the $`b(s)`$’s in $`P^n`$ such that $`p(s):={\displaystyle \frac{b(s)}{a(s)}}`$SPR is said to be the SPR region associated with $`a(s)`$,denoted as $`\mathrm{\Omega }_a.`$

Definition 4 Given $`a(s)H^n,`$ the set of the coefficients (in $`R^n`$) of all the $`b(s)`$’s in $`P^{n1}`$ such that $`p(s):={\displaystyle \frac{b(s)}{a(s)}}`$WSPR is said to be the WSPR region associated with $`a(s)`$,denoted as $`\mathrm{\Omega }_a^W.`$

For notational convenience, $`\mathrm{\Omega }_a(\mathrm{\Omega }_a^W)`$ sometimes also stands for the set of all the polynomials $`b(s)`$ in $`P^n(P^{n1})`$ such that $`p(s):={\displaystyle \frac{b(s)}{a(s)}}`$SPR(WSPR).

From the definitions above, it is easy to get the following properties:

Property 1 If $`p(s)`$SPR(WSPR), then $`|arg(b(j\omega ))arg(a(j\omega ))|<{\displaystyle \frac{\pi }{2}}`$ , $`\omega R,`$ where $`arg()`$ stands for the argument of the complex number,and the difference of two arguments can differ by an integer number of $`2\pi `$.

Property 2 Given $`a(s)H^n,`$ $`\mathrm{\Omega }_a`$ is a non-empty, open, convex cone in $`R^{n+1}`$.

Property 3 Given $`a(s)H^n,`$ we have $`\mathrm{\Omega }_aH^n,\mathrm{\Omega }_a^WH^{n1}.`$

The problem we are interested in is: Given a family of Hurwitz stable polynomials, how can find a fixed polynomial such that their ratios will be SPR-invariant? In what follows, we will first give some characterization of WSPR regions, and then propose an efficient design procedure for this problem.

## 3 Geometric Characterization of SPR Regions

By definition, an SPR (WSPR) transfer function times a positive integer is still SPR (WSPR). Thus, without loss of generality, let

$$a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n$$
(1)
Denote as $`\mathrm{\Omega }_{1a}`$ the set of the coefficients of all the $`b(s)=s^n+x_1s^{n1}+\mathrm{}+x_nP^n,`$ i.e., $`(x_1,x_2,\mathrm{},x_n)`$ in $`R^n`$ , such that $`p(s)={\displaystyle \frac{b(s)}{a(s)}}`$SPR; and denote as $`\mathrm{\Omega }_{1a}^W`$ the set of the coefficients of all the $`b(s)=s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}P^{n1},`$ i.e., $`(x_1,x_2,\mathrm{},x_{n1})`$ in $`R^{n1}`$ , such that $`p(s)={\displaystyle \frac{b(s)}{a(s)}}`$WSPR. Obviously, we have

$`\{1\}\times \mathrm{\Omega }_{1a}=\left\{(1,{\displaystyle \frac{b_1}{b_0}},{\displaystyle \frac{b_2}{b_0}},\mathrm{},{\displaystyle \frac{b_n}{b_0}})\right|`$
$`(b_0,b_1,b_2,\mathrm{},b_n)\mathrm{\Omega }_a\}`$ (2)
$`\{1\}\times \mathrm{\Omega }_{1a}^W=\left\{(1,{\displaystyle \frac{b_1}{b_0}},{\displaystyle \frac{b_2}{b_0}},\mathrm{},{\displaystyle \frac{b_{n1}}{b_0}})\right|`$
$`(b_0,b_1,b_2,\mathrm{},b_{n1})\mathrm{\Omega }_a^W\}`$ (3)
For notational convenience, $`\mathrm{\Omega }_{1a}(\mathrm{\Omega }_{1a}^W)`$ sometimes also stands for the corresponding polynomial set.

As we know, $`\mathrm{\Omega }_a`$ is a non-empty, open, convex cone in $`R^{n+1}`$. Thus, $`\mathrm{\Omega }_a`$ is an unbounded set in $`R^{n+1}`$. In what follows, we will show that $`\mathrm{\Omega }_{1a}`$ is also an unbounded set in $`R^n`$.

Theorem 1 Given $`a(s)H^n,`$ $`\mathrm{\Omega }_{1a}`$ is a non-empty, open, unbounded convex set in $`R^n`$.

Proof Obviously, we have $`a(s)\mathrm{\Omega }_a.`$ If the leading coefficient of $`a(s)`$ is $`a_0,`$ then

$`({\displaystyle \frac{a_1}{a_0}},{\displaystyle \frac{a_2}{a_0}},\mathrm{},{\displaystyle \frac{a_n}{a_0}})\mathrm{\Omega }_{1a}.`$

Thus, $`\mathrm{\Omega }_{1a}`$ is not empty.

Moreover, $`(1,x_1,x_2,\mathrm{},x_n)\mathrm{\Omega }_a,`$ $`(x_1,x_2,\mathrm{},x_n)\mathrm{\Omega }_{1a}.`$ By Property 2, $`\mathrm{\Omega }_a`$ is open. Thus, these exists $`\delta >0,`$ such that, when $`\sqrt{(1y_0)^2+(x_1y_1)^2+\mathrm{}+(x_ny_n)^2}<\delta `$, we have $`(y_0,y_1,y_2,\mathrm{},y_n)\mathrm{\Omega }_a.`$ For this $`\delta ,`$ if $`(z_1,z_2,\mathrm{},z_n)R^n`$ satisfies

$`\sqrt{(x_1z_1)^2+\mathrm{}+(x_nz_n)^2}<\delta `$.

Then, we have

$$\sqrt{(11)^2+(x_1z_1)^2+\mathrm{}+(x_nz_n)^2}<\delta .$$

Hence $`(1,z_1,z_2,\mathrm{},z_n)\mathrm{\Omega }_a.`$ Thus, we have $`(z_1,z_2,\mathrm{},z_n)\mathrm{\Omega }_{1a}.`$ Namely, $`\mathrm{\Omega }_{1a}`$ is an open set.

By definition, $`\mathrm{\Omega }_{1a}`$ is convex.

In what follows, we will prove that $`\mathrm{\Omega }_{1a}`$ is unbounded. For this purpose, we first introduce some notations, which are needed in other proofs as well. Let

$$a(s)=s^n+a_1s^{n1}+\mathrm{}+a_nH^n$$
(4)
$$b(s)=x_0s^n+x_1s^{n1}+\mathrm{}+x_nP^nP^{n1}$$
(5)
Then $`\omega R,`$ we have

$`\begin{array}{cc}\text{Re}[{\displaystyle \frac{b(j\omega )}{a(j\omega )}}]\hfill & ={\displaystyle \frac{1}{|a(j\omega )|^2}}\text{Re}[b(j\omega )a(j\omega )]\hfill \\ & ={\displaystyle \frac{1}{|a(j\omega )|^2}}{\displaystyle \underset{l=0}{\overset{n}{}}}({\displaystyle \underset{k=0}{\overset{n}{}}}a_kx_{2lk}(1)^{l+k})\omega ^{2(nl)}\hfill \\ & ={\displaystyle \frac{1}{|a(j\omega )|^2}}{\displaystyle \underset{l=0}{\overset{n}{}}}c_l\omega ^{2(nl)}\hfill \end{array}`$

where $`c_l:=\underset{k=0}{\overset{n}{}}a_kx_{2lk}(1)^{l+k},`$ where $`a_0=1,`$ and let $`a_i=0,x_i=0,`$ when $`i<0`$ or $`i>n,`$ $`l=0,1,\mathrm{},n.`$

Introducing the matrices

$$H_a:=\left[\begin{array}{ccccccc}a_1& 1& 0& 0& 0& \mathrm{}& 0\\ a_3& a_2& a_1& 1& 0& \mathrm{}& 0\\ a_5& a_4& a_3& a_2& a_1& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ a_{2n1}& a_{2n2}& a_{2n3}& a_{2n4}& a_{2n5}& \mathrm{}& a_n\end{array}\right],$$

$$E_n:=\left[\begin{array}{ccccc}1& & & & \\ & 1& & & \\ & & 1& & \\ & & & 1& \\ & & & & \mathrm{}\end{array}\right],$$

$$A:=\left[\begin{array}{cccc}1& 0& \mathrm{}& 0\\ a_2& & & \\ a_4& & E_nH_aE_n& \\ \mathrm{}& & & \end{array}\right],$$

$$b:=\left[\begin{array}{c}x_0\\ x_1\\ \mathrm{}\\ x_{n1}\end{array}\right],$$

$$c:=\left[\begin{array}{c}c_0\\ c_1\\ \mathrm{}\\ c_n\end{array}\right]$$

where $`a_i=0`$ when $`i>n.`$ Then, it is easy to verify that

$$c=Ab$$
(6)
Since $`a(s)H^n,`$ we know that $`E_nH_aE_n`$ is invertible. For any $`d=[d_1,d_2,\mathrm{},d_n]^TR^n`$ such that all elements of $`d`$ are positive. Denote $`\overline{a}=[a_2,a_4,a_6,a_8,\mathrm{},(1)^na_{2n}]^T,`$ where $`a_i=0`$ when $`i>n`$. Let $`\overline{b}=(E_nH_aE_n)^1(d\overline{a}):=[b_1,b_2,\mathrm{},b_n]^T,`$ Then obviously, we have $`[1,d_1,d_2,\mathrm{},d_n]^T=A[1,b_1,b_2,\mathrm{},b_n]^T,`$ and $`b_n={\displaystyle \frac{d_n}{a_n}}.`$ By $`c=Ab`$, we know that $`(b_1,b_2,\mathrm{},b_n)\mathrm{\Omega }_{1a}.`$

On the other hand, due to the arbitrariness of $`d`$, $`d_n`$ can be taken arbitrarily large. Therefore, $`b_n`$ can also be arbitrarily large. Namely, $`\mathrm{\Omega }_{1a}`$ is unbounded. This completes the proof.

Remark 1 Given two stable polynomials $`a_1(s)`$ and $`a_2(s)`$ , existence of a polynomial $`b(s)`$ such that $`\frac{b(s)}{a_1(s)}`$ and $`\frac{b(s)}{a_2(s)}`$ are both SPR is tantamount to non-emptiness of the intersection of the two SPR regions associated with $`a_1(s)`$ and $`a_2(s)`$ . Since $`\mathrm{\Omega }_a`$ and $`\mathrm{\Omega }_{1a}`$ are both unbounded sets. When dealing with robust SPR design problem, we must find the intersection of several unbounded sets (i.e., SPR regions), which is intractable. This is the reason that we introduce the concept of WSPR regions, which are bounded as shown below.

Theorem 2 Given $`a(s)H^n,`$ $`\mathrm{\Omega }_{1a}^W`$ is a non-empty, bounded convex set in $`R^{n1}`$.

Proof

$`H_a,E_n`$ and $`A`$ were defined in the proof of Theorem 1.

Denote $`B`$ as the $`(n1)\times (n1)`$ matrix formed by the first $`n1`$ row and last $`n1`$ column of the matrix $`E_nH_aE_n`$. Obviously, $`B`$ is also invertible.

Denote $`\overline{a}:=[a_1,a_3,a_5,a_7,\mathrm{},(1)^{n1}a_{2(n1)+1}]^T`$ ( $`a_i=0`$ when $`i>n`$).

Let $`\overline{b}:=B^1\overline{a}=[b_1,b_2,\mathrm{},b_{n1}]^T.`$ Since $`a(s)H^n,`$ it is easy to verify that $`b_{n1}>0.`$ Denote $`b=[0,1,b_1,b_2,\mathrm{},b_{n1}]^T.`$ Let $`c_0=c_1=\mathrm{}=c_{n1}=0,`$ $`c_n=a_nb_{n1},`$ in $`c:=[c_0,c_1,\mathrm{},c_n]`$. Then, it is easy to verify that $`c=Ab`$ is true. Thus, we have $`(b_1,b_2,\mathrm{},b_{n1})\mathrm{\Omega }_{1a}^W,`$ Namely, $`\mathrm{\Omega }_{1a}^W`$ is not empty.

By definition, $`\mathrm{\Omega }_{1a}^W`$ is convex.

We now prove that $`\mathrm{\Omega }_{1a}^W`$ is bounded.

For any $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }_{1a}^W,`$ we have

$$\frac{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}{s^n+a_1s^{n1}+\mathrm{}+a_n}WSPR,$$

By Property 3, $`s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}H^{n1}.`$ Moreover, $`\omega R,`$ we have

$$\text{Re}(\frac{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega })>0.$$
(7)
Thus

$$\text{Re}(\frac{s^n+a_1s^{n1}+\mathrm{}+a_n}{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}|_{s=j\omega })>0,\omega R.$$
(8)
Obviously

$`{\displaystyle \frac{s^n+a_1s^{n1}+\mathrm{}+a_n}{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}}=s+`$
$`{\displaystyle \frac{\left(a_1x_1\right)s^{n1}+\left(a_2x_2\right)s^{n2}+\mathrm{}+\left(a_{n1}x_{n1}\right)s+a_n}{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}}`$
Therefore

$`Re\left({\displaystyle \frac{\left(a_1x_1\right)s^{n1}+\left(a_2x_2\right)s^{n2}+\mathrm{}+\left(a_{n1}x_{n1}\right)s+a_n}{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}}|_{s=j\omega }\right)`$
$`=Re\left({\displaystyle \frac{s^n+a_1s^{n1}+\mathrm{}+a_n}{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}}|_{s=j\omega }\right)Re\left(j\omega \right)`$
$`=\text{Re}\left({\displaystyle \frac{s^n+a_1s^{n1}+\mathrm{}+a_n}{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}}|_{s=j\omega }\right)>0,\omega R.`$
It is easy to see that

$$\frac{\left(a_1x_1\right)s^{n1}+\left(a_2x_2\right)s^{n2}+\mathrm{}+\left(a_{n1}x_{n1}\right)s+a_n}{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}$$

$$\left\{\text{SPR}\right\}\left\{\text{WSPR}\right\}.$$

Again, by Property 3, we have

$`(a_1x_1)s^{n1}+(a_2x_2)s^{n2}+\mathrm{}`$
$`+(a_{n1}x_{n1})s+a_nH^{n1}H^{n2}`$ (9)
Hence

$$0<x_1a_1,0<x_2<a_2,\mathrm{},0<x_{n1}<a_{n1}$$
(10)
Namely

$$\mathrm{\Omega }_{1a}^W\{(x_1,x_2,\mathrm{},x_{n1})|\alpha (s):=$$

$$\underset{i=1}{\overset{n}{}}(a_ix_i)s^{ni}H^{n1}H^{n2},\text{where}x_n=0\}$$

$$\{(x_1,x_2,\mathrm{},x_{n1})|0<x_1a_1,0<x_2<a_2,\mathrm{},$$

$$0<x_{n1}<a_{n1}\}.$$

Thus, $`\mathrm{\Omega }_{1a}^W`$ is bounded. This completes the proof.

Remark 2 It should be pointed out that $`\mathrm{\Omega }_{1a}^W`$ is not an open set in $`R^{n1}`$. In fact, from the proof of Theorem 2, we know that $`\mathrm{\Omega }_{1a}^W`$ is tangent to the hyperplane $`x_1=a_1`$ in $`R^{n1}`$. And there exist some points of $`\mathrm{\Omega }_{1a}^W`$ in this hyperplane. Thus, $`\mathrm{\Omega }_{1a}^W`$ can not be an open set. Obviously,$`\mathrm{\Omega }_a^W`$ is a non-empty, convex cone in $`R^{n1}`$. Thus, $`\mathrm{\Omega }_a^W`$ is also unbounded. One may be tempted to believe that $`\mathrm{\Omega }_a^W`$ is not an open set either.

Though $`\mathrm{\Omega }_{1a}^W`$ is not an open set. The following theorem guarantees such a fact: when the intersection of two or more WSPR regions is not empty, then the intersection must be a region, not a single point. This means that Ackermann’s counterexample (that the unstable region is an isolated point so that gridding the parameter space can lead to erroneous conclusions no matter how dense the gridding is) does not happen in this case.

Theorem 3 Given $`a(s)H^n,`$ if $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }_{1a}^W,`$ then for sufficiently small $`\epsilon >0,`$ we have $`(x_1\epsilon ,x_2\epsilon ,\mathrm{},x_{n1}\epsilon )\mathrm{\Omega }_{1a}^W.`$

Proof $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }_{1a}^W,`$ and $`\omega R,`$ we have

$$\text{Re}(\frac{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega })>0.$$
(11)
$`\epsilon >0,`$ since

$$Re(\frac{s^{n1}+(x_1\epsilon )s^{n2}+\mathrm{}+(x_{n1}\epsilon )}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega })$$

$$=Re(\frac{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega })$$

$$+Re(\frac{(\epsilon )(s^{n2}+s^{n3}+\mathrm{}+1)}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega })$$

$$=Re(\frac{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega })$$

$$+\frac{(\epsilon )}{|a(j\omega )|^2}(\omega ^{2(n1)}+\stackrel{~}{c}(\omega )),$$

where $`\stackrel{~}{c}(\omega )`$ is a real polynomial of order less or equal to $`2(n2)`$. Thus, when $`\omega `$ is sufficiently large, the sign of $`(\epsilon )(\omega ^{2(n1)}+\stackrel{~}{c}(\omega ))`$ will be positive. Namely, there exists $`\omega _1>0`$ such that, for all $`\omega \omega _1,`$

$$\text{Re}(\frac{s^{n1}+(x_1\epsilon )s^{n2}+\mathrm{}+(x_{n1}\epsilon )}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega })>0.$$
(12)
Denote

$$M_1=\underset{\omega \omega _1}{inf}\text{Re}(\frac{s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1}}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega }),$$
(13)
$$N_1=\underset{\omega \omega _1}{sup}\text{Re}(\frac{1}{|a(j\omega )|^2}(\omega ^{2(n1)}\stackrel{~}{c}(\omega ))$$
(14)
Then $`M_1>0`$ and $`N_1>0.`$ Choosing $`0<\epsilon <{\displaystyle \frac{M_1}{N_1}},`$ by simple computation, we have

$$\text{Re}\left(\frac{s^{n1}+\left(x_1\epsilon \right)s^{n2}+\mathrm{}+\left(x_{n1}\epsilon \right)}{s^n+a_1s^{n1}+\mathrm{}+a_n}|_{s=j\omega }\right)>0,\omega R$$
(15)
Therefore,

$$\frac{s^{n1}+(x_1\epsilon )s^{n2}+\mathrm{}+(x_{n1}\epsilon )}{s^n+a_1s^{n1}+\mathrm{}+a_n}WSPR,$$

namely $`(x_1\epsilon ,x_2\epsilon ,\mathrm{},x_{n1}\epsilon )\mathrm{\Omega }_{1a}^W.`$ This completes the proof.

The following theorem shows the relationship between $`\mathrm{\Omega }_{1a}^W`$ and $`\mathrm{\Omega }_a`$, and plays an important role in robust SPR design.

Theorem 4 Given $`a(s)H^n,`$ if $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }_{1a}^W,`$ then $`(1,\alpha _1,\alpha _2,\mathrm{},\alpha _n)R^{n+1},`$ we can take sufficiently small $`\epsilon >0`$ such that $`(0,1,x_1,x_2,\mathrm{},x_{n1})+\epsilon (1,\alpha _1,\alpha _2,\mathrm{},\alpha _n)\mathrm{\Omega }_a.`$

Proof Denote

$$b(s)=s^{n1}+x_1s^{n2}+\mathrm{}+x_{n1},$$
(16)
$$\alpha (s)=s^n+\alpha _1s^{n1}+\mathrm{}+\alpha _n,$$
(17)
$$\stackrel{~}{b}(s)=b(s)+\epsilon \alpha (s).$$
(18)
Since $`(x_1,x_2,\mathrm{},x_{n1})\mathrm{\Omega }_{1a}^W,`$ we have

$$\text{Re}(\frac{b(j\omega )}{a(j\omega )})>0,\omega R.$$
(19)
We only need to show that, for sufficiently small $`\epsilon >0,`$

$$\text{Re}(\frac{\stackrel{~}{b}(j\omega )}{a(j\omega )})>0,\omega R$$
(20)
Obviously, $`\stackrel{~}{b}(s)`$ and $`a(s)`$ have same order $`n`$. Thus, there exists $`\omega _2>0`$ such that, for all $`\omega \omega _2,`$ we have Re$`({\displaystyle \frac{\stackrel{~}{b}(j\omega )}{a(j\omega )}})>0.`$

Denote

$$M_2=\underset{\omega \omega _2}{inf}\text{Re}(\frac{b(j\omega )}{a(j\omega )}),N_2=\underset{\omega \omega _2}{sup}\text{Re}(\frac{\alpha (j\omega )}{a(j\omega )})$$
(21)
Then $`M_2>0`$ and $`N_2>0.`$ Choosing $`0<\epsilon <{\displaystyle \frac{M_2}{N_2}},`$ by simple computation, we have

$$\text{Re}(\frac{\stackrel{~}{b}(j\omega )}{a(j\omega )})>0,\omega R$$
(22)
This completes the proof.

## 4 Applications in Robust Design of SPR Transfer Functions

Generally speaking, the design problem is more difficult than analysis problem, since it is usually constructive, i.e., it not only shows the existence of the solution, but also provides a constructive procedure to find it. In this section, we will propose an algorithm for robust design of SPR transfer functions. This algorithm works well for both low order and high order polynomial families. Illustrative examples are provided to show the effectiveness of this algorithm.

Suppose

$$F=\{a_i(s)=s^n+\underset{l=1}{\overset{n}{}}a_l^{(i)}s^{nl},i=1,2,\mathrm{},m.\}$$
(23)
How can we find a polynomial $`b(s),`$ such that $`p_i(s):={\displaystyle \frac{b(s)}{a_i(s)}}`$SPR,$`i=1,2,\mathrm{},m`$?

As observed earlier, existence of such a polynomial $`b(s)`$ boils down to the condition that the intersection of the SPR regions associated with $`a_i(s)`$ is not empty. From the results in the previous section, we know that SPR regions are unbounded, whereas WSPR regions are bounded. Thus, by a computational consideration, we first consider the intersection of WSPR regions, and then construct a polynomial $`b(s)`$ by using the technique presented in the previous section.

Since SPR(WSPR) transfer functions with fixed numerator (or denominator) enjoy convexity property, namely, if there exists a polynomial $`c(s)`$, such that $`{\displaystyle \frac{c(s)}{a(s)}}`$ and $`{\displaystyle \frac{c(s)}{b(s)}}`$ are both SPR(WSPR),then, it is easy to verify that, for any $`\alpha 0,\beta 0`$ and $`(\alpha ,\beta )(0,0),`$ we have $`{\displaystyle \frac{c(s)}{\alpha a(s)+\beta b(s)}}`$SPR(WSPR). Therefore, the assumptions made on $`F`$ do not lose any generality. Actually, the method proposed in our paper also applies to convex combination of polynomials, interval polynomials, and more generally, polytopic polynomials and multilinearly perturbed polynomials.

By the results presented in the previous section, we propose the following design procedure:

Step 1. Test the robust stability of the convex hull of $`F`$, i.e., $`\overline{F}`$. If $`\overline{F}`$ is robustly stable, then go to Step 2; otherwise, print ”there does not exist such a $`b(s)`$ ”; (by Definitions 1 and 2)

Step 2. Let $`\alpha _l=\mathrm{min}\{a_l^{(i)},i=1,2,\mathrm{},m\},l=1,2,\mathrm{},n1.`$ Gridding the hyperrectangle

D:=

$$\{(x_1,x_2,\mathrm{},x_{n1})|0<x_l<\alpha _l,l=1,2,\mathrm{},n1\}$$
(24)
according to the precision required; (by Theorem 2 and its proof)

Step 3. Take $`b:=(b_1,b_2,\mathrm{},b_{n1})`$ at each gridding point. Test if $`b`$ belongs to $`_{i=1}^m\mathrm{\Omega }_{1a_i}^W`$ by the following steps:

1) Test if the $`(n1)`$-th order polynomial with coefficients $`(1,b_1,b_2,\mathrm{},b_{n1})`$ belongs to $`H^{n1}`$ (by Property 3),

2) For $`i=1,2,\mathrm{},m,`$ test if the polynomial with coefficients $`(a_1^{(i)}b_1,a_2^{(i)}b_2,\mathrm{},a_{n1}^{(i)}b_{n1},a_n^{(i)})`$ belongs to $`H^{n1}H^{n2}`$, respectively (by Theorem 2 and its proof),

3) Test if $`b`$ belongs to $`\mathrm{\Omega }_{1a_i}^W,i=1,2,\mathrm{},m`$ .

If all 1), 2), 3) above are satisfied, go to Step 4; otherwise, move to the next gridding point and test 1), 2), 3) again; (If all gridding points have been tested, then print ”there does not exist such a $`b`$ in $`_{i=1}^m\mathrm{\Omega }_{1a_i}^W`$ with the given precision”).

Step 4. Take a sufficiently small $`\epsilon >0`$ such that $`(\epsilon ,1,b_1,b_2,\mathrm{},b_{n1})_{i=1}^m\mathrm{\Omega }_{a_i}.`$ Hence, the $`n`$-th order polynomial with coefficients $`(\epsilon ,1,b_1,b_2,\mathrm{},b_{n1})`$ satisfies the design requirement. (by Theorem 4).

For the low order stable interval polynomial family or low order stable convex combination, existence of the solution to the design problem is always guaranteed. Given adequate precision, our method will surely find a polynomial that satisfies the design requirement. As shown by numerous examples below, our method is also effective for higher order polynomial families.

Example 1

Let

$$\begin{array}{cc}F=\{\hfill & a_1(s)=s^4+11s^3+56s^2+88s+1,\hfill \\ & a_2(s)=s^4+11s^3+56s^2+88s+50,\hfill \\ & a_3(s)=s^4+89s^3+56s^2+88s+1,\hfill \\ & a_4(s)=s^4+89s^3+56s^2+88s+50\}\hfill \end{array}$$
(25)
the methods proposed in \[8-10,14,15\] do not work here. Using our method, it is easy to get $`b(s)=s^3+3.3s^2+2.24s+1.76_{i=1}^4\mathrm{\Omega }_{1a_i}^W`$. Then let $`c(s):=\epsilon s^4+b(s)`$, where $`\epsilon >0`$ is sufficiently small, e.g., let $`\epsilon 0.3`$ (which is determined by Theorem 4),it is easy to check that the design requirement has been met.

Note that the example above is constructed by overbounding the line segment in  by an interval polynomial family. Thus, instead of dealing with two vertex polynomials as in , we must now deal with four Kharitonov’s vertex polynomials.

In what follows, we will give two more examples of higher order polynomial families.

Example 2

Let

$`F=\{`$
$`a_1\left(s\right)=s^7+9s^6+31s^5+71.5s^4+111.5s^3+109s^2+76s+12.5,`$
$`a_2\left(s\right)=s^7+9.5s^6+31s^5+71s^4+111.5s^3+109.5s^2+76s+12,`$
$`a_3\left(s\right)=s^7+9s^6+31.5s^5+71.5s^4+111s^3+109s^2+76.5s+12.5,`$
$`a_4\left(s\right)=s^7+9.5s^6+31.5s^5+71s^4+111s^3+109.5s^2+76.5s+12\}`$ (26)
It is easy to see that the convex hull $`\overline{F}`$ of $`F`$ is robust stable. Using our method, it is easy to get $`b(s)=s^6+7.2s^5+18.6s^4+42.6s^3+44.4s^2+43.6s+15.2_{i=1}^4\mathrm{\Omega }_{1a_i}^W.`$

Then let $`c(s):=\epsilon s^7+b(s)`$, where $`\epsilon >0`$ is sufficiently small, e.g., let $`\epsilon 0.1,`$ it is easy to check that the design requirement has been met.

Example 3

Let

$$\begin{array}{cc}F=\{\hfill & a_1\left(s\right)=s^9+11s^8+52s^7+145s^6+266s^5+331s^4\hfill \\ & +280s^3+155s^2+49s+6,\hfill \\ & a_2\left(s\right)=s^9+11s^8+52s^7+146s^6+265.5s^5+332s^4\hfill \\ & +278.5s^3+151s^2+48s+2\}\hfill \end{array}$$
(27)
It can be verified that $`a_2(s)a_1(s)=s^60.5s^5+s^41.5s^34s^2s4`$ satisfies the extended Alternating Hurwitz Minor Condition namely, it is a convex direction for Hurwitz stability. Moreover, it is easy to see that $`a_1(s)`$ and $`a_2(s)`$ are both Hurwitz stable polynomials. Thus, the convex hull $`\overline{F}`$ of $`F`$ is robust stable. Using our method, it is easy to get $`b(s)=s^8+8.8s^7+41.6s^6+87s^5+159.3s^4+132.4s^3+111.4s^2+30.2s+9.6\mathrm{\Omega }_{1a_1}^W\mathrm{\Omega }_{1a_2}^W.`$ Thus, let $`c(s):=\epsilon s^9+b(s),\epsilon >0,\epsilon `$ sufficiently small, e.g., take $`\epsilon 0.07,`$ then the design requirement has been met.

Note that our design method is also effective when dealing with discrete time systems. Note also that, in the Examples 1-3, $`b(s)`$ is not unique. Using our method, we can get all such $`b(s)`$’s with given precision.

It should be pointed out that there is hardly any example with order higher than 4 in the literature. Recently, a sixth-order example of interval family was given in  as follows. Unfortunately, this example is incorrect.

Example 4

Suppose

$$\begin{array}{cc}F=\{\hfill & \\ a_1\left(s\right)=s^6+0.8s^5+58.06s^4+50.9s^3+1028.5s^2+163.82s+1042.5,\hfill & \\ a_2\left(s\right)=s^6+1.5s^5+58.06s^4+28.3s^3+1028.5s^2+376.36s+1042.5,\hfill & \\ a_3\left(s\right)=s^6+0.8s^5+68.62s^4+50.9s^3+755.47s^2+163.82s+3286.7,\hfill & \\ a_4\left(s\right)=s^6+1.5s^5+68.62s^4+28.3s^3+755.47s^2+376.36s+3286.7\},\hfill & \end{array}$$
(28)
find a polynomial $`b(s),`$ such that

$$p_i(s):=\frac{b(s)}{a_i(s)}SPR,i=1,2,3,4.$$

By Definition 1 , Definition 2 and Property 3, a prerequisite of the robust SPR design problem is that the convex hull $`\overline{F}`$ of $`F`$ is robustly stable. But it is easy to check that $`\overline{F}`$ is not robustly stable. In fact, $`(1.0446\pm 5.8969i)`$ are roots of $`a_1(s)`$ with positive real part; $`(1.037\pm 4.9835i)`$ are roots of $`a_2(s)`$ with positive real part; $`(0.03291\pm 7.5026i)`$ and $`(0.68089\pm 2.4933i)`$ are roots of $`a_3(s)`$ with positive real part; $`(0.87123\pm 2.867i)`$ are roots of $`a_4(s)`$ with positive real part. Thus, it does not make sense to consider the robust SPR design in this case.

It should also be pointed out that, for the vertex set of a general polytopic polynomial family $`F=\{a_i(s)=s^n+_{l=1}^na_l^{(i)}s^{nl},i=1,2,\mathrm{},m.\}`$, even if $`\overline{F}`$ is robustly stable,it is still possible that there does not exist a polynomial $`c(s),`$ such that, $`c(s)/f(s)`$WSPR, for all $`f(s)\overline{F}`$. Namely $`_{i=1}^m\mathrm{\Omega }_{1a_i}^W=\varphi .`$

To see this, let us look at an example of a third order triangle polynomial family.

Proposition 1

Let $`a(s)=s^3+a_1s^2+a_2s+a_3H^3,`$ then

$`\begin{array}{cc}\mathrm{\Omega }_{1a}^W=\hfill & \{(x,y)|a_2^2x^2+2(2a_3a_1a_2)xy+a_1^2y^2\hfill \\ & 2a_2a_3x2a_1a_3y+a_3^2<0\}\hfill \\ & \{(x,y)|xa_1,a_2xa_1ya_30,y>0\}.\hfill \end{array}`$

Example 5

Let

$$\begin{array}{cc}F=\{\hfill & a_1(s)=s^3+2.6s^2+37s+64,\hfill \\ & a_2(s)=s^3+17s^2+83s+978,\hfill \\ & a_3(s)=s^3+15s^2+28s+415\}\hfill \end{array}$$
(29)
it is easy to verify that $`a_i(s),i=1,2,3,`$ are Hurwitz stable. Moreover, all edges of $`\overline{F}`$, i.e., $`\lambda a_i(s)+(1\lambda )a_j(s),\lambda [0,1],i,j=1,2,3,`$ are also Hurwitz stable. Therefore, by Edge Theorem, $`\overline{F}`$ is robustly stable. On the other hand, by a direct computation using Proposition 1, we have $`\mathrm{\Omega }_{1a_1}^W\mathrm{\Omega }_{1a_2}^W\mathrm{\Omega }_{1a_3}^W=\varphi .`$ Henceforth, there does not exist a polynomial $`c(s)`$ such that $`c(s)/a_i(s)`$ WSPR, $`i=1,2,3`$ (although $`\mathrm{\Omega }_{1a_1}^W\mathrm{\Omega }_{1a_2}^W\varphi ,\mathrm{\Omega }_{1a_1}^W\mathrm{\Omega }_{1a_3}^W\varphi ,`$ and $`\mathrm{\Omega }_{1a_2}^W\mathrm{\Omega }_{1a_3}^W\varphi `$).

Note that, in this example, though we have $`\mathrm{\Omega }_{1a_1}^W\mathrm{\Omega }_{1a_2}^W\mathrm{\Omega }_{1a_3}^W=\varphi ,`$ we do not know whether $`\mathrm{\Omega }_{a_1}\mathrm{\Omega }_{a_2}\mathrm{\Omega }_{a_3}`$ is an empty set or not. This is a problem deserving further study.

For a fourth ( or higher ) order stable interval polynomial family ( or stable convex combination of two polynomials), does there exist a polynomial such that their ratios are SPR-invariant? This is still an open problem. From our numerous examples, it seems that such a polynomial can always be found. Thus, we conjecture that this problem has a positive answer.

References

  L. Wang and L. Huang, Extreme Point Results for Strict Positive Realness of Transfer Function Families, Systems Science and Mathematical Sciences, 7: 430-439, 1994.

  Y. D. Landau, Adaptive Control: The Model Reference Approach. Marcel Dekker, New York, 1979.

  V. M. Popov, Hyperstability of Automatic Control Systems.  Springer-Verlag, New York, 1973.

  L. Wang and L. Huang, Robust Stability of Polynomial Families and Robust Strict Positive Realness of Rational Function Families, Int. J. Systems Science, 23: 235-247, 1992.

  L. Huang, Stability Theory, Peking University Press, Beijing, 1989.

  S. P. Bhattacharyya, H. Chapellat and L . H . Keel, Robust Control: The Parametric Approach. Prentice Hall , New York, 1995.

  B. R. Barmish, New Tools for Robustness of Linear Systems. MacMillan Publishing Company, New York, 1994.

  S. Dasgupta and A. S. Bhagwat, Conditions for designing strictly positive real transfer functions for adaptive output error identification. IEEE Trans. on Circuits and Systems, CAS-34: 731-737, 1987.

  L. Huang, C. V. Hollot, Robust analysis of strictly positive real function set. Preprints of The Second Japan-China Joint Symposium on Systems Control Theory and its Applications, 210-220, 1990.

  C. V. Hollot, L. Huang, Designing strictly positive real transfer function families: a necessary and sufficient condition for low degree and structured families. Proc. of the International Conference on Mathematical Theory of Network and Systems, 215-227, 1989.

  H. Chapellat and S. P. Bhattacharyya, On robust nonlinear stability of interval control systems. IEEE Trans. on Automatic Control, AC-36: 59-67, 1991.

  L. Wang and L. Huang, Finite Verification of Strict Positive Realness of Interval Transfer Functions, Science Bulletin, 4: 262-268, 1991.

  A. Betser and E. Zeheb, Design of robust strictly positive real transfer functions. IEEE Trans. Circuits and Systems, CAS-40: 573-580, 1993.

  B. D. O. Anderson, S. Dasgupta, P. Khargonekar, F. J. Kraus and M. Mansour, Robust strict positive realness: characterization and construction. IEEE Trans. Circuits and Systems, CAS-37: 869-876, 1990.

  W. Yu and L. Huang, Necessary and Sufficient Conditions for Strictly Positive Real Stabilization of Low-Order Systems, Science Bulletin, 43: 2275-2279, 1998.

  W. Yu, Robust SPR Synthesis and Robust Stability Analysis. PhD Thesis, Peking University, Beijing, 1998.

  W. S. Yu and L. Wang, Some Remarks on the Definition of Strict Positive Realness of Transfer Functions. Proceedings of the Conference on Control and Decision, Nanjing, 1999.

  W. S. Yu, L. Wang and M. Tan, Complete characterization of strict positive realness regions in Coefficient space. Proceedings of IEEE Hong Kong Symposium on Robotics and Control, Hong Kong, 1999.

  J. Ackermann, Robust Control: Systems with Uncertain Physical Parameters. Springer-Verlag, London, 1993.

  A. C. Bartlett, C. V. Hollot and L. Huang, Root locations for an entire polytope of polynomial: it suffices to check the edges. Mathematics of Control, Signals, and Systems, 1: 61–71, 1988.

  Barmish, B. R. and Kang, I. H., Extreme point results for robust stability of interval plants: Beyond first order compensators. Automatica,  1992, 28: 1169-1180.

  Rantzer, A., Stability conditions for polytopes of polynomials. IEEE Trans. on Automatic Control,  1992, AC-37: 79-89.

  Marquez, H. J. and Agathoklis, P., On the existence of robust strictly positive real rational functions. IEEE Trans. on Circuits and Systems, Part I, 1998, CAS-45: 962-967.


text:
# 1 Factorization for the decay 𝐵→𝑀𝜋.

UdeM-GPP-TH-02-94

UTPT-02-01, TAUP 2696-02, WIS/07/02-JAN-DPP

February 2002

Charmless B Decays to Final States with Radially Excited Vector Mesons
Alakabha Datta111email: datta@lps.umontreal.ca a, Harry J. Lipkin 222email: harry.lipkin@weizmann.ac.il b and Patrick J. O’Donnell 333email: pat@medb.physics.utoronto.cac

a Laboratoire René J.-A. Lévesque, Université de Montréal,

C.P. 6128, succ. centre-ville, Montréal, QC, Canada H3C 3J7

b) Department of Particle Physics,
Weizmann Institute,
Rehovot 76100, Israel
and
School of Physics and Astronomy,
Tel-Aviv University,
Tel-Aviv 69978, Israel
c Department of Physics and Astronomy,
University of Toronto, Toronto, Canada.

## Abstract

We consider the weak decays of a B meson to final states that contain a S-wave radially excited vector meson. We consider vector-pseudoscalar final states and calculate ratios of the type $`B\rho ^{}\pi /B\rho \pi `$, $`B\omega ^{}\pi /B\omega \pi `$ and $`B\varphi ^{}\pi /B\varphi \pi `$ where $`\rho ^{}`$, $`\omega ^{}`$ and $`\varphi ^{}`$ are higher $`\rho `$, $`\omega `$ and $`\varphi `$ S-wave radial excitations. We find such decays to have larger or similar branching ratios compared to decays where the final state $`\rho `$, $`\omega `$ and $`\varphi `$ are in the ground state. We also study the effect of radial mixing in the vector system generated from hyperfine interaction and the annihilation term.

The new data accumulating from B factories and other accelerators will include transitions to many new final states which have not been previously studied in detail; e.g. radially excited meson states. Many decays involve a transition from a low momentum spectator quark to a high momentum relativistic meson. The form factors for such transitions are expected to be sensitive to the high momentum components of the final meson wave function, and therefore to favor radially excited states. The data on B decays to these states will thus provide important new information, particularly for the form factors to the radially excited states and probe the high-momentum tails of their wave functions.

In this paper we calculate predictions for the ratios

$`R_{\rho ^+}`$ $`=`$ $`BR(\overline{B}^0\pi ^{}\rho ^+)/BR(\overline{B}^0\pi ^{}\rho ^+)`$ (1)
$`R_{\rho ^0}`$ $`=`$ $`BR(\overline{B}^{}\pi ^{}\rho ^0)/BR(\overline{B}^{}\pi ^{}\rho ^0)`$ (2)
$`R_\omega `$ $`=`$ $`BR(B^{}\pi ^{}\omega ^{})/BR(B^{}\pi ^{}\omega )`$ (3)
$`R_\varphi `$ $`=`$ $`BR(B_s\pi ^0\varphi ^{})/BR(B_s\pi ^0\varphi )`$ (4)
where $`\rho ^{},\omega ^{}`$ and $`\varphi ^{}`$ are the radially excited states. Most studies of two body nonleptonic $`B`$ decays have concentrated on processes of the type $`BM_1M_2`$ where both $`M_1`$ and $`M_2`$ are mesons in the ground state configuration. Nonleptonic decays, where one of the mesons in the final state containing the spectator quark is a radially excited state, are expected to have larger or similar branching ratios compared to decays where the final state contains the same meson in the ground state. This is easily seen in a simple model in which $`B\pi M`$ and M is a simple flavor eigenstate with no flavor mixing beyond isospin. We follow the inactive spectator approach used  to treat B decays to charmonium in which the spectator quark does not participate in a flavor-changing interaction and later combines with a light antiquark to make the final light meson as shown in Fig. 1. The decay amplitudes are then described as the product of a b-quark decay amplitude and a hadronization function describing the combination of a quark-antiquark pair to make the final meson. Neglecting the relative Fermi momentum of the b quark and the spectator quark in the $`B`$ meson, the quark transition for the processes in Eq. (1-3) is

$$b\pi ^{}(\stackrel{}{p})u(\stackrel{}{p})$$
(5)
where the b quark is at rest and $`\stackrel{}{p}`$ denotes the final momentum of the $`\pi ^{}`$. For the process in Eq.4, the quark transition is essentially similar to the one above

$$b\pi ^0(\stackrel{}{p})s(\stackrel{}{p})$$
(6)
where now $`\stackrel{}{p}`$ denotes the final momentum of the $`\pi ^0`$.

Concentrating on the processes in Eq. (1-3) the transition matrix for the full decay has the form

$$\pi ^{}(\stackrel{}{p})M(\stackrel{}{p})\left|T\right|B=M(\stackrel{}{p})\left|F\right|u(\stackrel{}{p})\overline{q}(0)\pi ^{}(\stackrel{}{p})u(\stackrel{}{p}))\left|W\right|b$$
(7)
where $`T`$ denotes the transition matrix for the hadronic decay which factors, as shown in Fig. 1, into a weak matrix element at the quark level denoted by $`W`$ and a recombination matrix element denoted by $`F`$. This latter matrix element describes the transition of a quark with momentum $`\stackrel{}{p}`$ and an antiquark with zero momentum to make a meson with momentum $`\stackrel{}{p}`$. Radial excitations could be favoured over ground states, if the final momentum $`\stackrel{}{p}`$ is large, since the radial excitations are expected to have higher kinetic energies.

An alternative but equivalent way of understanding this effect is to note that for a heavy $`b`$ quark one can write

$$\pi ^{}(\stackrel{}{p})M(\stackrel{}{p})\left|T\right|B=M(\stackrel{}{p})\left|J_{1\mu }\right|B\pi ^{}(\stackrel{}{p})\left|J_2^\mu \right|0$$
(8)
where $`J_{1,2}`$ are currents that occur in $`W=J_1\times J_2`$. The transition matrix element for the hadronic decay can then be written in terms of $`BM`$ form factors and the pion decay constant. The form factors can be expressed as overlap integrals of the $`B`$ and the $`M`$ meson wavefunctions. When $`M`$ is a light meson, with a mass much smaller than the $`B`$ meson, the overlap integrals get contributions mainly from the high momentum components of the meson wavefunctions. It is now clear that for a radially excited meson $`M^{}`$, which has more higher momentum components, the overlap integrals will be enhanced compared to the overlap integral with a ground state meson $`M`$. Consequently the $`BM^{}`$ form factors are likely to be enhanced relative to the $`BM`$ form factors which would then translate into higher branching ratio for $`BM^{}\pi `$ relative to $`BM\pi `$.

Our discussion above assumed the physical states to be pure radial excitations. However, additional interactions can mix the various radial excited components. For instance hyperfine interactions can mix radial excitations with the same flavor structure and so in general, in the $`\rho `$, $`\omega `$ and $`\varphi `$ systems, the various physical states will be admixtures of radial excitations . Flavor mixing in the vector system is known to be small but is important in the pseudoscalar sector. We will consider the pseudoscalar case in a different publication. To make quantitative predictions we use constituent quark wave functions with several potentials to test the dependence of the results on the confining potential. We shall see that the effects of the potential dependence and mixing are small so that the results are reasonably robust and are not seriously dependent on the fine details of the model.

Even though our discussion has so far only included vector-pseudoscalar final states we can also consider vector-vector final states such as $`B\rho \rho `$ or $`BJ/\psi K^{}`$. However vector- vector final states are complicated since different partial waves are present. Our purpose here is to demonstrate the effects of radial mixing in the simpler physical system of the vector-pseudoscalar final state. If the effects of radial enhancements are observed in the vector-pseudoscalar case we would expect them to be also present in the vector-vector final state.

We will first review the study the masses and mixing in the vector meson sector following the simple nonrelativistic approach in Ref. To obtain the eigenstates and eigenvalues in the vector meson system we diagonalize the mass matrix

$`<q_a^{}\overline{q}_b^{},n^{}|M|q_a\overline{q}_b,n>`$ $`=`$ $`\delta _{aa^{}}\delta _{bb^{}}\delta _{nn^{}}(m_a+m_b+E_n)+\delta _{aa^{}}\delta _{bb^{}}{\displaystyle \frac{B}{m_am_b}}\stackrel{}{s}_a\stackrel{}{s}_b\psi _n(0)\psi _n^{}(0)`$ (9)
where $`\stackrel{}{s}_{a,b}`$ and $`m_{a,b}`$ are the quark spin operators and masses. Here $`n=0,1,2`$ and the basis states for the isovector mesons are chosen as $`|N,I=1,I_3=1=|u\overline{d}`$, $`|N,I=1,I_3=0=|u\overline{u}d\overline{d}/\sqrt{2}`$ and $`|N,I=1,I_3=1=|d\overline{u}`$ where $`I,I_3`$ stand for the isospin quantum numbers. In the above equation $`E_n`$ is the excitation energy of the $`n^{th}`$ radially excited state and $`B`$ is the strength of the hyperfine interaction.

To begin with, we use the same harmonic confining potential as well as the other parameters used in Ref to obtain the eigenstates and eigenvalues for the mass matrix in Eqn. 9. The various parameters used in the calculation are the constituent masses, $`m_u=m_d=`$ 0.350 GeV, $`m_s=`$ 0.503 GeV, the angular frequency, $`\omega =`$0.365 GeV and $`b=B/m_u^2=`$0.09. The eigenvalues and eigenstates for the $`\rho `$ system with a harmonic potential are shown in in Table. 1.

To see how this result changes with a different confining potential we use a power law potential $`V(r)=\lambda r^n`$  . We will use a linear and a quartic confining potential and compare the spectrum with that obtained with a harmonic oscillator potential. To fix the coefficient $`\lambda `$ we require that the energy eigenvalues of the Schrödinger equation are similar in the least square sense with the energy eigenvalues used in Ref. So for example, for the linear potential, we demand that

$$F=\underset{n}{}(E_n(harmonic)E_n(linear))^2$$

is a minimum. This fixes the constant $`\lambda `$ in $`V(r)=\lambda r`$ and we obtain the eigenvalues and eigenstates in Table. 2.

We follow the same procedure for the quartic potential and obtain the eigenvalues and eigenstates in Table. 3

We observe from Tables. (1-3) that the mass eigenstates and eigenvalues of the $`\rho `$ system are not very sensitive to the confining potential and the radial mixing effects are small.

To obtain the eigenstates and eigenvalues in the $`\omega \varphi `$ system we diagonalize the mass matrix

$`<q_a^{}\overline{q}_b^{},n^{}|M|q_a\overline{q}_b,n>`$ $`=`$ $`\delta _{aa^{}}\delta _{bb^{}}\delta _{nn^{}}(m_a+m_b+E_n)+\delta _{aa^{}}\delta _{bb^{}}{\displaystyle \frac{B}{m_am_b}}\stackrel{}{s}_a\stackrel{}{s}_b\psi _n(0)\psi _n^{}(0)`$ (10)
$`+`$ $`\delta _{ab}\delta _{a^{}b^{}}{\displaystyle \frac{A}{m_am_b}}\psi _n(0)\psi _n^{}(0)`$
This has a similar structure as the $`\rho `$ system but now we have the additional annihilation interaction with strength $`A`$ that causes flavor mixing. .

Diagonalizing the mass matrix in Eqn. 10, with the basis states $`|N=|u\overline{u}+d\overline{d}/\sqrt{2}`$ and $`|S=|s\overline{s}`$, we obtain the eigenvalues and the eigenstates of the $`\omega \varphi `$ system. We use the same value for the hyperfine interaction as used for the $`\rho `$ system. For the linear potential we obtain with $`B=0.09m_u^2`$ and $`A=0.005m_u^2`$ the eigenstates and eigenvalues in Table. 4.

For the harmonic potential we obtain with $`B=0.09m_u^2`$ and $`A=0.015m_u^2`$ the eigenstates and eigenvalues in Table. 5.

For the quartic potential we obtain with $`B=0.09m_u^2`$ and $`A=0.023m_u^2`$ the eigenstates and eigenvalues in Table. 6.

As in the $`\rho `$ system we find the mixing to be insensitive to the confining potential and the effects of radial mixing to be small. We also find, as expected, a small value for the annihilation term in the fits to the masses.

We now use these wavefunctions to predict the ratios in Eq. (1-4). These decays are dominated by diagrams which satisfy the inactive spectator approach  and are treated with Eq.(7-8). Some of the diagrams which violate this assumption; e.g. penguin and annihilation contributions, may not be as negligible here as in the charmonium case treated in Ref for the decays to the ground state configurations. But they are expected to have much smaller form factors for radial excitations. Therefore it is reasonable to neglect them for this preliminary investigation of the order of magnitude of these ratios. Note that for the $`B_s\pi ^0\varphi `$ decay the QCD penguin is isospin forbidden, the annihilation contribution is OZI forbidden and the electroweak penguin is also described by Eq. 7.

We obtain, using factorization for the nonleptonic amplitude,

$`R_{\rho ^+}`$ $`=`$ $`\left|{\displaystyle \frac{\rho ^+\left|\overline{u}\gamma ^\mu (1\gamma _5)b\right|\overline{B}^0\pi ^{}\left|\overline{d}\gamma _\mu (1\gamma _5)u\right|0}{\rho ^+\left|\overline{u}\gamma ^\mu (1\gamma _5)b\right|\overline{B}^0\pi ^{}\left|\overline{d}\gamma _\mu (1\gamma _5)u\right|0}}\right|^2`$ (11)
$`=`$ $`\left|{\displaystyle \frac{A_0^{\rho ^+}}{A_0^{\rho ^+}}}\right|^2{\displaystyle \frac{P_{\rho ^+}^3}{P_\rho ^3}}`$
where $`P`$ is the magnitude of the three momentum of the final states and the form factor $`A^0`$ is defined through

$`V_f\left|A_\mu \right|P_i`$ $`=`$ $`(M_i+M_f)A_1\left[ϵ_\mu ^{}{\displaystyle \frac{ϵ^{}.q}{q^2}}q_\mu \right]A_2{\displaystyle \frac{ϵ^{}.q}{M_i+M_f}}\left[(P_i+P_f)_\mu {\displaystyle \frac{M_i^2M_f^2}{q^2}}q_\mu \right]`$ (12)
$`+`$ $`2M_fA_0{\displaystyle \frac{ϵ^{}.q}{q^2}}q_\mu `$
where $`A_\mu `$ is the axial vector current. Similarly we obtain

$`R_{\rho ^0}`$ $``$ $`\left|{\displaystyle \frac{\rho ^0\left|\overline{u}\gamma ^\mu (1\gamma _5)b\right|\overline{B}^0\pi ^{}\left|\overline{d}\gamma _\mu (1\gamma _5)u\right|0}{\rho ^0\left|\overline{u}\gamma ^\mu (1\gamma _5)b\right|\overline{B}^0\pi ^{}\left|\overline{d}\gamma _\mu (1\gamma _5)u\right|0}}\right|^2`$ (13)
$`=`$ $`\left|{\displaystyle \frac{A_0^{\rho ^0}}{A_0^{\rho ^0}}}\right|^2{\displaystyle \frac{P_{\rho ^0}^3}{P_{\rho ^0}^3}}`$
$`R_\omega `$ $``$ $`\left|{\displaystyle \frac{\omega ^{}\left|\overline{u}\gamma ^\mu (1\gamma _5)b\right|\overline{B}^0\pi ^{}\left|\overline{d}\gamma _\mu (1\gamma _5)u\right|0}{\omega \left|\overline{u}\gamma ^\mu (1\gamma _5)b\right|\overline{B}^0\pi ^{}\left|\overline{d}\gamma _\mu (1\gamma _5)u\right|0}}\right|^2`$ (14)
$`=`$ $`\left|{\displaystyle \frac{A_0^\omega }{A_0^\omega }}\right|^2{\displaystyle \frac{P_\omega ^{}^3}{P_\omega ^3}}`$
Finally,

$`R_\varphi `$ $`=`$ $`\left|{\displaystyle \frac{\varphi \left|\overline{s}\gamma ^\mu (1\gamma _5)b\right|\overline{B}_s\pi ^0\left|\overline{u}\gamma _\mu (1\gamma _5)u\right|0}{\varphi \left|\overline{s}\gamma ^\mu (1\gamma _5)b\right|\overline{B}^0\pi ^0\left|\overline{u}\gamma _\mu (1\gamma _5)u\right|0}}\right|^2`$ (15)
$`=`$ $`\left|{\displaystyle \frac{A_0^\varphi }{A_0^\varphi }}\right|^2{\displaystyle \frac{P_\varphi ^3}{P_\varphi ^3}}`$
To calculate the ratios we need the form factor $`A_0`$. Note that the wavefunctions for the various vector meson states are not enough to calculate non leptonic decay amplitudes. In particular, with the factorization assumption for non leptonic decays, the calculations of decay amplitudes require the matrix elements of current operators between the physical states. These matrix elements can be expressed in terms of form factors and decay constants. In this work we use a constituent quark model(CQM) model for the form factors  that incorporates some relativistic features and is relatively simple to adapt to the case of radially excited states. In this model the form factor $`A_0`$ is given by

$`A_0`$ $`=`$ $`Z[I_1{\displaystyle \frac{M_{}}{M_+}}I_2]`$
where

$`Z`$ $`=`$ $`{\displaystyle \frac{\sqrt{4M_iM_f}M_+}{M_+^2q^2}}`$
$`I_1`$ $`=`$ $`{\displaystyle d^3p\varphi _f^{}(\stackrel{}{p}+\stackrel{}{a})\varphi _i(\stackrel{}{p})}`$
$`I_2`$ $`=`$ $`m_s{\displaystyle d^3p\varphi _f^{}(\stackrel{}{p}+\stackrel{}{a})\varphi _i(\stackrel{}{p})[\frac{\stackrel{}{p}.\stackrel{}{a}}{\mu a^2}+\frac{1}{m_f}]}`$
$`M_\pm `$ $`=`$ $`M_f\pm M_i`$
$`\stackrel{}{a}`$ $`=`$ $`2m_s\stackrel{}{\beta }=2{\displaystyle \frac{m_s\stackrel{~}{q}}{M_+}}`$
$`\stackrel{~}{q}^2`$ $`=`$ $`M_+^2{\displaystyle \frac{M_{}^2q^2}{M_+^2q^2}}`$
$`\mu `$ $`=`$ $`{\displaystyle \frac{m_im_f}{m_i+m_f}}`$ (16)
and $`\varphi _f`$ and $`\varphi _i`$ represent the momentum space wave functions while $`\stackrel{}{\beta }`$ is the velocity of the mesons in the equal velocity frame ( also called the Breit frame or the brick wall frame) and $`m_{i,f}`$ are the non spectator quark masses of the initial and final meson. The equal velocity frame in a convenient frame to calculate the Lorentz invariant form factors where the velocities, $`\stackrel{}{\beta _i}`$ and $`\stackrel{}{\beta }_f`$ of the mesons with masses $`M_i`$ and $`M_f`$ are equal in magnitude but opposite in direction. We use the momentum wavefunction $`\varphi _f`$ obtained from spectroscopy in section. (2) while for $`\varphi _i`$ we use the wave function

$`\varphi _i`$ $`=`$ $`\varphi _B=N_Be^{p^2/p_F^2}`$ (17)
where $`p_F`$ is the Fermi momentum of the B meson. In our calculations we will take $`p_F=300`$MeV. Note that in the analysis presented in the introduction we have neglected the Fermi momentum of the b quark, since $`p_F/m_b`$ is small.

For transitions to higher resonant states, we use the same quark masses as those used in the transition of the $`B`$ meson to the lowest resonant state. This is reasonable, as the spectator quark still comes from the $`B`$ meson and therefore has the same value for its mass irrespective of whether the final state is in the lowest or the first excited state. The values for the masses of the non spectator masses are taken to be the same as those used for spectroscopy. However for the calculation of the velocity $`\stackrel{}{\beta }`$ and hence $`\stackrel{}{a}`$ defined in Eqn. 16 we use the physical mass of the higher resonant state.

In Table. 7 we give our predictions for the various ratios defined above. We find that the transitions to higher excited states can be comparable or enhanced relative to the transitions to the ground state. From Table. 7 we see that the ratios of branching ratios are slightly sensitive to the confining potential and the ratios of branching ratios increase as we go from the quartic to the linear potential. This is because the wavefunction for the linear potential has a longer tail and hence more high momentum components than the wavefunction for the quadratic and the quartic potentials. The wavefunction for the quadratic potential, has in turn, a longer tail and hence more high momentum components than the wavefunction for the quartic potential. Hence we would expect the hierarchy $`(A_0)_{linear}>(A_0)_{quadratic}>(A_0)_{quartic}`$ and a similar one for the radially excited states $`(A_0^{})_{linear}>(A_0^{})_{quadratic}>(A_0^{})_{quartic}`$ where $`A_0`$ and $`A_0^{}`$ are the form factors for the transition of $`B`$ to the ground state and the first radially excited state of the meson $`M`$. We see from from Table. 7 that this hierarchy is maintained for the ratios of form factors and so we have $`(A_0^{}/A_0)_{linear}>(A_0^{}/A_0)_{quadratic}>(A_0^{}/A_0)_{quartic}`$. Note that the ratio of form factors also depend on the choice of the Fermi momentum of the $`B`$ meson, as a smaller(larger) Fermi momentum would make the form factors more(less) sensitive to the tail of the wavefunction of $`M`$, as well as mixing effects in the wavefunction of the meson $`M`$. The effect of mixing between the various radially excited states and the ground state is generally small.

We observe in Table. 7 that there can be a large enhancement for $`R_\varphi `$. This decay is suppressed in the standard model. One can get a rough estimate of the branching ratio for $`B_s\varphi \pi ^0`$ using factorization as

$`{\displaystyle \frac{BR[B_s\varphi \pi ^0]}{BR[\overline{B}^0\rho ^+\pi ^{}]}}`$ $``$ $`{\displaystyle \frac{1}{2}}\left|{\displaystyle \frac{V_{ub}V_{us}^{}(c_1+c_2/N_c)V_{tb}V_{ts}^{}3c_9/2}{V_{ub}V_{ud}^{}(c_2+c_1/N_c)}}\right|0.02`$
where we have neglected form factor and phase space differences between $`B_s\varphi \pi ^0`$ and $`\overline{B}^0\rho ^+\pi ^{}`$. The Wilson coefficients $`c_i`$ can be found in Ref while $`V_{ub},V_{us},V_{ud},V_{tb}`$ and $`V_{ts}`$ are the various CKM elements . Using the measured $`BR[B\rho ^+\pi ^{}]28\times 10^6`$  we get $`BR[B_s\varphi \pi ^0]6\times 10^7`$. Hence the large enhancement for $`R_\varphi `$ indicates that $`BR[B_s\varphi ^{}\pi ^0]`$ can be around $`4\times 10^6`$.

Note that in the $`\rho (\omega )`$ system there are two resonances, $`\rho (1450)[\omega (1420)]`$ and $`\rho (1700)[\omega [1650]`$, which can be identified a S-wave radial excitation(2S) and a D wave orbital excitation in the quark model. However recent studies of the decays of these resonances show that it is possible that these states are mixtures of $`q\overline{q}`$ and hybrid states Ref. Hence the state $`\rho (1450)[\omega (1420)]`$ is interpreted as a 2S state with a small mixture of a hybrid state. We do not take into account such possible mixing with a hybrid state in our calculation and the meson masses for these excited states used in our calculation are the ones we predict in section 2. For the $`\varphi `$ system there is only state at $`\varphi (1680)`$ which we interpret as a 2S state in the absence of mixing effects.

In conclusion we have considered the weak decays of a B meson to final states that are mixtures of S-wave radially excited components. We calculated nonleptonic decays of the type $`B\rho ^{}\pi /B\rho \pi `$, $`B\omega ^{}\pi /B\omega \pi `$ and $`B\varphi ^{}\pi /B\varphi \pi `$ where $`\rho ^{}`$, $`\omega ^{}`$ and $`\varphi ^{}`$ are higher $`\rho `$, $`\omega `$ and $`\varphi `$ resonances. We found that the transitions to the excited states can be comparable or enhanced relative to transitions to the ground state. It would, therefore, be extremely interesting to test these predictions. We also studied the effect of radial mixing in the vector system generated from hyperfine interaction and the annihilation term; these turn out to be generally small.

Acknowledgment

This work was supported by the US-Israel Bi-National Science Foundation and by the Natural Sciences and Engineering Research Council of Canada.