First created: Nov. 18, 2003 Updated: January 28, 2008
Elasticity With Variable Temperature: (SPH,0) Case

The continuum elastic problem of an embedded sphere is extended to the situation where temperature varies spatially and temporally

Mon. Nov 24, 2003:
For each of three examples for which u_r was calculated from FDTD yesterday, I decompose u_r(R_p,t) into components by subjective visual curve fitting. The function is assumed to be of the form

u_r(R_p,t) =

3
Σ
j=0

A_j e^-t/τ_j cos( 2 π f_j t + φ_j )

In what follows, pay attention only to the relative values of the A_j. The magnitude of the A_j has not be systematically standardized here.

Figure 1124(a) This is for initial electron temperature 2000 K, τ_e-ph = 850 fs. CFM results are shown in red. Calculated using vart18.cpp	j:	0	1	2	3
	A_j	.057	.006	.148	-5.2
	τ_j	7.2417 7.25	6.87 7.24	2.46 2.22	0.91	ps
	f_j	.1354 .1344	.288 .284	0 0	0	THz
	φ_j	1.11	1.32	0	0	rad

Figure 1124(b) This is for initial electron temperature 1300 K, τ_e-ph = 850 fs. CFM results are shown in red. Calculated using vart19.cpp	j:	0	1	2	3
	A_j	.040	.003	.048	-1.7
	τ_j	7.19 7.25	6.96 7.24	2.37 2.22	0.82	ps
	f_j	.135 .1344	.288 .284	0 0	0	THz
	φ_j	1.47	1.09	0	0	rad

Figure 1124(c) This is for initial electron temperature 400 K, τ_e-ph = 850 fs. CFM results are shown in red. Calculated using vart20.cpp	j:	0	1	2	3
	A_j	.235	.032	.082	-.16
	τ_j	7.263 7.25	6.742 7.24	2.617 2.22	0.728	ps
	f_j	.1353 .1344	.292 .284	0 0	0	THz
	φ_j	2.05	1.35	0	0	rad

Sun. Nov. 23, 2003: Use τ_e-ph = 850 fs [ from PRL90_177401_2003 ] for silver to calculate the variation of B(t) with time. Note that the calculation is not detailed enough to find the actual value of B. But the scale of B doesn't make a difference. The red colored scale for B(t) is not the same for each plot. Each plot is for a different value of the initial electronic temperature.

   These examples all involve silver in BaO-P2O5 with Rp = 13 nm. This is chosen to try to match the experimental pump-probe results in Nelet et al.
   Now the power spectrum of the Fourier transform is plotted. The window function is now square, and extends all the way to t=0. In other words, no window function is used for the Fourier transform.
   The only thing among these top six plots is the way that B(t) increases with time. The final value of B is not important. Increasing B by 10 or 100 times makes no difference. What is important is the rate at which B turns on. The parameter "Btime" is the time it takes B to increase from 0 to its maximum value, which is always 0.0002, appropriate for silver increased in temperature by 10 degrees. In one case an exponential turn-on pattern was tried. The dark blue curve shows B(t).

1. Opening Comments
   In femtosecond pump probe experiments, a short laser pulse (about 20 fs long) hitting a matrix embedded nanoparticle excites the electrons in the nanoparticle (the dipole plasmon mode in the case of a metal nanoparticle). It takes about 200 fs for this energy to be equilibrated among the conduction band electrons through e-e scattering [Hodak et al. JPCB104_9954_2000 - on page 9957]. The electrons may reach a temperature of several thousand degrees Kelvin. Energy from the hot electrons is transferred (with a decay time on the order of 850 fs) to the crystal lattice, raising its temperature by approximately 10 degrees K [Hodak et al. JPCB104_9954_2000 - Fig 7 on page 9960]. However, the exact time depends on the initial electronic temperature [Hodak et al. JPCB104_9954_2000 - on page 9958] The e-ph coupling time in bulk Au is 650 ps [Hodak et al. JPCB104_9954_2000]. The coupling between electrons and phonons is primarily in the bulk material, rather than at the nanoparticle surface, unless the nanoparticle is smaller than 1 nm in diameter [Hodak et al. JPCB104_9954_2000].
   Heating of the phonons makes the material want to expand, causing an impulsive wave to propagate away from the nanoparticle. The nanoparticle remains significantly heated for about 200 ps before it cools off [Hodak et al. JPCB104_9954_2000 - at bottom of page 9958]. The breathing mode of the nanoparticle (SPH,0,0) can continue to oscillate several times before dying out. The breathing mode modulates the dipole plasmon resonance, so that a "probe" laser pulse will have time varied absorption, depending on the delay of its arrival after the first "pump" pulse.

