Appendix/Ramblings/NumericallyDeterminedEigenvectors

From JETohlineWiki
Jump to navigation Jump to search

Numerically Determined Eigenvectors of a Zero-Zero Bipolytrope[edit]

Here we build on the analytic foundation summarized in an accompanying chapter and attempt to numerically construct a variety of eigenvectors that describe radial oscillations of bipolytropes for which, (nc,ne)=(0,0).

Setup[edit]

We'll begin with the linear-adiabatic wave equations that describe oscillations of the core and envelope, separately. We also will immediately restrict our investigation to configurations for which,

g2=                 g2=1+8q3(1+2q3)2,         and,         q3=𝒟=[ρe/ρc2(1ρe/ρc)]                 ρeρc=2q31+2q3.

For the core we have,

0

=

(1η2)d2xdη2+(46η2)1ηdxdη+𝔉corex,

where,

ηξg,         and         𝔉core3ωcore22πGγcρc2αc.

And, for the envelope we have,

0

=

(1q3ξ3)d2xdξ2+(36q3ξ3)1ξdxdξ+[q3𝔉envξ3αe]xξ2,

where,

𝔉env

3ωenv22πGγeρe2αe.


Initial Focus[edit]

Properties of 21Analytic Solution[edit]

Figure 1

reference analytic solution

Same analytic, (,j)=(2,1) eigenfunction as here, but renormalized to unity at the center; numerical label provides value of σc2.

Evidently, one analytic solution with quantum numbers, (,j)=(2,1), shown again here on the right, is available for a zero-zero bipolytrope that has the following properties:

q

0.6840119

ρeρc=2q31+2q3

0.3902664

γe=43+0.35

1.1940299

γc

1.845579

σc23ω22πGρc=20γc8

28.91158.

This means, as well, that,

c01+αe1

=

0.6510.1937742

g21+8q3(1+2q3)2

1.3236092

𝔉coreσc2+8γc6

=

14

𝔉env1γe[σc2(ρcρe)+8]6

=

(c02+17c0+66)=62.743385


In the envelope, the analytically defined eigenfunction is given by the expression,

x=2|env

=

ξc0[1+q3A21ξ3+q6A21B21ξ61+q3A21+q6A21B21],

where,

A21

(4c0+222c0+5)4.6016533,

B21

(c0+72c0+8)0.8940912;

and in the core, it is,

xj=1|core

=

5(1+8q3)7(1+2q3)2ξ25(1+8q3)7(1+2q3)2.

More succinctly we have,

xcore

=

17.326820+18.326820ξ2;

and,

axenv

=

ξ0.1937742[11.4726681ξ3+0.4213836ξ6],

where,

a[1+q3A21+q6A21B21]0.05128445.

Demonstrate Core Solution[edit]

This means that,

dxcoredξ

=

36.65364ξ,

and,

d2xcoredξ2

=

36.65364.

Therefore, the LAWE for the core becomes,

[LAWE]core

=

(1η2)d2xcoredη2+(46η2)1ηdxcoredη+𝔉corexcore

 

=

(g2ξ2)d2xcoredξ2+(4g26ξ2)1ξdxcoredξ+𝔉corexcore

 

=

36.65364(1.3236092ξ2)+36.65364(5.29443686ξ2)+14(17.326820+18.326820ξ2)

 

=

36.65364(1.3236092)+36.65364(5.2944368)+14(17.326820)+[36.65364(1)+36.65364(6)+14(18.326820)]ξ2

 

=

36.65364(6.618046)14(17.326820)+36.65364[16+7]ξ2

 

=

0.

Q.E.D.


Demonstrate Envelope Solution[edit]

Given that,

axenv

=

ξ0.1937742[11.4726681ξ3+0.4213836ξ6],

we deduce that,

adxenvdξ

=

0.1937742ξ1.1937742[11.4726681ξ3+0.4213836ξ6]+ξ0.1937742[4.4180043ξ2+2.5283016ξ5],

and,

ad2xenvdξ2

=

0.2313226ξ2.1937742[11.4726681ξ3+0.4213836ξ6]2×0.1937742ξ1.1937742[4.4180043ξ2+2.5283016ξ5]

 

 

+ξ0.1937742[8.8360086ξ+12.641508ξ4].

