Sedimentation of particles acted upon by a vertical, time-oscillating force

We analyse a spatially restricted one-dimensional diffusion process occurring in a half-space under the influence of a harmonically oscillating and space-homogeneous driving force. The force possesses a dc component which points against the reflecting boundary. Our approach is based on an exact time-asymptotic solution of the underlying non-convolution integral equation. It yields the full time- and space-resolved density profile in the asymptotic non-equilibrium regime. We identify scaling laws for the complex amplitudes which control the nonlinear response. The time-averaged density profile exhibits surprising features. At the boundary, it equals the equilibrium value in the corresponding problem without modulation. However, farther from the boundary, it is first smaller and then bigger in comparison with its static counterpart. It can even develop a well-pronounced partially depleted region in the vicinity of the boundary which is then followed by a region of increased concentration. The time-averaged mean position always exceeds the equilibrium mean position. Their difference, viewed as a function of the temperature, displays a resonance-like maximum.


3
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT The known solutions mostly concern the spatially unlimited diffusion, i.e. they do not incorporate absorbing and/or reflecting barriers. On the other hand, spatially restricted diffusion processes are important in many situations in physics, chemistry and biology. Moreover, the problems with a static boundary and a time-dependent potential can be mapped on to the settings with a static potential and a moving boundary [28]. The most interesting regime occurs if the boundary motion competes with the diffusive spreading of the probability density. The moving-boundary formulation has been addressed in the mathematical literature [29]- [31].
Our main interest here concerns the one-dimensional (1D) diffusion problem with a static reflecting boundary and with a harmonically oscillating external driving force. In our previous paper [32], the analysis of the same setting has been based on the numerical solution of the emerging integral equation (cf equation (3) below). The equation is not of the convolution type and it has a weakly singular kernel. Employing the product-integration method [33,34], the approach was oriented towards the description of the initial transitory stage of the evolution. On the other hand, the transitory effects are seldom interesting-the physically relevant features emerge in the time-asymptotic (still periodically time dependent, but after the decay of transitory effects) regime. We now propose a new method which is directly focused on the long-time dynamics per se and which reveals several of its non-intuitive features.

