text:
# Candidate Pre-Mainsequence F Stars with Circumstellar Dust Identified Using Combined 2MASS and uvby Data

## 1 Introduction

The detection of circumstellar dust disks provides us with an unprecedented opportunity to study the formation and evolution of exo-solar systems (see, e.g., Beckwith 1999) as well as the PMS stellar evolution (e.g., van den Ancker et al. 1997), given that many of these systems are associated with PMS stars. The bulk of the known PMS stars are comprised by Herbig Ae and T Tauri stars. The presence of large amounts of dust about many of these stars has been inferred from IRAS observations (Aumann et al. 1984; Aumann 1985; Gillett 1986; Backman & Gillett 1987; Walker & Wolstencroft 1988; Backman & Paresce 1993). The IRAS survey still remains a major source to search for candidate dust-disk systems (e.g., Mannings & Barlow 1998). Some 400 IRAS candidates exhibit an unusually large far-infrared (FIR) excess (e.g., Stencel & Backman 1991, Mannings & Barlow 1998), interpreted as being due to thermal emission of hot circumstellar dust ($`30<T<1000`$ K). In contrast to IRAS systems whose large FIR excess results from a high-density concentration of warm dust, there are young stellar objects (YSO) that show relatively small FIR excess. They include T Tauri stars exhibiting spectral energy distribution (SED) that declines longward of $`2`$ µm (YSOs with SED of Class II and Class III); these stars are interpreted as advanced evolutionary phases of T Tauri systems, with only a fraction of the original dust mass remaining in the disk (e.g., Lada and Wilking 1984; André & Montmerle 1994). Less massive dust disks or envelopes may not be detectable in the far infrared. Disk clearing may have started for these stars, but it may still be possible to detect the disk remnant material in the near infrared through reflected light.

In this paper we present a new criterion that uses a combined set of 2MASS and uvby data to identify F stars with circumstellar dust. The candidate stars found with this method appear to possess dusty disks or envelopes radiating the reflected and/or thermal near-infrared (NIR), causing a star to look in the NIR brighter than expected for a reddened normal F star at the same visual magnitude $`V`$. Many of them are not IRAS sources and may therefore represent a dusty environment substantially different from that found among the IRAS sources.

## 2 Data

Starting with a sample described in Suchkov & McMaster (1999), which includes F stars with uvby colors from Hauck & Mermilliod (1998) and observed by the Hipparcos, we found $`5230`$ counterparts in the 2MASS Point Source Catalog. Some of the stars among those with high reddening are, in fact, of the spectral type A rather than F. The reason is that the selection for the original sample was made based on the $`(by)`$ color, specifically, $`0.222<by<0.390`$. For unreddened stars, the latter criterion effectively selects stars later than spectral type A. However, if reddening is present, stars with intrinsically bluer colors, including A stars, can also be picked up by this criterion. The overall fraction of A stars in our sample is very small, and for the sake of simplicity, we’ll be referring to the sample stars as F stars. The criterion also selects a small fraction of early type G stars.

The uvby data were used to derive reddening $`E_{by}`$, effective temperature, and metallicity as described in Suchkov & McMaster (1999) and Suchkov (2001), based on the algorithm in Moon (1985) as well as the Schuster & Nissen (1989) metallicity calibration. Extinction was estimated as

$$A_V=\frac{E_{by}}{0.74}(3.354+\frac{0.0406E_{by}}{0.74}).$$
(1)
This equation follows from the results in Straižis (1977), Crawford (1975), and Sudžius (1974). Absolute magnitude was computed from the $`V`$ magnitude corrected for extinction, using the Hipparcos parallax; $`V`$ magnitudes were taken from the Hipparcos catalog.

## 3  Circumstellar dust from intrinsic near-infrared excess

The sample stars exhibit a tight temperature-induced correlation between the temperature-sensitive uvby color index $`(by)_0`$ and $`VK`$. In Figure 1, this correlation is shown for stars within 80 pc with small or no reddening, $`E_{by}<0.015`$. The scatter about the regression line, $`0.05`$ after 3-sigma correction, is due to both photometric errors and real physical effects, mainly related to absolute magnitude, duplicity, and remaining reddening. It reduces to $`\sigma _{VK}0.035`$ if the range for these parameters is narrowed down to $`2<M_V<3`$, $`E_{by}<0.005`$, and binaries are excluded. This should be close to the value resulting from photometric errors only. The linear regression shown in Figure 1 remains, however, the same, within uncertainties, as for the constrained sample. It was used to compute the deviation from the regression line for all sample stars,

$$\delta (VK)_{by}=(VK)(VK)_{by}.$$
(2)
A star having an excessive NIR flux with respect to the one expected for an unreddened standard star at the same effective temperature would deviate upward from this line,

$$\delta (VK)_{by}>0.$$
(3)
Circumstellar dust can contribute to observed infrared excess in two quite different ways. First, this material would absorb stellar light stronger at shorter wavelengths, which adds up with the effect of interstellar extinction and results in flux deficiency in the blue and visual with respect to the flux in the infrared (the reddening effect). The part of the infrared excess due to this effect will be referred to as apparent infrared excess. Second, a dust disk is a source of reflected and thermal radiation. Unlike circumstellar and interstellar extinction, this radiation actually adds to the obscured stellar flux in the infrared. The part of the near-infrared excess due to this extra flux will be referred to as intrinsic near-infrared excess. Isolating the latter part in the observed total infrared excess would thus amount to the detection of circumstellar dust.

The two parts of the infrared excess can be separated as follows. The reddening effect in the total near-infrared excess as defined by equations (2) and (3) is clearly seen in the $`E_{by}\delta (VK)_{by}`$ diagram in Figure 2. In the diagram, the bulk of the reddened stars tightly envelope the line representing the relationship between $`\delta (VK)_{by}`$ and $`E_{by}`$ expected from reddening only,

$$\delta (VK)_{by}=A_V(E_{by}),$$
(4)
where $`A_V`$ is defined by equation (1), and extinction in $`K`$ is neglected in comparison with $`A_V`$. Suppose a star has a circumstellar dust disk reradiating stellar optical and UV light in the infrared as well as reflecting unabsorbed stellar infrared radiation. In the $`E_{by}\delta (VK)_{by}`$ diagram, such a star would deviate upward from the line given by equation (4), the amount of deviation representing a measure of intrinsic infrared excess. Therefore, we can use this deviation as a criterion identifying candidate dusty stars:

$$\mathrm{\Delta }(VK)_{by}\delta (VK)_{by}A_V(E_{by})>0.$$
(5)
We have applied this criterion to our sample of F stars. Given that large masses of dust are expected only in very young systems, it makes sense to start a search for dusty stars by examining metal rich stars because high metallicity typically implies a young age. Therefore, we have restricted our search sample to stars with \[Fe/H\] $`>0`$. Also, since we expect a dusty system to be reddened, only stars with $`E_{by}>0`$ have been considered. These constraints reduce the working sample to $`900`$ stars. These stars are shown in Figure 2; those above the reddening line are predicted, in accordance with our criterion, to possess circumstellar dust.

The most conspicuous candidate dusty stars are shown in Figures 2 and 3 as filled circles; they are also listed in Table 1. These 21 stars have $`E_{by}>0.01`$ and $`\mathrm{\Delta }(VK)_{by}>0.1`$. About half of them are spatially associated with the giant molecular complex centered at the site of ongoing star formation in $`\rho `$ Oph.

IRAS fluxes in Table 1, $`F(12\mu )`$, $`F(25\mu )`$, and $`F(60\mu )`$, are from the IRAS Faint Source Catalog, Version 2. In total, we found seven out of the 21 selected candidate PMS dusty F stars to be IRAS sources (marked in Figure 2 by overlaid asterisk), which suggests the presence of heated circumstellar dust about them. These 7 stars are brighter than $`V7.5`$ mag (with a notable exception for HD 144432, which is 8.188 in $`V`$), which may hint that the thermal emission from possible dust disks around fainter F stars is simply below the detection limit of the IRAS; however, an intrinsic difference between the IRAS and non-IRAS stars cannot be ruled out. We have also searched through the ISO archive (using SIMBAD interface) to see if any of the stars in Table 1 was among the ISO pointed observations, but found none.

To gain insight into the nature of the predicted stars with circumstellar dust, we have examined the $`\mathrm{log}T_eM_V`$ diagrams for more extreme, 2 $`\sigma `$ candidates, $`\mathrm{\Delta }(VK)_{(by)_0}>0.07`$, and compared it with the one for “noncandidates” defined by $`\mathrm{\Delta }(VK)_{(by)_0}<0.03`$, selecting both groups from the sample shown in Figure 2. Inspection of Figure 3 reveals a large difference between the two diagrams. While the stars in the lower panel scatter across the entire temperature range displayed in the diagram, the stars in the upper panel exhibit a heavy concentration at the hot end of the temperature range where only very young stars can reside. Thus, the stars predicted to possess dusty disks (upper panel), are mostly very young whereas the non-candidates (lower panel) are dominated by a much older, on average, population. This is exactly what one would expect comparing stars with and without circumstellar dust.

The prominent clustering of the stars in the upper panel in Figure 3 within a narrow temperature range of $`3.86<\mathrm{log}T_e<3.83`$ (see also Figure 4), which span about 3 mag in $`M_V`$, offers a further clue to the nature of these objects. The large spread of absolute magnitude they exhibit is obviously inconsistent with the effects of stellar evolution away from the main sequence; if it were, the dominant population in the $`\mathrm{log}T_eM_V`$ diagram would have been at cooler rather than hotter temperatures, simply because the initial mass function increases toward lower masses. To put it differently, once we ask the question: “Where are stars in that diagram that must be numerous downward of $`\mathrm{log}T_e3.83`$?”, we would immediately realize that most of the hot stars above the ZAMS in Figure 3 are not post-ZAMS stars but rather are at various phases of PMS evolution—which perfectly matches the inferred presence of dust around these stars.

One may notice that also in the lower panel in Figure 3 and the respective middle panel in Figure 4, there are many hot stars separated from the bulk of cool stars by a poorly populated temperature interval around $`\mathrm{log}T_e3.832`$. The same arguments as above make us believe that these are also mostly PMS stars. Additional argument is that their reddening is about twice as large as that of the cool stars, $`E(by)=0.028\pm 0.003`$ versus $`0.012\pm 0.001`$; this should be expected if a significant part of it, if not all, occurs in a dense circumstellar environment of a PMS star. The absence of excessive NIR emission apparently indicates that the dust is colder than in the PMS candidates in the upper panel in Figures 3 and 4. Incidentally, these latter candidates are less reddened, $`E(by)=0.021\pm 0.001`$, then the stars with presumably cold dust, which should be expected if warm dust becomes more dominant when much of the original dust content has been dispersed.

Note that a large fraction of dusty PMS candidates inferred for the sample shown in Figure 3 is not at all typical for the general population of the field F stars. The temperature distribution of a sample that is more or less representative of a general population (the sample shown in Figure 1) is displayed in the lower panel of Figure 4. It is obvious that for this sample the fraction of stars in the temperature range dominated by PMS candidates is quite low and consistent with the expectations for an IMF declining toward higher stellar masses.

The dusty PMS candidates are contiguous in the $`\mathrm{log}T_eM_V`$ diagram with A stars that are known to contain a population of PMS emission-line, dusty stars usually referred to as a group as Herbig Ae/Be (HAeBe) stars. The properties of the candidates and the fact that they are adjacent to A stars suggest that they represent a continuation of the population of HAeBe stars into the spectral class F. Thé et al. (1994) extended the criteria for membership of the class of HAeBe stars to include emission line objects of spectral types earlier than F8. However, the number of stars later than F2 in their catalog of HAeBe objects is very small, only about half a dozen. The candidate PMS dusty F stars identified in this paper would represent an enormous increase in the number of F stars later than F2 that appear to be akin to HAeBe stars.

