How to create a scree plot for factor analysis given that eigenvalues depend on the number of extracted factors?

Question

I understand how Kaiser rule works for PCA, as no matter how many components I extract I always get the same eigenvalues.

For example, with 3 components I get

principal(r = corMat, nfactors = 3, rotate = "oblique")
....
....
                        PC1  PC2  PC3
SS loadings           12.62 3.99 2.57
Proportion Var         0.37 0.12 0.08
Cumulative Var         0.37 0.49 0.56
Proportion Explained   0.66 0.21 0.13
Cumulative Proportion  0.66 0.87 1.00

For six components, I basically just get additional columns

principal(r = corMat, nfactors = 6, rotate = "oblique")
....
....
                        PC1  PC2  PC3  PC4  PC5  PC6
SS loadings           12.62 3.99 2.57 1.39 1.13 1.09
Proportion Var         0.37 0.12 0.08 0.04 0.03 0.03
Cumulative Var         0.37 0.49 0.56 0.60 0.64 0.67
Proportion Explained   0.55 0.18 0.11 0.06 0.05 0.05
Cumulative Proportion  0.55 0.73 0.84 0.90 0.95 1.00

Thus, first you extract all of them, plot & pick, or use Kaiser rule > 1 or something like that..that can be done using nFactors package, it is giving the same numbers as principal function

principal(r = corMat, nfactors = 34, rotate = "oblique")
....
....
                        PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29 PC30 PC31 PC32 PC33 PC34
SS loadings           12.62 3.99 2.57 1.39 1.13 1.09 0.84 0.74 0.72 0.69 0.59 0.56 0.55 0.51 0.47 0.42 0.40 0.39 0.36 0.35 0.34 0.33 0.32 0.31 0.30 0.28 0.26 0.25 0.24 0.23 0.23 0.22 0.19 0.11
Proportion Var         0.37 0.12 0.08 0.04 0.03 0.03 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.00
Cumulative Var         0.37 0.49 0.56 0.60 0.64 0.67 0.69 0.72 0.74 0.76 0.78 0.79 0.81 0.82 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.96 0.97 0.98 0.98 0.99 1.00 1.00
Proportion Explained   0.37 0.12 0.08 0.04 0.03 0.03 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.00
Cumulative Proportion  0.37 0.49 0.56 0.60 0.64 0.67 0.69 0.72 0.74 0.76 0.78 0.79 0.81 0.82 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.96 0.97 0.98 0.98 0.99 1.00 1.00

However, with factor analysis (for example principal axis factoring), every time I specify a different number of factors, the eigenvalues which I get back are slightly different.

For example, if I specify 3 factors I get the following:

> fa(r = corMat, nfactors = 3, rotate = "oblique", fm="pa")
.....
.....

                        PA1  PA2  PA3
SS loadings           12.16 3.57 2.10
Proportion Var         0.36 0.10 0.06
Cumulative Var         0.36 0.46 0.52
Proportion Explained   0.68 0.20 0.12
Cumulative Proportion  0.68 0.88 1.00

but for six factors I get slightly different values

> fa(r = corMat, nfactors = 6, rotate = "oblique", fm="pa")
.....
.....    
                        PA1  PA2  PA3  PA4  PA5  PA6
SS loadings           12.23 3.65 2.17 1.04 0.76 0.73
Proportion Var         0.36 0.11 0.06 0.03 0.02 0.02
Cumulative Var         0.36 0.47 0.53 0.56 0.58 0.61
Proportion Explained   0.59 0.18 0.11 0.05 0.04 0.04
Cumulative Proportion  0.59 0.77 0.88 0.93 0.96 1.00

In many papers, I've seen that they selected the number of factors based on kaiser rule, and then show a scree plot. But I'm confused where the values for the scree plot come from?

principal axis factoring with Oblimin rotations was carried out. We attempted four and three-factor solutions. Both the Kaiser rule of eigenvalues greater than 1 and the scree plot (see Fig. 1) indicated that three-factor solution would fit the data the best

and then they show a typical scree plot. Do they first run a PCA to pick the number of factors or something else?