Therefore, the LAWE for the envelope becomes,

a[LAWE]env

=

a(1q3ξ3)d2xdξ2+a(36q3ξ3)1ξdxdξ+a[q3𝔉envξ3αe]xξ2,

 

=

a{d2xdξ2+3ξdxdξαe(xξ2)}aq3ξ3{d2xdξ2+6ξdxdξ𝔉env(xξ2)}.

Now, the first of these sub-expressions gives,

a{d2xdξ2+3ξdxdξαe(xξ2)}

=

0.2313226ξ2.1937742[11.4726681ξ3+0.4213836ξ6]2×0.1937742ξ1.1937742[4.4180043ξ2+2.5283016ξ5]

 

 

+ξ0.1937742[8.8360086ξ+12.641508ξ4]

 

 

0.5813226ξ2.1937742[11.4726681ξ3+0.4213836ξ6]+3ξ1.1937742[4.4180043ξ2+2.5283016ξ5]

 

 

+0.35ξ2.1937742[11.4726681ξ3+0.4213836ξ6]

 

=

(0.23132260.5813226+0.35)ξ2.1937742[11.4726681ξ3+0.4213836ξ6]

 

 

+(32×0.1937742)ξ1.1937742[4.4180043ξ2+2.5283016ξ5]+ξ0.1937742[8.8360086ξ+12.641508ξ4]

 

=

(32×0.1937742)ξ2.1937742[4.4180043ξ3+2.5283016ξ6]+ξ2.1937742[8.8360086ξ3+12.641508ξ6]

 

=

ξ2.1937742[20.37783ξ3+19.24657ξ6]

And the sub-expression inside the second set of curly braces gives,

a{d2xdξ2+6ξdxdξ𝔉env(xξ2)}

=

0.2313226ξ2.1937742[11.4726681ξ3+0.4213836ξ6]2×0.1937742ξ1.1937742[4.4180043ξ2+2.5283016ξ5]

 

 

+ξ0.1937742[8.8360086ξ+12.641508ξ4]

 

 

1.1626452ξ2.1937742[11.4726681ξ3+0.4213836ξ6]+6ξ1.1937742[4.4180043ξ2+2.5283016ξ5]

 

 

62.74339ξ2.1937742[11.4726681ξ3+0.4213836ξ6]

 

=

(0.23132261.162645262.74339)ξ2.1937742[11.4726681ξ3+0.4213836ξ6]

 

 

+(62×0.1937742)ξ1.1937742[4.4180043ξ2+2.5283016ξ5]+ξ0.1937742[8.8360086ξ+12.641508ξ4]

 

=

63.67471ξ2.1937742[11.4726681ξ3+0.4213836ξ6]

 

 

+5.612452ξ2.1937742[4.4180043ξ3+2.5283016ξ6]+ξ2.1937742[8.8360086ξ3+12.641508ξ6]

 

=

ξ2.1937742[63.67471+93.77171ξ326.83148ξ6]

 

 

+ξ2.1937742[24.79584ξ3+14.18997ξ6]+ξ2.1937742[8.8360086ξ3+12.641508ξ6]

 

=

ξ2.1937742[63.67471+(93.7717124.795848.8360086)ξ3+(26.83148+14.18997+12.641508)ξ6]

 

=

ξ2.1937742[63.67471+60.13986ξ3]

a(q3ξ3){d2xdξ2+6ξdxdξ𝔉env(xξ2)}

=

ξ2.1937742[20.37783+19.24657ξ6].

But these two sub-expressions cancel precisely, which means that our eigenfunction satisfies the LAWE! Q.E.D.

Boundary Conditions[edit]

Notice that for this particular eigenfunction solution, the value and first radial derivative at the center (ξ=0) of the configuration is,

xcore

=

17.326820+18.326820ξ20=17.326820;

and,

dxcoredξ

=

36.65364ξ0=0.

And, at the surface (ξ=q1) the value and first radial derivative are,

axenv

=

{ξ0.1937742[11.4726681ξ3+0.4213836ξ6]}ξ=1/q

 

0.47627246,

where,

a0.05128445;

and,

dlnxenvdlnξ

=

ξaxenv[adxenvdξ]

 

=