As discussed above, the stars in the upper panel of Figure 3 are concentrated within a narrow temperature range, $`3.860<\mathrm{log}T_e<3.832`$ (spectral range $``$ F2 to F5), their number plummeting downward of $`\mathrm{log}T_e3.84`$. This may have interesting implications pertinent to the nature of the class of HAeBe objects extended into the spectral range F, suggesting that such objects get extremely rare redward of F5. This phenomenon is likely related to fundamental properties of the PMS evolution of $`12M_{}`$ stars. So the detection of the cut-off in the temperature (spectral type) distribution of HAeBe+F objects at $`\mathrm{log}T_e3.84`$ (Sp $``$ F5) should be expected to provide new insight into the physics of PMS stars.

In summary, we have proposed a new criterion to identify F stars with warm circumstellar dust emitting in the near infrared. Application of this criterion to a sample of $`900`$ metal rich, moderately reddened F stars resulted in detection of $`70`$ 2-sigma candidates. They appear to be PMS stars, representing a continuation of the class of Herbig Ae/Be stars into the spectral range F and exhibiting a cutoff redward of $``$ F5. There is also evidence for a population of PMS stars with cold dust undetectable in the near infrared.


format: text

text:
# Characters of 𝒰_𝑞(𝑔𝑙(𝑛))-reflection equation algebra

## 1. Introduction

Let $`V`$ be a complex vector space and $`S\mathrm{End}^2(V)`$ a Yang-Baxter operator, satisfying the braid identity

$$S_{12}S_{23}S_{12}=S_{23}S_{12}S_{23}.$$

Here, $`S_{12}=S1`$ and $`S_{23}=1S`$. A matrix $`A\mathrm{End}(V)`$ is called a solution to the numerical reflection equation (RE) or a numerical RE matrix if the equality

(1) $`SA_2SA_2=A_2SA_2S,`$
where $`A_2=1A`$, holds in $`\mathrm{End}^2(V)`$. Numerical RE matrices define one-dimensional representations (characters) of the RE algebra, the quotient of the tensor algebra $`\mathrm{T}\left(\mathrm{End}^{}(V)\right)`$ by quadratic relations (1), \[KSkl, KS\]. RE and related algebraic structures appeared in the theory of integrable models \[Cher, Skl, AFS\]; they find applications in 3-dimensional topology \[K\] and problems of covariant quantization on G-spaces \[Do1, Do2\]. In the latter case, they are related to representations of quantum groups on $`V`$; then the operator $`S`$ is equal to $`PR`$, the product of the permutation $`P`$ in $`VV`$ and the image $`R`$ of the universal R-matrix, \[Dr\]. Numerical solutions to equation (1) are the key ingredient in the character method of covariant quantization on homogeneous spaces developed in \[DoM1, DoM2\].

In the present paper we study equation (1) associated with the fundamental vector representation of the quantum group $`𝒰_q\left(gl(n)\right)`$. An example of solution to this equation is the matrix $`D_n=\lambda _{i=1}^ne_{n+1i}^i`$, $`\lambda `$, as well as its similarity transformation by a diagonal invertible matrix, according to \[KSS\]. Therein, all the non-degenerate RE matrices were presented in dimensions $`n=2,3,4`$:

$$A^{1,1}=\left(\begin{array}{cc}\lambda +\mu & y_1\\ y_2& \end{array}\right),A^{2,1}=\left(\begin{array}{ccc}\lambda +\mu & & y_1\\ & \lambda & \\ y_3& & \end{array}\right),$$

$$A^{2,2}=\left(\begin{array}{cccc}\lambda +\mu & & & y_1\\ & \lambda +\mu & y_2& \\ & y_3& & \\ y_4& & & \end{array}\right),A^{3,1}=\left(\begin{array}{cccc}\lambda +\mu & & & y_1\\ & \lambda & & \\ & & \lambda & \\ y_4& & & \end{array}\right),$$

$$\text{where}y_iy_{n+1i}=\lambda \mu 0.$$

Here, the parameterization is adapted to our exposition. Besides, in \[DoM1\], there were found diagonal numerical RE matrices $`P_k=\lambda _{i=1}^ke_i^i`$, $`kn`$. The above mentioned matrices exhaust all (as far as we can judge) known solutions to (1) associated to the standard quantum group $`𝒰_q\left(gl(n)\right)`$. In the present paper, we list all possible solutions.

## 2. Result

The image of the universal R-matrix of the Drinfeld-Jimbo quantum group $`𝒰_q\left(gl(n)\right)`$ in the fundamental vector representation on $`^n`$ is, \[FRT\],

(2) $`R`$ $`=`$ $`q{\displaystyle \underset{i=1,\mathrm{},n}{}}e_i^ie_i^i+{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,j=1,\mathrm{},n}{ij}}{}}e_i^ie_j^j+\omega {\displaystyle \underset{\genfrac{}{}{0pt}{}{i,k=1,\mathrm{},n}{i<k}}{}}e_i^ke_k^i.`$
Here, $`\{e_j^i\}`$ is the standard basis in $`\mathrm{End}(^n)`$ with the multiplication $`e_j^ie_k^l=\delta _j^le_k^i`$ expressed through the Kronecker symbols $`\delta _j^l`$; $`\omega `$ stands for $`=qq^1`$. It is convenient to represent the corresponding braid matrix in the form

(6) $`S`$ $`=`$ $`{\displaystyle \underset{i,k=1,\mathrm{},n}{}}s_{ik}e_i^ie_k^k+{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,j=1,\mathrm{},n}{ij}}{}}e_j^ie_i^j,\text{where}s_{ik}=\{\begin{array}{cc}\omega ,\hfill & i<k,\hfill \\ q,\hfill & i=k,\hfill \\ 0,\hfill & i>k,\hfill \end{array}`$
We shall deal with subsets in $``$ and denote by $`[a,b]`$ the intervals $`\{k|akb\}`$. Instead of the square brackets, we use parentheses for the intervals defined by strict inequalities. To formulate the classification theorem, we need the following data.

###### Definition 2.1.

An admissible pair $`(Y,\sigma )`$ consists of an ordered subset $`YI=[1,n]`$ and a decreasing injective map $`\sigma :YI`$ without stable points.

Clearly, the map $`\sigma `$ is determined by its image $`\sigma (Y)`$. In the set $`Y`$, we distinguish two non-intersecting subsets $`Y_+=\{iY|i>\sigma (i)\}`$ and $`Y_{}=\{iY|i<\sigma (i)\}`$; obviously $`Y=Y_{}Y_+`$. Denote $`b_{}=\mathrm{max}\{Y_{}\sigma (Y_+)\}`$ and $`b_+=\mathrm{min}\{Y_+\sigma (Y_{})\}`$. Because the injection $`\sigma `$ is decreasing, $`b_{}<b_+`$. We adopt the convention $`b_{}=0`$ and $`b_+=n+1`$ if $`Y=\mathrm{}`$.

###### Theorem 2.2 (Classification).

General solution to equation (1) with the braid matrix (6) has the form

(7) $`A`$ $`=`$ $`{\displaystyle \underset{iI}{}}x_ie_i^i+{\displaystyle \underset{jY}{}}y_je_j^{\sigma (j)},`$
where

* $`(Y,\sigma )`$ is an admissible pair,

(8) $`Y`$ $`=`$ $`[1,b_{}][b_+,b_++b_{}1],`$
(9) $`\sigma (i)`$ $`=`$ $`b_++b_{}i`$
for $`b_{},b_+`$ subject to the conditions $`b_{}<b_+`$, $`b_{}+b_+n+1`$ and

(13) $`x_i`$ $`=`$ $`\{\begin{array}{cc}\lambda +\mu ,\hfill & i[1,b_{}],\hfill \\ \lambda ,\hfill & i(b_{},b_+),\hfill \\ 0,\hfill & i[b_+,n],\hfill \end{array}`$
(14) $`y_iy_{\sigma (i)}`$ $`=`$ $`\lambda \mu 0`$
for $`\lambda ,\mu `$;
* $`(Y,\sigma )`$ is an admissible pair such that $`Y\sigma (Y)=\mathrm{}`$ and

(17) $`x_i`$ $`=`$ $`\{\begin{array}{cc}\lambda ,& i[1,b],\\ 0,& i(b,n],\end{array}`$
(18) $`y_i`$ $``$ $`0`$
for $`b[b_{},b_+)`$ and $`\lambda `$.

###### Remark 2.3.

As follows from the theorem, there are two classes of numerical RE matrices, those corresponding to $`\sigma (Y)=Y`$ and $`\sigma (Y)Y=\mathrm{}`$. We call them solutions of Type 1 and 2, respectively.

## 3. Proof of the classification theorem

The proof of Theorem (2.2) formulated in the previous section is combinatorial. It is a result of the direct analysis of equation (1) organized into a sequence of lemmas. Putting $`A=_{\alpha ,\beta =1}^nA_\beta ^\alpha e_\alpha ^\beta `$, we compute

$$A_2SA_2=\underset{i,\beta ,\alpha ,\nu }{}s_{i\nu }A_\beta ^\nu A_\nu ^\alpha e_i^ie_\alpha ^\beta +\underset{\genfrac{}{}{0pt}{}{i,j,\beta ,\alpha }{ij}}{}A_\alpha ^jA_i^\beta e_j^ie_\beta ^\alpha .$$

Substituting this into the left- and right-hand sides of (1), we rewrite it in a more explicit form:

$`SA_2SA_2`$ $`=`$ $`{\displaystyle \underset{i,\beta ,\alpha ,\nu }{}}s_{i\beta }s_{i\nu }A_\beta ^\nu A_\nu ^\alpha e_i^ie_\alpha ^\beta +{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,j,\beta ,\alpha }{ji}}{}}s_{i\alpha }A_\alpha ^jA_i^\beta e_j^ie_\beta ^\alpha `$
$`+`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,\beta ,\alpha ,\nu }{i\beta }}{}}s_{i\nu }A_\beta ^\nu A_\nu ^\alpha e_i^\beta e_\alpha ^i+{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,j,\beta ,\alpha }{ji\alpha }}{}}A_\alpha ^jA_i^\beta e_j^\alpha e_\beta ^i,`$
$`A_2SA_2S`$ $`=`$ $`{\displaystyle \underset{i,\beta ,\alpha ,\nu }{}}s_{i\alpha }s_{i\nu }A_\beta ^\nu A_\nu ^\alpha e_i^ie_\alpha ^\beta +{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,j,\beta ,\alpha }{ji}}{}}s_{j\beta }A_\alpha ^jA_i^\beta e_j^ie_\beta ^\alpha `$
$`+`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,\beta ,\alpha ,\nu }{i\alpha }}{}}s_{i\nu }A_\beta ^\nu A_\nu ^\alpha e_\alpha ^ie_i^\beta +{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,j,\beta ,\alpha }{ij\beta }}{}}A_\alpha ^jA_i^\beta e_\beta ^ie_j^\alpha .`$
Comparison of (3) with (3) gives rise to the system of quadratic equations on the matrix elements $`A_j^i`$.

###### Lemma 3.1.

Equation (1) is equivalent to the following system of equations:

