Is it always better to extract more factors when they exist? Unlike principal components analysis, the solutions to factor analysis models are not necessarily nested.  That is, the loadings (for example) for the first factor won't necessarily be identical when only the first factor is extracted vs. when the first two factors are.  
With that in mind, consider a case where you have a set of manifest variables that are highly correlated and (by theoretical knowledge of their content) should be driven by a single factor.  Imagine that exploratory factor analyses (by whichever metric you prefer: parallel analysis, scree plot, eigen values >1, etc.) strongly suggests that there are $2$ factors: A large primary factor, and a small secondary factor.  You are interested in using the manifest variables and the factor solution to estimate (i.e., get factor scores) participants' values for the first factor.  In this scenario, would it be better to: 


*

*Fit a factor model to extract only $1$ factor, and get factor scores (etc.), or

*fit a factor model to extract both factors, get factor scores for the factors, but throw away / ignore the scores for the second factor?  


For whichever is the better practice, why?  Is there any research on this issue?  
 A: The issue you're alluding to is the 'approximate unidimensionality' topic when building psychological testing instruments, which has been discussed in the liturature quite a bit in the 80's. The inspiration existed in the past because practitioners wanted to use traditional item response theory (IRT) models for their items, and at the time these IRT models were exclusively limited to measuring unidimensional traits. So, test multidimensionality was hoped to be a nuisance that (hopefully) could be avoided or ignored. This is also what led to the creation of the parallel analysis techniques in factor analysis (Drasgow and Parsons, 1983) and the DETECT methods. These methods were --- and still are --- useful because linear factor analysis (what you are referring to) can be a decent limited-information proxy to full-information factor analysis for categorical data (which is what IRT is at its core), and in some cases (e.g., when polychoric correlations are used with a weighted estimator, such as WLSMV or DWLS) can even be asymptotically equivalent for a small selection of ordinal IRT models. 
The consequences of ignoring additional traits/factors, other than obviously fitting the wrong model to the data (i.e., ignoring information about potential model misfit; though it may of course be trivial), is that trait estimates on the dominant factor will become biased and therefore less efficient. These conclusions are of course dependent on how the properties of the additional traits (e.g., are they correlated with the primary dimension, do they have strong loadings, how many cross-loadings are there, etc), but the general theme is that secondary estimates for obtaining primary trait scores will be less effective. See the technical report here for a comparison between a miss-fitted unidimensional model and a bi-factor model; the technical report appears to be exactly what you are after.
From a practical perspective, using information criteria can be helpful when selecting the most optimal model, as well as model-fit statistics in general (RMSEA, CFI, etc) because the consequences of ignoring multidimensional information will negatively affect the overall fit to the data. But of course, overall model fit is only one indication of using an inappropriate model for the data at hand; it's entirely possible that improper functional forms are used, such as non-linearity or lack of monotonicity, so the respective items/variables should always be inspected as well.   
See also:
Drasgow, F. and Parsons, C. K. (1983). Application of Unidimensional Item Response Theory Models to Multidimensional Data. Applied Psychological Measurement, 7 (2), 189-199.
Drasgow, F. & Lissak, R. I. (1983). Modified parallel analysis: A procedure for examining the latent-dimensionality of dichotomously scored item responses. Journal of Applied Psychology, 68, 363-373.
Levent Kirisci, Tse-chi Hsu, and Lifa Yu (2001). Robustness of Item Parameter Estimation Programs to Assumptions of Unidimensionality and Normality. Applied Psychological Measurement, 25 (2), 146-162.
A: If you truly do not want to use the second factor, you should just use a one-factor model. But I am puzzled by your remark that the loadings for the first factor will change if you use a second factor.
Let's deal with that statement first. If you use principal components to extract the factors and do not use factor rotation, then the loadings will not change -- subject perhaps to scaling (or complete flipping: If $x$ is a factor, then $-x$ is a legitimate way to express it as well). If you use maximum likelihood extraction and/or factor rotations, then the loadings may depend on the number of factors you extracted.
Next, for the explanation of the effects of rotations. I am not good at drawing, so I will try to convince you using words. I will assume that your data are (approximately) normal, so that the factor scores are approximately normal also. If you extract one factor, you get a one-dimensional normal distribution, if you extract two factors, you get a bivariate normal distribution. 
The density of a bivariate distribution looks roughly speaking like a hat, but the exact shape depends on scaling as well as the correlation coefficient. So let's assume that the two components each have unit variance. In the uncorrelated case, you get a nice sombrero, with level curves that look like circles. A picture is here. Correlation "squashes" the hat, so that it looks more like a Napoleon hat. 
Let's assume that your original data set had three dimensions and yu want to extract two factors out of that. Let's also stick with normality. In this case the density is a four-dimensional object, but the level curves are three-dimensional and can at least be visualized. In the uncorrelated case the level curves are spherical (like a soccer ball). In the presence of correlation, the level curves will again be distorted, into a football, probably an underinflated one, so that the thickness at the seams is smaller than the thickness in the other directions.
If you extract two factors using PCA, you completely flatten the football into an ellipse (and you project every data point onto the plane of the ellipse). The unrotated first factor corresponds to the long axis of the ellipse, the second factor is perpendicular to it (i.e., the short axis). Rotation then chooses a coordinate system within this ellipse in order to satisfy some other handy criteria.
If you extract just a single factor, rotation is impossible, but you are guaranteed that the extracted PCA factor corresponds to the long axis of the ellipse.  
A: Why would you not use something like lavaan or MPlus to run two models (unidimensional model and a two dimension model aligned to your EFA results) and compare the relative and absolute fit indices of the different models (i.e., information criteria - AIC and BIC, RMSEA, SRMR, CFI/TLI)? Note that if you go down this road you would not want to use PCA for the EFA, but rather principal factors. Somebody really concerned with measurement would embed the CFA into a full structural equation model.
Edit: The approach I'm asking you to consider is more about figuring out how many latent variables actually explain the set of items. If you want to get the best estimate of the larger factor, I would vote for using the factor scores from the CFA model with the better fit, whichever that is.