0.1937742ξ0.1937742[11.4726681ξ3+0.4213836ξ6]+ξ0.1937742[4.4180043ξ3+2.5283016ξ6]ξ0.1937742[11.4726681ξ3+0.4213836ξ6]

 

=

0.1937742+[4.4180043ξ3+2.5283016ξ6][11.4726681ξ3+0.4213836ξ6]

{dlnxenvdlnξ}ξ=1/q

=

0.1937742+{[4.4180043ξ3+2.5283016ξ6][11.4726681ξ3+0.4213836ξ6]}ξ=1/q

 

=

0.1937742+21.22492=21.03115.

Finite-Difference Representation[edit]

General Approach[edit]

Working with the Taylor series expansion, we can write,

x(ξ)

x(a)+(ξa)xa+12(ξa)2xa,

and letting ξ±=a±Δ, we have,

x+

x(a)+Δxa+12Δ2xa,

and,

x

x(a)Δxa+12Δ2xa.

Subtracting the second of these two expressions from the first gives,

x+x

2Δxa

xa

x+x2Δ;

while, adding the two expressions together gives,

x+2xa+xΔ2

xa.

Integrating Outward Through the Core[edit]

From the LAWE for the core, we have,

a(g2a2)xa

=

(4g26a2)xaa𝔉corexa.

So, putting these last three expressions together gives an approximate relation between x+ and the previous two values of the function, x and xa, namely,

a(g2a2)[x+2xa+xΔ2]

(4g26a2)[x+x2Δ]a𝔉corexa

a(g2a2)[x+Δ2]+(4g26a2)[x+2Δ]

a(g2a2)[2xaxΔ2]+(4g26a2)[x2Δ]a𝔉corexa

x+[2a(g2a2)+Δ(4g26a2)]

2a(g2a2)[2xax]+(4g26a2)[Δx]2Δ2a𝔉corexa

 

[4a(g2a2)2Δ2a𝔉core]xa+[Δ(4g26a2)2a(g2a2)]x.

Now, at the very center of the configuration, (a=0), we expect the function, x(ξ), to be symmetric; that is, we expect x=x+. So for this case alone, we have,

x+[2a(g2a2)+Δ(4g26a2)Δ(4g26a2)+2a(g2a2)]

=

[4a(g2a2)2Δ2a𝔉core]xa

x+[2(g2a20)+2(g2a20)]

=

[4(g2a20)2Δ2𝔉core]xa

x+

=

[1Δ2𝔉core2g2]xa.

For all other coordinate locations, a=ξ, in the range, 0<ξ<1, we will use the general expression, namely,

x+

[4a(g2a2)2Δ2a𝔉core]xa+[Δ(4g26a2)2a(g2a2)]x[2a(g2a2)+Δ(4g26a2)].

Keep in mind that, when we move across the interface at a=1, we want both the value of the function, xq, and its first derivative, xq, to be the same as viewed from both the envelope and the core. In a numerical integration algorithm, it will be very straightforward to set the value of the eigenfunction at the interface. In order to properly handle the first derivative, I propose that we extend the core solution and evaluate the eigenfunction at one zone beyond the interface, and identify the values of the eigenfunction that straddles the interface as,

(x)q       and       (x+)q.

Then define the slope of the eigenfunction at the interface by the expression,

Slope at the Interface

xq

(x+)q(x)q2Δ.


Integrating Outward Through the Envelope[edit]

From the LAWE for the envelope, we have,

a2(1q3a3)xa

=

(36q3a3)axa[q3𝔉enva3αe]xa.

Inserting the same finite-difference expressions for the first and second derivatives, we therefore have,

a2(1q3a3)[x+2xa+xΔ2]

=

(36q3a3)a[x+x2Δ][q3𝔉enva3αe]xa

a2(1q3a3)[x+Δ2]+(36q3a3)a[x+2Δ]

=

(36q3a3)a[x2Δ]a2(1q3a3)[x2xaΔ2][q3𝔉enva3αe]xa

x+[2a2(1q3a3)+Δ(36q3a3)a]

=

[Δ(36q3a3)a2a2(1q3a3)]x+[4a2(1q3a3)2Δ2(q3𝔉enva3αe)]xa

Now, at the interface (only), we need to relate x to x+ in such a way that the slope gives the proper value at the interface. Specifically, we need to set,