(23) $`\{\begin{array}{ccc}\hfill A_i^mA_i^n& \hfill =& 0,\hfill \\ \hfill A_m^iA_n^i& \hfill =& 0,\hfill \end{array}mnim,`$
(26) $`A_i^nA_m^j=0,\{\begin{array}{c}jmni,\hfill \\ (mi)(nj)<0,\hfill \end{array}`$
(30) $`\{\begin{array}{cccc}\hfill (qs_{im})A_i^iA_i^m& \hfill =& _\nu s_{i\nu }A_i^\nu A_\nu ^m,\hfill & im,\hfill \\ \hfill (qs_{im})A_i^iA_m^i& \hfill =& _\nu s_{i\nu }A_m^\nu A_\nu ^i,\hfill & im,\hfill \\ \hfill 0& \hfill =& _\nu s_{i\nu }A_m^\nu A_\nu ^n,\hfill & (mi)(ni)<0,\hfill \end{array}`$
(31) $`A_i^nA_m^i{\displaystyle \underset{\nu }{}}s_{i\nu }A_m^\nu A_\nu ^n=(s_{ni}s_{im})A_i^iA_m^n,minm,`$
(32) $`\omega A_m^mA_i^i={\displaystyle \underset{\nu }{}}s_{i\nu }A_m^\nu A_\nu ^m{\displaystyle \underset{\nu }{}}s_{m\nu }A_i^\nu A_\nu ^i,i<m`$
###### Proof.

This statement is verified by the direct analysis. ∎

The next lemma accounts for equations (23) and (26).

###### Lemma 3.2.

If $`A`$ a solution to equation (1), then it can be represented in the form (7), where $`(Y,\sigma )`$ is an admissible pair, $`x_i=A_i^i`$ for $`iI`$, and $`y_i=A_i^{\sigma (i)}0`$ for $`iY`$.

###### Proof.

By virtue of equations (23), the matrix $`A`$ has at most one non-zero off-diagonal entry in every row and every column. Indices of such rows form a subset $`Y`$ in $`I`$, and the off-diagonal entries may be written as $`y_i=A_i^{\sigma (i)}0`$, $`iY`$, for some bijection $`\sigma :YI`$ with no stable points. This is equivalent to equation (23). The map $`\sigma `$ is decreasing; that is encoded in equation (26). ∎

###### Lemma 3.3.

The equations of system (30) are equivalent to

(33) $`{\displaystyle \underset{\alpha \mathrm{max}(i,m)}{}}A_i^\alpha A_\alpha ^m=0,im,`$
by virtue of (26). They imply $`x_i=0`$ for $`i=b_+`$.

###### Proof.

Setting $`i=j`$ in (26) and substituting it into the equations of system (30), we reduce them to (33). Setting $`m=\sigma (i)<iY_+`$ in (33) leads to $`x_iy_i=0`$ and therefore $`x_i=0`$. Similarly, the assumption $`m=\sigma (i)>iY_{}`$ reduces (33) to $`y_ix_{\sigma (i)}=0`$ and hence $`x_{\sigma (i)}=0`$. So $`x_i=0`$ for all $`iY_+\sigma (Y_{})`$ and $`i=b_+`$ in particular ∎

###### Lemma 3.4.

Lemma 3.3 taken into account, equation (31) is equivalent to the following two assertions.

1. For any $`mY`$, either $`\sigma ^2(m)=m`$ or $`\sigma (m)Y`$.
2. $`x_i=x_b_{}`$ and $`x_j=0`$ whenever $`ib_{}`$ and $`jb_+`$.

###### Proof.

First note that if $`mY`$, equation (31) holds identically. So we can assume $`mY`$ and rewrite (31) as

(34) $`A_i^nA_m^is_{im}A_m^mA_m^ns_{i\sigma (m)}A_m^{\sigma (m)}A_{\sigma (m)}^n=(s_{ni}s_{im})A_i^iA_m^n,`$
$$minm.$$

Supposing $`n\sigma (m)`$ we find $`A_i^nA_m^is_{i\sigma (m)}A_m^{\sigma (m)}A_{\sigma (m)}^n=0.`$ If $`i=\sigma (m)`$, then, having in mind $`s_{ii}=q1`$ and $`A_m^{\sigma (m)}=y_m0`$, we obtain $`A_{\sigma (m)}^n=0`$. Since $`n\sigma (m)=i`$ and $`nm`$, this means either $`\sigma (m)Y`$ or $`\sigma ^2(m)=m`$. Assuming $`i\sigma (m)`$, we come to the equation $`s_{i\sigma (m)}A_m^{\sigma (m)}A_{\sigma (m)}^n=0`$, which is fulfilled as well, due to $`A_{\sigma (m)}^n=0`$.

It remains to study the case $`n=\sigma (m)`$. Under this hypothesis, the term $`A_i^nA_m^i`$ vanishes. Indeed, since $`im`$, one has $`A_m^i0i=\sigma (m)`$. But this contradicts the condition $`in=\sigma (m)`$. In terms of the variables $`x_i`$ and $`y_i`$, equation (34) reads

$$s_{im}x_my_ms_{i\sigma (m)}y_mx_{\sigma (m)}=(s_{\sigma (m)i}s_{im})x_iy_m,$$

$$mi\sigma (m).$$

Depending on allocation of the indices $`i`$, $`m`$, and $`\sigma (m)`$, this equation splits into the following four implications.

(35) $`i<m\text{and}i<\sigma (m)`$ $``$ $`x_i=x_m+x_{\sigma (m)},`$
(36) $`i<m\text{and}i>\sigma (m)`$ $``$ $`x_m=0,`$
(37) $`i>m\text{and}i<\sigma (m)`$ $``$ $`x_{\sigma (m)}=0,`$
(38) $`i>m\text{and}i>\sigma (m)`$ $``$ $`x_i=0.`$
Recall that the index $`m`$ is assumed to be from $`Y`$. Equation (38) is equivalent to $`x_i=0`$ for $`i>b_+`$. By Lemma 3.3, this is also true for $`ib_+`$. Equations (36) and (37) are thus satisfied as well. Equation (35) states $`x_i=x_b_{}+x_{\sigma (b_{})}`$ if $`i<b_{}Y_{}`$ or $`x_i=x_b_{}+x_{\sigma ^1(b_{})}`$ if $`i<b_{}\sigma (Y_+)`$. Applying Lemma 3.3, we find $`x_i=x_b_{}`$, in either cases. ∎

###### Remark 3.5.

We used Lemma 3.3 in the proof of Lemma 3.4 in order to find the values $`x_{b_\pm }`$ on the boundaries of the intervals $`[1,b_{}]`$ and $`[b_+,n]`$. If $`(b_{},b_+)\mathrm{}`$, equations (36) and (37) are enough for that purpose, because we can put $`i(b_{},b_+)`$ there and come to the same result.

###### Lemma 3.6.

Equation (33) is fulfilled by virtue of Lemma 3.4.

###### Proof.

Consider the case $`m<i`$. Equation (33) holds if $`iY`$, because the sum turns into $`A_i^iA_i^m`$. So we may assume $`iY`$ and distinguish two cases: $`iY_{}`$ and $`iY_+`$. Assumption $`i<\sigma (i)`$ leads to $`A_i^iA_i^m+A_i^{\sigma (i)}A_{\sigma (i)}^m=0`$. The equality $`m=\sigma (i)`$ is impossible, since otherwise $`m=\sigma (i)<i<\sigma (i)`$. Therefore the first term vanishes and we come to $`A_{\sigma (i)}^m=0`$, $`mi`$. This condition is satisfied, by Lemma 3.4, Statement 1. The case $`iY_+`$ results in $`x_iA_i^m=0`$. Setting $`m=\sigma (i)`$, we come to $`x_i=0`$. Once $`iY_+ib_+`$, we encounter a particular case of Lemma 3.4, Statement 2.

We should study the situation $`i<m`$. Equation (33) evidently holds if $`\sigma (i)<m`$. We may think that $`m\sigma (i)`$; then

$$\begin{array}{ccccc}i<m=\sigma (i)\hfill & \hfill & \hfill A_m^m=0& \hfill & x_{\sigma (i)}=0,\hfill \\ i<m<\sigma (i)\hfill & \hfill & \hfill A_{\sigma (i)}^m=0& \hfill & \sigma (i)Y.\hfill \end{array}$$

These requirements are fulfilled, by Lemma 3.4. ∎

It remains to satisfy equation (32), to complete the proof of Theorem 2.2. Let $`Y_0`$ be the subset in $`Y`$, such that $`\sigma `$ restricted to $`Y_0`$ is involutive. By Lemma 3.4, either $`\sigma (i)Y`$ or $`i`$ and $`\sigma (i)`$ belong to $`Y_0`$ simultaneously. Equation (32) falls into the four equations

(39) $`x_m(x_ix_m)`$ $`=`$ $`0,iY_0,mY_0,`$
(40) $`\omega x_m(x_ix_m)`$ $`=`$ $`s_{i\sigma (m)}y_my_{\sigma (m)},iY_0,mY_0,`$
(41) $`\omega x_m(x_ix_m)`$ $`=`$ $`s_{m\sigma (i)}y_iy_{\sigma (i)},iY_0,mY_0,`$
(42) $`\omega x_m(x_ix_m)`$ $`=`$ $`s_{i\sigma (m)}y_my_{\sigma (m)}s_{m\sigma (i)}y_iy_{\sigma (i)},iY_0,mY_0.`$
Everywhere $`i<m`$, see (32).

###### Lemma 3.7.

Suppose $`Y_0\mathrm{}`$. Then $`Y=Y_0`$ and, moreover,

(43) $`Y_{}=\{1,\mathrm{},b_{}\},Y_+=\{b_+,\mathrm{},b_++b_{}1\},`$
(44) $`\sigma (i)=b_++b_{}i.`$
###### Proof.

Let $`kY_0Y_{}`$ and $`\sigma (k)Y_0Y_+`$ . Suppose either $`lY_{}\backslash Y_0`$ or $`l<\mathrm{min}(Y_{})`$. The assumption $`l<k`$ contradicts equation (40) if one sets $`i=l`$, $`m=k`$. The inequality $`k<l`$ does not agree with (41) if one sets $`i=k`$, $`m=l`$. In both cases one uses Lemma 3.4, Statement 2, and gets $`0=y_ky_{\sigma (k)}`$; that is impossible since $`y_k0`$ for all $`kY`$, by definition of $`Y`$. Therefore $`Y_{}Y_0`$ and $`\mathrm{min}(Y_{})=1`$. Assume now $`lY_+\backslash Y_0`$. Then, $`l`$ cannot exceed $`\sigma (k)`$, because otherwise $`\sigma (l)<\sigma ^2(k)=k\sigma (l)Y_0lY_0`$, the absurdity. The only possibility is $`l<\sigma (k)`$. But this again contradicts equation (41) if one sets $`i=k`$ and $`m=l`$. ∎

###### Lemma 3.8.

Suppose $`Y_0\mathrm{}`$. Then, there is $`a`$ such that $`y_iy_{\sigma (i)}=a`$ for all $`iY`$. If $`(b_{},b_+)\mathrm{}`$, then $`x_m=\lambda `$ for all $`m(b_{},b_+)`$, where $`\lambda `$ is a solution to the quadratic equation $`\lambda (\lambda x_b_{})=a`$.

###### Proof.

By Lemmas 3.4 and 3.7, equation (42) is fulfilled if and only if the product $`y_iy_{\sigma (i)}`$ does not depend on $`iY`$ and is equal to some $`a`$. If $`(b_{},b_+)\mathrm{}`$, equation (41) suggests $`x_m(x_mx_b_{})=a`$ as soon as $`m(b_{},b_+)`$. It is easy to see that equation (39) is equivalent to $`x_m=x_i=\lambda `$ for $`i,m(b_{},b_+)`$ and some $`\lambda `$. ∎

It is convenient to introduce the parameterization $`x_b_{}=\lambda +\mu `$, $`a=\mu \lambda `$. Then $`x_m=\lambda +\mu `$ for $`mb_{}`$, $`x_m=\lambda `$ for $`b_{}<m<b_+`$, and $`y_iy_{\sigma (i)}=\lambda \mu `$. This parameterization makes sense even if $`(b_{},b_+)=\mathrm{}`$. The scalars $`\lambda `$ and $`\mu `$ will have the meaning of eigenvalues of the matrix $`A`$.