2. Coordinates
  Let r denote coordinates in real space. Let R denote material coordinates. In this problem temperature can vary, spatially and temporally. The default equilibrium situation is where stress is zero everywhere and the temperature everywhere is 20^oC. In this situation, r = R everywhere. Note that increasing temperature will make the material want to expand. What we have in mind is a laser "pump" pulse a few fs in duration absorbed by a metal or semiconductor nanoparticle, making the nanoparticle about 10 K hotter than the surrounding material. The localized temperature increase will lead to displacement of the material. Position r now becomes dependent both on material coordinate and time as r(R,t). It is convenient to work with a displacement field u = r - R, which is also a function u(R,t).
   It will sometimes be necessary in what follows to use Cartesian components of u (denoted u_x, u_y, and u_z), and R (denoted X, Y and Z). Cartesian coordinates will be numbered x=1, y=2 and z=3 and sometimes t=0. A Roman index such as i or j will be used to denote 1, 2 or 3. In this way, the Cartesian components of u will be denoted by u_j and are functions of {X_j}, j=1,2,3. Partial derivatives will be represented with commas, such as

u_i, j =

∂ u_i

∂ X_j

u_i,jk =

∂² u_i

∂ X_j ∂ X_k

for second derivatives. Spatial derivatives will be understood to be made with respect to material coordinates {X_j}.
   It is also necessary to use spherical coordinates here because of the symmetry of the problem. The origin is located at the center of the nanoparticle. In spherical coordinates, the components of u are u_r, u_θ and u_φ. Material coordinates are R, θ and φ. Real space coordinates are r, θ and φ. Note that we don't make a distinction between real space and material coordinates for the angular coordinates.
   Because of the symmetry of the problem, u_θ = u_φ = 0, and u_r is a function of R and t only. So let the displacement field u(R,t) be (u_r(R,t),0,0).
   Later on, we will need the first and second spatial derivatives of u in Cartesian coordinates. The first step is to express u(R,t) in Cartesian coordinates after having assuming spherical symmetry:

u_x(X,Y,Z,t) = u_r(R,t) sin(θ)cos(φ) = u_r(R,t) (X/R)

u_y(X,Y,Z,t) = u_r(R,t) sin(θ)sin(φ) = u_r(R,t) (Y/R)

u_z(X,Y,Z,t) = u_r(R,t) cos(θ) = u_r(R,t) (Z/R)

where

R = sqrt(X² + Y² + Z²)

φ = arctan(Y/X)

3. Euler Lagrange Equation
The total kinetic energy of the system is T. The total potential energy of the system is U. The Lagrangian is L = T - U. Let

denote the Lagrangian density. So

. The Euler-Lagrange equation for an elastic continuum is:

dX_j

(

∂

∂ u_i, j

) -

∂

∂ u_i

= 0

[3.1]

where a summation over j from 0 to 3 is implied. The equation holds separately for i = 1, 2 and 3. This is adapted from equation (12-23) on page 552 of H. Goldstein "Classical Mechanics" 2nd edition Addison Wesley (1980).
The kinetic energy density is

(R,t) = (1/2) ρ(R) u_i,t u_i,t

[3.2]

where there is an implied summation over i from 1 to 3.
The potential energy density of an elastic medium is normally (but not here) given as

(R,t) = (1/2) c_{ijk l}(R) u_i, j u_{k, l}

[3.3]

where there are implied summations over i, j, k and l from 1 to 3. But we will introduce a modified version of

