There is no mathematical proof for a general form. In the elasticity context see discussion in the books by Barber, Constantinescu & Korsunsky, etc.

Thanks. I have a big problem from here. I do not have access to many literatures in particular to these books. I have found on the net something I am reading at the moment on 3D twistors and biharmonic equation but I have to understand first before saying anything else. Thanks for above reference.

Bi in biharmonic function indicates '2'. So, perhaps the name should be extended too. Anyway, In the sense of modeling  , perhaps the form of the pde of harmonic 3d can be as my attachment here. 

Yes this is the biharmonic equation in Cartesian coordinates. We are interested in a 3-dimensional function f, of general form, that satisfies above equation.

If you allow integrals in your sought "most general available expression", can one use an electrostatic analogy, applied twice, in order to obtain a rather general expression for the biharmonic functions in the given domain D?

Let's pick an arbitrary scalar function rho(r) ("charge density"), where r=(x,y,z), such that (1) rho(r) = 0 inside D and (2) the integral ("electric field potential")  over the whole space U1(r) = \int rho(r') dr' / |r - r'| exists for any r.

Note that  ∆U1(r) = 4*pi*rho(r) (Poisson equation in electrostatic) and therefore U1(r) is a harmonic function inside D.

Repeat the same procedure again, by calculating the following integral over the whole space U2(r) = \int U1(r') dr' / |r - r'| . If this integral exists, ∆U2 = 4*pi*U1 and ∆∆U2 = 0 inside D.

As a result, U2(r) is a biharmonic function inside D, parameterized by a broad(?) set of "charge density" functions rho(r), for which both above integrals exist.

In particular, to simplify the calculations, rho(r) can be selected to be a combination of charge distributions, for which the first integral can be calculated analytically - e.g. a combination of point charges, point dipoles, point multipoles, uniformly charged straight lines, planes, etc.  

The domain D in my study is the whole unbounded space .

Excuse me, why " If this integral exists, ∆U2 = 4pi*U1"?

Ok, I understand my last question because ∆(1/ |r - r'|) behaves as a Dirac delta function.

My study involves infinite dislocation and interface in elasticity.

Yes, ∆U2 = 4*pi*U1 , because U2(r) = \int U1(r') dr' / |r - r'| . One way to explain this is to refer to the delta function relation, as you did. A less formal way is to notice that  U2(r) = \int U1(r') dr' / |r - r'| means that U2 is an "electric field potential" corresponding to the "electric charges" distributed with the "charge density" U1 . As a result, the "electric field potential" and "charge density" satisfy the Gauss Law  ∆U2 = 4*pi*U1 (in Gauss units).

More importantly, a similar, but much simpler, method to parameterize the biharmonic functions in the bounded 3D domain D can be obtained using the known expression for the fundamental solution for the biharmonic equation in 3D:  F(r, r') = | r - r' | / (8*pi)   ( http://en.wikipedia.org/wiki/Fundamental_solution#Biharmonic_equation )  

Unfortunately, the methods to parameterize the biharmonic functions based on "double electrostatic analogy" or the known fundamental solution are not useful for parameterizing the biharmonic functions in the entire space. The parameterizing "charge density" rho(r) needs to be zero inside D, so if D is the whole space, the only biharmonic function yielded by this parameterization can be identical zero.

For the entire 3D space, one broad family of biharmonic functions are the functions of just 2 Cartesian coordinates derived from a pair of complex analytic functions - see http://en.wikipedia.org/wiki/Biharmonic_equation#2-dimensional_space  

Again, I am reading at the moment on 3D twistors and biharmonic equation but I have to understand first before saying anything else. Thanks for above comments

Please let me add another point. Could someone suggest a possible extension of biharmonic functions based on the form I have provided in the Q&A, i.e. (... + a/(k.x) + a0 + a1(k.x) +a2 (k.x)2 +...) exp(i k.x)? This form seems convenient according to my first calculations.

Note that: (k.x)2 reads (k.x) to the power two.

Dear Colleague,

I have a file of book with general form of solution of polyharmonic equation. I hope it will be useful (especially, if you have difficulties with the literature access). Unfortunately, it is in Russian (there are English editions but I have no English file ), but formulas are international. See page 160, formulas (5.31)-(5.34). It is general solution for domain D with smooth boundary S,   \sigma_i  are, in general, arbitrary (it can be determined from the boundary conditions (i.e. from integral equations (5.34)).

Patrick, what does "k.k=0" in your original question mean?

Also, since your formula includes the complex numbers ai and also exp(ik.x), what makes its value real-valued? Do you also take Re or Im of the value given by this formula? 

In addition to my previous answer, I can add file of my paper (in English) with the corresponding theorem  (p. 19, Lemma 3)

Regarding the reference to the book in Russian, posted above by I.Yu.Popov, equations (5.31)-(5.34) on p.160 deal with the boundary value  problem (5.34) for the polyharmonic equation (5.32), which becomes biharmonic when m=2. The text suggests a set of functions (5.31) (potentials of the dipoles/multipoles placed at the surface of D?), satisfying (5.32) inside D. Each of these functions can be parameterized by a "density" function \sigma_i defined on the boundary of D. The functions (5.31) can be combined linearly as (5.33), resulting in a polyharmonic function parameterized by "densities" \sigma_i, with the latter defined on the surface of D. The boundary value problem then reduces to a set of integral equations for these unknown "densities".

This reads like an interesting generalization of the well-known approach reducing the boundary-value problem for the Laplace equation to integral equations for the functions defined on the surface of the domain.

However, it is not clear how can this help with the original question, seeking the general expression for biharmonic functions in the entire 3D space.    

Thank you all for above comments.

1) To Yu. Popov: I will first read the given files and then give my comments

2) To Azat Latypov: The form given in the Q&A will be first integrated with respect to some components of k and then equated to a real-valued expression with the aim of matching boundary conditions. k.k = 0 means that the three components of k are not independent.

Let f be a spatial function depending on vector position x = (x1, x2, x3). We assume that f is biharmonic : ∆∆f = 0, ∆ Laplace operator. We seek a general form for f. We express f in its Fourier form: f(x) = ∫α(k) exp(ik.x) dk, where k = (k1, k2, k3), so that ∆∆f(x) = ∫α(k) ∆∆exp(ik.x) dk. It becomes evident that ∆∆exp(ik.x) =0 implies that f is biharmonic; this occurs when k12 + k22 + k32 = 0.We assume that, for the problem under consideration, it is convenient to express k2 as a function of k1 and k3. Consequently, we can manage f to read

f(x) = ∫α(k) exp(ik.x) dk1 dk3

where k = (k1, k2, k3) with k2 =± i√(k12 + k32) .

      This is a "general form" for a three-dimensional biharmonic function that can be used as components of Galerkin vectors. This form is very effective in problems where boundary conditions have to be satisfied (boundaries with infinite sizes). By writing the boundary conditions in Fourier forms, α(k) may be determined.

      We are using successfully the form above in work with title: A study of the elastic fields of interfacial dislocations (straight and sinusoidal) using Galerkin vectors with three-dimensional biharmonic functions in Fourier forms (see our contributions in Research Gate).

     I have not investigated, at the moment, whether the forms provided in the book by Yu Popov apply in the entire 3-D space.

Again the form of three-dimensional biharmonic function above permitted to complete a recent work, entitled “Non-planar interface crack under general loading III. Dislocation, crack-tip stress and crack extension force” , in our contributions in ResearchGate. In this work, expressions for crack-tip stresses and crack extension force are provided for a model of non-planar interface crack in bi-elastic solids with oscillatory front ξ= ξ (x1, x3) spreading in x2x3-plane perpendicular to fracture propagation x1-direction. The average crack plane is Ox1x3, perpendicular to the applied tension x2- direction. The shearing stresses are applied in the x1 and x3 directions, parallel to Ox1x3. Also considered are induced shearing and normal stresses that originate from Poisson contractions in the media. Any contribution is welcome through comments.