###### Lemma 3.9.

Let $`A`$ be a matrix satisfying equations (23)-(31) and $`(Y,\sigma )`$ the corresponding admissible pair. Then, equation (32) gives rise to the alternative

* $`Y_0=Y`$. $`A`$ is a solution of Type 1 from Theorem 2.2.
* $`Y_0=\mathrm{}`$. $`A`$ is a solution of Type 2 from Theorem 2.2.

###### Proof.

If $`Y_0=\mathrm{}`$, then $`Y\sigma (Y)=\mathrm{}`$ by Lemma 3.4, Statement 1. Equation (32) is reduced to (39). It is satisfied if and only if $`x_m=\lambda `$, $`ib`$, and $`x_m=0`$, $`b<i`$, for some $`\lambda `$ and $`b`$. Comparing this with Lemma 3.4, Statement 2, we come to a solution of Type 2 with $`b[b_{},b_+)`$.

Suppose $`Y_0\mathrm{}`$. Then, equation (32) implies Lemmas 3.7 and 3.8. Conversely, these lemmas ensure (39)–(42), which are equivalent to (32). Altogether, this gives a solution of Type 1. ∎

Thus we accomplish the proof of Theorem 2.2 and proceed to the analysis of the solutions obtained.

## 4. Structure of solutions

The first question about the structure of the numerical RE matrices listed in Theorem 2.2 is what pairs $`(Y,\sigma )`$ may participate in the classification. The solutions of Type 1 involve two integers $`b_{}`$, $`b_+`$ subject to the conditions of the theorem. They determine the admissible pair $`(Y,\sigma )`$. The solutions of Type 2 are labeled with subsets $`YI`$ and injective maps $`\sigma :YI\backslash Y`$. Clearly, $`K=\mathrm{card}(Y)\frac{n}{2}`$ and, with $`Y`$ given, there are $`C_{nK}^K`$ possibilities for $`\sigma `$, i.e. the number of subsets in $`I\backslash Y`$ with $`\mathrm{card}=K`$. Further we describe the properties of solutions to equation (1). We consider the standard basis $`\{e^i\}`$ in $`^n`$; the action of $`\mathrm{End}(^n)`$ on $`^n`$ being given by $`e_j^ie^k=\delta _j^ke^i`$.

###### Proposition 4.1.

Let $`(Y,\sigma )`$ be an admissible pair corresponding to a numerical RE matrix $`A`$. The subspaces $`V_i^2=e^ie^{\sigma (i)}`$, $`iY`$, and $`V_i^1=e^i`$, $`iY\sigma (Y)`$, are $`A`$-invariant.

###### Proof.

This follows from the general form (7) of the RE matrices and from Lemma 3.4, Statement 1. ∎

The spectral properties of numerical RE matrices are described by the following proposition.

###### Proposition 4.2.

An RE matrix of Type 1 from Theorem 2.2 has eigenvalues $`\mu `$, $`\lambda `$, and $`0`$ of multiplicities $`b_{}`$, $`b_+b_{}1`$, and $`n(b_{}+b_+)+1`$, respectively. It is semisimple if and only if $`\lambda \mu `$.

An RE matrix of Type 2 from Theorem 2.2 has eigenvalues $`\lambda `$ and $`0`$ of multiplicities $`b`$ and $`bK`$ respectively. It is semisimple if and only if $`\lambda 0`$.

###### Proof.

Restricted to $`V_i^2`$, $`i[1,b_{}]`$, an RE matrix $`A`$ of of Type 1 is equal to $`\begin{array}{cc}\lambda +\mu & y_i\\ y_{\sigma (i)}& 0\end{array}`$. Because $`y_iy_{\sigma (i)}=\lambda \mu `$, this $`2\times 2`$ matrix has eigenvalues $`\mu `$ and $`\lambda `$. On the subspaces $`V_i^1`$, it acts as multiplication by $`\lambda `$ when $`i(b_{},b_+)`$, and $`0`$ when $`i[b_{}+b_+,n]`$.

Restricted to $`V_i^2`$, $`iY`$, an RE matrix $`A`$ of of Type 2 is equal to $`\begin{array}{cc}\lambda & y_i\\ 0& 0\end{array}`$ if $`iY_{}`$ and $`\begin{array}{cc}\lambda & 0\\ y_i& 0\end{array}`$ if $`iY_+`$. In either cases, it has eigenvalues $`\lambda `$ and $`0`$. Observe that $`Y_{}\sigma (Y_+)[1,b]`$ and $`Y_{}\sigma (Y_+)=\mathrm{}`$; thus $`\mathrm{card}\left(Y_{}\sigma (Y_+)\right)=\mathrm{card}(Y)=K`$. Therefore there are $`bK`$ one-dimensional subspaces $`V_i^1`$ where $`A`$ acts as multiplication by $`\lambda `$. So the eigenvalue $`\lambda `$ has multiplicity $`b`$. The zero eigenvalue enters with multiplicity $`nb`$. ∎

Note that when one of the eigenvalues $`\lambda ,\mu `$ tends to zero, the RE matrices of Type 1 turns into RE matrices of Type 2. Non-degenerate numerical RE matrices belong to the first class. This is the case when $`b_{}+b_+1=n`$ and $`\lambda 0\mu `$. The RE matrices from \[KSS\] written out in Introduction are of this kind, and the Type 1 solutions from Theorem 2.2 are their generalization to higher dimensions. Let us illustrate Theorem 2.2 on those examples.

1. The matrix $`D_n`$ is a solution of Type 1. Here, $`\lambda =\mu =y_i`$, $`iY`$, where the set $`Y`$ coincides with $`I`$ if $`n=0mod2`$ and $`Y=I\backslash \{\frac{n+1}{2}\}`$ if $`n=1mod2`$. The similarity transformation by a diagonal matrix changes the elements $`y_i`$ however preserving the condition $`y_iy_{n+1i}=\lambda ^2`$.
2. The matrix $`A^{1,1}`$ has $`b_{}=1`$, $`b_+=2`$, while for the matrix $`A^{2,2}`$ these parameters take the values $`b_{}=2`$ and $`b_+=3`$. These solutions have in common $`(b_{},b_+)=\mathrm{}`$. On the contrary, the interval $`(b_{},b_+)`$ is not empty in the case of matrices $`A^{2,1}`$ and $`A^{3,1}`$. For them, one has $`b_{}=1`$, $`b_+=3`$ and $`b_{}=1`$, $`b_+=4`$.
3. The diagonal matrix $`P_k=\lambda _{i=1}^ke_i^i`$ is an RE matrix of Type 2 with $`Y=\mathrm{}`$ and $`k=bn`$.

###### Remark 4.3.

Numerical RE matrices play an important role in $`𝒰_q\left(gl(n)\right)`$-covariant quantization on the adjoint orbits in $`\mathrm{End}^2(^n)`$, \[DoM1, DoM2\]. They allow to explicitly represent the quantized algebras of functions on orbits as quotients of the RE algebra and, simultaneously, as subalgebras of functions on the quantum group. Theorem 2.2 classify all such realizations. It implies, in particular, that all the symmetric and bisymmetric orbits (consisting of matrices with two and three different eigenvalues) can be represented in this way.


text:
# 1 Introduction

## 1 Introduction

One of the many pioneer contributions by A.C. Newell has been the study of solitary waves (which we will simply call “solitons” here, in line with currently adopted terminology). These occur not only in integrable models, but also in many nonintegrable nonlinear-wave systems . It is toward an improved understanding of solitons in nonintegrable systems, that we present this work, dedicated to Prof. Newell.

In nonintegrable models, more complex forms of solitons are found as one considers higher-order systems. An issue of fundamental importance is to find where, in the space of the soliton parameters, such solutions could exist. Where this could be is dominantly determined by the model’s linear dispersion relation, $`\omega (k)`$. Typically, dispersion curves have gaps or partial gaps, which are values of $`\omega `$ at which $`k`$ would take imaginary or complex values. For example, for the nonlinear Schrödinger (NLS) equation, one has $`\omega =Dk^2`$, where $`D`$ is a positive dispersion coefficient, hence $`k`$ must be imaginary if $`\omega <0`$. This is precisely the region where one finds solitons of the NLS equation. For more complex systems, with several branches of the dispersion relation, one can have regions of $`\omega `$ where $`k`$ can be complex as well as real. Here one typically encounters “delocalized solitons” (see the book  and relevant examples from nonlinear-optical models in Refs. ), which are quasi-solitary waves with nonvanishing oscillatory tails. Such objects, obviously, have an infinite energy, and therefore are unphysical except, possibly, in finite-size systems. However, it may occur that the amplitude of the tail vanishes at some special values of $`\omega `$. Then, one has a truly localized object at an isolated (discrete) value of $`\omega `$ at which a real value of $`k`$ does exist. Since this value of $`\omega `$ lies in a continuous part of the spectrum, these objects are called embedded solitons (ESs). A number of examples of ESs in physically meaningful models are now known (see examples of such solitons found in the hydrodynamical models in Refs. , and a short review in Ref. ), their stability (which turns out to be semi-stability) being quite distinct from the stability of ordinary solitary waves .

Currently, the only method for locating ESs is to search for them numerically. One does know that they may be found inside a partial gap of the dispersion relation, but otherwise one has no analytical tool for locating them (in exceptional cases, exact ES solutions can be found by guess ). The first objective of this paper is to develop an approximate analytical method for locating ESs. The method is based on the variational approximation (VA; see a recent review of the application of this technique to solitons in Ref. ), which incorporates an amplitude of the infinitesimal tail of the related delocalized soliton as a key variational parameter. The approach will be tested on two versions of an ES-generating model introduced in Ref. , which combines the second-harmonic generation (SHG) through quadratic \[$`\chi ^{(2)}`$\] nonlinearity, and the usual cubic \[$`\chi ^{(3)}`$\] nonlinear terms. The difference between the two versions of the model is that one is a full model (it has the usual Lagrangian representation), in terms of the $`\chi ^{(3)}`$ terms, while the one model is a *truncated model*, wherein, upon assuming that the fundamental-harmonic (FH) field is much stronger than the second harmonic (SH), some of the SH $`\chi ^{(3)}`$ terms are omitted. The truncated model does not admit the usual Lagrangian representation. However, it was this truncated version of the $`\chi ^{(2)}:\chi ^{(3)}`$ model in which the above-mentioned exact analytic ES solution was found in Ref. . On the other hand, in the full version of the model, ES solutions could only be found by means of numerical methods. We will obtain a natural variational criterion which makes it possible to distinguish ESs from delocalized solitons in each of these systems. Furthermore, we will demonstrate that, in the case of the truncated model, this method yields an exact result, and in the full model, a relative error is $`1\%`$ in the prediction of the location of ES (in comparison with numerical results).

A distinctive feature of the truncated model is that, due to a missing (omitted) term, it does not have a “complete” Lagrangian representation. Namely, it can be obtained from a Lagrangian, but only if one term in it is not subjected to the variation with respect to the FH field, when deriving the system of the FH and SH equations. Accordingly, the truncated model does not conserve any Hamiltonian. But it is, nevertheless, a conservative system, as it conserves the norm of the solution (which is usually called “energy” in nonlinear optics and is different from the Hamiltonian).