for our purposes later on.
In this problem the Lagrangian density does not depend directly on u, but only on its derivatives. The Euler-Lagrange equation therefore becomes:

(

∂

∂ u_i,t

) +

(

∂

∂ u_i,X

) +

(

∂

∂ u_i,Y

) +

(

∂

∂ u_i,Z

) = 0

[3.4]

separately holding for i = 1, 2 and 3.

4. Variable Temperature Elasticity
Let T denote the temperature in Kelvin degrees. Ordinarily, in elasticity theory, the potential energy density

is locally minimized when u_i, j = 0 for all i and j. This implies that u is a constant vector independent of material coordinate R. As a result real space coordinate r = u + R. In other words, apart from an offset, real space coordinates agree with material coordinates. [one objection seems to be that a rigid rotation should be allowed -- I don't know the resolution to this yet]. But the important point here is that the "density" of the two coordinate systems is the same. In order to allow the possibility of the density of the medium varying with temperature T, we introduce the modified potential energy function:

(R,t) = (1/2) c_{ijk l}(R) (u_i, j - B(R,t) δ_ij) (u_{k, l} - B(R,t) δ_{k l})

[4.1]

where δ_ij is the Kronecker delta.
   It is important to emphasize that this is our proposed form for the potential energy, and not something derived from fundamental principles. So the following should be regarded as a possible model of elasticity with temperature variation. There may be other ways to do this.
   This new potential energy density will be locally minimized when u_i, j = B(R,t) δ_ij for all i and j. One solution for this is: (x - X) = BX, (y - Y) = BY and (z- Z) = BZ or x = (B+1)X, y = (B+1)Y and z = (B+1)Z. That is, a positive B corresponds to increasing the volume of material per unit mass by a factor of (B+1)³. The density is decreased when B is positive. The choice of material coordinates R to label atoms is arbitrary. So let us fix that R is defined at T = 293 K (20 C). Then B will be a function of (T-293) such that B = 0 when T = 293 K. The simplest approximation would be that B = (const.)(T-293) where (const.) is related to the thermal expansion coefficient of the material. Since temperature can be a function T(R,t) of material position and time, B will also be a function B(R,t).

5. Equation of Motion
   Adopting the modified potential energy density, the Lagrangian density is now:

= (1/2) ρ u_i,t u_i,t - (1/2) c_{ijk l} (u_i, j - B δ_ij) (u_{k, l} - B δ_{k l})

[5.1]

where ρ depends on material position R but not on the time. This is a consequence of conservation of mass, and not the way that the material is being modelled. c_{ijk l} and B all may have R and t dependence. It is certainly the case that elastic constants of materials are temperature dependent, as is the density. But in the simplest version of our calculation we will let B be temperature dependent, but will let c_{ijk l} be functions of material coordinate only.
First evaluate

∂

∂u_i,t

= ρ u_i,t

[5.2]

where "t" denotes time and i could be 1, 2 or 3. And next

(

∂

∂u_i,t

) = ρ u_i,tt = ρ

..
u_i

[5.3]

Note that the time-independence of ρ is important here. This illustrates the nature of ρ(R) as we define it in constrast to ρ_rsc(r,t), the mass density in real space coordinates. The latter quantity fluctuates with time as the material vibrates. Material coordinate density ρ(R) does not. Material coordinate density does not vary with temperature either.
Next, for any i and j

∂

∂u_i, j

= -

Σ
k l

c_{ijk l}(R,t) ( u_{k, l} - B(R,t) δ_{k l} )

[5.4]

This allows for possible temperature dependence both of the elastic constants c_{ijk l} and B.
Now, take the spatial derivatives for each j, and summing:

dX_j

(

∂

∂u_i, j

) = -

Σ
j k l

dX_j

(

c_{ijk l}(R,t) ( u_{k, l} - B(R,t) δ_{k l} )

)

[5.5]

where j, k and l are each summed from 1 to 3. This equation holds separately for i = 1, 2 and 3.
The equation of motion resulting from the Euler-Lagrange equation may now be written in the form:

ρ(R)

..
u_i

Σ
j k l

dX_j

(

c_{ijk l}(R,t) ( u_{k, l} - B(R,t) δ_{k l} )

)

[5.6]

Taking the derivatives, the right side breaks into four terms:

ρ(R)

..
u_i

Σ
j k l

c_{ijk l} u_{k, lj}

Σ
j k l

c_{ijk l, j} u_{k, l}

Σ
j k l

c_{ijk l, j} B δ_{k l}

Σ
j k l

c_{ijk l} B_, j δ_{k l}

[5.7]

or, after removing Kronecker deltas,

ρ(R)

..
u_i

Σ
j k l

c_{ijk l} u_{k, lj} +

Σ
j k l

c_{ijk l, j} u_{k, l} -

Σ
j k

c_ijkk, j B -

Σ
j k

c_ijkk B_, j

[5.8]

6. Isotropic but Inhomogeneous Elasticity
We want to assume that the elasticity is locally isotropic, i.e. invariant at a point under any rotation that keeps that point fixed. On the other hand, the material may be inhomogeneous, meaning that the elastic constants c_{ijk l} are functions of material coordinate. Even though there are 3×3×3×3 = 81 elements in the fourth elasticity tensor, an isotropic material has only two independent elastic constants, C₁₁ and C₄₄. For compactness it will be also be useful to use C₁₂ = C₁₁ - 2 C₄₄. The following tablee shows what is what:

c₁₁₁₁ c₂₂₂₂ c₃₃₃₃	C₁₁
c₁₂₁₂ c₁₂₂₁ c₂₁₁₂ c₂₁₂₁ c₁₃₁₃ c₃₁₃₁ c₃₁₁₃ c₁₃₃₁ c₂₃₂₃ c₃₂₃₂ c₃₂₂₃ c₂₃₃₂	C₄₄
c₁₁₂₂ c₂₂₁₁ c₁₁₃₃ c₃₃₁₁ c₂₂₃₃ c₃₃₂₂	C₁₂
c₁₁₁₂ c₁₁₂₁ c₁₁₁₃ c₁₁₃₁ c₁₁₂₃ c₁₁₃₂ c₂₂₁₂ c₂₂₂₁ c₂₂₁₃ c₂₂₃₁ c₂₂₂₃ c₂₂₃₂ c₃₃₁₂ c₃₃₂₁ c₃₃₁₃ c₃₃₃₂ c₃₃₂₃ c₃₃₃₁ c₁₂₁₁ c₂₁₁₁ c₁₃₁₁ c₃₁₁₁ c₂₃₁₁ c₃₂₁₁ c₁₂₂₂ c₂₁₂₂ c₁₃₂₂ c₃₁₂₂ c₂₃₂₂ c₃₂₂₂ c₁₂₃₃ c₂₁₃₃ c₁₃₃₃ c₃₁₃₃ c₂₃₃₃ c₃₂₃₃ c₁₂₁₃ c₁₂₃₁ c₂₁₃₁ c₂₁₁₃ c₁₃₁₂ c₃₁₁₂ c₃₁₂₁ c₁₃₂₁ c₁₂₂₃ c₁₂₃₂ c₂₁₃₂ c₂₁₂₃ c₂₃₁₂ c₃₂₁₂ c₃₂₂₁ c₂₃₂₁ c₁₃₂₃ c₁₃₃₂ c₃₁₃₂ c₃₁₂₃ c₂₃₁₃ c₃₂₁₃ c₃₂₃₁ c₂₃₃₁	0

7. Spherical Symmetry
We assume that the following depend on R = ||R|| and time t but not on θ or φ: ρ, C₁₁, C₄₄, B and u_r. But note that ρ depends on R but not t. u_θ and u_φ are both 0. We assume that the origins of the Cartesian and spherical coordinate systems coincide. In that case, the propagation of an elastic wave can be understood by looking at the positive X-axis only.
Take equation of motion [5.8] and specialize to the positive X-axis, at X=R with Y=Z=0:

..
u₁

Σ
j k l

c_{1jk l} u_{k, lj} +

Σ
j k l

c_{1jk l, j} u_{k, l} -

Σ
j k

c_1jkk, j B -

Σ
j k

c_1jkk B_, j

[7.1]

Since c_{ijk l} and B are both scalar fields depending on R and not θ or φ, only their X partial derivatives are nonzero on the X-axis. Partial derivatives with respect to Y and Z are zero. Therefore the equation of motion simplifies to:

..
u₁

Σ
j k l

c_{1jk l} u_{k, lj} +

Σ
k l

c_{11k l,1} u_{k, l} -

Σ
k

c_11kk,1 B -

Σ
k

c_11kk B_,1

[7.2]

In addition, note that for an isotropic material c_{11k l} is nonzero only if k = l. So this becomes slightly simpler:

..
u₁

Σ
j k l

c_{1jk l} u_{k, lj} +

Σ
k

c_11kk,1 u_k,k -

Σ
k

c_11kk,1 B -

Σ
k

c_11kk B_,1

[7.3]

It is necessary to translate each Cartesian component and partial derivative in the above into corresponding quantities in spherical coordinates in terms of the radial displacement field u_r(R,t) and its derivatives with respect to R: u_r'(R,t) and u_r"(R,t):

..
u₁(R,0,0,t)

..
u_r(R,t)

[7.4]

u_1,1(R,0,0,t) = u_r'(R,t)

[7.5]

u_2,2(R,0,0,t) = u_r(R,t) / R

[7.6]

u_3,3(R,0,0,t) = u_r(R,t) / R

[7.7]

u_i, j(R,0,0,t) = 0 if i ≠ j

[7.8]

u_1,11(R,0,0,t) = u_r"(R,t)

[7.9]

u_1,22(R,0,0,t) =

u_r'(R,t)

u_r

R²

[7.10]

and on the X-axis, u_1,22 = u_1,33 = u_2,12 = u_2,21 = u_3,13 = u_3,31. Otherwise u_i, jk = 0. To be sure, I verified all of these identities numerically using single or double numerical differentiation in the C++ computer program vart3.cpp.
The right side of equation [7.3] has the following terms:

Table 2:
i j k l	Cartesian	Spherical
1111	C₁₁ u_1,11	C₁₁ u_r"
1122	C₁₂ u_2,21	C₁₂(u_r' / R - u_r / R² )
1133	C₁₂ u_3,31	C₁₂(u_r' / R - u_r / R² )
1212	C₄₄ u_1,22	C₄₄(u_r' / R - u_r / R² )
1221	C₄₄ u_2,12	C₄₄(u_r' / R - u_r / R² )
1313	C₄₄ u_1,33	C₄₄(u_r' / R - u_r / R² )
1331	C₄₄ u_3,13	C₄₄(u_r' / R - u_r / R² )
- -1-	C_11,1 u_1,1	dC₁₁/dR u_r'
- -2-	C_12,1 u_2,2	dC₁₂/dR u_r / R
- -3-	C_12,1 u_3,3	dC₁₂/dR u_r / R
- -1-	-B C_11,1	-B dC₁₁/dR
- -2-	-B C_12,1	-B dC₁₂/dR
- -3-	-B C_12,1	-B dC₁₂/dR
- -1-	-C₁₁ B_,1 &nbssp;	-C₁₁ dB/dR
- -2-	-C₁₂ B_,1 &nbssp;	-C₁₂ dB/dR
- -3-	-C₁₂ B_,1 &nbssp;	-C₁₂ dB/dR

All of these terms in the right column, when added up, are equal to ρ(R) ∂²u_r /∂t². This provides a second order differential equation for u_r(R,t) which can be numerically integrated starting with initial conditions.

References

J. H. Hodak, A. Henglein and G. V. Hartland, "Photophysics of Nanometer Sized Particles: Electron-Phonon Coupling and Coherent Excitation of Breathing Vibrational Modes" J. Phys. Chem. B 104 9954-9965 (2000).

Daniel Murray
Associate Professor
Math, Stats & Physics Unit
University of British Columbia - Okanagan
Kelowna, BC, Canada
daniel "dot" murray "at" ubc "dot" ca

For a list of related articles click here.

Hosted by www.Geocities.ws