An efficient spectral method for computing dynamics of rotating two-component Bose–Einstein condensates via coordinate transformation

# An efficient spectral method for computing dynamics of rotating two-component Bose–Einstein condensates via coordinate transformation

Ming Ju Beijing Computational Science Research Center, No. 3 He-Qing Road, Hai-Dian District, Beijing, P.R. China 100084 Qinglin Tang Department of Mathematics and Center for Computational Science and Engineering, National University of Singapore, Singapore 119076 Yanzhi Zhang Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409-0020, USA
###### Abstract

In this paper, we propose an efficient and accurate numerical method for computing the dynamics of rotating two-component Bose–Einstein condensates (BECs) which is described by coupled Gross–Pitaevskii equations (CGPEs) with an angular momentum rotation term and an external driving field. By introducing rotating Lagrangian coordinates, we eliminate the angular momentum rotation term from the CGPEs, which allows us to develop an efficient numerical method. Our method has spectral accuracy in all spatial dimensions and moreover it can be easily implemented in practice. To examine its performance, we compare our method with those reported in literature. Numerical results show that to achieve the same accuracy, our method needs much shorter computing time. We also applied our method to study the dynamic properties of rotating two-component BECs. Furthermore, we generalize our method to solve the vector Gross–Pitaevskii equations (VGPEs) which is used to study rotating multi-component BECs.

###### keywords:
Rotating two-component BECs, coupled/vector Gross–Pitaevskii equations, angular momentum rotation, rotating Lagrangian coordinates, time-splitting.

## 1 Introduction

The Bose–Einstein condensation (BEC), which affords an astonishing glimpse into the macroscopic quantum world, has been extensively studied since its first realization in 1995 Anderson1995 (); Bradley (); Davis1995 (). Later, with the observation of quantized vortices in BECs Matthews1999 (); Madison2000 (), attention has been broaden to explore vortex states and their dynamics associated with superfluidity. Rotating BECs which are known to exhibit highly regular vortex lattices have been heavily studied both experimentally and theoretically Abo2001 (); Madison2001 (); Lieb1 (); Fetter2009 (). On the other hand, multi-component BECs admit numerous interesting phenomena absent from single-component condensates, for example, domain walls, vortons, square vortex lattices and so on; see Hall1998 (); Ho1996 (); Chui2000 (); Jezek2001 (); Kasamatsu2004 (); Kasamatsu2005 (); Kobyakov2011 () and references therein. As the simplest cases, two-component BECs provide a good opportunity to investigate the properties of multi-component condensation.

The first experiment of two-component BECs was carried out in and hyperfine states of MBGCW (). At temperatures much smaller than the critical temperature , a rotating two-component BEC with an external driving field (or an internal Josephson junction) can be well described by two self-consistent nonlinear Schrödinger equations (NLSEs), also known as the coupled Gross–Pitaevskii equations (CGPEs). The dimensionless CGPEs has the following form Kasamatsu2003 (); Kasamatsu2005 (); Wang2007 (); Zhang2007 (); Bao2008 (); Bao2011 (); Kobyakov2011 ():

 i∂ψ1(x,t)∂t=[−12∇2+V1(x)+(β11|ψ1|2+β12|ψ2|2)−ΩLz]ψ1−λψ2, (1.1) i∂ψ2(x,t)∂t=[−12∇2+V2(x)+(β21|ψ1|2+β22|ψ2|2)−ΩLz]ψ2−λψ1,x∈Rd, t>0. (1.2)