As the truncated model does not admit the full Lagrangian representation, it is necessary to work out a special version of the VA for it, which is another objective of this paper. In fact, this can be done in a very simple way: after inserting the variational ansatz into the Lagrangian and performing the integration in order to obtain the corresponding effective Lagrangian (the one which is an explicit function of variational parameters, rather than a functional depending on the field variables), the above-mentioned special term in the effective Lagrangian should not be varied with respect to variational parameters that belong to the FH component of the ansatz. Thus, other models similar to the above-mentioned truncated one, even though they cannot be represented in the usual Lagrangian form, can still be handled by means of VA. Since one cannot freely vary all the fields in the Lagrangian of these systems, we call them semi-Lagrangian systems.

The rest of the paper is organized as follows. In section 2, we introduce the full and truncated models and the Lagrangian representations of each. In section 3, the VA for the semi-Lagrangian case is developed. The general VA-based analytical condition for identifying ESs (in an approximate form) is obtained in section 4. Section 5 concludes the paper.

## 2 The full and truncated models

Following Refs.  and , we first introduce the full $`\chi ^{(2)}:\chi ^{(3)}`$ model:

$$iu_z+(1/2)u_{tt}+u^{}v+\gamma _1\left(|u|^2+2|v|^2\right)u=0,$$
(1)
$$iv_z(1/2)\delta v_{tt}+qv+(1/2)u^2+2\gamma _2\left(|v|^2+2|u|^2\right)v=0,$$
(2)
which is written in the usual “optical” notation, so that $`z`$ and $`t`$ are the propagation distance and reduced time, $`u`$ and $`v`$ are the FH and SH fields, $`\delta `$ is the relative dispersion coefficient at SH, $`q`$ is the SHG mismatch, and $`\gamma _{1,2}`$ are Kerr coefficients. In fact, the ratio of the self-phase-modulation (SPM) and cross-phase-modulation (XPM) coefficients is not necessarily $`1:2`$ simultaneously in both equations, as it is written in Eqs. (1) and (2), but this feature of the model is not a crucially important one.

The form of Eqs. (1) and (2) implies that the group-velocity dispersion is anomalous at FH, while at SH it may be both normal, if $`\delta >0`$, and anomalous, if $`\delta <0`$ (both cases are physically possible). As for the Kerr coefficients $`\gamma _1`$ and $`\gamma _2`$, they always have one sign. Most typically, they are positive (corresponding to the self-focusing nonlinearity), but may be negative too, see a detailed discussion of this point in Ref. .

In many cases, the SH field is much weaker than the FH field - for instance, if the mismatch is large. Then, assuming that $`\left|v\right|^2\left|u\right|^2`$, one may neglect the XPM term in comparison with the SPM one in Eq. (1), and the SPM term in comparison with its XPM counterpart in Eq. (2), which leads to the truncated model,

$$iu_z+(1/2)u_{tt}+u^{}v+\gamma _1|u|^2u=0,$$
(3)
$$iv_z(1/2)\delta v_{tt}+qv+(1/2)u^2+4\gamma _2|u|^2v=0.$$
(4)
In this work, we are interested in stationary fundamental-soliton solutions, which are looked in the form

$$u(z,t)=e^{ikz}U(t),\text{ }v(z,t)=e^{2ikz}V(t),$$
(5)
where $`k`$ is the FH wavenumber, and real even functions $`U(t)`$ and $`V(t)`$ with a single maximum at $`t=0`$ exponentially decay at $`t\mathrm{}`$. Note that ordinary (non-embedded) solitons may exist in the regions

$$0<k<q/2,\mathrm{if}\text{ }\delta >0;k>\mathrm{max}\{0,q/2\},\mathrm{if}\delta <0$$
(6)
(which implies that $`q`$ must be positive for the existence of ordinary solitons if $`\delta `$ is positive). On the other hand, ESs may exist in a range of $`k`$ which does not overlap with the continuous spectrum in the FH equation, (1) or (3), but falls into the continuous spectrum of the SH equation, (2) or (4). Thus, ES may exist in the regions \[cf. Eqs. (6)\]

$$k>\mathrm{max}\{0,q/2\},\mathrm{if}\text{ }\delta >0;\mathrm{\hspace{0.17em}0}<k<q/2,\mathrm{if}\delta <0.$$
(7)
This implies that, in the case $`\delta <0`$, the existence of ESs makes it necessary to have positive $`q`$.

Substituting the expressions (5) into Eqs. (1) and (2), we arrive at a system of ordinary differential equations,

$$kU+(1/2)U^{\prime \prime }+U^{}V+\gamma _1\left(U^2+2V^2\right)U=0,$$
(8)
$$2kV(1/2)\delta V^{\prime \prime }+qv+(1/2)U^2+2\gamma _2\left(V^2+2U^2\right)V=0,$$
(9)
where the prime stands for $`d/dt`$. Accordingly, the system of stationary equations corresponding to the truncated model is

$$kU+(1/2)U^{\prime \prime }+U^{}V+\gamma _1U^3=0,$$
(10)
$$2kV(1/2)\delta V^{\prime \prime }+qV+(1/2)U^2+4\gamma _2U^2V=0.$$
(11)
It is obvious that Eqs. (8) and (9) can be derived from the Lagrangian

$$L=\frac{1}{2}_{\mathrm{}}^+\mathrm{}\left[kU^2\left(2k+q\right)V^2\frac{1}{2}\left(U^{}\right)^2+\frac{\delta }{2}\left(V^{}\right)^2+U^2V+\frac{\gamma _1}{2}U^4+4\gamma _2U^2V^2+\frac{\gamma _2}{2}V^4\right]𝑑t.$$
(12)
In the case of the truncated system of stationary equations, (10) and (11), one may still obtain this system from the variational principle, provided we do the following. First, the last term in the integrand in Eq. (12) should be dropped. Second, the next to the last term should only be subjected to varying in $`V`$, but not in the $`U`$ field.

## 3 The variational approximation for the semi-Lagrangian system

In this section, we focus on the application of VA to the semi-Lagrangian (truncated) system (10), (11). The most natural variational ansatz to search for embedded solitons proper (without the tail, which will be considered in the next section) in this system is based on the following expressions:

$$U=A\mathrm{sech}\left(\sqrt{2k}x\right)x,V=B\mathrm{sech}^2\left(\sqrt{2k}x\right),$$
(13)
where the amplitudes $`A`$ and $`B`$ are variational parameters, while the inverse width $`\sqrt{2k}`$ has been fixed to match to the linearized form of Eqs. (8) and (9) or (10) and (11) at $`t\mathrm{}`$. Substituting the ansatz (13) into the Lagrangian (12), dropping the last term, and then performing the integration, we find the effective Lagrangian of the truncated system:

$$3\sqrt{2k}L_{\mathrm{eff}}=4kA^22\left[2\left(1\frac{2\delta }{5}\right)kq\right]B^2+2A^2B+\gamma _1A^4+\frac{35}{5}\gamma _2A^2B^2.$$
(14)
As it was said above, in order to go from the full system to its truncated counterpart, one not only has to omit the last term in Eq. (12), but also must take care to not vary the term $`4\gamma _2U^2V^2`$ with respect to $`U`$. In terms of the effective Lagrangian (14), this means that, when deriving the variational equations for $`A`$ and $`B`$, one should not vary the last term of the effective Lagrangian with respect to $`A`$. In this case, the variation with respect to $`A`$ yields a simple equation which allows one to eliminate $`B`$,

$$B=\gamma _1A^2+2k.$$
(15)
Using this result, the equation produced by varying the Lagrangian with respect to $`B`$ can be cast into the final form of a biquadratic equation for $`A`$,

$$\frac{32}{5}\gamma _1\gamma _2A^4\left[1+\frac{64}{5}k\gamma _2+2\gamma _1\left(2\left(1\frac{2\delta }{5}\right)kq\right)\right]A^2+4k\left[2\left(1\frac{2\delta }{5}\right)kq\right]=0.$$
(16)
Thus, depending on values of the parameters, Eq. (16) may give up to two different physical solutions. The simplest nontrivial case (which still allows an ES to exist) is with $`\gamma _1=0`$, but $`\gamma _20`$. Then, Eq. (16) yields a single solution,

$$A^2=\frac{4k\left[2\left(52\delta \right)k5q\right]}{5+64k\gamma _2},$$
(17)
which is physical if it gives $`A^2>0`$; the amplitude $`B`$ can then be obtained from Eq. (15).

It should be stressed that, as it was demonstrated in Ref. , the ansatz (13) yields exact soliton solutions to Eqs. (10) and (11) at some uniquely selected value of $`k`$ \[an expression for it is given below in Eq. (25)\]. Depending on the value of the mismatch parameter $`q`$, this exact solution may be either ES or an ordinary soliton. Comparing Eqs. (15) and (16) with that solution, one can easily verify that expressions (15) and (16) are precisely parts of the exact solution.

Although the exact solution for the ES in the truncated model is available at the *single* value of $`k`$ at which it may exist, in the region (6), where a *continuous* family of ordinary solitons is expected to exist, no general exact solution is known. So, to illustrate the accuracy and reliability of the modified VA for producing approximate solutions to the truncated equations (10) and (11), in Fig. 1 we display a typical example of a numerically found ordinary soliton in the region $`0<k<q/2`$ (in this example, $`\delta =1`$), together with the analytical approximation generated by Eqs. (13), (15), and (16). As one looks at the numerical solution in Fig. 1, one observes that the SH component goes slightly negative along the shoulder of the soliton, and then appears to oscillate as it decays. Note that in the region where this behavior of the SH component is observed, $`U`$ is not small, and, according to Eq. (4), the solution for $`V`$ in this region should oscillate indeed. Thus, ordinary solitons in the truncated system may have fine features which are not found in the ESs.

## 4 An analytical criterion to identify embedded solitons

### 4.1 General analysis

In the case when the wavenumber $`k`$ falls into the regions (7), ESs can exist in both the full and truncated models, i.e., (1), (2) and (3), (4) . However, except for using the exact ES solutions in the truncated model to guess where ESs might exist in its full counterpart, the only known method for locating ESs in the full model was to search for them by means of direct numerical computations.

Thus, there is a need for an analytical approach to the quest for ESs. Such an approach can be based on VA, if one assumes that, at values of $`k`$ close to $`k_{\mathrm{ES}}`$, at which an ES exists, there also exists a family of delocalized solitons, with small-amplitude oscillating tails in the SH component that vanish when $`k`$ become equal to $`k_{\mathrm{ES}}`$. As it immediately follows from Eqs. (9) or (11), the free oscillating tail of a delocalized soliton, which has an infinitesimal amplitude $`b`$ and arbitrary phase shift $`\psi `$, is given by the expression

$$V_{\mathrm{tail}}=b\mathrm{cos}\left(\sqrt{\left(2/\delta \right)\left(2kq\right)}t+\psi \right),$$
(18)
that neglects the nonlinear terms. Now, one may add this tail to the ansatz (13), to have a more general tailed ansatz,

$$V(t)=V_{\mathrm{sol}}(t)+b\mathrm{cos}\left(\sqrt{\left(2/\delta \right)\left(2kq\right)}t+\psi \right),$$
(19)
where $`V_{\mathrm{sol}}(t)`$ corresponds to the ansatz for the core of the delocalized soliton.

As $`b`$ is an extra variational parameter, one should add the variational equation

$$L_{\mathrm{eff}}/b=0$$
(20)
to the set of equations obtained by varying the effective Lagrangian with respect to all the other free parameters (irrespective of the fact if the system is complete Lagrangian or semi-Lagrangian). As we are interested in the location of an ES which, by itself, has $`b=0`$ (no tail), one should set $`b=0`$ *after* completing the differentiation in Eq. (20). This means that, prior to varying in $`b`$, one should only keep terms in $`L_{\mathrm{eff}}`$ which are linear in $`b`$, hence Eq. (20) takes the general form,