x

=

x+2Δ(xq),

where, xq takes the value that was determined for the core. Hence, at the interface (a=1), the first step into the envelope is special and demands that,

x+[2a2(1q3a3)+Δ(36q3a3)a]

=

[Δ(36q3a3)a2a2(1q3a3)][x+2Δ(xq)]+[4a2(1q3a3)2Δ2(q3𝔉enva3αe)]xa

x+[4a2(1q3a3)]

=

[Δ(36q3a3)a2a2(1q3a3)][2Δ(xq)]+[4a2(1q3a3)2Δ2(q3𝔉enva3αe)]xa

 

=

2Δ[Δ(36q3a3)a2a2(1q3a3)]xq+[4a2(1q3a3)2Δ2(q3𝔉enva3αe)]xa

and, setting,       a=1x+

=

2Δ[2(1q3)Δ(36q3)]xq+[4(1q3)2Δ2(q3𝔉envαe)]xa4(1q3).

Varying the Oscillation Frequency[edit]

Approach[edit]

First, we fix q, γe, and γc; in the example, here (as above) we choose: (q,γe,γc)=(0.6840119,1.1940299,1.845579). For this example, we will also retain the constraint, g2=, in which case,

ρeρc=2q31+2q3

0.3902664.

Next, we pick various values of the (square of the) dimensionless oscillation frequency, σc2, and, in order to ensure that the dimensional frequency in the envelope matches the dimensional frequency of the core, from each value we set,

𝔉core

=

σc2+8γc6,

𝔉env

=

1γe[σc2(ρeρc)1+8]6.

For the finite-difference algorithm, we divide the core — radial coordinate range, 0ξ1 — into Ncore zones, and the envelope — radial coordinate range, 1ξ1/q — into Nenv zones. This means that the spacing between successive radial zones in the core and envelope is, respectively,

Δc1Ncore

      and      

Δeq11Nenv.

Starting at the center of the configuration (ξ=0), where we arbitrarily set the value of the eigenfunction to x0=1, the value of the eigenfunction at the first grid point away from the center (ξ=Δc) is,

xk=1

=

[1Δc2𝔉core2g2]x0.

Thereafter — moving out toward and just beyond the interface location (ξ=1), the radial coordinate of each successive grid point is ξk=kΔc, and the numerically determined value of the eigenfunction at each successive grid point (k=1Ncore) is,

xk+1

[4ξk(g2ξk2)2Δc2ξk𝔉core]xk+[Δc(4g26ξk2)2ξk(g2ξk2)]xk1[2ξk(g2ξk2)+Δc(4g26ξk2)].

Then, at the interface, which is associated with k=Ncore, we define the reference slope as,

xq

=

xk+1xk12Δc.

Next, we move outward into the envelope, using the integer index, n=1Nenv, to label successive radial grid locations (ξn=1+nΔe). Letting the value of the eigenfunction at the interface be represented by xq, at the first grid location outside the interface (ξ=1+Δe), the value of the eigenfunction is,

xn=1

=

2Δe[2(1q3)Δe(36q3)]xq+[4(1q3)2Δe2(q3𝔉envαe)]xq4(1q3).

Thereafter, moving outward through the envelope to the surface, the value of the eigenfunction at each successive grid location is,

xn+1

=

[Δe(36q3ξn3)ξn2ξn2(1q3ξn3)]xn1+[4ξn2(1q3ξn3)2Δe2(q3𝔉envξn3αe)]xn[2ξn2(1q3ξn3)+Δe(36q3ξn3)ξn].

TEST:   We tested this finite-difference algorithm on a grid of resolution, Ncore=Ncore=50, by first setting σc2=28.91158. The resulting, numerically constructed eigenfunction matched to high accuracy the analytically generated eigenfunction shown, above, as Figure 1; see also, the middle image in the top panel of Figure 2. Representative values of the numerically determined eigenfunction, x(ξ) at various discrete grid locations are provided in Table 1, along with the numerically determined value of the slope at the interface, xq. At each grid location, the associated value of the dimensionless radius, r/R, may be obtained by simply multiplying each tabulated value of ξ by the parameter, q.

Table 1:
Example Numerical Determination of Eigenfunction

(q,γe,γc)

=

