Double beta decays comprise decays of radioactive nuclei where two neutrons decay into two protons, or vice versa. Depending on whether neutrinos are emitted concurrently, the decays can be classified as two-neutrino (2νββ) or neutrinoless (0νββ). As second-order processes, the only situation in which double beta decays represent a significant decay branch is the decay of even-even nuclei[1], where a single beta decay of an even-even parent nucleus, stabilised by pairing, into a less stable odd daughter nucleus is forbidden[1]. This means that there is a small number of experimentally verified nuclides that decay through 2νββ decay and a small number of candidates for the experimentally unobserved 0νββ decay. Once considered of little significance[2], double beta decays have become an object of considerable theoretical interest in recent years[3, 4]. The revival of interest in the topic of double beta decays is chiefly due to the connection of the 0νββ decay mode to physics beyond the standard model, particularly to the question of the nature of the neutrino (namely whether the neutrino is a Majorana or a Dirac particle)[5]. This connection is proven to hold regardless of the underlying physics by the Schechter-Valle "black box" theorem [6]. Despite the importance of neutrinoless decays, the 2νββ remains important as the only double beta decay mode we have experimental evidence of, despite significant experimentalist effort to find a signal of neutrinoless double beta decay[7, 8, 9, 10, 11, 12, 13], and as such is indispensable for developing and benchmarking theoretical frameworks for double beta decay. The present work addresses both modes of double beta decay, in the framework of a protonneutron (or charge exchange) quasiparticle random-phase approximation (PN-RQRPA)[14] based on relativistic nuclear energy density functional[15, 16]. This work presents the first calculation of the nuclear properties of double beta decay based on the PN-RQRPA and as the first implementation of our theoretical framework represents a starting point for further studies within the relativistic framework. The most important quantity to be calculated, from a nuclear structure standpoint, is the nuclear matrix element. Within approaches based on the QRPA, this is achieved by using data (transition matrix elements, QRPA amplitudes, occupation numbers and energies) from two single beta decay calculations, starting from the initial and from the final nucleus. The difference in the final states for both, representing the intermediate nucleus, is taken into account with an overlap factor calculated from occupation factors and QRPA amplitudes for both single beta calculations. In the present work we use both the most widespread form for the overlap factor[17] and a more involved form due to Simkovic[18], however by comparing the overall values of the NME we find that the differences are small. The PN-RQRPA calculation is built on top of a relativistic mean-field theory, which in the present work is realised as either relativistic Hartree-BCS or relativistic Hartree-Bogoliubov model. The relativistic Lagrangians associated with the nuclear energy density functional include three density dependent interactions: the DD-ME2[19] meson-exchange interaction based on quantum hadrodynamics[20] and the point coupling interactions DD-PC1[21] and DD-PCX[22]. The use of different interactions allows us to estimate systematically the theoretical spread of values for the matrix elements in our framework. All of the calculations were maximally self-consistent, including quantities such as nuclear radii and Q-values derived directly from our relativistic mean-field calculations. The implementation of the relativistic mean field theory in the description of ground state properties is based on the DIRHB code for the relativistic Hartree-Bogoliubov model[23] and alternately the relativistic Hartree-BCS model[24] and the code for QRPA already in use at the University of Zagreb[25]. The relativistic charge-exchange QRPA (PN-RQRPA) allows us to calculate the states in the intermediate nuclei for both kinds of double beta decay as excitations from an even-even nucleus. For the calculation of the final nuclear matrix elements, custom code has been written that connects the outputs of the two QRPA calculations. The strength of the isoscalar pairing in the residual interaction of the PN-RQRPA, which we describe using a dimensionless parameter V_0pp, is a crucial parameter of our model. For the 2νββ decay mode, we show both the dependence of the nuclear matrix element on V_0pp and the values of the NME at optimal T=0 pairing strength. The optimal values of V_0pp were determined following the approach of Ref. [26, 27, 28], based on a global fit of the half-lives of single β decay to the experimental data. We use even-even nuclides spanning 8 ≤ Z ≤ 82, for which experimental data on single β decays are available, and where the effect of changing V_0pp on the relative change of T_1/2 is larger than 20%, and where T_1/2 < 103 s. The ansatz we work with depends on (N −Z), introducing an isotopic dependence on V_0pp. The optimal values of V_0pp obtained in this way are mostly clustered around 1.0, showing a very weak breaking of the isospin symmetry. As this is close to the so-called breaking point of the PN-RQRPA[29], where the lowest root of the PN-RQRPA equation becomes imaginary and a new mean field description is necessary, for most nuclides, it contributes to an additional and considerable lowering of the NME from the values at V_0pp = 0. The problem of the value of the axial coupling constant, gA, that is essential in weak processes, is briefly discussed in the work; more information can be found, e.g., in Ref.[3]. Concerning the value of the gA parameter in our work, we have decided not to use excessively quenched values below gA = 1.0, as sometimes used in QRPA calculations[3, 30]. We either use the vacuum value, gA = 1.27, or a reasonably quenched value of gA = 1.0. An argument in favour of the latter can be made as it has proven to be the case that some quenching is necessary in order to reproduce single beta decay half-lives[31]. In the case of 2νββ decay we show the results for both values, although we use gA = 1.0 in tables and in comparisons with the results of other theoretical approaches. In the case of 0νββ decays, our calculations use gA = 1.27 in order to facilitate comparison to the results of calculations of 0νββ nuclear matrix elements in other theoretical frameworks. As is known, the closure approximation is not appropriate for the treatment of 2νββ decays, so in all cases except when explicitly noted, the 2νββ nuclear matrix elements have been calculated using the full energy denominator. The only cases where the closure approximation was used was in comparison with the Interacting Boson Model (IBM) results[32] which were themselves calculated in the closure approximation. The energy denominator is taken from Ref. [33] and differs for each nucleus. The 0νββ results were calculated in the closure approximation, which is appropriate for this mode of decay[34]. The average closure energies were the same as those used in the 2νββ case. Our results for the nuclear matrix elements for the 2νββ mode are in general smaller than the results of other approaches[3, 30, 35, 36, 37, 38, 39]. However we argue they are just as good if not improved when compared to the experimental NMEs[40] for known double beta decays. The two nuclides whose NMEs diverge from the experimental results the most, 76Ge and 82Se, can be shown to display marked triaxiality[41], which means that a more successful description of their 2νββ decays might have to wait for the development of a deformed version of the PN-RQRPA code. In addition to the nuclear matrix elements, we calculate the phase space factors and the halflives for those 2νββ nuclides where the data about the optimal V_0pp is available. The calculation of the phase space factors follows the approach of Ref.[42, 33], although we report new values of the PSF that do not correspond entirely to the results in the works cited. The calculated halflives show variance from the experimental results, although most of them are within an order of magnitude from the experimental T_1/2. Concerning the calculations of 0νββ nuclear matrix elements, the main feature that distinguishes this work from previous studies of neutrinoless double beta decay is that, following recent work mostly in the framework of the generator coordinate method[43, 44, 45], we extend our calculations to also use the full relativistic transition operator in addition to the transition operator that arises from a nonrelativistic reduction of that operator. As a result, instead of grouping the contributions to our nuclear matrix element into distinct Fermi, Gamow-Teller and tensor terms, we can directly show the contributions of the operators containing vector-vector (VV), axial-vector-axial-vector (AA), axial-vector-pseudovector (AP), pseudoscalar-pseudoscalar (PP) and weak magnetic (MP) couplings to the final result. We obtain results consistent with the fact, already reported[44], that the AA operator contributes the most to the final value of the nuclear matrix element, with the other couplings contributing at most around 1% of the total absolute value to the 0νββ nuclear matrix element, whether the fully relativistic operator is used or not, except in 116^Cd where the contributions go up to 7%. Our calculations of the 2νββ nuclear matrix elements converged at a maximum pn pair energy of around 100 MeV. Unfortunately we find that the 0νββ matrix element converges much slower, and calculations need to be extended to higher energies. Therefore the calculation of 0νββ NMEs requires considerable computing power, which is exacerbated by the fact that we have to sum over a large number of JΠ states of the intermediate nucleus. Since this is the first application of the PN-RQRPA-based model, we have decided to limit ourselves to a maximum value of the multipole of the transition to the intermediate nucleus of 6+, and to disregard the weak magnetic contribution (which we expect to be low[44]) from the fully relativistic calculation. Furthermore we have been unable to implement an easy way to take short-range correlations[34] into account, particularly since our calculation of the 0νββ matrix elements proceeds in a novel way. These and other limitations still leave our work as proof of concept that 0νββ NMEs can be calculated by directly calculating the contributions of various coupling channels, without recoupling the result into a matrix element between two-neutron and two-proton states, as is commonly done in the relevant literature[46]. Finally, the nuclear matrix elements we obtain are almost systematically lower than those calculated in other theoretical approaches, even those based on the pn-QRPA. This was already the case for the 2νββ nuclear matrix elements we obtain, but is more pronounced in the 0νββ case. Of course, in the latter case we also do not have any experimental data to compare our results to. Calculating the 0νββ matrix elements for the case in which isoscalar pairing vanishes does not help much as the dependence of the nuclear matrix elements on the strength parameter for isoscalar pairing V_0pp is rather slight, as expected[47]. Only for the DD-ME2 results we obtain an appreciable enhancement of the 0νββ NMEs. In fact, our results for the 0νββ nuclear matrix elements are closer to those of the interacting shell model than most of those reported from the nonrelativistic QRPA approaches. Curiously, a similar Skyrme-based (the similarity of Skyrme and point coupling interactions has been noted[48]) QRPA calculation by Terasaki reported even lower values of the 0νββ NME, which was blamed on excessive correlations within the QRPA[48] and the high-momentum components of the interaction. However, we believe that the low value of our results for the 0νββ nuclear matrix elements points to a need to extend the calculation to cover transitions of higher multipolarity, and further developments of the model going beyond this study, especially the inclusion of nuclear deformation effects. The ratio of our results for the nuclear matrix elements obtained using the nonrelativistic reduction of the transition operator and the results obtained by using the full relativistic transition operator also needs to be noted. Previous work within the GCM reported these two as being very similar, and therefore the nonrelativistic approximation as justified[45, 43]. In contrast, for some nuclides we find a ratio of the result obtained with the nonrelativistically reduced operator to the one obtained using a fully relativistic transition operator as low as 5%, closer to what was reported in Ref. [44] as an anomalous result for 150Nd, raising the question of how anomalous this result truly is and if the nonrelativistic reduction of the transition operator is as justified as it first appears to be. As stated, this work represents the first implementation of the theoretical framework based on the PN-RQRPA in the study of double beta decays. The perspectives for future research as continuation of the present work are numerous. In the immediate aftermath, we will focus on extending our calculations of the 0νββ nuclear matrix elements in the sense of extending it to transitions of higher multipolarity and inclusion of nuclear deformation effects. In addition, we can consider other ways of finding the optimal value of the isoscalar pairing strength parameter V_0pp, for decays proceeding through both modes of double beta decay. Further extensions are possible. In order to implement a viable treatment of short-range correlations within our 0νββ decay study, we can extend our treatment of the ground state wavefunction to include a so-called unitary correlation operator. This procedure, consequently referred to as the Unitary Correlation Operator Method (UCOM)[49, 50], does not violate the normalisation of the wave functions as the simpler Jastrow function approach does[51] and might be more appropriate for our approach. The PN-RQRPA code we use is based on a code that is able to treat charge-exchange excitations at finite temperature. Double beta decay at finite temperature has only been sporadically mentioned in the literature, which means that the effects of finite temperature, while possibly small, mostly remain unexplored. By opening new charge-exchange transitions in nuclei at finite temperature, the NMEs could become larger. A further possibility is to extend our calculation to the case of deformed nuclei. As mentioned earlier, already in the case of 2νββ decays we have knowledge about the importance of triaxiality to some of the nuclei being studied, while others have deformations that could be approximated with axial symmetry. The theory framework for a deformed PN-RQRPA is under development at the University of Zagreb; once complete we will be able to use it to extend our study to the treatment of deformed nuclei as well, although this represents an increase in the computational complexity of the problem, particularly for 0νββ decays [47]. Such an extension would then represent the first QRPA calculation of deformed nuclei using the full relativistic operator.