$$_{\mathrm{}}^+\mathrm{}\left[\left(\frac{\delta L}{\delta V}\right)|_{U=U_{\mathrm{sol}}(t),V=V_{\mathrm{sol}}(t)}\right]\mathrm{cos}\left(\sqrt{\left(2/\delta \right)\left(2kq\right)}t+\psi \right)𝑑t=0,$$
(21)
with $`\delta /\delta V`$ standing for the variational derivative of the underlying Lagrangian \[the one given by Eq. (12)\].

Let us make a couple of observations here. First, the above applies equally well to the full system as well as to the truncated system, since only variations in the SH field are involved, while the differences between the two types of systems involve solely the variation in $`U`$. Furthermore, since the solitons sought for are even, the expression $`\delta L/\delta V`$ with $`V`$ substituted by $`V_{\mathrm{sol}}(t)`$ is also even. Hence Eq. (21) amounts to an *orthogonality* condition between the infinitesimal tail and the soliton,

$$_{\mathrm{}}^+\mathrm{}\left[\left(\frac{\delta L}{\delta V}\right)|_{U=U_{\mathrm{sol}}(t),V=V_{\mathrm{sol}}(t)}\right]\mathrm{cos}\left(\sqrt{\left(2/\delta \right)\left(2kq\right)}t\right)𝑑t=0.$$
(22)
This consideration also shows that the phase parameter $`\psi `$ in the expression (18) is not important in the limit of $`b=0`$.

We notice that the above derivation circumvents the formal problem of the divergence of the integral expression (12) for the Lagrangian, when there is a tail which does not vanish as $`\left|t\right|\mathrm{}`$. The divergence did not appear since the tail was designed as a solution to the linearized version of Eqs. (9) or (11). One can readily verify that the divergence can also be eliminated in the fully nonlinear case, with respect to the tail (provided one suitably adjusts the frequency of the tail according to the nonlinearity).

Thus, Eq. (22) is a general criterion that can be used to locate ES solutions within the framework of VA. In Eq. (22), we recognize that the variational derivative $`\left(\delta L/\delta V\right)|_{U=U_{\mathrm{sol}},V=V_{\mathrm{sol}}(t)}`$ is just the left-hand side (l.h.s.) of the stationary equation for SH, with $`U`$ and $`V`$ taken as per the chosen ansatz for the core of the delocalized soliton. Another point is that upon applying integration by parts to the second-derivative term in Eq. (11), one sees that all contributions from the linear terms cancel in Eq. (22), leaving only the nonlinear terms in Eq. (11) to determine this condition.

### 4.2 Embedded soliton in the truncated system

To test the efficiency of the criterion (22), we first apply it to the truncated system. Then, the variational derivative in Eq. (22) is replaced by the nonlinear part of l.h.s. of Eq. (11),

$$_{\mathrm{}}^+\mathrm{}\left[\frac{1}{2}U^2(t)+4\gamma _2U^2(t)V(t)\right]\mathrm{cos}\left(\sqrt{\left(2/\delta \right)\left(2kq\right)}t\right)𝑑t=0.$$
(23)
Using the ansatz (13), it is easy to explicitly perform the integration in Eq. (23), which finally yields a simple result (note that the FH amplitude drops out),

$$4\gamma _2B=\frac{3\delta k}{2k\left(1+2\delta \right)q}.$$
(24)
Now, combining the above results (16) and (15), which were obtained by means of the VA, with the relation (24) that locates where the ES must be (and is also based on VA), one can easily verify that this set of three relations is *precisely tantamount* to the exact analytical solution for the ES which was found by guess in Ref. . In particular, an eventual expression for the wavenumber of ES is

$$k_{\mathrm{ES}}=\frac{1}{2}\left(1+2\delta \right)^1\left[q\frac{3}{2}\delta \left(4\gamma _2+3\delta \gamma _1\right)^1\right].$$
(25)
The fact that, in the case of the semi-Lagrangian truncated system, VA reproduces the *exact* ES is remarkable, although the reason for this occurrence is not fully understood. We also note that this exact soliton is not always an embedded one, as the wavenumber (25) is not necessarily restricted to the region (7) in which ES may exist: the wavenumber may instead fall into the region (6), in which ordinary solitons are to be found. A condition (inequality) showing whether the exact soliton is embedded or ordinary has already been given in Ref. .

### 4.3 Embedded soliton in the full system

The next step is to apply the general ES-selecting criterion (22) to the full system which is based on Eqs. (8) and (9). To this end, we assume that the soliton proper may again be approximated by the ansatz (13). Then, substituting the nonlinear part of l.h.s. of Eq. (9) into Eq. (22), we arrive, instead of Eq. (24), at a more complicated relation,

$$4\gamma _2B\left[1+\frac{B^2}{A^2}\frac{2\left(8\delta +1\right)k_{\mathrm{ES}}q}{40\delta k_{\mathrm{ES}}}\right]=\frac{3\delta k_{\mathrm{ES}}}{2k_{\mathrm{ES}}\left(1+2\delta \right)q}.$$
(26)
We test the validity of this relation in the following way: take a particular example of the ES in the full system that was found in a numerical form in Ref. , for which $`\delta =1`$, $`q=1`$, and $`\gamma _1=\gamma _2=0.05`$. Borrowing values of the amplitudes $`A`$ and $`B`$ directly from the numerical data, we find that $`A=3.794`$ and $`B=2.735`$. Substituting these values into Eq. (26) yields $`k_{\mathrm{ES}}=0.688`$, while the numerical value found in Ref.  was $`k_{\mathrm{ES}}=0.696`$. Thus,the relative error of the criterion (22) for this case is $`1.1\%`$.

## 5 Conclusion

In this work, we have put forward two modifications to the technique of the variational approximation for solitons. First, it may happen that a physical model does not admit the full Lagrangian representation, as some terms may be missing due to various reasons. In the case of the $`\chi ^{(2)}:\chi ^{(3)}`$ model considered in this work, this means that one term in the Lagrangian should not be varied when deriving the equation for the fundamental wave. We demonstrate that the VA can be applied to such  semi-Lagrangian systems as efficiently as to their full Lagrangian counterparts. Second, we have shown that, by the addition of an infinitesimal tail, which does not vanish at infinity, to the usual soliton ansatz, and demanding, after performing the variation, that the amplitude of the tail be zero, we obtain an approximate analytical criterion for locating embedded solitons inside a family of delocalized ones, i.e., isolated truly localized solutions existing inside the continuous spectrum of radiation modes. The criterion takes the form of orthogonality between the radiation mode contained in the infinitesimal tail and the core of the delocalized soliton. To test the criterion, we have applied it to both the semi-Lagrangian truncated version of the $`\chi ^{(2)}:\chi ^{(3)}`$ model and to the same model in its full form. In the former model, the criterion, combined with the VA for the soliton proper, yields a result which completely coincides with the previously found exact solution for the embedded soliton. In the latter model, the criterion predicts the wavenumber corresponding to the embedded soliton with a relative error $`1\%`$.

## Acknowledgement

B.A.M. appreciates hospitality of the Department of Mathematics at the University of Central Florida. This research has been supported in part by NSF grant # DMS0129714.

## Figure Captions

Fig. 1. A comparison between the numerical shape (solid curves) of an ordinary (non-embedded) soliton of the semi-Lagrangian (truncated) system of equations (10) and (11), as found by the shooting method, and the result (dashed curves) obtained from the modified variational approximation based on Eqs. (13), (15), and (16). The parameters are $`q=1`$, $`\delta =1`$, $`\gamma _1=0`$, and $`\gamma _2=1/4`$, and both the numerical and analytical solutions are taken for $`k=1/4`$.


text:
# Configurational entropy of network-forming materials

## Abstract

We present a computationally efficient method to calculate the configurational entropy of network-forming materials. The method requires only the atomic coordinates and bonds of a single well-relaxed configuration. This is in contrast to the multiple simulations that are required for other methods to determine entropy, such as thermodynamic integration. We use our method to obtain the configurational entropy of well-relaxed networks of amorphous silicon and vitreous silica. For these materials we find configurational entropies of $`1.02k_B`$ and $`0.97k_B`$ per silicon atom, respectively.

In materials such as vitreous silica, amorphous silicon or vitreous ice, the structure is determined by the set of bonds (covalent or hydrogen) between particles. In this manuscript, we will refer to these materials as network-forming materials. While the local environment of each particle is usually governed by strict rules, the bonded network can show a wide variety of different topologies. The focus of this work is to present a method to estimate the number of topologies that a network-forming material can take, or more precisely, its configurational entropy.

The common computational procedure to estimate the entropy $`S`$ at temperature $`T_2`$ is to measure the energy $`E`$ as a function of temperature and then to integrate from a temperature $`T_1`$ at which the entropy is known:

$$S(T_2)S(T_1)=_{T_1}^{T_2}\frac{1}{T}\frac{E}{T}𝑑T.$$
(1)
This requires sampling a large number of different networks and can therefore only be applied to systems with fast dynamics. In this work, we present a completely different approach to estimate the entropy, based on information theory shannon , and related to the work of Schlijper et al. who determined the entropy of the Ising and three-states Potts models schlijper .

One important concept in information theory is the Shannon entropy. It is commonly explained in the context of a string of $`n`$ bits. In this case, the Shannon entropy $`H(n)`$ is defined as:

$$H(n)=\underset{i}{}p(i)\mathrm{log}_2p(i),$$
(2)
where the index $`i`$ runs over all possible bit sequences of length $`n`$ and $`p(i)`$ is the probability of sequence $`i`$ occurring. A related concept is the entropy density $`s`$ of a large string of $`N`$ bits:

$$s=\underset{n\mathrm{}}{lim}\left[H(n+1)H(n)\right].$$
(3)
A practical procedure to estimate $`s`$ for such a long string is to extract from it a large number $`m`$ sub-sequences, each containing $`n`$ bits with $`nN`$. An estimate for the probabilities $`p(i)`$ is then given by:

$$p(i)\frac{f_i}{m},$$
(4)
where $`f_i`$ equals the number of times sub-sequence $`i`$ was observed. The estimates for $`p(i)`$ in combination with Eq. (2) yield $`H(n)`$. The entropy density $`s`$ is then obtained using Eq. (3). Usually, the limit $`n\mathrm{}`$ converges rapidly and even moderate values of $`n`$ are sufficient to predict $`s`$ accurately.

For systems in equilibrium it is easily shown that the Shannon entropy and the thermodynamic entropy of Eq. (1) are equivalent, apart from a factor of $`\mathrm{ln}(2)`$. In this case, the probabilities $`p(i)`$ are simply the Boltzmann weights: $`p(i)=\mathrm{exp}(\beta E_i)/Z`$, where $`\beta `$ is the inverse temperature, $`Z`$ the partition function and $`E_i`$ the energy of state $`i`$. In the present work we show how information theory can also be used to obtain the configurational entropy of network forming materials. The method we present can be applied to any network provided (1) the atomic coordinates are known and (2) a list of bonds is supplied or can be constructed — for instance based on a distance criterion — which uniquely determines the network.

To determine the configurational entropy of a network-forming material we choose a large number $`m`$ of random positions in the simulation cell. For each position we find the nearest $`n`$ particles and identify the graph formed by the bonds connecting these particles. We then assign a label to this graph, based on the graph automorphism nauty . We count the number of times a graph is observed and use Eq. (4) to estimate its corresponding probability of occurrence. These probabilities are fed into Eq. (2) to obtain $`H(n)`$.

A small displacement of one of these random positions will usually result in exactly the same list of $`n`$ nearest particles and thus the same graph. Consequently, there is an upper bound to the number $`m`$ of random positions that one should choose in a simulation cell containing $`N`$ particles. An estimate of this upper bound is obtained from the typical distance over which a random position can be displaced without altering the selected graph. From this we obtained as upper bounds $`m=1.6nN`$ and $`m=3.4n^2N`$ for two-dimensional and three-dimensional networks, respectively.