(0.684012,1.194030,1.845579)

and

σc2=28.91158

Core   Envelope

g2

=

1.323609

𝔉core

=

14

Δc

=

0.02

αe

=

0.35

𝔉env

=

62.74338

Δe

=

0.00923926

k ξ x n ξ x
0 0.00 1.000000 0 1.00 -0.057649
1 0.02 0.997885 1 1.0092393 -0.076955
2 0.04 0.997182 2 1.0184785 -0.095792
49 0.98 -0.015811 49 1.452724 0.466484
50 1.00 -0.057649 50 1.461963 0.535957
xq=2.113043  

Results[edit]

Motivated by Analytic21[edit]

Continuing with our analysis of the equilibrium model that is defined by the parameters, (q,γe,γc)=(0.6840119,1.1940299,1.845579), we have used the above-described numerical algorithm, to construct 26 different eigenfunctions that simultaneously satisfy the LAWE of the core and the LAWE of the envelope for 26 different values of σc2 in the range, 300σc20. The curve traced by a sequence of small circular markers (red = core; green = envelope) in the bottom panel of Figure 2 displays each of these numerically constructed eigenfunctions in succession — in order of decreasing values of σc2 — in the form of a looped animation sequence. Also displayed in each frame of the animation, for reference, is the relevant value of σc2, as well as an unchanging, smooth, thin red/green curve that traces the analytically derived eigenfunction shown in Figure 1, for which σc2=28.91158.

Figure 2:

(q,γe,γc)=(0.6840119,1.1940299,1.845579)

Three movie frames
Three movie frames
Eigenfunction movie
Eigenfunction movie

Three frames from the animation sequence have been displayed side-by-side in the top panel of Figure 2. This image montage is presented, in part, to illustrate the degree to which our numerically generated eigenfunction matches the analytically generated eigenfunction in the specific case (σc2=28.9) for which we have been able to obtain an analytic solution to the combined/matched, core/envelope LAWEs.

Motivated by Analytic22[edit]

We have also numerically constructed an eigenfunction that matches our accompanying analytic Illustration22. In Figure 3, the numerically derived solution has been plotted on top of the analytically derived solution.

Figure 3:

(q,γe,γc)=(0.886575,1.798817,1.021798)

Numerically generated eigenfunction plotted on top of the analytically derived, Illustration22
Numerically generated eigenfunction plotted on top of the analytically derived, Illustration22

Motivated by Analytic31[edit]

We have also numerically constructed an eigenfunction that matches our accompanying analytic Illustration31. In Figure 4, the numerically derived solution has been plotted on top of the analytically derived solution.

Figure 4:

(q,γe,γc)=(0.4059596,1.180462,1.008887)

Numerically generated eigenfunction plotted on top of the analytically derived, Illustration31
Numerically generated eigenfunction plotted on top of the analytically derived, Illustration31

Unconstrained LAWEs[edit]

Here we use the most general expressions for the pair of governing LAWEs; that is, we will not force g2=. Drawing from our most general summary discussion, the LAWE for the core is,

0

=

(1η2)d2xdη2+(46η2)1ηdxdη+𝔉corex.

where,

ηξg,         and         𝔉core3ωcore22πGγcρc2αc=(σc2+8)γc6,

and,

g2

1+(ρeρc)[2(1ρeρc)(1q)+ρeρc(1q21)].

Hence, we may also write the core's LAWE as,

0

=

(g2ξ2)d2xdξ2+(4g26ξ2)1ξdxdξ+𝔉corex,

and we should be able to numerically integrate from the center, outward through the core, exactly as described above.

Separately, the LAWE for the envelope is,

0

=

[1+(g2)ξ𝒜𝒟ξ3]d2xdξ2+{3+4(g2)ξ𝒜6𝒟ξ3}1ξdxdξ+[𝒟𝔉envξ3αe]xξ2,

where,

𝒜

2(ρeρc)(1ρeρc);

1+2(ρeρc)3(ρeρc)2,

𝒟

1𝒜(ρeρc)2=[ρe/ρc2(1ρe/ρc)],

𝔉env

3ωenv22πGγeρe2αe=1γe[σc2(ρe/ρc)+8]6.

After defining the new parameter,

g2𝒜,

this LAWE for the envelope may be written as,