Here, ( or ) is the Cartesian coordinate vector, is the time and is the complex-valued macroscopic wave function of the th () component. The interaction constants (for ), where is the total number of atoms in two-component BECs, is the dimensionless spatial unit and represents the -wave scattering lengths between the th and th components (positive for repulsive interaction and negative for attractive interaction). The constant describes the effective Rabi frequency to realize the internal atomic Josephson junction by a Raman transition, represents the speed of angular momentum rotation and is the -component of the angular momentum operator. The real-valued function () represents the external trapping potential imposed on the th component. In most BEC experiments, a harmonic potential is used, i.e.,

 Vj(x)=12{γ2x,jx2+γ2y,jy2,d=2,γ2x,jx2+γ2y,jy2+γ2z,jz2,d=3,  j=1,2. (1.3)

The initial conditions of (1.1)–(1.2) are given by

 ψj(x,0)=ψ0j(x), x∈Rd,  j=1,2. (1.4)

There are two important invariants associated with the CGPEs in (1.1)–(1.2): the total mass (or normalization), i.e.,

 N(t):=∥Ψ(⋅,t)∥2=N1(t)+N2(t)≡∥Ψ(⋅,0)∥2=1, t≥0, (1.5)

where and (t) is the mass of the th component at time , which is defined by

 Nj(t):=∥ψj(⋅,t)∥2=∫Rd|ψj(x,t)|2dx, t≥0,  j=1,2, (1.6)

and the energy

 E(Ψ(⋅,t)) = ∫Rd[2∑j=1(12|∇ψj|2+Vj(x)|ψj|2+βjj2|ψj|4−ΩRe(ψ∗jLzψj)) (1.7) +β12|ψ1|2|ψ2|2−2λRe(ψ1ψ∗2)]dx=E(Ψ(⋅,0)), t≥0,

where and represent the real part and the conjugate of a function , respectively. In fact, if there is no external driving filed (i.e., in (1.1)–(1.2)), the mass of each component is also conserved, i.e., () for . These invariants can be used, in particular, as benchmarks and validation of numerical algorithms for solving the CGPEs (1.1)–(1.2).

Many numerical methods have been proposed to study the dynamics of the non-rotating two-component BECs, i.e., when , with/without external driving field Bao2004 (); Sepulveda2008 (); JGSGZ (); CW (). Compared to other methods, the time-splitting pseudo-spectral method in Bao2004 () is one of the most successful methods. It has spectral order of accuracy in space and can be easily implemented, i.e., they can achieve both the accuracy and efficiency. However, the appearance of the angular rotational term hinders the direct application of those methods to study the rotating two-component BECs. Recently, several numerical methods were proposed for simulating the dynamics of rotating two-component BECs Zhang2007 (); Wang2007 (); Bao2008 (); Corro2009 (); Hsueh2011 (); Jin2013 (). For example, in Zhang2007 (), a pseudo-spectral type method was proposed by reformulating the problem in two-dimensional polar coordinates or three-dimensional cylindrical coordinates. While in Wang2007 (), the authors designed a time-splitting alternating direction implicit method, where the angular rotation term is treated in - and -directions separately. Although these methods have higher spatial accuracy compared to those finite difference/element methods, they have their own limitations. The method in Zhang2007 () is only of second-order or fourth-order in the radial direction, while the implementation of the method in Wang2007 () could become quite involved. One possible approach to overcome those limitations is to relax the constrain of the rotational term, which is the main aim of this paper. In this paper, we propose a simple and efficient numerical method to solve the CGPEs (1.1)–(1.2). The main merits of our method are: (i) Using a rotating Lagrangian coordinate transform, we reformulate the original CGPEs in (1.1)–(1.2) to one without angular momentum rotation term. Then, the time-splitting pseudo-spectral method designed for the non-rotating BECs, which are of spectral order accuracy in space and easy to implemented, can be directly applied to solve the CGPEs in new coordinates. Moreover, (ii) our method solves the CGPEs in two splitting steps instead of three steps in literature Zhang2007 (); Wang2007 (), which makes our method more efficient.

The paper is organized as follows. In Section 2, we introduce a rotating Lagrangian coordinate and then cast the CGPEs (1.1)–(1.4) in the new coordinate system. A simple and efficient numerical method is introduced to discretize the CGPEs under a rotating Lagrangian coordinate in Section 3, which is subsequently generalized in Section 4 to solve the VGPEs for multi-component BECs. To test its performance, we compare our method with those reported in literature and apply it to study the dynamics of rotating two-component BECs in Section 5. In Section 6, we make some concluding remarks.

## 2 CGPEs under a rotating Lagrangian coordinate

In this section, we first introduce a rotating Lagrangian coordinate and then reformulate the CGPEs (1.1)–(1.4) in the new coordinate system. In the following, we will always refer the Cartesian coordinates as the Eulerian coordinates. For any time , let be an orthogonal rotational matrix defined as Garcia2001 (); Antonelli2012 (); Bao2013 ()

 A(t)=(cos(Ωt)sin(Ωt)−sin(Ωt)cos(Ωt)),  if \ \ d=2, (2.1)

and

 A(t)=⎛⎜⎝cos(Ωt)sin(Ωt)0−sin(Ωt)cos(Ωt)0001⎞⎟⎠,  if \ \ d=3. (2.2)

It is easy to verify that for any and with the identity matrix. Referring the Cartesian coordinates as the Eulerian coordinates, we introduce the rotating Lagrangian coordinates as

 ˜x=A−1(t)x=AT(t)x⇔x=A(t)˜x, x∈Rd,t≥0. (2.3)

and reformulate the wave functions in the new coordinates as

 ϕj(˜x,t):=ψj(x,t)=ψj(A(t)˜x,t), ˜x∈Rd,t≥0,  j=1,2. (2.4)

We see that when , the transformation in (2.3) does not change the coordinate in -direction, that is, and the coordinate transformation essentially occurs only in -plane for any . Fig. 1 illustrates the geometrical relation between the -plane in the Eulerian coordinates and -plane in the rotating Lagrangian coordinates for .

Furthermore, it is easy to see that when the rotating Lagrangian coordinates become exactly the same as the Eulerian coordinates , i.e., when . Note that we assume that in this paper.

From (2.3)–(2.4), we obtain that

 ∂tϕj(˜x,t)=∂tψj(x,t)+∇ψj(x,t)⋅(˙A(t)˜x)=∂tψj(x,t)−Ω(x∂y−y∂x)ψj(x,t), ∇ϕj(˜x,t)=A−1(t)∇ψj(x,t),∇2ϕj(˜x,t)=∇2ψj(x,t), x∈Rd,t≥0, j=1,2.

Substituting the above derivatives into (1.1)–(1.2) gives the following -dimensional CGPEs in the rotating Lagrangian coordinates :

 i∂ϕ1(˜x,t)∂t=[−12∇2+W1(˜x,t)+(β11|ϕ1|2+β12|ϕ2|2)]ϕ1−λϕ2, (2.5) i∂ϕ2(˜x,t)∂t=[−12∇2+W2(˜x,t)+(β21|ϕ1|2+β22|ϕ2|2)]ϕ2−λϕ1, ˜x∈Rd,t>0. (2.6)

The corresponding initial conditions are

 ϕj(˜x,0):=ϕ0j(˜x)=ψj(x,0)=ψ0j(x), ˜x=x∈Rd. (2.7)

In (2.5)–(2.6), () denotes the effective potential of the th component, which is obtained from

 Wj(˜x,t)=Vj(A(t)˜x), ˜x∈Rd,t≥0,  j=1,2. (2.8)

In particular, if is a harmonic potential as defined in (1.3), then has the form

 Wj(˜x,t)

for . Hence, when the external harmonic potentials are radially symmetric in two dimensions (2D) or cylindrically symmetric in three dimensions (3D), i.e., , the potential

 Wj(˜x,t)=Vj(˜x), ˜x∈Rd,t≥0, j=1,2, (2.9)

become time-independent.

In rotating Lagrangian coordinates, the wave functions satisfy the normalization

 ˜N(t):=∥Φ(⋅,t)∥2=˜N1(t)+˜N2(t)≡∥Φ(⋅,0)∥2=1, t≥0, (2.10)

where is the mass of the th component at time , i.e.,

 ˜Nj(t):=∥ϕj(⋅,t)∥2=∫Rd|ϕj(˜x,t)|2d˜x=Nj(t),t≥0,  j=1,2. (2.11)

Similarly, when the mass of each component is also conserved, i.e., for and . The energy associated with the CGPEs (2.5)–(2.6) is

 ˜E(Φ(⋅,t))=∫Rd[2∑j=1(12|∇ϕj|2+Wj(˜x,t)|ϕj|2−∫t0|ϕj|2∂τWj(˜x,τ)dτ+βjj2|ϕj|4) +β12|ϕ1|2|ϕ2|2−2λRe(ϕ1ϕ∗2)]d˜x=˜E(Φ(⋅,0)), t≥0. (2.12)

As we see in (2.9), when (), the potential is time-independent, which implies that the term of in this case.

Compared to (1.1)–(1.2), the CGPEs (2.5)–(2.6) in rotating Lagrangian coordinates does not have the angular momentum rotational term, which eliminates the difficulties in discretizing the CGPEs and allows us to develop an efficient spectral method to solve (2.5)–(2.6).

## 3 Numerical method

In this section, we present a time-splitting spectral method to study the dynamics of rotating two-component BECs. To the best of our knowledge, so far all numerical methods computing dynamics of rotating two-component BECs in literature are based on discretizing the CGPEs (1.1)–(1.2) in Eulerian coordinates Zhang2007 (); Wang2007 (); Bao2008 (). However, the appearance of rotating angular momentum term in Eulerian coordinates makes it very challenging to develop an efficient methods with higher accuracy but less computational efforts. In the following, instead of simulating (1.1)–(1.2) in Eulerian coordinates, we solve the CGPEs (2.5)–(2.6) in rotating Lagrangian coordinates. Hence, we avoid the discretization of the angular rotational term and it makes numerical method simpler and more efficient than those reported in Wang2007 (); Zhang2007 ().

In practical computations, we truncate the problem (2.5)–(2.6) into a bounded computational domain and consider

 i∂tϕ1(˜x,t)=[−12∇2+W1(˜x,t)+(β11|ϕ1|2+β12|ϕ2|2)]ϕ1−λϕ2, (3.1) i∂tϕ2(˜x,t)=[−12∇2+W2(˜x,t)+(β21|ϕ1|2+β22|ϕ2|2)]ϕ2−λϕ1, ˜x∈D,t>0, (3.2)

along with the initial conditions

 ϕj(˜x,t)=ϕ0j(˜x), ˜x∈¯¯¯¯D,with∫¯¯¯¯D(|ϕ01(˜x)|2+|ϕ02(˜x)|2)d˜x=1. (3.3)

The following homogeneous Dirichlet boundary conditions are considered here, i.e.,

 ϕj(˜x,t)=0, ˜x∈∂D,t>0,  j=1,2. (3.4)

Due to the confinement of the external potential and conservation of the normalization (2.10) and energy (2), the wave function vanishes as . Hence, it is natural to impose homogeneous Dirichlet boundary conditions to the truncated problem (3.1)–(3.3). The use of more sophisticated boundary conditions for more generalized cases, e.g., absence of trapping potential, is an interesting topic that remains to be examined in the future ABK (); BT (). In practical simulations, the computational domain is chosen as if and if . Moreover, we use sufficiently large domain to ensure the homogeneous Dirichlet boundary conditions do not introduce aliasing error. Usually, the diameter of the bounded computational domain depends on the problem. In general, it should be larger than the “Thomas-Fermi radius” Bao2005 (); Bao2006 ().

### 3.1 Time-splitting method

In the following, we use the time-splitting method to discretize the problem (3.1)–(3.4) in time. To do it, we choose a time step and define time sequence for . Then from time to , we numerically solve the CGPEs (3.1)–(3.2) in two steps, i.e., solving

 i∂ϕj(˜x,t)∂t=−12∇2ϕj(˜x,t)−λϕ(3−j)(˜x,t), ˜x∈D,tn≤t≤tn+1,  j=1,2, (3.5)

and

 (3.6)

In fact, Eq. (3.5) is coupled linear Schrödinger equations and its discretization will be discussed later.

We notice that in (3.6), both and are invariants in time , i.e., () for any . Thus, for time , (3.6) is equivalent to

 i∂tϕj(˜x,t)=(Wj(˜x,t)+2∑k=1βjk|ϕk(˜x,tn)|2)ϕj(˜x,t), ˜x∈D,  j=1,2. (3.7)

Integrating (3.7) exactly in time leads to the solution of (3.6), i.e.,

 ϕj(˜x,t)=ϕj(˜x,tn)exp[−i((t−tn)2∑k=1βjk|ϕk(˜x,tn)|2+∫ttnWj(˜x,τ)dτ)],  j=1,2, (3.8)

for and .

###### Remark 3.1.

If ( is a harmonic potential as defined in (1.3), then the integral in (3.8) can be evaluated analytically, i.e.,

 ∫ttnWj(˜x,τ)dτ=(γ2x,j+γ2y,j)(˜x2+˜y2)4(t−tn)+U(˜x,t)+{0,d=2,12γ2z,j˜z2(t−tn),d=3, (3.9)

where

 U(˜x,t) = (γ2x,j−γ2y,j)4∫ttn[(˜x2−˜y2)cos(2Ωτ)+2˜x˜ysin(2Ωτ)]dτ = (γ2x,j−γ2y,j)(˜x2−˜y2)8Ω[sin(2Ωt)−sin(2Ωtn)]−(γ2x,j−γ2y,j)˜x˜y4Ω[cos(2Ωt)−cos(2Ωtn)].

Typically, when (), we have .

For a general potential , if the integral in (3.8) can not be found analytically, numerical quadratures such as Trapezoidal rule or Simpson’s rule can be used to calculate its approximation Bao2006 (); Bao2013 ().

### 3.2 Discretization of coupled linear Schrödinger equations

In the following, we first introduce a linear transformation of the wave functions () such that the coupled linear Schrödinger equations become independent of each other. Then we describe the sine pseudospectral discretization in two-dimensional case. Its generalization to three dimensions is straightforward.

Let the matrix

 P=(111−1), (3.10)

and denote

 (φ1(˜x,t)φ2(˜x,t))=P(ϕ1(˜x,t)ϕ2(˜x,t))=(ϕ1(˜x,t)+ϕ2(˜x,t)ϕ1(˜x,t)−ϕ2(˜x,t)), ˜x∈Rd,t≥0. (3.11)

Combining (3.5) and (3.11), we obtain the following equations for ():

 i∂tφ1(˜x,t)=−12∇2φ1(˜x,t)−λφ1(˜x,t), (3.12) i∂tφ2(˜x,t)=−12∇2φ2(˜x,t)+λφ2(˜x,t), ˜x∈D,   tn≤t≤tn+1. (3.13)

It is easy to see that the functions and are independent in (3.12)–(3.13), which allows us to solve them separately.

Choose two even integers and denote the index set

 TJK={(p,q)|1≤p≤J−1,  1≤q≤K−1}.

Define the function

 Upq(˜x)=sin(μxp(˜x−a))sin(μyq(˜y−c)), ˜x=(˜x,˜y)T∈D,  (p,q)∈TJK.

Assume that

 φj(˜x,t)=J−1∑p=1K−1∑q=1ˆφj,pq(t)Upq(˜x), ˜x∈D,t∈[tn,tn+1],  j=1,2, (3.14)

where is the discrete sine transform of corresponding to frequencies and

 μxp=pπb−a, μyq=qπe−c,(p,q)∈TJK.

Substituting (3.14) into (3.12)–(3.13) leads to

 ˆφ1,pq(t)=ˆφ1,pq(tn)exp[−i((μxp)2+(μyq)22−λ)(t−tn)], (3.15) ˆφ2,pq(t)=ˆφ2,pq(tn)exp[−i((μxp)2+(μyq)22+λ)(t−tn)],(p,q)∈TJK,  t∈[tn,tn+1]. (3.16)

Combining (3.15)–(3.16) and (3.14) and noticing the linear transformation in (3.11), we obtain a sine pseudospectral approximation to (3.5), i.e.,

 ϕj(˜x,t)=J−1∑p=1K−1∑q=1[cos(λ(t−tn))ˆϕj,pq(tn)+isin(λ(t−tn))ˆϕ(3−j),pq(tn)]ηpq(t)Upq(˜x), (3.17)

for , where

 ηpq(t)=exp[−i2((μxp)2+(μyq)2)(t−tn)],(p,q)∈TJK.

In (3.17), is the discrete sine transform of () corresponding to the frequency . We remark here that although the solution (3.17) is found via (3.15)–(3.16), (3.14) and (3.11), in practice we only need to compute and to obtain (3.17).

### 3.3 Implementation of the method

For convenience of the readers, in the following we will summarize our method and describe its implementation. For simplicity of notations, the method will only be presented in two-dimensional case. Choose spatial mesh sizes and in - and -directions, respectively. Define

 ˜xs=a+sh˜x,0≤s≤J;˜yl=c+lh˜y,0≤l≤K.

Let denote the numerical approximation to . From to , we use the second-order Strang splitting method Strang1968 (); Glowinski1989 (); Bao2005 () to combine the two steps in (3.5) and (3.6), i.e,

 ϕ(1)j,sl=ϕnj,slexp[−i(Δt22∑k=1βjk|ϕnk,sl|2+∫t+Δt/2tnWj(˜xs,˜yl,τ)dτ)], (3.18) ϕ(2)j,sl=J−1∑p=1K−1∑q=1e−iΔt2[(μxp)2+(μyq)2][cos(λΔt)ˆϕ(1)j,pq+isin(λΔt)ˆϕ(1)(3−j),pq]sin(spπJ)sin(lqπK), (3.19) ϕn+1j,sl=ϕ(2)j,slexp[−i(Δt22∑k=1βjk|ϕ(2)k,sl|2+∫tn+1t+Δt/2Wj(˜xs,˜yl,τ)dτ)],j=1,2, (3.20)

for , and . At , the initial conditions (3.3) are discretized as

 ϕ0j,sl=ϕ0j(˜xs,˜yl), 0≤s≤J,0≤l≤K,j=1,2, (3.21)

Our method described in (3.18)–(3.21) is explicit and it is easy to implement. Furthermore, the memory cost is and the computational cost per time step is if a 2D CGPEs is solved. In 3D case, the memory cost and the computational cost per time step are and , respectively, where the even integer and is the number of grid points in -direction in 3D.

###### Remark 3.2.

The solutions and obtained from (3.18)–(3.20) are grid functions on the bounded computational domain in rotating Lagrangian coordinates. To obtain the wave functions and satisfying the CGPEs (1.1)–(1.2) over a set of fixed grid points in the Eulerian coordinates , we can use the standard Fourier/sine interpolation operators from the discrete numerical solution to construct an interpolation continuous function over Boyd1992 (); Shen ().

###### Remark 3.3.

If the potential in (1.3) is replaced by a time-dependent potential, e.g., , the rotating Lagrangian coordinate transformation and the numerical method are still valid provided that we replace in (2.8) by for and .

## 4 Extension to rotating multi-component BECs

In Sections 23, we presented an efficient and accurate numerical method to compute the dynamics of rotating two-component BECs with internal Josephson junction. In fact, this method can be easily generalized to solve the vector Gross–Pitaevskii equations (VGPEs) with an angular momentum rotation term and an external driving field, which describes the dynamics of rotating multi-component BECs Bao2004 (); CLL (); Liu (); Lieb1 ().

Suppose there are species in multi-component BECs. Denote the complex-valued macroscopic wave function for the th component as for . Let . Then the evolution of the wave function is governed by the following self-consistent VGPEs Bao2004 (); CLL (); Liu (); Lieb1 (); LW ():

 i∂Ψ(x,t)∂t=[−12∇2+V(x)+F(Ψ)−ΩLz+g(t)B]Ψ, x∈Rd,t>0 (4.1)

with the initial conditions

 Ψ(x,0)=Ψ0(x)=(ψ01(x),…,ψ0M(x))T,  x∈Rd. (4.2)

The matrix represents the external traping potentials and with

 Fj(Ψ)=M∑k=1βjk|ψk(x,t)|2,  j=1,…,M,

where the constant describes the interaction strength between the th and th components. is a real-valued scalar function and is a real-valued diagonalizable constant matrix.

To solve (4.1)–(4.2), similarly we introduce the rotating Lagrangian coordinates as defined in (2.3)–(2.4) and cast the VGPEs in the new coordinates. Then we truncate it into a bounded computational domain and consider the following VGPEs with homogenous Dirichlet boundary conditions:

 i∂Φ(˜x,t)∂t=[−12∇2+W(˜x,t)+F(Φ)+g(t)B]Φ, ˜x∈D,t>0, (4.3) Φ(˜x,0)=Φ0(˜x), ˜x∈¯¯¯¯DandΦ(˜x,t)=0, ˜x∈∂D,t>0, (4.4)

where and with for . From time to , we split the VGPEs (4.3) into two subproblems and solve

 i∂tΦ(˜x,t)=[−12∇2+g(t)B]Φ, ˜x∈D,tn≤t≤tn+1 (4.5)

for a time step of length , followed by solving

 i∂tΦ(˜x,t)=[W(˜x,t)+F(Φ)]Φ, ˜x∈D,tn≤t≤tn+1 (4.6)

for the same time step.

Equation (4.6) can be integrated exactly in time and the solution is

 ϕj(˜x,t)=ϕj(˜x,tn)exp[−i((t−tn)Fj(Φ(˜x,tn))+∫ttnWj(˜x,τ)d