A related side-effect of choosing random positions is that the number of different graphs observed gets multiplied with a factor proportional to $`n^{d1}`$. This results in a correction to $`H(n)`$ of the form $`g(n)=(d1)\mathrm{log}(n)`$, with $`d`$ the spatial dimension of the network. The latter can be verified in crystalline networks where the configurational entropy is zero.

Graphs with a probability $`p(i)`$ smaller than $`1/m`$ will likely be observed only once, if at all. This finite-size effect grows with $`n`$, when the selected graphs become very complex. To monitor the impact of this effect, we record the quantity $`H_1(n)`$, defined as the contribution to $`H(n)`$ of the topologies observed once. We reject the measurements for which $`H_1(n)`$ exceeds one percent of $`H(n)`$. Furthermore, we observed that the quantity $`H(n)+cH_1(n)/H(n)`$ converges much faster with increasing $`m`$ than $`H(n)`$ itself, with a suitable choice of $`c`$. We therefore use this extrapolated value as our best estimate for $`H(n)`$.

In Fig. 1 we show the typical behavior of the corrected Shannon entropy $`H_c(n)H(n)g(n)`$ as a function of graphsize $`n`$. The entropy follows from Eq. (3). The rapid convergence of the limit is demonstrated by the linear behavior of $`H_c(n)`$ for intermediate $`n`$ shown in Fig. 1 by the dashed line. The entropy per atom equals the slope of this line. In Fig. 1 we have marked with circles the points for which $`H_1(n)`$ exceeds one percent of $`H(n)`$. These points suffer from finite-size effects and should not be used.

To test the usefulness of the above procedure, we apply it to the sillium model www , one of the prototype models to study network-forming materials. In this model, tuned for amorphous silicon, an explicit list of covalent bonds between pairs of Si atoms is kept, with the property that each Si atom is bonded to four neighboring atoms. The energy is described by the Keating potential, which contains a quadratic penalty for bond-length deviations from the crystalline distance of 2.35 Å, and a quadratic penalty for bond-angle deviations from the tetrahedral angle $`\mathrm{\Theta }_0=\mathrm{arccos}(1/3)`$. The list of bonds determines the atomic positions uniquely, since in this model the energy is minimized at all times. As a result, the phase space of this model is limited to a finite number of $`3N`$-dimensional points. The evolution of the network consists of a large number of random bond transpositions, each accepted with the Metropolis probability:

$$P=\mathrm{min}[1,\mathrm{exp}\left(\frac{E_bE_f}{k_BT}\right)],$$
(5)
where $`T`$ is the temperature and $`E_b`$ and $`E_f`$ are the total (minimized) energies of the system before and after the proposed bond transposition. A typical sillium network is shown in Fig. 2. Our networks are generated as described in Refs. barkemahq, ; vinkhq, , in which a number of algorithmic improvements have been proposed as compared to the original algorithm of Wooten, Winer and Weaire.

To verify the validity of the information theory approach we first turn to a two-dimensional version of the sillium model. Contrary to the three-dimensional sillium model, the dynamics here does not lead to glassy states. This allows us to determine the entropy in equilibrium using both the standard approach — via Eq. (1) — as well as information theory. In this two-dimensional model, atoms are three-fold coordinated and the ideal bond angle is 120 degrees; a typical configuration is shown in Fig. 3.

We simulate a two-dimensional sillium network containing 1008 atoms for a number of (physically interesting) temperatures ranging from $`k_BT=0.15`$ eV to $`k_BT=0.60`$ eV using periodic boundary conditions in both directions. At each temperature we bring the system to equilibrium with 100 attempted bond-transpositions per atom. Next, 20 snapshots of the network are stored; each separated by five attempted bond transpositions per atom. From these snapshots the average energy and the Shannon entropy are obtained. In determining the Shannon entropy we draw clusters from each snapshot simultaneously. After the simulation has covered the entire temperature range, the entropy is also obtained from the average energy measurements and Eq. (1).

In Fig. 4 we show the configurational entropy as a function of temperature. The solid curve shows the entropy per atom as obtained using the standard thermodynamic approach; crosses show the entropy per atom as obtained using information theory. There appears to be a phase transition near $`k_BT=0.37`$ eV shown by the sudden jump in entropy. The remarkable feature of the displayed results is the overall good agreement between the two methods, even near the phase transition.

Next, we use information theory to determine the configurational entropy of a well-relaxed three-dimensional sillium network containing 20,000 atoms. This network was generated using the improved Wooten, Winer and Weaire algorithm barkemahq ; vinkhq . The Keating energy of this network is 0.286 eV per atom; the standard deviation in the mean bond-angle is 9.63 degrees. Structural and electronic properties of this network are in excellent agreement with experimental properties of amorphous silicon vinkhq . Applying our method to this model, we obtain a configurational entropy of $`1.02k_B`$ per atom.

The (classical) vibrational entropy is obtained from the eigenvalues of the dynamical matrix and for this network was found to be $`3.07k_B`$ per atom. For the crystalline phase the configurational entropy is zero and only the vibrational entropy contributes to the entropy, in this case $`3.18k_B`$ per atom. The difference in entropy between the crystalline and amorphous phases of silicon is therefore estimated to be $`0.91k_B`$ per atom.

The energy difference between the crystalline and amorphous phases of silicon was recently determined by Biswas biswas  using tight-binding, who found an energy difference of 0.18 eV per atom.

With the above estimates for the differences in energy and entropy we can calculate the transition temperature $`T_c`$, defined as the temperature where the free energy difference between the crystalline phase and the amorphous phase changes sign. It is given by:

$$T_c=\frac{\mathrm{\Delta }E}{\mathrm{\Delta }S},$$
(6)
where $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }S`$ are the energy difference and the entropy difference between the crystalline and the amorphous phase, respectively. Substitution of our estimates for $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }S`$ into Eq. (6) yields $`T_c2300`$K. This temperature compares well with the value of $`T_c2500`$K as inferred from calorimetric experiments sinke , and is well above the melting point of silicon. It confirms that the amorphous phase is not thermodynamically stable at any temperature.

As a final application of our method we determine the configurational entropy of vitreous silica. The structure of this material is formed by covalent bonds between silicon atoms and oxygen atoms. Barring rare defects, each silicon atom is bonded to four oxygen atoms and each oxygen atom is bonded to two silicon atoms. We generate a silica network containing 3000 atoms (with periodic boundary conditions) in a spirit similar to amorphous silicon, but we replace the Keating potential by the one of Tu and Tersoff tu1 ; tu2 . The resulting network is then quenched with the more realistic BKS potential beest , using parameters as described in Ref. vollmayr, . The BKS energy per silicon atom of the quenched network is 0.13 eV above that of the $`\beta `$-cristobalite structure. For comparison, samples prepared by molecular dynamics typically yield a much larger energy difference of 0.30 eV. This clearly demonstrates that our network is well relaxed. After quenching with the BKS potential, the atoms remain perfectly coordinated: the silicon-oxygen bonds are easily reconstructed from a distance criterion and the silicon-oxygen radial distribution function. For our network, the average O-Si-O bond-angle is 109.43 degrees with a standard deviation of 4.40 degrees; the average Si-O-Si bond-angle was found to be 150.85 degrees with a standard deviation of 12.09 degrees.

To estimate the configurational entropy per silicon atom in well-relaxed vitreous silica, we used the above network and replaced each oxygen atom plus its two bonds by a single silicon-silicon bond. Next, we applied the information-theoretic method described earlier on this network and found for the configurational entropy $`0.97k_B`$ per silicon atom. Given the limited size of this sample we estimate that the actual entropy might be up to 10% larger. Ignoring the difference in vibrational entropy, as well as the entropy contribution of possible multiple oxygen positions, we obtain for the transition temperature of vitreous silica $`T_c1600`$K.

In summary, we have developed a computationally efficient method to determine the configurational entropy of network forming materials. For well-relaxed samples of amorphous silicon and vitreous silica, we find for the entropy per silicon atom $`1.02k_B`$ and $`0.97k_B`$, respectively. In future research, we hope to extend the applicability of this method to other disordered materials such as colloidal systems and metallic glasses.

We thank Henk van Beijeren and Normand Mousseau for stimulating discussion, and Normand Mousseau also for the computation of the BKS energies.


text:
# Role of spin-orbit coupling in the far infrared absorption of lateral semiconductor dots

## Abstract

We investigate the relevance of the spin-orbit coupling to the far-infrared absorption of two-dimensional semiconductor dots. Varying the strength of the Dresselhaus term, a mechanism feasible in experiment by changing the dot width, distinctive splittings of the Kohn peak as well as additional low energy modes are predicted in a non-interacting model. Each mode has a spatial distribution of charge, perpendicular and in-plane spin densities that correlate in a peculiar way with the frequency and polarization of the external field. We study the robustness of these features against electron-electron interactions as well as the appearance of interaction-induced additional characteristics.

The role of the electron spin-orbit (SO) coupling in the general properties of lateral semiconductor dots has attracted much interest in recent years. To a large extent this is motivated by the technological impact attributed to the so-called spintronic devices, whose operation principle exploits precisely the SO mechanism to control the spin transport . Actually, spintronics is an incipient field whose fundamental principles need a deeper understanding . In this general context, experimental and theoretical investigations have very recently shown that the SO coupling affects the quantum transport and, more specifically, the conductance fluctuations of chaotic quantum dots in a parallel magnetic field .

The different SO contributions in semiconductor dots have been reviewed by Voskoboynikov et al. . In general, SO coupling has two contributions: the Rashba term due to the nanostructure confinement potential; and the Dresselhaus term originating from the crystalline inversion asymmetry. In lateral dots the strength of the latter depends on the vertical width, which can be experimentally varied with vertical electric fields. For an embedding GaAs medium both Rashba and Dresselhaus parameters are reasonably determined .

In this Letter we study the influence of SO couplings on the dipole absorption of GaAs dots with a parabolic confinement. While the Rashba term provides a negligible effect we show that the Dresselhaus contribution leads to characteristic splittings of the Kohn modes and to the appearance of quasi-spin flip modes in the far-infrared absorption (FIR). For each mode we find a correlation between the spatial patterns of oscillating electron density, horizontal (in-plane) spin density, and vertical spin density. As in most of the literature treating the SO coupling , we shall first consider a non-interacting model, using both analytical and numerical treatments. Nevertheless, the role of the interaction will be also considered within self-consistent mean-field for Hartree and exchange-correlation contributions, using a spinorial formalism of local-spin-density theory (LSDA) .

Within the effective Hamiltonian formalism and for a quasi two-dimensional system the Rashba and Dresselhaus \[for the standard (001) plane of GaAs\] contributions read, respectively,

$`_R`$ $`=`$ $`{\displaystyle \frac{\lambda _R}{\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left[\stackrel{}{\sigma }\left(V(𝐫)\times 𝐏\right)\right]_i`$ (1)
$`_D`$ $`=`$ $`{\displaystyle \frac{\lambda _D}{\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left[P_x\sigma _xP_y\sigma _y\right]_i,`$ (2)
where the $`\lambda `$’s control the coupling strengths, the $`\sigma `$’s refer to the Pauli matrices and $`V(𝐫)`$ stands for the total electrostatic potential. In Eq. (1), $`𝐏`$ gives the kinetic momentum. Although we shall not consider in this work an external magnetic field (or vector potential $`𝐀`$) it can be easily included using the canonical momentum $`𝐏=i\mathrm{}+\frac{e}{c}𝐀`$ in Eq. (1).