a2(1+a𝒟a3)xa

=

(3+4a6𝒟a3)axa[𝒟𝔉enva3αe]xa.

Inserting the same finite-difference expressions for the first and second derivatives, we therefore have,

a2(1+a𝒟a3)[x+2xa+xΔ2]

=

(3+4a6𝒟a3)a[x+x2Δ][𝒟𝔉enva3αe]xa

a2(1+a𝒟a3)[x+Δ2]+(3+4a6𝒟a3)a[x+2Δ]

=

(3+4a6𝒟a3)a[x2Δ]a2(1+a𝒟a3)[xΔ2][𝒟𝔉enva3αe]xa+a2(1+a𝒟a3)[2xaΔ2]

x+[2a2(1+a𝒟a3)+Δ(3+4a6𝒟a3)a]

=

x[Δ(3+4a6𝒟a3)a2a2(1+a𝒟a3)]+xa[4a2(1+a𝒟a3)2Δ2(𝒟𝔉enva3αe)].

Now, at the interface (only), we need to relate x to x+ in such a way that the slope gives the proper value at the interface. Specifically, we need to set,

x

=

x+2Δ(xq),

where, xq takes the value that was determined for the core. Hence, at the interface (a=1), the first step into the envelope is special and demands that,

x+[2a2(1+a𝒟a3)+Δ(3+4a6𝒟a3)a]

=

[x+2Δ(xq)][Δ(3+4a6𝒟a3)a2a2(1+a𝒟a3)]+xa[4a2(1+a𝒟a3)2Δ2(𝒟𝔉enva3αe)]

x+[4a2(1+a𝒟a3)]

=

2Δ[Δ(3+4a6𝒟a3)a2a2(1+a𝒟a3)]xq+xa[4a2(1+a𝒟a3)2Δ2(𝒟𝔉enva3αe)]

and, setting,       a=1x+

=

Δ[2(1+𝒟)Δ(3+46𝒟)]xq+xa[2(1+𝒟)Δ2(𝒟𝔉envαe)]2(1+𝒟).

Discussion[edit]

Conclusions[edit]

Initially, I was surprised to find that, by employing our above-described numerical algorithm, we were able to solve the combined/matched LAWEs for a continuum set — rather than a discrete set — of dimensionless oscillation frequencies, σc2. After all, I have been taught to believe that radial oscillation modes are obtained by solving an eigenvalue problem. After a bit of thought, I recognized that the continuum set of solutions has been obtained in the absence of a specified surface boundary condition. I suspect that the continuum of solutions can only be relevant to a real astrophysical problem after a physically meaningful surface boundary condition has been imposed; for example, a specification of the slope of the eigenfunction at the equilibrium configuration's surface. This should naturally reduce the continuum set to a discrete set of eigenvectors.

Watching the animation sequence reveals, for example, that as the value of σc2 is reduced, the number of nodes inside the configuration is reduced, in a predictable, quantized fashion. At the same time — between each drop in the integer number of nodes — the slope of the eigenfunction at the surface (r/R=1) varies between large positive, and large negative values. Hence, we should be able to find a matched solution whose slope at the surface also matches any reasonably specified boundary condition.

Additional Possibilities[edit]

  • We should be able to numerically identify a wide range of quantized radial modes of oscillation by specifying a physically reasonable surface slope then, for each successive quantum node count, tuning the choice of σc2 until the desired slope has been encountered.
  • It seems unlikely that any of our analytically derived eigenfunctions will happen to satisfy the specified boundary condition. The analytic functions have nevertheless proven to be useful in the sense that they provide a terrific check for the computational algorithm that, then, can be used to identify physically meaningful solutions numerically.
  • Tie this new (numerical) technique into the associated discussion of the relationship between stability determinations via LAWEs and global free-energy considerations.
  • The Figure 2 animation displays numerically determined solutions, x(r/R), to the combined/matching LAWEs only for non-negative values of σc2; these should be relevant only to stable radial oscillations. We should also see what can be learned from solutions associated with negative values of σc2, which may be relevant to unstable radial modes of oscillation.
Tiled Menu

Appendices: | VisTrailsEquations | VisTrailsVariables | References | Ramblings | VisTrailsImages | myphys.lsu | ADS |