Green function for the restricted diffusion
We are interested in the solution of the Langevin equation for an overdamped Brownian particle moving in the time-dependent potential V(x, t), where V(x, t) = −xF(t), if x 0, and V(x, t) = ∞, for x < 0. Here F L (t) is the δ-correlated Langevin force, and equals the particle mass times the viscous friction coefficient. Differently speaking, while localised on the positive half-line, the particle is acted upon by the Langevin force F L (t) and by the spatially-homogenous, time-dependent force F(t). Additionally, we assume a reflecting barrier at the origin, i.e. the diffusion is restricted on the positive half/line. As an auxiliary problem, consider first the spatially unrestricted 1D diffusion in the field of a spatially-homogenous and time-dependent force F(t). The time-evolution of the probability density for the position of the diffusing particle reads [21]- [23] G(x, t; x , t ) The Green function yields the solution of the Smoluchowski diffusion equation for the initial condition π(x) = δ(x − x ) imposed at time t . Qualitatively, it represents the gradually spreading Gaussian curve whose centre moves in time, the drift being controlled by the protocol (timedependent scenario) of the external force. The momentary value of the mean particle position is given by the expression x + t t ds v(s), where v(t) = F(t)/ is the time-dependent drift velocity. The spreading of the Gaussian curve is controlled by the thermal-noise strength parameter Returning to our central problem, we now assume that the particle is initially (i.e. at time zero) fully localized at a fixed point x > 0, and we place at the origin of the coordinate system the reflecting boundary. The Green function which corresponds to the problem with the reflecting boundary, say U(x, t; x , 0), can be constructed in two steps (cf the detailed derivation in our previous paper [32]). First, one has to solve the Volterra integral equation of the first kind Here both the kernel and the right-hand side follow directly from equation (2). The integral equation should be solved for the unknown function U(0, t; x , 0) which represents, as the designation suggests, the time evolution of the probability density for the restricted diffusion at the boundary. Therefore, the solution must be a non-negative function of time. It represents the first orientation when considering the full space-resolved probability density in the restricted diffusion problem. As a matter of fact, having obtained the probability density at the boundary, the final space-resolved solution emerges after performing just one additional quadrature. We have proved the formula [32] U(x, t; The resulting function is properly normalized, i.e. we have ∞ 0 U(x, t; x , 0) dx = 1 for any t 0 and for any fixed initial position x > 0.
Up to now, our reasoning is valid for any form of the external driving force F(t) = v(t). A negative instantaneous force pushes the particle to the left, i.e. against the reflecting boundary at the origin. In this case, the force acts against the general spreading tendency stemming from the thermal Langevin force. A positive instantaneous force amplifies the diffusion in driving the particle to the right. In the present paper, we focus on the case of a harmonically oscillating driving force. Specifically, we insert into equation (2) The three parameters, v 0 , v 1 , and ω occurring in this formula together with the diffusion constant D represent the full description of our setting and they form a firm frame for the complete discussion of the results. Having a static force, i.e. if v 1 = 0, the explicit solution of the integral equation (3) is well known, cf the formula (29) in [32]. U(0, t; x , 0) approaches in this case either zero, if v 0 0, or the value |v 0 |/D, if v 0 < 0. Having the oscillating force, the most interesting physics emerges if the symmetrically oscillating component superposes with a negative static force, i.e. if v 1 > 0, and v 0 < 0. This case is treated in the rest of the paper.

Probability density at the boundary
Considering the integral equation (3), the basic difficulty is related with the non-convolution structure of the integral on the left-hand side. It may appear that any attempt to perform the Laplace transformation must fail. But we shall show that this need not be the case. We start with an appropriate scaling of the time variable. After any such scaling, the four model parameters will form certain dimensionless groups. However, there are just two 'master' combinations of the four model parameters which control the substantial features of the long-time asymptotic solution. These combinations emerge after we introduce the dimensionless time τ = [v 2 0 /(4D)] t (we assume D > 0). If we insert the scaled time into equation (2), the exponent will include solely the combinations κ = |v 0 |v 1 /(2ωD) and θ = 4ωD/v 2 0 . The first of them measures the scaled amplitude of the oscillating force, the second one its scaled frequency.
Another remark concerns the square in the exponent on the right-hand side of equation (2). We can get rid of the square at the price of introducing the Fourier transformation of the whole expression. Accomplishing the transformation, the Green function (2) reads where τ is defined in the text above. In the next step, we make use of the well known formula exp [β cos (α)] = +∞ k=−∞ I k (β) exp (ikα) [35], where I k (β) is the modified Bessel function of order k with argument β. In other words, we expand the expressions exp [−iuκ cos (θτ)] and exp [+iuκ cos (θτ )] into Fourier series. As a result, we acquire a double infinite summation involving products of the modified Bessel functions I m (−iuκ) I n (+iuκ), and also two exponential factors exp (imθτ), exp (inθτ), m, n = 0, ±1, ±2, . . ..
After these preparatory steps, we insert the function G(0, t; 0, t ) into the integral equation (3). Two important observations are in order. Firstly, for any v 0 < 0, the right-hand side of the integral equation (3) approaches unity. Hence we can represent it as 1 +r(x , t), with lim t→∞r (x , t) = 0. Secondly, after the transitory components die out, the probability density at the boundary will lose its dependence on x and it will oscillate with the fundamental frequency ω. Therefore we can write U(0, t; x , 0) = (|v 0 |/D)[f a (t) +f (x , t)], where again lim t→∞f (x , t) = 0. The reduced time-asymptotic form of the probability density at the boundary f a (t) can be represented by the Fourier series Here and below, the Laplace-transformation variable z already corresponds to the scaled-time variable τ. The complex numbers f k , k = 0, ±1, ±2, . . . , will be referred to as the complex amplitudes. If κ = 0, we must get f k = δ 0k , where δ mn is the Kronecker delta. We now perform the Laplace transformation of the integral equation (3). The exponentials exp (imθτ) and exp (inθτ ) will induce shifts in the Laplace variable z. The resulting equation reads where z m = z − imθ for an integer m. On the left-hand side, the integrated function possesses first-order singularities at u ± (z m ) = ±i[ √ 1 + z m ∓ 1]. Using the residue theorem, we can now carry out the u-integration. We close the integration path in the upper half of the complex 6 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT u-plane and take into account the singularities within the integration contour, i.e. at u + (z m ). Accomplishing the integration, we arrive at the relation We are now able to recognize the terms which survive in the time-asymptotic limit. We multiply the above equation by the Laplace variable z and take the limit z → 0. As noted above, the functionf (x , t) relaxes to zero. Hence we have lim z→0 zf (x , z m+n ) = 0 for any two integers m and n. Similarly, on the right-hand side, we apply lim z→0 zr(x , z) = 0. The limit lim z→0 zf a (z m+n ) yields the complex amplitude f k=m+n . Otherwise, e.g. in the arguments of the modified Bessel functions, we put z = 0. Summing up, in the time-asymptotic limit, the integral equation (3) implies the difference equation for the unknown complex amplitudes f k . In order to get further insight into the structure of this difference equation, we introduce three infinite matrices L +− , D, and R −+ , their matrix elements being where m and n are integers (in these expressions and below we use the bra-ket notation, to make the algebraic operations more transparent). Notice that the matrix elements depend solely on the dimensionless combinations κ and θ. Invoking the substitution m + n = k in the summation index n and using the above matrices, equation (10) assumes the compact form The matrix D is diagonal and its diagonal matrix element which corresponds to the zeroth row equals one. As for the matrices L +− and R −+ , the modified Bessel functions along a diagonal which is parallel with the main diagonal have the same index. Notice, however, that the arguments of the modified Bessel functions in the matrix L +− are specified solely by the column index. Further, the argument of the modified Bessel functions in the zeroth column of the matrix L +− is zero and hence all matrix elements in this column vanish except of the zeroth diagonal element which equals to one, k|L +− |0 = δ 0k . In other words, we can write L +− |0 = |0 . Similarly, the arguments of the modified Bessel functions in the matrix elements of the matrix R −+ are specified by the row index. The zeroth row is simply 0|R −+ |k = δ 0k , i.e. we can write 0|R −+ = 0|.

7
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT We are now approaching the last step of our heuristic procedure. Assuming the uniqueness of the solution, the complex amplitudes must necessarily be given by the expressions Actually, substituting the proposed solution back into equation (14), the proof can be completed by invoking the resolution of the identity matrix in the form +∞ k=−∞ |k k| = I. Let us now consider the zeroth column of the inverse matrix L −1 +− . Forming the fractions of the algebraical complements and the determinant of the matrix L +− , we readily ascertain that the zeroth column of the inverse matrix L −1 +− equals that of the original matrix L +− . A similar observation holds true for the inverse matrix R −1 −+ . Its zeroth row is simply 0|R −1 −+ = 0|. Using these properties in equation (15), the complex amplitudes are given by the expressions that is, they form the zeroth column of the matrix R −1 −+ . Moreover, the zeroth complex amplitude f 0 equals one. Analysing the matrix elements of the inverse matrix, one can prove that the amplitudes f k and f −k are complex conjugated numbers. Generally, their absolute value |f k | decreases with increasing the index k. The even (odd) amplitudes are even (odd) in the parameter κ.
Let us summarize the main results of the present section. Asymptotically, the probability density at the origin reads Here the amplitudes of the higher harmonics and their phase shifts relate to the complex amplitudes through the standard formulae We have emphasized their dependence on the parameters κ and θ. Except for the multiplicative factor |v 0 |/D, the asymptotic form of the probability density at the boundary is controlled solely by these two parameters. For example, changing the diffusion constant D and, at the same time, keeping a constant value of the product Dω, the time-asymptotic form of the reduced function f a (t) = (D/|v 0 |)U a (0, t) will not change. For any value of the parameters κ and θ, the time average of the probability density at the boundary equals the equilibrium value of this quantity in the problem without driving force. In symbols, we have Turning to the numerical analysis of the results, the complex amplitudes (16) have been calculated via a direct numerical inversion of the matrix R −+ defined in (13). Of course, the infinite-order matrix R −+ must be first reduced onto its finite-order central block. The matrix elements of the reduced matrix are again given by equation (13); presently, however, m, n = 0, ±1, ±2, . . . , ±N. The integer N has been taken large enough such that its further increase does not change the results, within a predefined precision. In this sense, all numerical results in the present paper represent the exact long-time solution of the problem in question. Figure 1 illustrates the reduced probability density at the boundary in the asymptotic regime together with the amplitudes and phase shifts as introduced in equation (17). We have chosen in the following section a particular set of parameters which yields interesting effects.

Space-resolved probability density
Up to now, we have only discussed the time dependence of the probability density at the boundary. As a matter of fact, the knowledge of the complex amplitudes f k allows for a rather detailed discussion of many other features of the emerging diffusion process. First of all, in this section, we focus on the time-and space-resolved probability density for the particle coordinate. We 9 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT recall that, regardless of the initial condition, a time-independent drift towards the origin induces the gradual constitution of the unique equilibrium density π eq (x) Assuming again the oscillating drift v(t) = v 0 + v 1 sin (ωt) with v 0 < 0 and v 1 > 0, we are primarily interested in the time-asymptotic dynamics. In the regime, the probability density U(x, t; x , 0) does not depend on the initial condition (as represented by the variable x ), and it exhibits at any fixed point x 0 oscillations with the fundamental frequency ω. Similarly as in the preceding section, we expand the asymptotic density into the Fourier series Presently, however, the Fourier coefficients u k (x) depend on the coordinate x. An interesting quantity will be the time-averaged value of the density in the asymptotic regime. This is simply the dc component of the above series and can be written as the time average We already know that the value of this function at the origin is u 0 (0) = |v 0 |/D, i.e. it equals the value of the equilibrium density at the origin, u 0 (0) = π eq (0). In other words, if we generically call the difference between the stationary value of a quantity and the corresponding equilibrium value of this quantity as 'dynamical shift', there is no dynamical shift of the density profile at the origin. But what happens for x > 0? Assume that the complex amplitudes f k are known. Hence we know also the time-asymptotic solution of the integral equation (3) and the subsequent analysis can be based on the expression (4). The first term on the right-hand side of equation (4) does not contribute to the stationary regime. Actually, for any fixed x, the function G(x, t; x , 0) exhibits damped oscillations and vanishes finally . As for the second term on the right-hand side of equation (4), we use the explicit form (6) and perform the required x-derivative. Further, we insert into this term the solution of the integral equation equation (3) in the form U(0, t; x , 0) = (|v 0 |/D)[f a (t) +f (x , t)], where f a (t) is given in equation (7). We now proceed along similar lines as in the preceding section. Namely, we introduce the modified Bessel functions and carry out the Laplace transformation. Afterwards, the second term on the right-hand side of equation (4) yields the expression where z m and u ± (z m ) have been introduced below equation (8 complex u-plane. Inside the integration contour, the singularities of the integrand are located at . At the same time, the above u-dependent factor must be taken into account when calculating the residues; it eventually produces the matrix elements (24) below. After the u-integration, we can already identify the terms which survive in the stationary regime.
Leaving out the details, the x-dependent Fourier coefficients in equation (20) are given by the expression Here |f is the column vector of the complex amplitudes, i.e. f k = k|f . Moreover, we have introduced the diagonal matrix E(x) with the complex, x-dependent diagonal elements and the two matrices L −− , R ++ with the matrix elements where m and n are integers.
Let us now focus on the numerical analysis of the main result (23). As described in the preceding section, we have to reduce both the vector of the complex amplitudes |f and the matrices L −− , E(x), and R ++ . Except for this controllable approximation, we can already reconstruct the full time-and space-resolved nonlinear 'waves' of the probability density, as represented by the function U a (x, t). We simply use the Fourier coefficients u k (x) in equation (20). Figure 2 illustrates the time-asymptotic density within two periods of the external driving force. Surprising features emerge provided both κ 1, and θ 1. Under these conditions, the timeaveraged probability density u 0 (x) exhibits in the vicinity of the boundary a strong dependence on the x-coordinate. It can even develop a well-pronounced minimum close to the boundary and, simultaneously, a well-pronounced maximum localized farther from the boundary. In between the two extreme values, there exists a spatial region where the time-averaged gradient of the concentration points against the time-averaged force. The situation is depicted in figure 3 where we have used the same parameters as in figures 1 and 2.

Dynamical shift of the mean position
describes the motion of the 'centre of mass' of the diffusing particles or, equivalently, the mean position of a representative particle within its probabilistic description. This function is intimately related to the mean energy of the particle, to the work done by the external agent, and to the heat released to the reservoir [15,32]. In the static case, i.e. if v 1 = 0, the function µ(t, x ) either diverges, if v 0 0, or approaches the value µ eq = D/|v 0 |, if v 0 < 0. The oscillatory driving force keeps the particle permanently in a non-equilibrium state and, of course, there is no time-asymptotic value of the mean position. In the asymptotic regime, we can again write and we can first inquire what is the stationary mean position In the following, the difference σ = µ 0 − µ eq will be referred to as the dynamical shift of the mean coordinate. It is depicted in figure 3 by the horizontal arrow and it is also illustrated in figure 4 for a wider range of the model parameters. The main result of the present section expresses the dynamical shift through the complex amplitude f 1 ; we shall show that σ = (v 1 /ω) Ref 1 . Except for a trivial multiplicative factor, the shift depends on the combinations κ and θ. We start the proof by inserting the Green function (4) into the definition (27). Hence we have The first term on the right-hand side is simply the (conditional) mean position for the unrestricted diffusion, under the condition of x ≥ 0. This term exhibits damped nonlinear oscillations and it finally vanishes. Actually, considering the unrestricted diffusion (without boundary) with a negative value of the parameter v 0 , the particle will eventually glide towards −∞. As for the second term on the right-hand side, it is tempting to carry out the x-integration in the curly brackets. However, this step would preclude any further time-asymptotic analysis of the ensuing non-convolution time integral. Instead, we shall treat the curly-bracketed expression in equation (30) using the steps delineated in section 3. Namely, we first use the expression (6) with x = 0 and then perform the x-derivation and x-integration in the curly brackets in equation (30). During the x-integration, we assume that the variable u in equation (6) has an infinitesimal negative imaginary part. Secondly, we perform the Laplace transformation. After these two steps, the second term in equation (30) yields the expression where z m and u ± (z m ) have been introduced below equation (8). Notice the additional singularity of the integrand at u = i . Precisely speaking, the variable u in the denominator represents a where the prime denotes the excluded summand with m = n = 0. The last expression forms the starting point for the time-resolved asymptotic analysis of the mean position (27). However, we are primarily interested in the time-averaged mean position (29). In other words, we need to identify the terms which survive if we multiply the expression (32) by the variable z and then take the limit z → ∞. A careful evaluation of the limit yields the result where f 1 is the complex amplitude given in equation (16). The proof is completed by subtracting the static equilibrium mean position µ eq = D/|v 0 |, and using the definition κ = |v 0 |v 1 /(2ωD).
Notice that the small-v 1 expansion of the dynamical shift σ = (v 1 /ω) Ref 1 starts with the term v 2 1 . In the case where we had used a variant of the linear response theory, the dynamical shift would escape our attention.

Concluding remarks
Returning to the physical setting as mentioned in the introduction, the approach elaborated in the present paper yields a rather complete picture of the time-asymptotic motion of the suspended particles. Depending on their distance from the impenetrable boundary, they exhibit non-harmonic oscillations which can be represented as a linear combination of several higher harmonics. The amplitudes and the phases of the harmonics are strongly sensitive to the distance from the boundary. Our approach is directly focused on the properties of the long-time regime and it is not based on any small-parameter expansion. In our numerical analysis, we have tested many combinations of the four governing parameters v 0 , v 1 , D, and ω. The most interesting dynamics emerges if, simultaneously, the following three condition apply: κ 1, θ 1, and κθ 1. Here again κ = |v 0 |v 1 /(2ωD), θ = 4ωD/v 2 0 , i.e. κθ = 2v 1 /|v 0 |. In other words, we assume the strong driving, low temperature and/or slow modulation, and a small regressive force F 0 = v 0 . In this domain of the model parameters, the formation of a spatial region exhibiting an increased particle concentration, as discussed in the fourth section, represents an interesting feature of the emerging stationary regime. Within our formulation, the feature is an exact, purely dynamical consequence of the governing Langevin equation. It is an interesting question whether the feature implies specific extremal properties of e.g. the heat dissipated over one period of the time-asymptotic non-equilibrium dynamics.
Although we have discussed only purely kinetic quantities, several words concerning the energetics are in order [32]. In the time-asymptotic regime, the system exhibits periodic changes of its state. With the time-dependent potential V(x, t) = −xF(t), the internal energy is also periodic. The work done on the system during one period must be equal to the heat dissipated during the period. But an interesting quantity is the stationary value of the internal energy It is obviously controlled by the first and the second Fourier coefficient in the expansion (28). We can show that E 0 is always bigger then the equilibrium internal energy E eq = D = k B T . In other words, in the time-averaged sense, the external driving enforces a permanent increase of the internal energy, as compared to its equilibrium value. The external agent does work on the system by increasing the potential V(x, t), while the position of the particle is fixed. Thus the work done at a given instant depends on the momentary position of the particle. Summing over all possible positions and over the time period, we get [32] This quantity must be positive. Otherwise, the system in contact with the single heat bath at one fixed temperature would produce positive work during the cyclic process in question. The average work equals the area enclosed by the hysteresis curve which represents the parametric plot of the oscillating force F 1 (t) versus the mean coordinate µ a (t). Hence the work done on the system during a definite time interval within the period can be both positive and negative. The work is positive (negative) within two quarter-periods when the slope of the potential increases (decreases). But its actual value depends on the localization of the particles, within a given quarter-period. The absolute value of the work is bigger if the particles are farther from the boundary.
Turning to a different issue, it is sometimes argued [28] that, in the slow-driving limit, one can invoke the so-called adiabatic approximation. One first calculates the equilibrium quantities for the a static force F . Afterwards, one performs the substitution F → F(t) directly in the results of the previous equilibrium calculation. Let us translate this prescription literally into our setting. In the first step, we consider the equilibrium density at the boundary: it is zero, if F 0, and it assumes the value |F |/( D) = |v|/D, if v = F/ < 0. Following the prescription for the adiabatic approximation, we should write The resulting function equals zero in time intervals where v 0 + v 1 sin (ωt) > 0. If |v 0 | v 1 , this case never occurs and U (adia) a (0, t) harmonically oscillates around the value |v 0 |/D. Hence its time average is again |v 0 |/D. On the other hand, if |v 0 | < v 1 , there exists a time interval within the period, where the instantaneous drift velocity assumes positive values. Within this interval, the function U (adia) a (0, t) vanishes. Its averaging over the period yields a value which exceeds |v 0 |/D. Recall that the exact treatment in the third section gives the value |v 0 |/D irrespective of the magnitude of the parameters v 0 and v 1 . Here we observe a qualitative failure of the adiabatic approximation. The disagreement exists for any value of the driving frequency ω. Formally speaking, the observation is related with the fact that the limits t → ∞ and ω → 0 + are not interchangeable. The proper slow-driving time-asymptotic dynamics emerges if we order the two limits as lim ω→0 + lim t→∞ . . . . On the other hand, the above described adiabatic 'results' are based on the double limit lim t→∞ lim ω→0 + . . . . More physically, if the frequency is small but still finite, the exact solution still contains all higher harmonics. Their summary effect can survive even in the limit ω → 0 + and it still influences the time-asymptotic periodic motion. For example, all higher harmonics are needed to guarantee, that probability density at the origin is positive for any time.
In our last remark, we would like to attract the attention of the mathematically oriented reader to the 'magic' matrixes (11), (13), (25), and (26). A retrospection of the present paper reveals a possible improvement of our procedure. Namely, instead of resorting to numerical inversion of the matrix R −+ (remind that the inverse is an important step as it yields, via equation (16), the complex amplitudes f k ), one could consider the problem of an analytic calculation of the