For GaAs one has $`\lambda _R=5.33`$ Å2 and a material-specific constant $`\gamma =27.5`$ eVÅ3 . The latter determines the Dresselhaus coefficient in terms of the dot vertical width $`z_0`$ as  $`\lambda _D\gamma (\pi /z_0)^2`$ . As a representative case, we consider an external parabolic confinement with energy $`\mathrm{}\omega _0=4.2`$ meV and a number of electrons ranging from 6 to 12. A strict two-dimensional motion is assumed and the dot width is simply taken as a parameter for the Dresselhaus coupling. The Rashba term turns out to be negligible in all cases we have studied. However, the Dresselhaus term increases its importance for decreasing doth widths, leading to relevant effects as we shall show below already for widths $`z_050`$ Å. The influence of the Rashba term has already been studied in the literature for InSb dots , where band non-parabolicity is important, and in phenomenologically fitted models .

We have implemented the numerical solution on a spatial grid of the full Hamiltonian $`_{\mathrm{f}\mathrm{u}\mathrm{l}\mathrm{l}}=_0+_R+_D`$, with

$$_0=\underset{i=1}{\overset{N}{}}\left[\frac{𝐏^2}{2m}+\frac{1}{2}m\omega _0^2r^2\right]_i.$$
(3)
However, before presenting numerical results it is worth to obtain an approximate analytical solution of the relevant contributions $`_0+_D`$ with an approach similar to that of Aleiner and Fal’ko . A diagonalization to second order in $`\lambda _D`$ is obtained by means of the transformation $`\stackrel{~}{}=U^+U`$, where

$$U=\mathrm{exp}\left[\text{}i\lambda _D\frac{m}{\mathrm{}^2}\underset{j=1}{\overset{N}{}}(x\sigma _xy\sigma _y)_j\text{}\right].$$
(4)
We find

$`\stackrel{~}{}`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}\left[\text{}{\displaystyle \frac{𝐏^2}{2m}}+{\displaystyle \frac{1}{2}}m\omega _0^2r^2+\lambda _D^2{\displaystyle \frac{m}{\mathrm{}^3}}(xP_yyP_x)\sigma _z\text{}\right]_j`$ (5)
$``$ $`N\lambda _D^2{\displaystyle \frac{m}{\mathrm{}^2}}+O(\lambda _D^3).`$ (6)
Identifying the orbital angular momentum $`\mathrm{}_z=xP_yyP_x`$ in (5), we realize that the Dresselhaus term introduces a correction depending on $`\mathrm{}_z\sigma _z`$ quite similar to that of the Rashba term for a circular system, but in the transformed reference frame. In fact, in the laboratory frame the $``$ eigenstates to first order in $`\lambda _D`$ are spinors that slightly deviate from the $`\sigma _z`$ eigenspinors. We shall call them quasi-up $`\chi _{qu}`$ and quasi-down $`\chi _{qd}`$ states

$`\chi _{qu}`$ $``$ $`\phi _n\mathrm{}(𝐫)\left({\displaystyle \genfrac{}{}{0pt}{}{1}{\lambda _D\frac{m}{\mathrm{}^2}(y+ix)}}\right),`$ (7)
$`\chi _{qd}`$ $``$ $`\phi _n\mathrm{}(𝐫)\left({\displaystyle \genfrac{}{}{0pt}{}{\lambda _D\frac{m}{\mathrm{}^2}(yix)}{1}}\right),`$ (8)
where $`\phi _n\mathrm{}(𝐫)`$ are the eigenfunctions of the circular harmonic oscillator $`_0`$. Since the dipole operator is invariant by the transformation (4), the analytical solution predicts dipole transitions with $`\mathrm{\Delta }\mathrm{}=\pm 1`$ for quasi-up and quasi-down subsystems (see Fig. 1) at the energies

$`\mathrm{}\omega _{qu,\pm 1}`$ $`=`$ $`\mathrm{}\omega _0\pm \lambda _D^2{\displaystyle \frac{m}{\mathrm{}^2}}`$ (9)
$`\mathrm{}\omega _{qd,\pm 1}`$ $`=`$ $`\mathrm{}\omega _0\lambda _D^2{\displaystyle \frac{m}{\mathrm{}^2}}.`$ (10)
An inmediate consequence is the splitting of the original Kohn mode at $`\mathrm{}\omega _0`$ in two peaks separated by an energy $`\mathrm{\Delta }\mathrm{}\omega =2\lambda _D^2m/\mathrm{}^2`$.

Left panels of Fig. 2 show the numerical result for the FIR absorption of $`_{\mathrm{f}\mathrm{u}\mathrm{l}\mathrm{l}}`$ obtained within a real-time simulation after an initial rigid displacement of the system as in Ref. . The results correspond to $`N=6`$ and 10 electrons with a Dresselhaus parameter for a dot width of $`z_0=50`$ Å. The splitting around $`\mathrm{}\omega _0`$ nicely agrees with the analytical prediction (arrows), thus proving the validity of the model. For the $`N=10`$ case we note the existence at low energy of a peak with much less strength than the dominant ones (notice the enlarged scale). This too can be explained within the analytical model realizing that $`N=10`$ has a partial occupation of the third oscillator major shell and thus the possibility of quasi-spin flip transitions $`\chi _{qu}\chi _{qd}`$ at an energy $`4\lambda _D^2m/\mathrm{}^2`$, indicated with dotted lines in Fig. 1 and with an arrow in Fig. 2. Note that in spite of having $`\mathrm{\Delta }\mathrm{}_z=0`$ these transitions can manifest in the dipole spectrum through the higher order corrections in Eq. (5).

Due to the spin-orbit coupling the charge dipole oscillation also induces a spin oscillation. This is proved in the right panels of Fig. 2 where the Fourier transform of $`S_x+S_y`$ is plotted. The same frequencies of the dipole absorption manifest in the spin channels, although with different relative strengths. It is also interesting to analyze the spatial patterns of induced density $`\delta \rho (𝐫,t)`$, parallel spin $`\delta 𝐒_{}(𝐫,t)`$ and vertical spin $`\delta S_z(𝐫,t)`$ for a given oscillation mode and a fixed time (Fig. 3), where we define the spin density as $`𝐒(𝐫)=_{i=1}^N\chi _i|\stackrel{}{\sigma }\delta (𝐫𝐫_i)|\chi _i`$. This local-signal analysis has been performed using the method of Ref. . Quite remarkably, the oscillation of $`\delta S_z(𝐫,t)`$ is always in a perpendicular direction to that of the charge density, which is explained in the analytical model as a coherent quasi-up and quasi-down excitation due to the SO term. On the contrary, the parallel spin patterns depend on whether the mode has quasi-spin flip character or not. For the split Kohn modes (peaks a and b) $`\delta 𝐒_{}(𝐫,t)`$ is localized in the two regions perpendicular to the charge oscillation axis. However, the quasi-spin flip peak (c) shows no angular dependence of the pattern. It can be shown that this behaviour is in complete agreement with the prediction of the analytical model using the spinors of Eq. (7).

We consider next how the features discussed above are modified by the addition of the Coulomb interaction

$$_C=\underset{i<j}{\overset{N}{}}\frac{e^2}{\kappa |𝐫_i𝐫_j|},$$
(11)
where $`\kappa `$ is the semiconductor dielectric constant ($`\kappa =12.4`$ for GaAs). We shall estimate interaction effects using a selfconsistent mean-field for the Hartree electrostatic potential

$$V_H(𝐫)=\frac{e^2}{\kappa }𝑑𝐫^{}\frac{\rho (𝐫^{})}{|𝐫^{}𝐫|},$$
(12)
and the LSDA for electronic exchange and correlation. Within the present spinor formalism the dots have a vectorial spin magnetization $`𝐒(𝐫)`$ and, therefore, a general exchange-correlation functional of the type $`E_{XC}[\rho ,𝐒]`$ is necessary. For this we have resorted to the non-collinear-spin formulation of LSDA theory, described for instance in Ref. , and employed the Monte-Carlo results of Ref.  for the bulk electron gas.

The modifications stemming from the Coulomb interaction are twofold; first, the electron-hole transitions move in energy due to the change in the effective mean field and; second, new collective peaks may appear as a result of coherent electron-hole excitations. Figure 4 shows the interacting results for $`N=6`$ with $`z_0=50`$ and 35 Å. Note that the splitting of the Kohn-like modes disappears for the larger width and, at the same time, a collective peak begins to appear at an energy below the electron-hole states. These results clearly show the formation of the collective plasmon and magnon states lying above and below the electron-hole transitions, respectively. We stress that in the absence of SO coupling the magnon state does not contribute to the FIR absorption. The lower left panel of Fig. 4 indicates that a great enhancement of the magnon contribution is obtained by reducing the dot width, which is accompanied by an important Landau damping of the plasmon state into a bundle of peaks. The Coulomb interaction thus suppresses the splitting of the Kohn-like modes although it introduces the mechanism of Landau damping that ultimately, and rather abruptly with $`z_0`$, dominates for decreasing widths.

The right panels of Fig. 4 show the horizontal spin spectra induced by the dipole shift. These results permit to ascertain the spin character of the different modes and, by comparison with the left panels, quantify their relevance to the FIR absorption. The relative spin or density character of each state is also seen from the amplitudes $`𝒜`$ of density and spin oscillation. For instance, when $`z_0=50`$ Å we find $`𝒜[\delta S_z(𝐫,t)]60𝒜[\delta \rho (𝐫,t)]`$ for the magnon and $`𝒜[\delta S_z(𝐫,t)]0.4𝒜[\delta \rho (𝐫,t)]`$ for the plasmon, which shows the dominant spin and density type for magnon and plasmon states, respectively.

The characterization of the interacting peaks is facilitated by the comparison of the spatial patterns with the corresponding ones in the non-interacting model. Figure 5 displays the density and spin patterns for two of the peaks of the $`z_0=35`$ Å spectrum (lower left panel of Fig. 4). As in the non-interacting model, the oscillation of $`\delta S_z(𝐫,t)`$ is orthogonal to that of $`\delta \rho (𝐫,t)`$ and the horizontal spin patterns show a priviledged spin direction, i.e., they have a net induced spin. Note in particular that peak b of Fig. 4 has no angular dependence in the horizontal spin pattern while for peak a there are two regions with enhanced $`\delta 𝐒_{}(𝐫,t)`$. We have checked that these features are not changed when the polarization of the mode, given by the direction of the density oscillation, is varied. We stress that the mechanism discussed here would permit a control of the horizontal spin oscillation by tuning the frequency of the applied electric field. This effect of spin-to-charge conversion is a major requisite for hybrid spintronic devices .

Finally, Fig. 6 presents the spectra with Coulomb interaction obtained for $`N=10`$, an open-shell system in the analytical model. At $`z_0=50`$ Å the open-shell character of this dot is clearly maintained, with low-energy electron hole transitions and an associated collective mode at $`\omega 0.7`$ meV similar to the quasi-spin flip excitations of Fig. 1. For narrower dots, however, the energy of quasi-spin flip transitions increases and as a result this mode merges with the magnon and other electron-hole states. Contrary to the quasi-spin flip state, both magnon and plasmon mean energies are not dependent on the Dresselhaus coupling.

To summarize, the role of the SO coupling in the FIR absorption of GaAs dots has been studied. In an analytical model without interaction the relevant Dresselhaus contribution induces characteristic splittings and quasi-spin flip modes. The Coulomb interaction favors the formation of collective plasmon and magnon states and it suppresses the splittings of the analytical model. Nevertheless, below a certain dot width Landau damping dominates the spectrum. Each mode has characteristic patterns of induced charge, horizontal spin and vertical spin, as well as a ratio of oscillating density and spin amplitudes. Therefore, by tuning the frequency of the applied electric field one can effectively control the dynamical spin properties and the relative intensities of the charge and spin oscillation.

This work was supported by Grant No. PB98-0124 from DGESeIC, Spain.