d orbitals PDOS analysis

I believe there could be two reasons for a single d-orbital to have multiple peaks in the density of states plot:

1. Hybridization: Because the metal's (Ti) d-orbitals are hybridizing with ligands' (most likely) p-orbitals.

   Just like mentioned @Andrey Poletayev's answer, this will result in the metal's d-orbitals contributing to both the bonding and anti-bonding states. However, I'd like to add that the degree of the contribution will differ depending on the polarity of the bond. For example, for pure covalent bonds, the contribution of the d-orbitals to the bonding and anti-bonding states shall be completely identical. For ionic (polar) bonds, these contributions will be largely different and biased.

   Your plot shows that $d_{z^2}$ orbital has two major peaks (up spin channel), one at ~0.8 eV and the other at ~2.5 eV, and their amplitudes are almost identical. Considering your system is most likely ionic (or with polar covalent bonds), this rules out the hybridization explanation.

2. Rotation: The metal's d-orbitals are rotated and are no longer aligned with the Cartesian axis.

   This point is more subtle as I rarely see anyone talking about it: The d-orbitals are not always aligned with the Cartesian axis. They can be rotated in a way that they are no longer aligned with the Cartesian axis due to the environment (tilting or rotated octahedral for example).

   Since most DFT codes project the calculated wavefunctions onto a set of atom-centered real spherical harmonics that are aligned with the Cartesian axis, the results of a rotated frame will become harder to interpret.

   For example, for p-orbitals, if the $p_z$ orbitals are rotated by 45 degrees along the polar angle by the y-axis, the projected wavefunctions will be a linear combination of the $p_x$ and $p_z$ orbitals. This will result in a density of states plot that has two $p_z$ peaks instead of one, and the contribution of $p_z$ (which is reflected by its amplitude) should be equally shared between these two peaks.

   For d-orbitals, the mixing would become more complex, but the idea is the same. If you want, you can create the transformation matrix for d- (or even f-) orbitals and apply that to the Cartesian axis projected DOS to get the rotated version.

One thing to note is that DFT codes usually output the square of the projection coefficients but for the rotation to work we need the original projection coefficients.

I actually wrote a code that can help with this: https://github.com/Chengcheng-Xiao/RotSph. Please see https://github.com/Chengcheng-Xiao/RotSph/tree/master/example/manual_rotation for an example of how this code works. Also, check out the automated version that automatically finds the best Euler angles for your atomic orbitals' local frame.

So the orbitals we are getting below 0 eV, we should not get above 0 eV. But in my plot, we are getting a $d_{z^2}$ orbital peak both below the Fermi level and above the Fermi level.

This question is not really about the DOS, but about chemical bonding. Grouping $d$ orbitals by symmetry into $t_{2g}$ and $e_g$ is a consequence of a crystal field of octahedral symmetry in crystal field theory. But crystal field theory is not sufficient to describe bonding of the (typically metal) center with its ligands. To answer this question, we need to take the next step from crystal field theory to ligand field theory.

Ligand field theory considers the bonding between the (metal) center and its ligands. To do so, the ligand orbitals (for example, oxygen $p$ orbitals in oxides, or similarly both $p$ and $\pi$ orbitals of ligands such as CO and CN) are also grouped by symmetry. Bonding will be possible between metal and ligand orbitals of matching symmetry. For the simplest example, $e_g$ orbitals can participate in $\sigma$ bonding, and $t_{2g}$ in $\pi$ bonding. When the metal orbitals and the ligand orbitals of matching symmetry hybridize, two sets of molecular orbitals are created. Now, crucially, both of those sets of molecular orbitals have partial contributions from both the "original" metal orbitals and the "original" ligand orbitals.

For example, if normal sigma bonding occurs, then metal $e_g$ orbitals will make partial contributions to both (usually filled) $\sigma$ and (usually empty) $\sigma^*$ molecular orbitals. See for example the second schematic here for a typical molecular orbital diagram. Stronger sigma bonding will change the energies of the molecular orbitals with (partial) $e_g$ character. Pi bonding will change the energies of molecular orbitals with (also partial) $t_{2g}$ character. See for example Figure 6.4.2 here. So the final energies at which contibutions of specific $d$ orbitals show up in the DOS depend not only on the crystal field symmetry, but also on the kind and strength of bonding between the metals and their ligands.

Then, these molecular orbitals are filled with electrons "originating" from both the ligand and the metal, which will together determine where the Fermi level falls in the extended solid.

This is in fact what happens in the DOS and PDOS in the original question. It looks like your $d_{z^2}$ is making partial contributions to filled bonding orbitals below the Fermi level and empty antibonding orbitals above the Fermi level. That is why the $d_{z^2}$ is showing up in two places on the DOS rather than in one place as crystal field theory would have predicted. The key piece of information missing from the PDOS as plotted is whether the states to which the $d$ orbitals are contributing are in fact 100% composed of those $d$ orbitals. Unless the $d$ orbitals are non-bonding, that will not be the case. You can overlay your ligand orbitals on the DOS plot together with the metal ones, and you should see some overlap for most of the orbitals you plotted.

But in short, crystal field theory does not correctly describe the overall system and its DOS and leads to erroneous intuition.

First of all, since you already have $d_{xy}$, $d_{yz}$, $d_{xz}$, $d_{z}^2$, and $d_{x^2-y^2}$ PDOS data, you can also plot $t_{2g}$ and $e_g$ resolved PDOS. \begin{equation} t_{2g} = d_{xy} + d_{yz} + d_{xz} \tag{1} \end{equation} \begin{equation} e_g = d_{z}^2 + d_{x^2-y^2} \tag{2} \end{equation} So, either create two new columns using the above definition and plot them, or you can just plot their sum and label them accordingly. Secondly,

Moreover the main thing you see here, as we know the below Fermi level all states are occupied,and above is unoccupied. So which orbital we are getting below ta 0eV, we should not get above the 0eV.

I think your this statement is not entirely correct. Indeed below Fermi levels, all the energy states are occupied and vice-versa. However, this does not mean that orbitals having allowed energy levels below the Fermi level should not have allowed energy levels above the Fermi level. They can have which is just not occupied. If you check PDOS plots from literature, you will notice that not only orbitals have energy levels both below and above Fermi levels, but they can have peaks both below and above.

I am not saying your PDOS plot is correct! There might be other issues and without knowing how you obtained this, its correctness cannot be affirmed but the issues you raised are not of concerns. So, if the methodology is correct, your PDOS plot is good I believe.

Yes we can get states(those states are corresponded to antibonding state) above the fermi level where electrons can be occupied when we will increase the temperature. Lastly, my PDOS is also correct according to my starting structure, but now the problem is my started strcture contains 'tilted octahedrals' with respect to unit cell 'x,y,z' axis (did calculation in VASP), Thats why I am not getting correct $m_l$ projected DOS. So I have to rotate these octahedrals locally along $z$ axis of unit cell to get corrected ligand positions.